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)

  1. 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.

  • $\begingroup$ 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. $\endgroup$
    – terdon
    Apr 17, 2023 at 22:13
  • $\begingroup$ @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. $\endgroup$
    – Ahsan
    Apr 18, 2023 at 18:52

1 Answer 1


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.