Phylogeny building in R from FASTA files:

Question

Im working on building a phylogeny from scratch with downloaded FASTA sequences from GeneBank. I think Im doing alright up until the multi sequence alignment in the msa package where I get an error Error in convertAlnRows(result$msa, type) : There is an invalid aln file! As always thanks in advance for taking a look: Here is the reproducible example:

library (ape)
library (rentrez)
library (msa)

#Download Species Data 
# B_terrestris
B_terrestris <- entrez_fetch(db = "nucleotide", 
                                  id = "NC_045179.1", 
                                  rettype = "fasta")

# B_hypocrita
B_hypocrita <- entrez_fetch(db = "nucleotide", 
                             id = "NC_011923.1", 
                             rettype = "fasta")

#Vespa
Vespa <- entrez_fetch(db = "nucleotide", 
                      id = " MT137096.1", 
                      rettype = "fasta")

seq<- c(B_terrestris, B_hypocrita, Vespa) # gotta clean it 

#FASTA Cleaner Function 
fasta_cleaner <- function(fasta_object, parse = TRUE)
{
  fasta_object <- sub("^(>)(.*?)(\\n)(.*)(\\n\\n)","\\4", fasta_object)
  fasta_object <- gsub("\n", "", fasta_object)
  
  if(parse == TRUE){
    fasta_object <- stringr::str_split(fasta_object,
                                       pattern = "",
                                       simplify = FALSE)
  }
  
  return(fasta_object[[1]])
}


#This should be ready to go! 
clean.seq<-fasta_cleaner (seq, parse = FALSE)

str (clean.seq)

msa (clean.seq,  type="dna", method="Muscle", cluster="neighborjoining")

Answer 1

The formal question I've been given is what aligner would I use? There is an issue between forming a data pipeline for ape and the latest and greatest. Its trade and the compromise would be msaClustalOmega(), but the rationale is complicated.

I strictly use muscle5 or specifically muscle -super5 option as a standalone. I am certain R msa has not introduced muscle5 and I don't even know if it's published yet. It's on bioRxiv.org. The super5 option must be specified and once done its amongst the most incredible algorithms per se notably in speed - no other aligner comes close - and personal benchmarks in accuracy. To give some idea, I was going to write an aligner as part of a package (its makes stuff simpler) and thought - whats the point with muscle5?

Its competitor was MAFFT and published benchmarks demonstrated it was the strongest, before muscle5 was released.

The dilemma is simply that ape is cool and if you are wanting to have a seamless pipeline then greatest aligners don't fit with that. MAFFT isn't in msa and seriously doubt muscle5 will be. HOWEVER Clustal Omega is in msa and that is the strongest (most recent) of the options it has. Omega does have published benchmarks and its neither quite as fast nor quite as accurate as MAFFT.

If you've got 50000 sequences - it's muscle5 no contest. If you've 50 sequences then keep the pipeline together via msaClustalOmega.

---

The original question and code appeared to involve converting fasta inside R onto msa. Thats how subprocessing works outside R (or Python) but it's not the OO approach.

My original answer was in comments and remains an unofficial answer because usually I run the code I'm recommending and I don't fancy subprocessing muscle via R.

The error was the msa algorithm was attempting to shift the object into a format that will subprocess via muscle and it doesn't understand fasta format (it wants the raw object). Thus, you appear to be wanting to pass "official fasta" into msa. It doesn't need it, it simply needs a sequence object. It's imported as fasta but then it should be a sequence object (it's no longer fasta). That should go straight into msa. Thus seq should load straight into msa with "fasta" option removed. Finally could you remove msa and used msaMuscle then simply msaMuscle(seq).