<p>The formal question I've been given is what aligner would I use? There is an issue between forming a data pipeline for <code>ape</code> and the latest and greatest. Its trade and the compromise would be <code>msaClustalOmega()</code>, but the rationale is complicated.</p>
<p>I strictly use <a href="https://drive5.com/muscle5/manual/cmd_super5.html" rel="nofollow noreferrer"><code>muscle5</code></a> or specifically <code>muscle -super5</code> option as a standalone. I am certain <code>R</code> <code>msa</code> has not introduced <code>muscle5</code> and I don't even know if it's published yet. It's on bioRxiv.org. The <code>super5</code> option must be specified and once done its amongst the most incredible algorithms <em>per se</em> 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 <code>muscle5</code>?</p>
<p>Its competitor was MAFFT and published benchmarks demonstrated it was the strongest, before <code>muscle5</code> was released.</p>
<p>The dilemma is simply that <code>ape</code> is cool and if you are wanting to have a seamless pipeline then greatest aligners don't fit with that. MAFFT isn't in <code>msa</code> and seriously doubt <code>muscle5</code> will be. HOWEVER Clustal Omega is in <code>msa</code> 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.</p>
<p>If you've got 50000 sequences - it's <code>muscle5</code> no contest. If you've 50 sequences then keep the pipeline together via <code>msaClustalOmega</code>.</p>
<hr />
<p>The original question and code appeared to involve converting fasta inside <code>R</code> onto <code>msa</code>. Thats how subprocessing works outside <code>R</code> (or <code>Python</code>) but it's not the OO approach.</p>
<p>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 <code>R</code>.</p>
<p>The error was the <code>msa</code> 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 &quot;official fasta&quot; 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 <code>seq</code> should load straight into msa with &quot;fasta&quot; option removed. Finally could you remove msa and used <code>msaMuscle</code> then simply <code>msaMuscle(seq)</code>.</p>

CC BY-SA 4.0
2023-02-20T20:31:13.907Z
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?
I Del Toro
CC BY-SA 4.0
2023-02-20T20:47:59.730Z
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.
M__
CC BY-SA 4.0
2023-02-20T17:11:34.430Z
2023-02-20T17:18:10.683Z
2023-02-20T17:18:10.683Z
M__
<p>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 <code>msa</code> package where I get an error <code>Error in convertAlnRows(result$msa, type) : There is an invalid aln file!</code> As always thanks in advance for taking a look: Here is the reproducible example:</p>
<pre><code>library (ape)
library (rentrez)
library (msa)

#Download Species Data 
# B_terrestris
B_terrestris &lt;- entrez_fetch(db = &quot;nucleotide&quot;, 
                                  id = &quot;NC_045179.1&quot;, 
                                  rettype = &quot;fasta&quot;)

# B_hypocrita
B_hypocrita &lt;- entrez_fetch(db = &quot;nucleotide&quot;, 
                             id = &quot;NC_011923.1&quot;, 
                             rettype = &quot;fasta&quot;)

#Vespa
Vespa &lt;- entrez_fetch(db = &quot;nucleotide&quot;, 
                      id = &quot; MT137096.1&quot;, 
                      rettype = &quot;fasta&quot;)

seq&lt;- c(B_terrestris, B_hypocrita, Vespa) # gotta clean it 

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


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

str (clean.seq)

msa (clean.seq,  type=&quot;dna&quot;, method=&quot;Muscle&quot;, cluster=&quot;neighborjoining&quot;)
</code></pre>

CC BY-SA 4.0
2023-02-20T16:38:58.750Z
@M__ Thank you! What sequence alignment method would you reccomend instead?
I Del Toro
<r><phylogenetics><phylogeny><fasta>
Phylogeny building in R from FASTA files:
# Phylogeny building in R from FASTA files:

(Asked by: username_0 on 2023-02-20)

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

**Tags:** r, phylogenetics, phylogeny, fasta
> **Comments:**
> **username_0**: @M__ Thank you! What sequence alignment method would you reccomend instead?  

---

## 1 answers

### Answer by username_1 on 2023-02-20. Score: 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`](https://drive5.com/muscle5/manual/cmd_super5.html) 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)`.

> **Comments:**
> **username_0**: 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?  
> **username_1**: 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.  

---

