2.1: New Dependencies Provide Up To Date Access, New Features, and Stability

This new version of AnnotationBustR is a complete rewrite of AnnotationBustR. Notably, AnnotationBustR 2.0+ now uses various Bioconductor dependencies to extract subsequences. Previously, AnnotationBustR accessed subsequences through the ACNUC via the package seqinr. There are several reasons for this change:

1. The GenBank version accessed by seqinr is quite out of date. We can see that the version of GenBank being used by ACNUC in seqinr is from 2021. Thus, newer accessions were not accessible in AnnotationBustR. We can see this by running the following:

ACNUC.GB.INFO <- seqinr::choosebank("genbank", infobank=T)
ACNUC.GB.INFO#return info on date

2. Several users reported connectivity issues, specifically with connecting to ACNUC via seqinr getting server connection issues or having their connection dropped with the ACNUC server. We now retrieve accessions from NCBI directly through the rentrez package, which is more stable.

3. ACNUC via seqinr did not always provide access to certain refseq sequences. These are now accessible in the current version written with rentrez and various Bioconductor packages.

4. ACNUC had its own predefined tags for feature subsequences. The new version of AnnotationBustR now parses GenBank feature files directly, allowing for a more flexible capture of potential feature types.

2.2: Other Minor New Features

1. Users no longer have to provide a translation code when requesting coding sequences be translated.

2. Users can now request coding sequences to be returned as DNA or the translated sequence or both.

3.1: Installation From CRAN

In order to install the stable CRAN version of the AnnotationBustR package:

install.packages("AnnotationBustR")

3.2: Installation of Development Version From GitHub

While we recommend use of the stable CRAN version of this package, you can use the R package devtools to install the development version of the package from GitHub if for any reason you wish to use it :

#1. Install 'devtools' if you do not already have it installed:
install.packages("devtools")

#2. Load the 'devtools' package and temporarily install the development version of
#'AnnotationBustR' from GitHub:
library(devtools)
dev_mode(on=T)
install_github("sborstein/AnnotationBustR")  # install the package from GitHub
library(AnnotationBustR)# load the package

#3. Leave developers mode after using the development version of 'AnnotationBustR' so it will not remain on 
#your system permanently.
dev_mode(on=F)

library(AnnotationBustR)

#install (if necessary) and load reutils
#install.packages("rentrez")#install if necessary
library(rentrez)

#search for Barbonymus mitogenome sequences
demo.search <- rentrez::entrez_search(db = "nuccore",  term = "Barbonymus[orgn] AND complete genome[title]", use_history = TRUE)#search
accessions <- rentrez::entrez_fetch(db = "nuccore", web_history = demo.search$web_history, rettype = "acc")#fetch accessions
accessions <- strsplit(accessions, "\n")[[1]]#split out accessions from new line characters

4.1:(Optional Step) Finding the Longest Available

AnnotationBustR’s FindLongestSeq function finds the longest available sequence for each species in a given set of GenBank accession numbers. All the user needs is to obtain a list of GenBank accession numbers they would like to input. The only function argument for FindLongestSeq is Accessions, which takes a vector of accession numbers as input. We can run the function below on the six Barbonymus accessions we found above by:

#Create a vector of GenBank nucleotide accession numbers. In order this contains accessions 
my.longest.seqs<-FindLongestSeq(accessions)#Run the FindLongestSeq function
my.longest.seqs#returns the longest seqs found per species

In this case we can see that the function worked as it only returned accession PP937076 (16577 bp) for Barbonymus schwanefeldii which was longer than accessions AP011317.1, KU498040.2, KU233186.2, KJ573467.1 (16577, 16570,16478, and 16576 bp respectively) as well as the single accessions for Barbonymus altus and Barbonymus gonionnotus. The table returns a three-column data frame with the species name, the corresponding accession number, and the length.

4.2: Load a Data Frame of Search Terms of Gene Synonyms to Search With:

AnnotationBustR works by searching through the annotation features table for a feature of interest using search terms for it (i.e. possible synonyms it may be listed under). These search terms are formatted to have three columns:

• Feature: The name you wish to refer to for a feature. This also controls the names of the FASTA files written for when that feature is extracted. It is important that you use names that will not confuse R, so don’t start these with numbers or include other characters like “.” or “-” that R uses for math.
• Type: The type of sequence to search for. For example, CDS, tRNA, rRNA, misc_RNA, D-Loop, or misc_feature to name a few. Features are parsed in AnnotationBustR dynamically, so any feature type should be extractable.
• Name: A possible synonym that the feature could be listed under that you wish to search for.

For extracting introns and exons, an additional fourth column is needed (which will be discussed in more detail later in the tutorial): -IntronExonNumber: The number of the intron or exon to extract (e.g. intron 2).

Below (Figure 2) is an example of where these corresponding items would be in the GenBank features table:

So, if we wanted to use AnnotationBustR to capture these and write them to a FASTA, we could set up a data frame that looks like the following.

##   Feature Type    Name
## 1    ATP8  CDS ATPase8
## 2    ATP6  CDS ATPase6

While AnnotationBustR will work with any data frame formatted as discussed above, we have included in it pre-made search terms for animal and plant mitochondrial DNA (mtDNA), chloroplast DNA (cpDNA), and ribosomal DNA (rDNA). These can be loaded from AnnotationBustR using:

#Load in pre-made data frames of search terms
data(mtDNAterms)#loads the mitochondrial DNA search terms for metazoans
data(mtDNAtermsPlants)#loads the mitochondrial DNA search terms for plants
data(cpDNAterms)#loads the chloroplast DNA search terms
data(rDNAterms)#loads the ribosomal DNA search terms

These data frames can also easily be manipulated to select only the loci of interest. For instance, if we were only interested in tRNAs from mitochondrial genomes, we could easily subset out the tRNAs from the premade mtDNAterms object by:

data(mtDNAterms)#load the data frame of mitochondrial DNA terms
tRNA.terms <- mtDNAterms[mtDNAterms$Type=="tRNA",]#subset out the tRNAs into a new data frame

4.3:(Optional Step) Merge Search Terms If Neccessary

While we have tried to cover as many synonyms for genes in our pre-made data frames, it is likely that some synonyms may not be represented due to the vast array of synonyms a single gene may be listed under in the features table. To solve this we have included the function MergeSearchTerms.

For example, let’s imagine that we found a completely new annotation for the gene cytochrome oxidase subunit 1 (COI) listed as CX1. The MergeSearchTerms function only has two arguments, ..., which takes two or more objects of class data.frame and the logical Sort.Genes, which We could easily add this to other mitochondrial gene terms by:

#Add imaginary gene synonym for cytochrome oxidase subunit 1, CX1
add.name <- data.frame("COI","CDS", "CX1", stringsAsFactors = FALSE)
colnames(add.name) <- colnames(mtDNAterms)#make the columnames for this synonym those needed for AnnotationBustR

#Run the merge search term function without sorting based on gene name.
new.terms <- MergeSearchTerms(add.name, mtDNAterms, SortGenes=FALSE)

#Run the merge search term function with sorting based on gene name.
new.terms <- MergeSearchTerms(add.name, mtDNAterms, SortGenes=TRUE)

We will use this function again in a more realistic example later in this vignette.

4.4 Extract sequences with AnnotationBust

The main function of AnnotationBustR is AnnotationBust. This function extracts the sub-sequence(s) of interest from the accessions and writes them to FASTA files in the current working directory. In addition to writing sub-sequences to FASTA files, AnnotationBust also generates an accession table for all found sub-sequences written to FASTA files which can then be written to a csv file using base R write.csv. AnnotationBustR requires at least two arguments, a vector of accessions for Accessions and a data frame of search terms formatted as discussed in 3.2 and 3.3 for Terms.

Additional arguments include the ability to specify duplicate genes you wish to recover as a vector of gene names using the Duplicates argument and specifying the number of duplicate instances to extract using a numeric vector (which must be the same length as Duplicates) using the DuplicateInstances argument.

AnnotationBustR also has arguments to translate coding sequences into the corresponding peptide sequence setting the TranslateSeqs argument to TRUE. If TranslateSeqs=TRUE, users should also specify the GenBank translation code number corresponding to their sequences using the TranslateCode argument. A list of GenBank translaton codes for taxa is available here: https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi

An additional argument of AnnotationBust is DuplicateSpecies which when set to DuplicateSpecies=TRUE adds the accession number to the species name for the FASTA file. This can be useful for later analyses involving FASTA files as duplicate names in FASTA headers can pose problems with some programs. It is important to note that if users select DuplicateSpecies=TRUE that while FASTA file will contain species names with their respected accession number, the corresponding accession table will either have a single row per species containing all the accession numbers for each subsequence locus found for that species seperated by a comma or a seperate row for each accession number dependent on what the user chooses for the TidyAccessions argument (see description below).

The argument Prefixcan either be NULL or a character vector of length 1. If Prefix is not NULL it will add the prefix specified to all FASTA files. If left as Prefix=NULL files names will just be the name of the locus.

The final argument TidyAccessions effects the format of the final accession table. If TidyAccessions=TRUE, all sequence for a single species will be collapsed into a single cell per locus and accessions will be seperated by commas. If TidyAccessions=FALSE it will leave each accession number in its own row in the accession table.

For the tutorial, we will use the accessions we created in examples 3.1 in the object my.longest.accessions. This is a vector that contains three accessions for mitogenomes of three species of Barbonymus. These accessions provide a good highlight as to the utility of this package.If we look at the annotations of these three accessions(https://www.ncbi.nlm.nih.gov/nuccore/AP011317.1,https://www.ncbi.nlm.nih.gov/nuccore/AP011181.1,https://www.ncbi.nlm.nih.gov/nuccore/AB238966.1), we can see that while AP011317 and AP011181 are annotated with similar nomenclature, AB238966 differs substantially (e.x. COI vs CO2, Cyt b vs Cb, etc.). For this example, we will use all the arguments of AnnotationBust to extract all 38 subsequences for the four accessions (22 tRNAs, 13 CDS, 2 rRNAs, and 1 D-loop). For this we will have to specify duplicates, in this case for tRNA-Leu and tRNA-Ser, which occur twice in vertebrate mitogenomes. We will translate the CDS using TranslateSeqs = "Both" to return both the DNA and peptide sequence. Because we have a single accession for each species, we will specify DuplicateSpecies=FALSE and have out FASTA files output with the prefix “Demo” by specifying Prefix="Demo" and with TiddyAccessions=TRUE. This will create 38 FASTA files for all mitochondrial loci in the mitochondrial genome in the working director with the prefix “Demo”.

#run AnnotationBust function for two duplicate tRNA genes occuring twice and translate CDS
my.seqs <- AnnotationBust(Accessions=my.longest.seqs$Accession, Terms=mtDNAterms,Duplicates=c("tRNA_Leu","tRNA_Ser"), DuplicateInstances=c(2,2), TranslateSeqs="Both", DuplicateSpecies=TRUE, Prefix="Demo", TidyAccessions=TRUE)

#We can return the accession table and write it to a CSV file.
my.seqs#retutn the accession table
write.csv(my.seqs, file="AccessionTable.csv")#Write the accession table to a csv file

We can also use AnnotationBustR to extract introns from sequences. One limitation to this is that the introns must be annotated within the subsequence, which is lacking from some accessions (e.x. not all chloroplasts have introns annotated). As stated above, this requires the use of a fourth column in the search terms, labelled IntronExonNumber, which is used to find which number intron is to be extracted. For instance, lets imagine we want to extract the coding sequence matK and exon 2 and intron 1 of trnK from the chloroplast genomes KX687911.1 and KX687910.1 We could set up the search terms as follows using the current cpDNAterms:

#Subset out matK from cpDNAterms
cds.terms <- subset(cpDNAterms,cpDNAterms$Feature=="matK")
#Create a vecotr of NA so we can merge with the search terms for introns and exons
cds.terms <- cbind(cds.terms,(rep(NA,length(cds.terms$Feature))))
colnames(cds.terms)[4] <- "IntronExonNumber"

#Prepare a search term table for the intron and exons to remove
#We can start with the cpDNAterms for trnK
IntronExon.terms <- subset(cpDNAterms,cpDNAterms$Feature=="trnK")
#As we want to go for two exons, we will want the synonyms repeated as we are doing and intron and an exon
IntronExon.terms <- rbind(IntronExon.terms,IntronExon.terms)#duplicate the terms
IntronExon.terms$Type <- rep(c("intron","intron","exon","exon"))#rep the sequence type we want to extract
IntronExon.terms$Feature <- rep(c("trnK_Intron","trnK_Exon2"),each=2)
IntronExon.terms <- cbind(IntronExon.terms,rep(c(1,1,2,2)))#Add intron/exon number info
colnames(IntronExon.terms)[4] <- "IntronExonNumber"#change column name for number info for IntronExon name

#We can then merge everything together with MergeSearchTerms terms
IntronExonExampleTerms <- MergeSearchTerms(IntronExon.terms,cds.terms)

#Run AnnotationBust
IntronExon.example <- AnnotationBust(Accessions=c("KX687911.1", "KX687910.1"), Terms=IntronExonExampleTerms, Prefix="DemoIntronExon")

4.5 Troubleshooting

As AnnotationBustR is dependent on access to online databases, it is important that you make sure you have a reliable internet connection while using the package, especially for the FindLongestSeq and AnnotationBust functions. Additionally, as these functions require access to NCBI, if you have internet access yet are still having issues, check that these sites and there servers are not down or undergoing maintenance, which they do occasionally. For AnnotationBust, you will get the followingCould not connect to NCBI message printed to the screen in the event that it cannot connect to the NCBI database.

The AnnotationBust function may also provide the following warnings and messages as well as report the following in the accessions table:

Working On Accession 1 of 1: ACCESSION# not found on NCBI. This indicates that the accession number was not found on NCBI. You should check that it is correct and exists. This will also fill in the accession table with Failed to find ACCESSION#. Here, ACCESSION# will be the actual accession.