Erin_TestScript.Rmd

---
title: "Erin_TestScript"
author: "Erin Grey"
date: "2024-02-13"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
rm(list=ls())
```

### DEFINE THE FOLLOWING VARIABLES FOR YOUR ANALYSIS
## locus choices are: (GENES) "ATP6" "ATP8"  "COI"  "COII"  "COIII"  "CYTB"  "ND1"  "ND2"  "ND3"  "ND4" "ND4L" "ND5"  "ND6"  and (OTHER FEATURES) "rRNA_12S" "rRNA_16S" "D_loop" "tRNA_Ala" "tRNA_Arg" "tRNA_Asn" "tRNA_Asp" "tRNA_Cys" "tRNA_Gln" "tRNA_Glu" "tRNA_Gly" "tRNA_His" "tRNA_Ile" "tRNA_Leu" "tRNA_Lys" "tRNA_Met" "tRNA_Phe" "tRNA_Pro" "tRNA_Ser" "tRNA_Thr" "tRNA_Trp" "tRNA_Tyr" "tRNA_Val"
## species_list must have "search_name" as the first column and it must be a binomial species name
```{r define_variables}
entrez_key <- "You_Entrez_Key_Here" #GET YOUR OWN ENTREZ KEY AND PUT IT HERE!
locus = "COI" #name of target locus, your choices are 
output_folder <- "Test_Erin-2024-02-13" #name of your output folder every time
species_list <- read.csv("workingfiles/MaineSpeciesList_Clean-EKG.csv") # your species list
colnames(species_list)<- c("search_name") 
order_list <- read.csv("workingfiles/MetazoaSpeciesByOrder_2023-09-19.csv") # list of vertebrate orders
```

### LOAD PACKAGES, DEFINE TERM, AND CREATE DATABASE SKELETONS. Note - you need to locally install an older version of genbankr because it's not up to date on Bioconductor - go here https://bioconductor.org/packages/3.16/bioc/html/genbankr.html and download the appropiate binary then go to your terminal and "R CMD INSTALL pathtobinary.tar.gz". There will be other dependencies that you'll have to install too. Sorry. The genbankr authors say they will fix this by April
```{r load_packages_and_terms}
library(taxizedb) #checks taxonomy
library(rentrez) #queries ENTREZ databases and downloads accessions
library(AnnotationBustR) #finds longest accessions, slice genes from mitogenomes
library(reutils) #other packages need it
library(ape) #convert fasta, fastq, etc.
library(ggplot2) #plots
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("Biostrings")
library(Biostrings)
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("genbankr")
library(genbankr) ##genbankr is not functioning on latest Bioconductor, install locally instead
set_entrez_key(entrez_key) #set the Entrez API key
dir.create(output_folder)

# create search terms for ENTREZ
data(mtDNAterms) #AnnotationBustR's list of synonyms for different loci
more_12Ssynonyms <- data.frame(Locus="rRNA_12S", Type="rRNA", Name= "small ribosomal RNA subunit RNA") # other synonyms that we find go here and get added to AnnotationBustR's list
mtDNAterms <- rbind(mtDNAterms, more_12Ssynonyms) #format
is_mtgene <- locus %in% c("ATP6", "ATP8",  "COI",  "COII",  "COIII",  "CYTB",  "ND1",  "ND2",  "ND3",  "ND4", "ND4L", "ND5",  "ND6") #check whether the locus is a gene or other feature, needed for the mitogenome scrape
target_locus_synonyms <- mtDNAterms[mtDNAterms$Locus==locus,] #the target synonyms
target_locus_synonyms$Terms <- paste0("OR ", target_locus_synonyms$Name, "[TITL]") # format for ENTREZ search terms
target_locus_synonyms$Terms[1] <- paste0("AND (", target_locus_synonyms$Name[1], "[TITL]") # first term starts with "AND ("
target_locus_synonyms$Terms[dim(target_locus_synonyms)[1]] <- paste0("OR ", target_locus_synonyms$Name[dim(target_locus_synonyms)[1]], "[TITL])") #last term ends with a ")"
target_locus_searchterm <- paste(as.vector(target_locus_synonyms$Terms), collapse=" ") # the big ENTREZ search term

# create empty dataframes to fill in during loops
tax_df <- as.data.frame(matrix(nrow=1, ncol=9, dimnames=list(NULL, c("tax_query", "superkingdom", "kingdom", "phylum", "class", "order", "family", "genus", "species")))) # taxonomy skeleton
order_seqs <- as.data.frame(matrix(nrow=1, ncol=8, dimnames=list(NULL, c("superkingdom", "kingdom", "phylum", "class", "order", "species_id", "ids_mitogenome","ids_target")))) #  missing orders skeleton
a02_REFDB <- data.frame(seq_header=NA, sequence=NA, seq_accession=NA, type=NA, species=NA) #create reference database skeleton
t = 1 # taxonomy for species list loop counter
u = 1  # taxonomy for order sequences loop counter
o = 1 # missing orders loop counter
j = 1 # missing orders species loop counter
i = 1 # ENTREZ search loop counter
s = 1 # sequence scrape loop counter
m = 1 # mitogenome scrape loop counter
z = 1 # output format loop counter
d = 1 # dada2 taxonomy counter
```

## Get taxonomy for each species in the list.
```{r taxonomy}
species_list_dedup <- unique(species_list$search_name) #de-duplicated the species list
db_tax_NCBI <- db_download_ncbi(verbose = TRUE, overwrite = FALSE) #download the NCBI taxonomy database
taxonomies_cls <- taxizedb::classification(species_list_dedup, db="ncbi") #get full taxonomies for each species, output is a list-like "classification" object that sucks; had to remove Rondeletia bicolor
#now loop through the sucky classification object to populate the empty dataframe
for (t in 1:length(taxonomies_cls)) {
  x <- as.data.frame(taxonomies_cls[t])
  tax_query <- strsplit(colnames(x)[1], split="[.]")
  tax_query <- paste(tax_query[[1]][1], tax_query[[1]][2], sep=" ")
  # update the taxonomy dataframe skeleton
  tax_df[nrow(tax_df)+1,] <- rep("NA", dim(tax_df)[2])
  tax_df$tax_query[t] <- tax_query
  if (dim(x)[1]>1) {
        c00 <- paste0(x[which(x[,2]=="superkingdom"),c(1,3)], collapse="_")
        c0 <- paste0(x[which(x[,2]=="kingdom"),c(1,3)], collapse="_")
        c1 <- paste0(x[which(x[,2]=="phylum"),c(1,3)], collapse="_")
        c2 <- paste0(x[which(x[,2]=="class"),c(1,3)], collapse="_")
        c3 <- paste0(x[which(x[,2]=="order"),c(1,3)], collapse="_")
        c4 <- paste0(x[which(x[,2]=="family"),c(1,3)], collapse="_")
        c5 <- paste0(x[which(x[,2]=="genus"),c(1,3)], collapse="_")
        c6 <- paste0(x[which(x[,2]=="species"),c(1,3)][1,], collapse="_")
  }
  
  if (exists("c00")){
      if (length(c1) > 0) {
        tax_df$superkingdom[t] <- c00
        } else {
        tax_df$superkingdom[t] <- "na"    
        }
  } 
  
  if (exists("c0")){
      if (length(c1) > 0) {
        tax_df$kingdom[t] <- c0
        } else {
        tax_df$kingdom[t] <- "na"    
        }
  } 
  
  if (exists("c1")){
      if (length(c1) > 0) {
        tax_df$phylum[t] <- c1
        } else {
        tax_df$phylum[t] <- "na"    
        }
    } 
  
  if (exists("c2")){
      if (length(c2) > 0) {
        tax_df$class[t] <- c2
        } else {
        tax_df$class[t] <- "na"    
        }
    } 
  
  if (exists("c3")){
      if (length(c3) > 0) {
        tax_df$order[t] <- c3
        } else {
        tax_df$order[t] <- "na"    
        }
    } 
  
   if (exists("c4")){
      if (length(c4) > 0) {
        tax_df$family[t] <- c4
        } else {
        tax_df$family[t] <- "na"    
        }
    } 
  
  if (exists("c5")){
      if (length(c5) > 0) {
        tax_df$genus[t] <- c5
        } else {
        tax_df$genus[t] <- "na"    
        }
    } 
  
   if (exists("c6")){
      if (length(c6) > 0) {
        tax_df$species[t] <- c6
        } else {
        tax_df$species[t] <- "na"    
        }
    } 
  
 suppressWarnings(rm(list = c("x", "tax_query", "c00", "c0", "c1","c2","c3","c4","c5","c6")))
}

a01_NAMES <- merge(species_list, tax_df, by.x="search_name", by.y="tax_query", all.x=TRUE)
```

## Remove species without taxonomies and identify metazoan orders that aren't represented in your list
```{r tax_subset_and_missing_orders}
a01_NAMES[a01_NAMES == "NA"] <- NA #format
a01_NAMES_missing <- droplevels(subset(a01_NAMES, is.na(a01_NAMES$species))) #names that are missing genbank taxonomy (could be misspelled, or just not in genbank)
a01_NAMES <- droplevels(subset(a01_NAMES, !is.na(a01_NAMES$species))) #use only the species with species names moving forward
a01_NAMES$species_id <- data.frame(do.call(rbind, strsplit(a01_NAMES$species, split = "_")))[,2] #format for future merging
orders_missing <- droplevels(subset(order_list, !(order_list$order %in% unique(factor(a01_NAMES$order))))) # orders that are missing from species list
row.names(orders_missing) = NULL # format
write.csv(a01_NAMES, file.path(output_folder, "a01_NAMES_wTaxonomy.csv"), row.names = FALSE)
write.csv(a01_NAMES_missing, file.path(output_folder, "a01_NAMES_withoutTaxonomy.csv"), row.names = FALSE)
#a01_NAMES <- droplevels(subset(a01_NAMES,kingdom == "Metazoa_33208")) #subset if desired
```

## Find up to 3 mitogenomes or target sequences for missing orders. If you run into an HTTP 500 error, just start the loop again and you should be good.
```{r missing_orders_seqs}
# find sequences for up to three species for each missing order; search mitogenomes first, then regular accessions
while (o <= dim(orders_missing)[1]) { # for every missing order
    cat("\r", "finding ref seqs for order", o, "of", dim(orders_missing)[1]) #counter
    STOPIT <- "NO" # set STOPIT to NO for this order
    while (STOPIT == "NO"){ # keep searching as long as STOPIT = "NO"
        species_ids_ls <- strsplit(orders_missing$spp_list[[o]], ";") # get the species ids for that order
        if (is.na(species_ids_ls[[1]][1])) {STOPIT <- "YES"} #IF THERE ARE NO GENBANK SPECIES FOR THIS ORDER, STOP SEARCHING FOR IT
        if (length(species_ids_ls[[1]])<=3) { 
            species_ids <- species_ids_ls[[1]] #format species ids
            for (j in 1:length(species_ids)){ # for each species
                search_name <- paste0("txid",species_ids[j],"[Organism]")
                mitogenomes <- tryCatch(entrez_search(db="nucleotide", term <- paste(search_name, "AND mitochondrion[TITL] AND complete genome[TITL]"), retmax=9999)) # search mitogenomes
                if(class(mitogenomes)[1] != "try-catch") { #if the search went through
                    if (length(mitogenomes$ids)>0) {
                    mito_id <- sample(mitogenomes$ids,1) #choose a random mitogenome for this species
                    row_info <- c(orders_missing$superkingdom[o], orders_missing$kingdom[o], orders_missing$phylum[o], orders_missing$class[o], orders_missing$order[o], species_ids[i], mito_id, "NA")
                    order_seqs <- rbind(order_seqs, row_info)
                    rm(list=c("mito_id","row_info"))
                    Sys.sleep(1) #slow down request to the Entrez server or you'll get kicked out
                    } else {
                        targets <- tryCatch(entrez_search(db="nucleotide", term <- paste(search_name, target_locus_searchterm, collapse=" "), retmax=999999))
                        if(class(targets)[1] != "try-catch") {
                            if (length(targets$ids)>0) {
                                target_id <- sample(targets$ids,1) #choose a random accession for this species
                                row_info <- c(orders_missing$superkingdom[o], orders_missing$kingdom[o], orders_missing$phylum[o], orders_missing$class[o], orders_missing$order[o],species_ids[j], "NA", target_id)
                                order_seqs <- rbind(order_seqs,row_info)
                                rm(list=c("target_id","row_info"))
                                Sys.sleep(1) #slow down request to the Entrez server or you'll get kicked out
                                }
                            }
                            rm(targets)
                        }
                        rm(search_name); rm(mitogenomes)
                    }
                    STOPIT <- "YES"
                }
            }   #If there are 3 or fewer species, search them all for a mitogenome or accession then STOP SEARCHING
        if (length(species_ids_ls[[1]])>3){ 
            species_ids <- sample(species_ids_ls[[1]]) #randomize the species for each order
            finds = 0 #start with zero finds
            for (j in 1:length(species_ids)){
                search_name <- paste0("txid",species_ids[j],"[Organism]")
                mitogenomes <- tryCatch(entrez_search(db="nucleotide", term <- paste(search_name, "AND mitochondrion[TITL] AND complete genome[TITL]"), retmax=9999)) # search mitogenomes
                if(class(mitogenomes)[1] != "try-catch") { #if the search went through
                    if (length(mitogenomes$ids)>0) {
                        mito_id <- sample(mitogenomes$ids,1) #choose a random mitogenome for this species
                        row_info <- c(orders_missing$superkingdom[o], orders_missing$kingdom[o], orders_missing$phylum[o], orders_missing$class[o], orders_missing$order[o], species_ids[j], mito_id, "NA")
                        order_seqs <- rbind(order_seqs,row_info)
                        finds=finds+1 #add to find count
                        rm(mito_id) 
                        Sys.sleep(1) #slow down request to the Entrez server or you'll get kicked out
                    } else {
                            targets <- tryCatch(entrez_search(db="nucleotide", term <- paste(search_name, target_locus_searchterm, collapse=" "), retmax=999999))
                            if(class(targets)[1] != "try-catch") {
                                if (length(targets$ids)>0) {
                                    target_id <- sample(targets$ids,1) #choose a random accession for this species
                                    row_info <- c(orders_missing$superkingdom[o], orders_missing$kingdom[o], orders_missing$phylum[o], orders_missing$class[o], orders_missing$order[o],species_ids[j], "NA", target_id)
                                    order_seqs <- rbind(order_seqs,row_info)
                                    finds=finds+1 #add to find count
                                    rm(target_id); rm(row_info)
                                    Sys.sleep(1) #slow down request to the Entrez server or you'll get kicked out
                                    }
                                }
                                rm(targets)
                            }
                            rm(search_name); rm(mitogenomes)
                            if (finds >= 3) break
                            STOPIT <- "YES"
                        }
                } #If there are >=3 species, randomize and search mitogenomes or accessions until you get 3 or have searched them all
            }
            rm(species_ids_ls)
            o=o+1
    }
}
order_seqs <- na.omit(order_seqs) #subtract the first row because its all NAs
order_seqs$search_name<-"NA"; order_seqs$n_mitogenome<-"NA"; order_seqs$n_target<-"NA"; order_seqs$species<-"NA"; order_seqs$family<-"NA"; order_seqs$genus<-"NA" #make column names same as a01_NAME for merging
order_seqs$n_mitogenome <- ifelse(order_seqs$ids_mitogenome!="NA",1,0) #update the n_mitogenome column
order_seqs$n_target <- ifelse(order_seqs$ids_target!="NA",1,0) #update the n_target column
```

## Get taxonomy for each species representative for the missing orders.
```{r taxonomy_for_order_seqs}
taxonomies_orderseqs <- taxizedb::classification(order_seqs$species_id, db="ncbi") #get full taxonomies for each species selected to represent missing orders, output is a list-like "classification" object that sucks
#now loop through the sucky classification object to populate the empty dataframe
for (u in 1:length(taxonomies_orderseqs)) {
  x <- as.data.frame(taxonomies_orderseqs[u])
  
  if (dim(x)[1]>1) {
        c4 <- paste0(x[which(x[,2]=="family"),c(1,3)], collapse="_")
        c5 <- paste0(x[which(x[,2]=="genus"),c(1,3)], collapse="_")
        c6 <- paste0(x[which(x[,2]=="species"),c(1,3)][1,], collapse="_")
  }

   if (exists("c4")){
      if (length(c4) > 0) {
        order_seqs$family[u] <- c4
        } else {
        order_seqs$family[u] <- NA    
        }
    } 
  
  if (exists("c5")){
      if (length(c5) > 0) {
         order_seqs$genus[u] <- c5
        } else {
         order_seqs$genus[u] <- NA    
        }
    } 
  
   if (exists("c6")){
      if (length(c6) > 0) {
        order_seqs$species[u] <- c6
        } else {
        order_seqs$species[u] <- NA    
        }
    } 
  
 suppressWarnings(rm(list = c("x", "tax_query", "c4","c5","c6")))
}
order_seqs$search_name <- data.frame(do.call(rbind, strsplit(order_seqs$species, split = "_")))[,1]
```

## Search Entrez for mitogenome and target sequences for each species. If you run into an HTTP 500 error, just start the loop again and you should be good.
```{r search_entrez}
while (i <= dim(a01_NAMES)[1]){
  cat("\r", "finding Genbank accessions & mitogenomes for species", i, "of", dim(a01_NAMES)[1]) #counter
  # define search terms for species
  search_name <- paste0(a01_NAMES$search_name[i],"[ORGN]") #format species name for ENTREZ search
  search_term <- paste(search_name, target_locus_searchterm, collapse=" ") #mash species & locus terms into one search term
  # search ENTREZ for mitogenomes and then for accessions
   tryCatch(entrez_search(db="nucleotide", term <- paste(search_name, "AND mitochondrion[TITL] AND complete genome[TITL]"), retmax=9999)) -> mitogenomes # search mitogenomes
      if(class(mitogenomes)[1] != "try-catch") {
        a01_NAMES$n_mitogenome[i] <-mitogenomes$count #add mitogenome count names dataframe
        a01_NAMES$ids_mitogenome[i] <-paste(mitogenomes$ids, collapse="|") #add mitogenome ids to name dataframe
            tryCatch(entrez_search(db="nucleotide", term <- search_term, retmax=9999)) -> targets # search accessions
                if(class(targets)[1] != "try-catch") {
                    a01_NAMES$n_target[i] <-targets$count
                    a01_NAMES$ids_target[i] <- paste(targets$ids, collapse="|")
                }
            }
  Sys.sleep(0.5) #slow down request to the Entrez server or you'll get kicked out
  # reset loop variables
  mitogenomes <- "na"
  targets <- "na"
  search_name <- "na"
  search_term <- "na"
  i <- i + 1
}

order_seqs <- order_seqs[,colnames(a01_NAMES)] # format the order_seqs dataframe to add to the a01_NAMES datafame
a01_NAMES<- rbind(a01_NAMES, order_seqs) # add the order_seqs dataframe to the a03BESTANAMES dataframe
a01_NAMES$n_target <- as.numeric(a01_NAMES$n_target) #format
a01_NAMES$n_mitogenome <- as.numeric(a01_NAMES$n_mitogenome) #format
row.names(a01_NAMES) = NULL # format
write.csv(a01_NAMES, file.path(output_folder, "a01_NAMES_v2.csv"), row.names = FALSE) # save the mitogenome and target accessions ids for scraping
```

## Scrape accessions identified above from GenBank.
```{r scrape_accessions}
while (s <= dim(a01_NAMES)[1]){ #for every good species name
    cat("\r","scraping accessions for species", s, "of",dim(a01_NAMES)[1])
    ids <- "na"
    seqs_target <- "na"
    if (a01_NAMES$n_target[s]>0 && a01_NAMES$n_target[s]<100) { # scrape GenBank target sequences if available, but skip if >= 100 targets
     ids <- c(unlist(strsplit(a01_NAMES$ids_target[s], split="\\|")))
      } else if  (a01_NAMES$n_target[s]>100) {   # if more than 200 accessions, randomly select 100  
     ids <- sample(c(unlist(strsplit(a01_NAMES$ids_target[s], split="\\|"))),100)
      }
    if (ids[1] !="na"){ # if there are accessions, fetch them from GenBank
        seqs_target <- tryCatch(entrez_fetch(db="nuccore", id=ids, rettype="fasta"))
        }
        if(class(seqs_target) != "try-catch"){
            if (seqs_target != "na"){
            write(seqs_target, file.path(output_folder, paste(a01_NAMES$search_name[s], paste0(locus, ".fasta")))) # formatting - write out the sequences
            fasta_target <- readDNAStringSet(file.path(output_folder, paste(a01_NAMES$search_name[s], paste0(locus, ".fasta"))), format="fasta") #formatting - read them back in as fasta
            seqs_target_accessions <- entrez_fetch(db="nuccore", id=ids, rettype="acc") # get accession numbers
            seq_header <- names(fasta_target) #formatting
            sequence <- paste(fasta_target) #formatting
            seq_accession <- unlist(strsplit(seqs_target_accessions, split="\n")) # formatting
            tempDB <- data.frame(seq_header, sequence, seq_accession, type="accession", species=a01_NAMES$search_name[s]) # make a temporary database with all sequences, their header, accession number, etc.
            a02_REFDB <- rbind(a02_REFDB, tempDB) # append temporary database to the full database
                }
            s=s+1
            }
    #reset loop variables
    fasta_target<-"na"; seqs_target_accessions<-"na"; seq_header<-"na"; sequence<-"na"; seq_accession<-"na";
    tempDB <-"na"
    #slow down request to the Entrez server or you'll get kicked out
    Sys.sleep(1)
 }

```

## Scrape your target locus out of mitogenomes identified in the Entrez search. Note - I need to fix some things here (the loop breaks on parsing errors). If you get stuck, just "m=m+1" and re-start loop
```{r scrape_mitogenomes}
while (m <= dim(a01_NAMES)[1]) { #for every good species name
    cat("\r","scraping mitogenomes for species", m, "of", dim(a01_NAMES)[1])
    mito_ids <- "na"
    if (a01_NAMES$n_mitogenome[m]>0 && a01_NAMES$n_mitogenome[m]<20) { #if mitogenomes available and <20
        mito_ids <- unlist(strsplit(a01_NAMES$ids_mitogenome[m], split="\\|")) #format ids
    } else if (a01_NAMES$n_mitogenome[m]>20) { #if >20 mitogenomes, subsample 20 mitogenomes randomly
        mito_ids <- sample(unlist(strsplit(a01_NAMES$ids_mitogenome[m], split="\\|")),20) # format ids  
    }
   
    if (mito_ids[1] != "na") { #if there are mitogenome ids for species m
        mito_accessions <- tryCatch(entrez_fetch(mito_ids, db="nuccore", rettype="acc")) #ask ENTREZ for id's accession
        if (class(mito_accessions) != "try-catch"){ #if no api error
            mito_accessions <- unlist(strsplit(mito_accessions, split="\n")) #format accession numbers
            for (n in 1:length(mito_accessions)){ # loop through and scrape each mitogenome accession
              new_row <- c(paste("Unparsed mitochondrion", mito_accessions[n], sep=" "), "na", mito_accessions[n], "scrape", species=a01_NAMES$search_name[m]) #blank row for each accession
              gb <- tryCatch(readGenBank(GBAccession(mito_accessions[n]))) # get the Genbank annotation for accession
                if (class(gb) != "try-catch"){ #if no error in getting gb
                  if (is_mtgene==FALSE) {# if target locus is not a gene, look in otherFeatures()
                  target_feature <- tryCatch(which(otherFeatures(gb)$product %in% as.character(target_locus_synonyms$Name))) # find target locus annotation metadata 
                    if (class(target_feature) != "try-catch"){
                      if(length(target_feature) > 0) { # if target feature is found
                        target_range <- tryCatch(otherFeatures(gb)@ranges[target_feature]) #extract the target range info
                          if(class(target_range) != "try-catch"){
                            target_seq <- tryCatch(subseq(getSeq(gb), start=target_range@start, width=target_range@width)) #scrape seq
                              if(class(target_seq) != "try-catch"){
                                scraped_seq <- paste(target_seq) #format
                                new_row <- c(paste(names(target_seq),"mitochondrion", mito_accessions[n], sep=" "), paste(target_seq), mito_accessions[n], "scrape", a01_NAMES$search_name[m]) #update information
                            }
                          }
                      }
                    }
                  }
                  else if (is_mtgene==TRUE){ # if target locus is a gene, look in gene()
                    target_feature <- tryCatch(which(genes(gb)$gene %in% as.character(target_locus_synonyms$Name))) # find target locus annotation metadata 
                    if (class(target_feature) != "try-catch"){
                      if(length(target_feature) > 0) { # if target feature is found
                        target_range <- tryCatch(genes(gb)@ranges[target_feature])  #extract the target range info
                          if(class(target_range) != "try-catch"){
                            target_seq <- tryCatch(subseq(getSeq(gb), start=target_range@start, width=target_range@width)) #scrape seq
                              if(class(target_seq) != "try-catch"){
                                scraped_seq <- paste(target_seq) #format
                                new_row <- c(paste(names(target_seq),"mitochondrion", mito_accessions[n], sep=" "), paste(target_seq), mito_accessions[n], "scrape", a01_NAMES$search_name[m]) #update information
                            }
                          }
                      }
                    }
                  }
                }
              a02_REFDB <- rbind(a02_REFDB, new_row) # update the database
              rm(gb, target_feature, target_range, target_seq, scraped_seq, new_row) # reset n loop variables
              Sys.sleep(0.5) #slow down request to the Entrez server or you'll get kicked out
            } # close n loop (each "n" mitogenome accession per "m" species)
        }
    }
    m=m+1 #update species m variable
    rm(mito_ids, mito_accessions) # reset loop variables
    Sys.sleep(0.5) #slow down request to the Entrez server or you'll get kicked out
} # close for each species m loop

a02_REFDB <- a02_REFDB[-1,] #format - remove the top row of NAs
a02_REFDB_unparsed <- subset(a02_REFDB, sequence == "na")
a02_REFDB_parsed <- subset(a02_REFDB, sequence != "na")
write.csv(a02_REFDB_parsed, file.path(output_folder, "a02_REFDB.csv"), row.names=FALSE)
```

## De-duplicate the reference database, but keep track of duplicate accessions/scrapes in Names
```{r remove_duplicates}
a03_UNIQUEDB <- a02_REFDB_parsed[!duplicated(a02_REFDB_parsed$sequence),] #remove duplicates
for (z in  1:dim(a03_UNIQUEDB)[1]){ # for every row in the unique db file
  dups <- subset(a02_REFDB_parsed, sequence == a03_UNIQUEDB$sequence[z]) # find identical sequences in the full db file 
  a03_UNIQUEDB$duplicate_accessions[z] <- paste(dups$seq_accession, collapse = "|") # paste all those accessions together into a new unique db field
  a03_UNIQUEDB$duplicate_species[z] <- paste(dups[!duplicated(dups$species),"species"], collapse = "|")
  z=z+1
}
row.names(a03_UNIQUEDB) <- NULL #format
write.csv(a03_UNIQUEDB, file.path(output_folder, "a03_UniqueRefDB.csv"), row.names=FALSE)
```

## Format the reference library for dada2 assignTaxonomy()
```{r format_dada2}
a04_UNIQUEDBdada <- merge(x=a03_UNIQUEDB, y=a01_NAMES, by.x="species", by.y="search_name", all.x=TRUE, all.y=FALSE) 
my_cols <- c("phylum", "class","order","family","genus","species.y")
a04_UNIQUEDBdada$header <- do.call(paste, c(a04_UNIQUEDBdada[my_cols], sep = ";"))
a04_UNIQUEDBdada$header2 <- paste(">",a04_UNIQUEDBdada$header, sep="")
REFDB_dada <- a04_UNIQUEDBdada[,c("header2","sequence")]
write.table(REFDB_dada, file.path(output_folder, "a04_REFDBdada.fasta"), sep="\n", col.names=FALSE, row.names=FALSE, quote=FALSE)
```

## Summarize taxonomic & reference sequence results for species in the original list
```{r summarize_species}
#Count reference sequences by species
counts_species <- as.data.frame(table(a03_UNIQUEDB$species)) #count the number of unique reference sequences per species
colnames(counts_species) <- c("species", "n_unique_seqs") #format
a01_NAMES <- merge(a01_NAMES, counts_species, by.x= "search_name", by.y="species", all.x=TRUE, all.y=TRUE) # update the summary database with reference sequence counts
a00_SPECIES_SUMMARY <- a01_NAMES[a01_NAMES$search_name %in% species_list$search_name,]

#add in duplicates information to BESTNAMES
a00_SPECIES_SUMMARY$dup_accessions <- "na"
a00_SPECIES_SUMMARY$dup_species <- "na"
a00_SPECIES_SUMMARY$dup_species_n <- "na"

for (b in 1:dim(a00_SPECIES_SUMMARY)[1]){
  dup_accessions <- paste(a03_UNIQUEDB[a03_UNIQUEDB$species == a00_SPECIES_SUMMARY$search_name[b], "duplicate_accessions"], collapse="|")
  dup_accessions <- unlist(strsplit(dup_accessions, split="\\|"))
  dup_accessions_unique <- paste(unique(dup_accessions), collapse="|")
  
  dup_species <- paste(a03_UNIQUEDB[a03_UNIQUEDB$species == a00_SPECIES_SUMMARY$search_name[b], "duplicate_species"], collapse="|")
  dup_species <- unlist(strsplit(dup_species, split="\\|"))
  dup_species_unique <- unique(dup_species)
  
    if (length(dup_accessions)>0){
  a00_SPECIES_SUMMARY$dup_accessions[b]<-dup_accessions_unique
  a00_SPECIES_SUMMARY$dup_species[b]<- paste(dup_species_unique, collapse="|")
  a00_SPECIES_SUMMARY$dup_species_n[b]<- length(dup_species_unique)
  }
  
  #reset variables
  rm(dup_accessions)
  rm(dup_species)
}

a00_SPECIES_SUMMARY$n_all_seqs <- a00_SPECIES_SUMMARY$n_mitogenome + a00_SPECIES_SUMMARY$n_target
n_species_withnames <- dim(a00_SPECIES_SUMMARY)[1]
n_species_withnamesANDtargets <- length(which(a00_SPECIES_SUMMARY$n_all_seqs>0))
n_species_nonames <- dim(a01_NAMES_missing)[1]
write.csv(a00_SPECIES_SUMMARY, file.path(output_folder, "a00_SPECIES_SUMMARY.csv"), row.names=FALSE)

n_species_withnames; n_species_withnamesANDtargets; n_species_nonames
```