8  Working with Sequences in R

8.1 ape

Introducing my favourite R package in this topic: ape. I simply use read.FASTA() and write.FASTA() for sequence I/O in R.

The sequence in ape is stored as a DNAbin class. There are a few other ways to create a DNAbin object, one of the most important being as.DNAbin(). But it takes sequences represented in vectors, meaning:

library(ape)
seq.char <- "AGTCGATAG"
seq <- as.DNAbin(seq.char) # This will fail, because `seq.char` is character of length 1. Sequence is not represented by vector, but by one single character string. 

seq 
# No error was shown but there are only NAs in the object. 

# We have to split character string into vector first.
seq.vec <- strsplit(seq.char, "")
seq <- as.DNAbin(seq.vec)

seq
# Now the actual sequence is recognized.

as.DNAbin can also take multiple seuqences at once in form of a list of vectors. Names of each list elements will become the sequence name, so as the sequence name in fasta file (those characters coming after >). You can extract or edit sequence names using names().

With these minimal functions from ape and data wranglings, you can do a lot of descent automations on managing sequence. For example, cleaning the sequence names exported by Geneious (all the “extraction”, “reversed” annotations):

library(tidyverse)
contig.actcp <- ape::read.FASTA("data/sequence/contig/COI_contigs.fasta") 

# Clean up sequence name
names(contig.actcp) <- names(contig.actcp) %>% # modify sequence names
    str_extract("^[0-9]{3}") %>% # extract leading 3-digit voucher using regex
    paste0("JN_Mu_", .) # add prefix "JN_Mu_" to the voucher
    
ape::write.FASTA(contig.actcp, "data/sequence/contig/COI_contigs_cleaned.fasta")

8.2 Multiple Sequence Alignment (MSA)

Using my zseq package, you can also call MAFFT from within R. The MAFFT needs to be installed locally. See ?zseq::mafft() for syntax.

contig.combined.path <- "data/sequence/contig/COI_combined.fasta"
aln.path <- "data/sequence/alignment/COI_combined.fasta"

zseq::mafft(contig.combined.path, aln.path)