<blockquote>
<p>I am not using any variant caller</p>
</blockquote>
<p>Use a variant caller. <code>bcftools</code>, at least, will output variant frequency per variant:</p>
<pre><code>bcftools mpileup -ugf ref.fa sample.bam | bcftools call -mv &gt; output.vcf
</code></pre>
<p>As an alternative approach, I wrote <a href="https://gitlab.com/gringer/bioinfscripts/-/blob/master/readstomper.pl" rel="nofollow noreferrer">some code</a> to parse the mpileup output and extract counts for different alleles; this output may be easier for you to parse and calculate your desired frequencies:</p>
<pre><code># readstomper.pl -- generate base-packed proportion statistics for
#   output from 'samtools mpileup'
#
# example use:
# samtools mpileup -B -Q 0 -f circ-Nb-ec3-mtDNA.fasta LAST_OL_132394_vs_mtDNA.bam | \
#   /bioinf/scripts/readstomper.pl &gt; LAST_OLmpileupProp_132394_vs_mtDNA.txt.gz
#
# example output:
# Assembly,Position,Coverage,ref,cR,pR,A,C,G,T,d,i,InsMode
# Nb_mtDNA,11456,86,T,79.00,0.92,0.0000,0.0465,0.0000,0.0000,0.03,0.01,A;1.00
# Nb_mtDNA,11457,86,A,81.00,0.94,0.0000,0.0000,0.0233,0.0000,0.03,0.02,TA;0.50
# Nb_mtDNA,11458,86,G,79.00,0.92,0.0349,0.0233,0.0000,0.0116,0.01,0.00
# Nb_mtDNA,11459,86,A,78.00,0.92,0.0000,0.0000,0.0235,0.0471,0.01,0.00
# Nb_mtDNA,11460,86,T,5.00,0.06,0.0233,0.8140,0.0233,0.0000,0.08,0.00
# Nb_mtDNA,11461,86,G,57.00,0.66,0.1512,0.1163,0.0000,0.0465,0.02,0.01,AATAACAACACGTAACCGA;1.00
# Nb_mtDNA,11462,86,G,78.00,0.91,0.0349,0.0349,0.0000,0.0116,0.01,0.00
# Nb_mtDNA,11463,86,G,66.00,0.77,0.0465,0.0698,0.0000,0.0116,0.10,0.01,GCAATA;1.00
# Nb_mtDNA,11464,86,T,74.00,0.86,0.0000,0.0814,0.0233,0.0000,0.03,0.00

</code></pre>

CC BY-SA 4.0
2023-04-18T19:31:21.230Z
2023-04-18T19:31:21.230Z
gringer
<p>I am working with PAR-CLIP data where experimental procedure induces T to C mismatches. The other types of mismatches we can consider coming from artifact such as PCR errors. It could also be an existing SNP. For simplicity, let's consider T to C mismatches are experimentally induced and rest of the mismatches are background mismatch rate.</p>
<p>From Illumina NGS seq after performing alignment, for each of he <code>read groups</code> (reads collapsed in a 30-70 base region with minimum 10 read coverage). Up to this point its already done.</p>
<p>I want to compare T to C mismatches vs background mismatches. To do so I need to calculate the rate of signal mismatches (T to C) and background mismatches.</p>
<p>Using <code>samtools mpileup</code> I can get total number of all kind of mismatches for a <code>read group</code>, however, takes some time.</p>
<ol>
<li>How do I calculate background mismatch rate from this?</li>
<li>How do I calculate total mismatches in a bam?</li>
</ol>
<p>I was thinking about something like<br />
for a read group:
<code>mismatch_rate_group=total_mismatches_group/(total_read_group*total_base_group)</code>
for a BAM file:
<code>mismatch_rate_global=total_mismatches_bam/(total_reads_aligned_bam*genome_size_base)</code></p>
<ol start="3">
<li>What important factor I might missing in my approach?</li>
</ol>
<p>I  am not using any variant caller (e.g. bcftools). Initially I tried with samtools mpileup and then used a shell script to figure out the mismatches. But I am not sure is this is the right thing to do. Or even there are any easier methods. If one read has a mutation, I consider that a mismatch and I presume it to be artifact and/or PCR errors.</p>

CC BY-SA 4.0
2023-04-17T22:13:32.807Z
Please [edit] your question and explain this. You really should use a dedicated variant caller. NGS is error prone, and you cannot just assume that >1 read means the variant is real. In any case, writing your own thing to do this would be far more work and would almost certainly result in worse data than if you use a dedicated variant caller which has been designed for the job. If [edit] your question and give some more details (what sequencing technology, what species, what coverage etc) we will be able to recommend some tools.
terdon
CC BY-SA 4.0
2023-04-18T18:52:21.270Z
@terdon , thanks for your suggestion. I tried to add some more details in my question. I hope this is better now.
Additionally, `NGS is error prone` : I am less interested in the variant, rather I want to know the rate of such errors so that I can differentiate my signal from all of such errors.
Ahsan
<ngs><sam><point-mutation>
calculate mismatch frequency/rate from a BAM file
# calculate mismatch frequency/rate from a BAM file

(Asked by: username_0 on 2023-04-17)

I am working with PAR-CLIP data where experimental procedure induces T to C mismatches. The other types of mismatches we can consider coming from artifact such as PCR errors. It could also be an existing SNP. For simplicity, let's consider T to C mismatches are experimentally induced and rest of the mismatches are background mismatch rate.

From Illumina NGS seq after performing alignment, for each of he `read groups` (reads collapsed in a 30-70 base region with minimum 10 read coverage). Up to this point its already done.

I want to compare T to C mismatches vs background mismatches. To do so I need to calculate the rate of signal mismatches (T to C) and background mismatches.

Using `samtools mpileup` I can get total number of all kind of mismatches for a `read group`, however, takes some time.

1. How do I calculate background mismatch rate from this?
2. How do I calculate total mismatches in a bam?

I was thinking about something like  
for a read group:
`mismatch_rate_group=total_mismatches_group/(total_read_group*total_base_group)`
for a BAM file:
`mismatch_rate_global=total_mismatches_bam/(total_reads_aligned_bam*genome_size_base)`

3. What important factor I might missing in my approach?

I am not using any variant caller (e.g. bcftools). Initially I tried with samtools mpileup and then used a shell script to figure out the mismatches. But I am not sure is this is the right thing to do. Or even there are any easier methods. If one read has a mutation, I consider that a mismatch and I presume it to be artifact and/or PCR errors.

**Tags:** ngs, sam, point-mutation
> **Comments:**
> **username_1**: Please [edit] your question and explain this. You really should use a dedicated variant caller. NGS is error prone, and you cannot just assume that >1 read means the variant is real. In any case, writing your own thing to do this would be far more work and would almost certainly result in worse data than if you use a dedicated variant caller which has been designed for the job. If [edit] your question and give some more details (what sequencing technology, what species, what coverage etc) we will be able to recommend some tools.  
> **username_0**: @terdon , thanks for your suggestion. I tried to add some more details in my question. I hope this is better now. Additionally, `NGS is error prone` : I am less interested in the variant, rather I want to know the rate of such errors so that I can differentiate my signal from all of such errors.  

---

## 1 answers

### Answer by username_2 on 2023-04-18. Score: 0

I am not using any variant caller

Use a variant caller. bcftools, at least, will output variant frequency per variant:
bcftools mpileup -ugf ref.fa sample.bam | bcftools call -mv > output.vcf

As an alternative approach, I wrote some code to parse the mpileup output and extract counts for different alleles; this output may be easier for you to parse and calculate your desired frequencies:
# readstomper.pl -- generate base-packed proportion statistics for
#   output from 'samtools mpileup'
#
# example use:
# samtools mpileup -B -Q 0 -f circ-Nb-ec3-mtDNA.fasta LAST_OL_132394_vs_mtDNA.bam | \
#   /bioinf/scripts/readstomper.pl > LAST_OLmpileupProp_132394_vs_mtDNA.txt.gz
#
# example output:
# Assembly,Position,Coverage,ref,cR,pR,A,C,G,T,d,i,InsMode
# Nb_mtDNA,11456,86,T,79.00,0.92,0.0000,0.0465,0.0000,0.0000,0.03,0.01,A;1.00
# Nb_mtDNA,11457,86,A,81.00,0.94,0.0000,0.0000,0.0233,0.0000,0.03,0.02,TA;0.50
# Nb_mtDNA,11458,86,G,79.00,0.92,0.0349,0.0233,0.0000,0.0116,0.01,0.00
# Nb_mtDNA,11459,86,A,78.00,0.92,0.0000,0.0000,0.0235,0.0471,0.01,0.00
# Nb_mtDNA,11460,86,T,5.00,0.06,0.0233,0.8140,0.0233,0.0000,0.08,0.00
# Nb_mtDNA,11461,86,G,57.00,0.66,0.1512,0.1163,0.0000,0.0465,0.02,0.01,AATAACAACACGTAACCGA;1.00
# Nb_mtDNA,11462,86,G,78.00,0.91,0.0349,0.0349,0.0000,0.0116,0.01,0.00
# Nb_mtDNA,11463,86,G,66.00,0.77,0.0465,0.0698,0.0000,0.0116,0.10,0.01,GCAATA;1.00
# Nb_mtDNA,11464,86,T,74.00,0.86,0.0000,0.0814,0.0233,0.0000,0.03,0.00


---

