1
$\begingroup$

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")
$\endgroup$
1
  • $\begingroup$ @M__ Thank you! What sequence alignment method would you reccomend instead? $\endgroup$
    – I Del Toro
    Feb 20, 2023 at 16:38

1 Answer 1

1
$\begingroup$

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).

$\endgroup$
2
  • 1
    $\begingroup$ I tied this bit with the set of sequences ` msaClustalW(c(B_terrestris, Vespa, B_hypocrita), type="dna", cluster="nj")` but now I get this error: "Error in .Call2("new_XStringSet_from_CHARACTER", class(x0), elementType(x0), : key 69 (char 'E') not in lookup table"- Any further suggestions? $\endgroup$
    – I Del Toro
    Feb 20, 2023 at 20:31
  • $\begingroup$ At least it can see the object. If this was Python I could debug it (well I could write the code to do it straight up), R will be more challenging. You need to print the seq object and check its integrity. It seems to be complaining that E isn't DNA, suggesting either this is protein (unlikely) or an E from the header has got caught in the sequence. There are alternative ways to construct the sequence object, so as a last resort could do that. The idea here is to split the object attributes into sequence and ids then reassemble it. $\endgroup$
    – M__
    Feb 20, 2023 at 20:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.