<p>I have clusters of DNA sequences that contain some mutations. All the sequences in a cluster should contain the same mutations.</p>
<p>As the mutations aren't known, those clusters of sequences were pairwise-aligned against the non-mutated reference with <code>Minimap2</code>. Then an interesting fact cropped-up: a few mutations close to each other will make the alignment go awry in that zone.</p>
<ul>
<li><p>Zone with a few SNPs:</p>
<pre class="lang-none prettyprint-override"><code>ref  GTATGCGGATACCG

sq1  GTATCGTGACCG
sq2  GTATCTCGTGACCG
sq3  GTGTCTATTCGTGACCG
sq4  GTATCTCGTGACCG
</code></pre>
</li>
<li><p>Independent pairwise alignments against the reference (I removed the insertions from the output):</p>
<pre class="lang-none prettyprint-override"><code>ref  GTATGCGGATACCG

sq1  GTAT-CGG--ACCG
sq2  GTATCTCGTGACCG
sq3  GTATTCGG--ACCG
sq4  GTATCTCGTGACCG
</code></pre>
</li>
</ul>
<p>So I made a MSA of this small cluster (adding the reference in it) with <code>Muscle5</code>; the difference is noticeable, as we can clearly &quot;see&quot; the mutations now:</p>
<ul>
<li><p>MSA alignment that includes the reference:</p>
<pre class="lang-none prettyprint-override"><code>ref  GTATG---CGGATACCG

sq1  GTAT-----CGTGACCG
sq2  GTATC---TCGTGACCG
sq3  GTGTCTATTCGTGACCG
sq4  GTATC---TCGTGACCG
</code></pre>
</li>
<li><p>Same MSA, but trimmed:</p>
<pre class="lang-none prettyprint-override"><code>ref  GTATGCGGATACCG

sq1  GTAT--CGTGACCG
sq2  GTATCTCGTGACCG
sq3  GTGTCTCGTGACCG
sq4  GTATCTCGTGACCG
</code></pre>
</li>
</ul>
<p>I would say calculating the MSA seems to be the way to go but the clusters sizes range from 2 to 1000 sequences, the sequences are 3000nt long, and there's 2,000,000 sequences in total; we cannot afford to do the MSA of each cluster on the full size of the sequences.
<br>
On the other hand, we could probably manage to do a MSA of a few small zones in the sequences.</p>
<p>So, given a cluster and the pairwise alignments of its sequences against the reference, <strong>how would you determine the zones that need to be re-aligned?</strong></p>
<p><strong>note:</strong> I'm using Muscle5 because there's a Python binding of it and the results seem good. I also tried FAMSA, which is very fast, but it didn't yield good alignments.</p>

CC BY-SA 4.0
2023-10-10T12:25:43.910Z
Are sequences 1-4 in your 3 example supposed to be the same sequence? It is not really clear to me what you are trying to do / which problem you are trying to solve. You seem to be dealing mainly with indels rather than SNPs? Its also not clear why the alignment is awry, it is probably the best alignment with the given penalty matrix..
Niklas
CC BY-SA 4.0
2023-10-10T13:13:25.733Z
@Niklas Thank you for the comment. Yes, they are the same sequences (I updated the examples to make it clearer). By "awry" I mean that I cannot  make a consensus using the pairwise alignments, as they are subjectively "wrong" (mostly because of the indels)
Fravadona
CC BY-SA 4.0
2023-10-10T17:01:42.900Z
OK. Generally, all alignments are equally "correct" or "wrong", since they minimize penalty score given the penalty matrix, which you can edit as you prefer. You would expect a msa to yield different results from pairwise alignments since it is a different problem.. regarding your actual question: I assume you could identify "bad" regions from a long stretch by the presence of many gaps in that region. But, keep in mind that when you re-align short sequences you may get different results than for the long sequences.
Niklas
CC BY-SA 4.0
2023-10-10T17:42:55.643Z
Maybe something like pileup is what you are looking for? https://samtools.github.io/bcftools/howtos/variant-calling.html
Niklas
CC BY-SA 4.0
2023-10-10T20:34:18.950Z
@Niklas Thank you, `bfctool` will definitely be useful
Fravadona
<python><phylogenetics><snp><multiple-sequence-alignment>
SNP How to detect regions of erratic pairwise alignment?
# SNP How to detect regions of erratic pairwise alignment?

(Asked by: username_0 on 2023-10-10)

I have clusters of DNA sequences that contain some mutations. All the sequences in a cluster should contain the same mutations.
As the mutations aren't known, those clusters of sequences were pairwise-aligned against the non-mutated reference with Minimap2. Then an interesting fact cropped-up: a few mutations close to each other will make the alignment go awry in that zone.

Zone with a few SNPs:
ref  GTATGCGGATACCG

sq1  GTATCGTGACCG
sq2  GTATCTCGTGACCG
sq3  GTGTCTATTCGTGACCG
sq4  GTATCTCGTGACCG


Independent pairwise alignments against the reference (I removed the insertions from the output):
ref  GTATGCGGATACCG

sq1  GTAT-CGG--ACCG
sq2  GTATCTCGTGACCG
sq3  GTATTCGG--ACCG
sq4  GTATCTCGTGACCG



So I made a MSA of this small cluster (adding the reference in it) with Muscle5; the difference is noticeable, as we can clearly "see" the mutations now:

MSA alignment that includes the reference:
ref  GTATG---CGGATACCG

sq1  GTAT-----CGTGACCG
sq2  GTATC---TCGTGACCG
sq3  GTGTCTATTCGTGACCG
sq4  GTATC---TCGTGACCG


Same MSA, but trimmed:
ref  GTATGCGGATACCG

sq1  GTAT--CGTGACCG
sq2  GTATCTCGTGACCG
sq3  GTGTCTCGTGACCG
sq4  GTATCTCGTGACCG



I would say calculating the MSA seems to be the way to go but the clusters sizes range from 2 to 1000 sequences, the sequences are 3000nt long, and there's 2,000,000 sequences in total; we cannot afford to do the MSA of each cluster on the full size of the sequences.

On the other hand, we could probably manage to do a MSA of a few small zones in the sequences.
So, given a cluster and the pairwise alignments of its sequences against the reference, how would you determine the zones that need to be re-aligned?
note: I'm using Muscle5 because there's a Python binding of it and the results seem good. I also tried FAMSA, which is very fast, but it didn't yield good alignments.

**Tags:** python, phylogenetics, snp, multiple-sequence-alignment
> **Comments:**
> **username_1**: Are sequences 1-4 in your 3 example supposed to be the same sequence? It is not really clear to me what you are trying to do / which problem you are trying to solve. You seem to be dealing mainly with indels rather than SNPs? Its also not clear why the alignment is awry, it is probably the best alignment with the given penalty matrix..  
> **username_0**: @Niklas Thank you for the comment. Yes, they are the same sequences (I updated the examples to make it clearer). By "awry" I mean that I cannot  make a consensus using the pairwise alignments, as they are subjectively "wrong" (mostly because of the indels)  
> **username_1**: OK. Generally, all alignments are equally "correct" or "wrong", since they minimize penalty score given the penalty matrix, which you can edit as you prefer. You would expect a msa to yield different results from pairwise alignments since it is a different problem.. regarding your actual question: I assume you could identify "bad" regions from a long stretch by the presence of many gaps in that region. But, keep in mind that when you re-align short sequences you may get different results than for the long sequences.  
> **username_1**: Maybe something like pileup is what you are looking for? https://samtools.github.io/bcftools/howtos/variant-calling.html  
> **username_0**: @Niklas Thank you, `bfctool` will definitely be useful  

---

## 0 answers

