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.
samtools mpileup I can get total number of all kind of mismatches for a
read group, however, takes some time.
- How do I calculate background mismatch rate from this?
- How do I calculate total mismatches in a bam?
I was thinking about something like
for a read group:
for a BAM file:
- 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.