text: user2425: I have a bam file aligned to a reference sequence. I am interested in finding the number of reads for read length 18 to 34, where the last base is mis-matched. I got this bam file from bowtie alignment allowing 3 mismatches (-v 0) option, so I am expecting some of the reads to have terminal bases mis-matched. How can we count the number of sequences of a specific length with terminal mismatch in bam file?
  **UPDATE:
  I have created a bowtie aligned .sam file with this command: 
  bowtie -S -p 18 /media/owner/ref.fa test1.fastq test1.sam
  which I then converted to .bam file like this: 
  samtools view -bS -o test1.bam test1.sam
  From this test1.bam, I was able to get the total number of plus and minus strand of reads aligned to the reference loci using bedtools coverage maping (for positives reads and negative reads) like this: 
  bedtools coverage -a my_region.bed -b test1.bam -bed -d -s | gzip > test1_region_pos.tsv.gz
  Since I used default bowtie option which allows for three mismatches, I want to know how many reads are aligned with terminal mismatches (i.e,last three, last two and last one base mismatched at 3' region of the aligned reads at their loci position). 
  
  user425: Could you provide some examples of lines that should be counted? I'm not sure how mismatches are encoded in the version of Bowtie you're using. As far as I remember, Bowtie doesn't distinguish between = and X in CIGARs (and uses M for both), therefore, mismatches must be detected from the MD tag. 
  You can use SAMsift for filtering the alignments. The following filtering criterion should keep only reads with  read length from 18 to 34 that have an MD tag with 0 as the last character, which must be preceded by a non-number (this should signalise a mismatch at the last position).
  samsift \
    -i test1.bam \
    -f '18 <= len(SEQ) <= 34 and MD[-1]=="0" and not MD[-2].isdigit()' \
    -m nonstop-remove
  
  The obtained alignments can be counted by adding | samtools view | wc -l.
  Explanation:
  
  -f ... – filtering criterion
  18 <= len(SEQ) <= 34 – the allowed read lengths
  MD[-1]=="0" – the last character of the MD tag is 0
  not MD[-2].isdigit() – the penultimate character of the MD tag is not a digit
  -m nonstop-remove – remove alignments without any MD tag (probably unaligned, they would cause an error).
  
  
original_text: 2425: <p>I have a bam file aligned to a reference sequence. I am interested in finding the number of reads for read length 18 to 34, where the last base is mis-matched. I got this bam file from bowtie alignment allowing 3 mismatches (-v 0) option, so I am expecting some of the reads to have terminal bases mis-matched. How can we count the number of sequences of a specific length with terminal mismatch in bam file?</p>
  
  <p>**UPDATE:</p>
  
  <p>I have created a bowtie aligned .sam file with this command: 
  <code>bowtie -S -p 18 /media/owner/ref.fa test1.fastq test1.sam</code></p>
  
  <p>which I then converted to .bam file like this: 
  <code>samtools view -bS -o test1.bam test1.sam</code></p>
  
  <p>From this test1.bam, I was able to get the total number of plus and minus strand of reads aligned to the reference loci using bedtools coverage maping (for positives reads and negative reads) like this: 
  <code>bedtools coverage -a my_region.bed -b test1.bam -bed -d -s | gzip &gt; test1_region_pos.tsv.gz</code></p>
  
  <p>Since I used default bowtie option which allows for three mismatches, I want to know how many reads are aligned with terminal mismatches (i.e,last three, last two and last one base mismatched at 3' region of the aligned reads at their loci position). </p>
   425: <p>Could you provide some examples of lines that should be counted? I'm not sure how mismatches are encoded in the version of Bowtie you're using. As far as I remember, Bowtie doesn't distinguish between <code>=</code> and <code>X</code> in CIGARs (and uses <code>M</code> for both), therefore, mismatches must be detected from the MD tag. </p>
  
  <p>You can use <a href="http://github.com/karel-brinda/samsift" rel="nofollow noreferrer">SAMsift</a> for filtering the alignments. The following filtering criterion should keep only reads with  read length from 18 to 34 that have an MD tag with <code>0</code> as the last character, which must be preceded by a non-number (this should signalise a mismatch at the last position).</p>
  
  <pre><code>samsift \
    -i test1.bam \
    -f '18 &lt;= len(SEQ) &lt;= 34 and MD[-1]=="0" and not MD[-2].isdigit()' \
    -m nonstop-remove
  </code></pre>
  
  <p>The obtained alignments can be counted by adding <code>| samtools view | wc -l</code>.</p>
  
  <p>Explanation:</p>
  
  <ul>
  <li><code>-f ...</code> – filtering criterion</li>
  <li><code>18 &lt;= len(SEQ) &lt;= 34</code> – the allowed read lengths</li>
  <li><code>MD[-1]=="0"</code> – the last character of the <code>MD</code> tag is <code>0</code></li>
  <li><code>not MD[-2].isdigit()</code> – the penultimate character of the <code>MD</code> tag is not a digit</li>
  <li><code>-m nonstop-remove</code> – remove alignments without any <code>MD</code> tag (probably unaligned, they would cause an error).</li>
  </ul>
  