Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-10-10T12:25:43.910Z
Comments.Text: 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..
Comments.UserDisplayName: Niklas
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-10-10T13:13:25.733Z
Comments.Text: @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)
Comments.UserDisplayName: Fravadona
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-10-10T17:01:42.900Z
Comments.Text: 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.
Comments.UserDisplayName: Niklas
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-10-10T17:42:55.643Z
Comments.Text: Maybe something like pileup is what you are looking for? https://samtools.github.io/bcftools/howtos/variant-calling.html
Comments.UserDisplayName: Niklas
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-10-10T20:34:18.950Z
Comments.Text: @Niklas Thank you, `bfctool` will definitely be useful
Comments.UserDisplayName: Fravadona
Tags: <python><phylogenetics><snp><multiple-sequence-alignment>
Title: SNP How to detect regions of erratic pairwise alignment?
ThreadText: # 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:
  
    ```none
    ref  GTATGCGGATACCG
  
    sq1  GTATCGTGACCG
    sq2  GTATCTCGTGACCG
    sq3  GTGTCTATTCGTGACCG
    sq4  GTATCTCGTGACCG
  
    ```
  * Independent pairwise alignments against the reference (I removed the insertions from the output):
  
    ```none
    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:
  
    ```none
    ref  GTATG---CGGATACCG
  
    sq1  GTAT-----CGTGACCG
    sq2  GTATC---TCGTGACCG
    sq3  GTGTCTATTCGTGACCG
    sq4  GTATC---TCGTGACCG
  
    ```
  * Same MSA, but trimmed:
  
    ```none
    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**: @username_1 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**: @username_1 Thank you, `bfctool` will definitely be useful  
  
  ---
  
  ## 0 answers
  
  