calculate mismatch frequency/rate from a BAM file

Question

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.

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