Introducing coil: an R package for DNA barcode data cleaning

Full analysis pipeline for a single sequence

An example execution of the complete coi5p analysis pipeline with default options is demonstrated below using an example COI-5P barcode DNA sequence.

output = coi5p_pipe(example_nt_string)
output

## coi5p barcode sequence
## raw sequence:
## ctctacttgatttttggtgcatgag...ggacccaattctctatcaacactta
## framed sequence:
## ---ctctacttgatttttggtgcat...ggacccaattctctatcaacactta
## Amino acid sequence:
## -LYLIFGAWAG?VG?ALSLLIRAEL...LTDRNLNTTFFDPAGGGDPILYQHL
## Raw sequence was trimmed: FALSE
## Stop codon present: FALSE, Amino acid PHMM score:-206.22045
## The sequence likely does not contain an insertion or deletion.
## Base pair 1 of the raw sequence is base pair 4 of the COI-5P region.

Executing the entire pipeline yields a coi5p object. Calling the variable name prints the coi5p object’s summary and shows important information about the sequence.

Individual components can be obtained from the object using the dollar sign notation.

#see the available components
names(output)

##  [1] "name"         "raw"          "data"         "framed"      
##  [5] "was_trimmed"  "align_report" "aaSeq"        "aaScore"     
##  [9] "indel_likely" "stop_codons"

#retrieve only the amino acid sequence from the object
output$aaScore

## [1] -206.2204

By default the pipeline conducts censored translation, avoiding translation of the codons that are known to code for different amino acids within different species. If taxonomic information is available for the sample (available ranks: family, order, class, phylum), in most cases you can use the helper function which_trans_table to determine the proper genetic code to use. If the taxonomic group contains species that have different genetic codes, a 0 is returned to indicate that it is a good idea to stick to censored translation.

ex_table_to_use = which_trans_table("Scyliorhinidae")
ex_table_to_use

## [1] 2

The analysis can then be run with a non-censored translation step. Note below that the amino acid sequence is now devoid of question marks and the PHMM score is lower.

output = coi5p_pipe(example_nt_string, trans_table = ex_table_to_use)
output

## coi5p barcode sequence
## raw sequence:
## ctctacttgatttttggtgcatgag...ggacccaattctctatcaacactta
## framed sequence:
## ---ctctacttgatttttggtgcat...ggacccaattctctatcaacactta
## Amino acid sequence:
## -LYLIFGAWAGMVGMALSLLIRAEL...LTDRNLNTTFFDPAGGGDPILYQHL
## Raw sequence was trimmed: FALSE
## Stop codon present: FALSE, Amino acid PHMM score:-103.65363
## The sequence likely does not contain an insertion or deletion.
## Base pair 1 of the raw sequence is base pair 4 of the COI-5P region.

Calling functions individually

There are four functions that constitute the coi5p analysis pipeline: coi5p, frame, translate and indel_check. These are available to the user individually, for instances where only part of the analysis pipeline is needed (i.e. if you wish to frame sequences but not waste resources translating them, you could run only coi5p and frame).

  #build the coi5p object
  dat = coi5p(example_nt_string, name="example_sequence_1")
  #frame the sequence
  dat = frame(dat)
  #since we determined the genetic code above, we can use
  #the proper transaltion table as opposed to conducting 
  #the default censored translation
  dat = translate(dat, trans_table = 2)
  #check to see if an insertion or deletion is likely
  dat = indel_check(dat)
  dat

## coi5p barcode sequence: example_sequence_1
## raw sequence:
## ctctacttgatttttggtgcatgag...ggacccaattctctatcaacactta
## framed sequence:
## ---ctctacttgatttttggtgcat...ggacccaattctctatcaacactta
## Amino acid sequence:
## -LYLIFGAWAGMVGMALSLLIRAEL...LTDRNLNTTFFDPAGGGDPILYQHL
## Raw sequence was trimmed: FALSE
## Stop codon present: FALSE, Amino acid PHMM score:-103.65363
## The sequence likely does not contain an insertion or deletion.
## Base pair 1 of the raw sequence is base pair 4 of the COI-5P region.

Example of a batch analysis of barcode sequences

Here we will be working with the example dataframe: example_barcode_data. Although loading and outputting DNA sequence data in R is outside of the scope of the coil package, the supplementary section at the end of this vignette includes an example of how one can load a fasta file into a dataframe with a structure matching that of example_barcode_data.

example_barcode_data contains 9 barcode sequences that demonstrate the different abilities of the coil package. Some sequences are longer than the barcode COI-5P barcode region, some are shorter and some have insertion or deletion errors.

#this is the example data 
dim(example_barcode_data)

## [1] 9 5

names(example_barcode_data)

## [1] "id"           "genetic_code" "taxa"         "sequence"    
## [5] "notes"

# to look at the full dataframe:
# example_barcode_data

The coi5p analysis pipeline can be applied to the a dataframe in a iterative fashion. Here the pipeline is implemented through the use of the lapply function, which lets the unique sequence, id and genetic code of each row in the dataframe be passed into the coi5p_pipe function.

example_barcode_data$coi_output = lapply(1:length(example_barcode_data$id), function(i){
  coi5p_pipe(example_barcode_data$sequence[i], 
             name = example_barcode_data$id[i], 
             trans_table = example_barcode_data$genetic_code[i])
})

example_barcode_data$coi_output[[1]] #example of the first output

## coi5p barcode sequence: ex_1_clean
## raw sequence:
## acgctttactttatttttggcatgt...taaccctattctttaccagcatttg
## framed sequence:
## acgctttactttatttttggcatgt...taaccctattctttaccagcatttg
## Amino acid sequence:
## TLYFIFGMWAGFIGLSMSLLIRMEL...LTDRNFNTSFFDPSGGGNPILYQHL
## Raw sequence was trimmed: FALSE
## Stop codon present: FALSE, Amino acid PHMM score:-184.28122
## The sequence likely does not contain an insertion or deletion.
## Base pair 1 of the raw sequence is base pair 1 of the COI-5P region.

The coi5p objects are nested within the dataframe. Individual components can be extracted from the object as needed using the dollar sign notation. Below lapply is used to extract the framed sequence from each coi5p object and turn it into its own column in the dataframe. As we can see from the output below, dashes have been added to the front of the short sequence, if we compare the framed sequences from the long inputs (rows 5 and 6) to their original sequence, we see that coi5p has trimmed the sequence outside of the barcode region.

example_barcode_data$framed_seq = unlist(lapply(example_barcode_data$coi_output, function(x){
  x$framed
}))

#has coi5p trimmed characters?
nchar(example_barcode_data$framed_seq[[5]]) < nchar(example_barcode_data$sequence[[5]])

## [1] TRUE

The lapply notation used above is rather clunky, so the coi5p package contains a helper function to aid the user in flattening a list of coi5p objects into a dataframe. By default, all of the available object components will be output to the dataframe, but the user can choose a subset of components they wish to extract from the object. Note: this function assumes that the coi5p objects in the list have been put through the same workflow and therefore have a matching set of components. As an example of how you could break it, if you’ve applied the translate function to only one member of the list and not the others then the coi5p objects will have non-matching sets of components and flatten_coi5p will not work properly.

#extract only a single column
col_df = flatten_coi5p(example_barcode_data$coi_output, keep_cols = 'aaSeq')
#extract multiple columns
multi_df = flatten_coi5p(example_barcode_data$coi_output, keep_cols = c('framed','aaSeq'))
#extract all columns
full_coi5p_df = flatten_coi5p(example_barcode_data$coi_output)
#full_coi5p_df

The memory requirements of the method demonstrated above are trivial because the example dataframe has only nine rows. If millions of sequences are being processed, then keeping all of the coi5p objects in memory at once may become prohibitive. This will depend on the amount of RAM available on your machine. The average coi5p_pipe output is ~6KB in size, so processing 1 million sequences at once would occupy ~6GB of RAM. If you are trying to limit RAM usage, the following workflow can help keep RAM requirements modest by instantiating only one coi5p object at a time, but as a tradeoff may take slightly more time to execute.

full_coi5p_df = data.frame(matrix(ncol = 7, nrow = 0),stringsAsFactors = FALSE )
colnames(full_coi5p_df) = c("name", "raw", "framed", 
                            "aaSeq", "aaScore", "indel_likely", "stop_codons")

for(i in 1:length(example_barcode_data$id)){
    out_data = coi5p_pipe(example_barcode_data$sequence[i], 
                            name = example_barcode_data$id[i], 
                            trans_table = example_barcode_data$genetic_code[i])
  #for extreme memory conservation - could write each line of output to a .csv
    #instead of binding it to an output dataframe.
    full_coi5p_df = rbind(full_coi5p_df, flatten_coi5p(list(out_data)))
}

As demonstrated here, the coil package allows for robust cleaning, contextualization and error assessment of novel COI-5P barcode data. The package’s analysis pipeline is designed in a modular fashion, allowing the user to run only the functions required for their given use case. The pipeline is designed with scalability in mind; each sequence is processed individually, allowing for parallelization to optimize analysis speed (i.e. via R’s mclapply function) when computational resources are abundant or for the sequential analysis of sequences when limited memory is available.

Loading and manipulating a fasta file

The example of processing batch DNA barcode data above begins with the data in a clean dataframe. Since barcode data is not always obtained in a tidy format, some preprocessing by the user will likely be required. The following is provided to aid the user in developing a workflow for loading their barcode sequence data into R and constructing a tidy dataframe before beginning sequence analysis.

The example presented below shows how one can load a fasta file containing DNA sequences into R and then convert the sequence and header data into a tidy dataframe object. More information on the read.fasta function can be found in the seqinr documentation.

Information found in the header lines of fasta files varies, so the reader will likely need to alter this code for novel data sources. In this example, the header contains four fields (id, genetic code, taxa and notes) separated by a pipe character (|). The code below turns this fasta file into a dataframe that matches the example_barcode_data file used above.

library(seqinr)

## 
## Attaching package: 'seqinr'

## The following object is masked from 'package:coil':
## 
##     translate

# load the example fasta file included with coil
# included in the file's header line:
# the name of the sample, its genetic code, taxonomic designation and some notes
ex_fasta_file = system.file("extdata/example_barcode_data.fasta", package = "coil")

#read in the example fasta file using seqinr
ex_data = seqinr::read.fasta(ex_fasta_file, as.string = TRUE)

#here is what the output from read.fasta looks like
#head(ex_data)

#parse the data in the header line by splitting the name on the | character
parsed_names_data = lapply(1:length(ex_data), function(i){
  unlist(strsplit(names(ex_data)[[i]],"\\|"))
})

# subset the components of the header line and build these and the sequence 
# into a dataframe matching the style used in the coi5p batch example
example_barcode_data_from_scratch = data.frame(
  id = sapply(parsed_names_data, function(x) x[[1]]),
  genetic_code = sapply(parsed_names_data, function(x) x[[2]]),
  taxa = sapply(parsed_names_data, function(x) x[[3]]),
  sequence = unname(unlist(ex_data)),
  notes = sapply(parsed_names_data, function(x) x[[4]])
)

#head(example_barcode_data_from_scratch)