Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-27T02:00:13.070Z
Answers.Text: Thanks for your answer. I have updated the question with how I created BAM file.
Answers.UserDisplayName: MAPK
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-27T02:01:07.917Z
Answers.Text: Could you explain this a bit please `MD[-1]=="0" and not MD[-2].isdigit()`?
Answers.UserDisplayName: MAPK
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-27T02:02:08.703Z
Answers.Text: Does this also consider for `+` and `-` strands?
Answers.UserDisplayName: MAPK
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-27T02:04:12.823Z
Answers.Text: Please, could you show at least one line that should pass the filter? I need it to verify the `MD` filtering step. Yes, it's possible to filter using `+`/`-` via `FLAG`.
Answers.UserDisplayName: Karel Břinda
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-27T02:09:06.320Z
Answers.Text: The bowtie -v 1 option allows one mismatch, -v 2, alows 2  mismatches anywhere in the sequences. The default option allows maximum of 3 mismatches.
Answers.UserDisplayName: MAPK
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-27T02:09:46.727Z
Answers.Text: I mean 1 line from the SAM file.
Answers.UserDisplayName: Karel Břinda
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-27T02:14:08.733Z
Answers.Text: Did you mean: `samtools view -H R1_001.sam       
  @HD VN:1.0 SO:unsorted
  @SQ SN:NC_013116 LN:2166
  @PG ID:Bowtie VN:1.1.2 CL:"bowtie --wrapper basic-0 -S -p 18 /m`
Answers.UserDisplayName: MAPK
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-27T02:48:58.020Z
Answers.Text: Why am I getting permission denied while using samsift? I did chmod u+x.
Answers.UserDisplayName: MAPK
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-27T02:55:24.983Z
Answers.Text: `/usr/bin/python: can't find '__main__' module in 'samsift'` Getting this error now?
Answers.UserDisplayName: MAPK
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-27T15:50:52.740Z
Answers.Text: If you have problems with installation, could you create an issue on Github, please? It might be related to python2 vs. python3. https://github.com/karel-brinda/samsift/issues
Answers.UserDisplayName: Karel Břinda
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-27T01:53:06.107Z
Answers.LastActivityDate: 2018-02-27T02:09:14.463Z
Answers.LastEditDate: 2018-02-27T02:09:14.463Z
Answers.UserDisplayName: Karel Břinda
Body: <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>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-02-27T13:54:21.573Z
Comments.Text: Have you looked into using the R Bioconductor Rsamtools::scanBam function? It can import a BAM file into R and then you can easily get the CIGAR string and check the CIGAR string for mismatches?
Comments.UserDisplayName: Biomagician
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-02-27T13:56:46.240Z
Comments.Text: Yes I tried using pileup package, but couldn't figure out how to get cigar string.
Comments.UserDisplayName: MAPK
Tags: <sam>
Title: How to you find read statistics for bam files with terminal mismatch?
ThreadText: # How to you find read statistics for bam files with terminal mismatch?
  
  (Asked by: username_0 on 2018-02-27)
  
  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).
  
  **Tags:** sam
  
  > **Comments:**
  > **username_1**: Have you looked into using the R Bioconductor Rsamtools::scanBam function? It can import a BAM file into R and then you can easily get the CIGAR string and check the CIGAR string for mismatches?  
  > **username_0**: Yes I tried using pileup package, but couldn't figure out how to get cigar string.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2018-02-27. Score: 2
  
  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](http://github.com/username_2-brinda/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).
  
  > **Comments:**
  > **username_0**: Thanks for your answer. I have updated the question with how I created BAM file.  
  > **username_0**: Could you explain this a bit please `MD[-1]=="0" and not MD[-2].isdigit()`?  
  > **username_0**: Does this also consider for `+` and `-` strands?  
  > **username_2**: Please, could you show at least one line that should pass the filter? I need it to verify the `MD` filtering step. Yes, it's possible to filter using `+`/`-` via `FLAG`.  
  > **username_0**: The bowtie -v 1 option allows one mismatch, -v 2, alows 2  mismatches anywhere in the sequences. The default option allows maximum of 3 mismatches.  
  > **username_2**: I mean 1 line from the SAM file.  
  > **username_0**: Did you mean: `samtools view -H R1_001.sam        @HD VN:1.0 SO:unsorted @SQ SN:NC_013116 LN:2166 @PG ID:Bowtie VN:1.1.2 CL:"bowtie --wrapper basic-0 -S -p 18 /m`  
  > **username_0**: Why am I getting permission denied while using samsift? I did chmod u+x.  
  > **username_0**: `/usr/bin/python: can't find '__main__' module in 'samsift'` Getting this error now?  
  > **username_2**: If you have problems with installation, could you create an issue on Github, please? It might be related to python2 vs. python3. https://github.com/username_2-brinda/samsift/issues  
  
  ---
  
  