1
$\begingroup$

I received (from manufacturer) several .bam files and I used four callers (samtools, freebayes, haplotypecaller, deepvariant) to find some sequence variants. In obtained .vcf files, I took a closer look to some calls. I found interesting, homozygous one rs477033 (C/G Ref/Alt) with flag 'COMMON=0' and very low MAF. I also checked the .bam file (in IGV) and this rs itself in ensembl database. It looks like this:

igv

The problem is: in this case ref and alt are switched. The G is reference (ancestral) allele and C is altered. Nevertheless in the .bam file it's interpreted in opposite way. Of course it's false positive variant - there's no mutation there at all. All variant callers returned it as high quality variant. My. bam files are full of them.

What went wrong in this case? Problem with mapping?

How can I fix this thing, to get only novel and rare true positive variants?

$\endgroup$
2
  • $\begingroup$ [link] drive.google.com/file/d/1eheG5MqS_a1xK5HSqsNJ0Tito_WdL_YV/… $\endgroup$ Commented Jun 26, 2022 at 1:34
  • $\begingroup$ A reference sequence can obviously not represent all populations, as shown above in the linked chapter (• page 30 •), a reference allele may be the minor allele in the population you are studying ( in some cases ) $\endgroup$ Commented Jun 26, 2022 at 1:36

1 Answer 1

3
$\begingroup$

The ancestral allele is not always the reference allele. Ref and alt aren't switched, ref is just rare in at least the population you're looking at (it appears to be rare in humans as a whole). So, the C is indeed reference and G is indeed altered, it's just that people homozygous for the G are the overwhelming majority. This is evidenced by the rs477033 variant being named NC_000001.11:g.165890393C>G. So, in short, nothing went wrong, all tools functioned correctly.

You'll notice a bunch of these over time, either make a database of them (or find one, someone may have made one) and filter them out or see if maybe the 1000 genomes reference genome has some of these modified (or maybe the T2T genome that was recently published).

$\endgroup$
1
  • 1
    $\begingroup$ I think that I found the solution and I'd like to add it here. It might be helpful for others. I just simply performed additional annotation of my .vcf files using 1KGenomes data. I've found it here: here Then each variant is described in INFO by AA which means ancestral allele. If this allele is different comparing to my REF I filter it out. $\endgroup$
    – Adamm
    Commented Jun 23, 2022 at 10:23

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.