Hello,

I'm working on variant calling to compare two strains and identify variants unique to strain 2. I performed RNA-seq on three replicates from each strain. I aligned all the samples to a reference genome using STAR, then used bedtools intersect to extract the genes positions I was interested in.

With the 6 new BAM files (2 strains X 3 replicates), I performed the following steps:

• samtools sort, samtools index, samtools rmdup, samtools index
• bcftools mpileup, bcftools call on the three replicates for each strain, which resulted in two VCF files (one for each strain).
• I applied tabix to both VCF files and used bcftools view.
• Finally, I filtered for variants that appear in the VCF for strain 2 but not in the VCF for strain 1, and applied filters for depth (DP) and quality (QUAL).

The issue I'm facing is that I’m seeing variants in strain 2 that don’t appear in strain 1. However, upon examining the data in IGV, it seems that these variants unique to strain 2 may be due to a lack of coverage (or very low coverage) in strain 1, not an actual difference between the two strains. In other words, I’m not sure if these variants in strain 2 are truly absent in strain 1, or if they just aren’t detectable in strain 1 due to poor or no coverage.

How can I determine which variants are truly unique to strain 2 and not just due to a lack of coverage in strain 1?

I’d greatly appreciate any help or suggestions.
Thank you!