Sorry. I've ignored issues report after leaving academia. I have absolutely no idea what is going on here and would love to have a look. Any chance you could share a BAM slice?

Hi Dario, to which email address did you send the file? Could you either upload it and share the link here or share with my private email address: [email protected]. Thanks

Thanks Dario. I can' t reproduce the problem. Which version are you using? With 2.1.4 (which should in this respect be identical to 2.1.3) gives the following:

$ lofreq call test_file_chr10.bam -r 10:123241496-123241496 -f human_g1k_v37.fasta | tail -n 1
10      123241496       .       T       C       3114    PASS    DP=123;AF=1.000000;SB=0;DP4=0,0,52,71

Sorry Andreas but you need the quality parameters -Q 30 -q 30 -m 30. Anyway I'm using version 2.1.3
My point is, irrespective of quality filters the AF shouldn't change when you have only reads covering the reference allele.
I am also confused by the fact that the raw depth (DP) is 123 (while it's 137 according to IGV) even without any filter. Am I missing something?

Got it. This is a bug. The culprit is -q 30. This quality filter was meant to only affect variant quality computation, and not AF or DP etc, however it does (note to self: through the later used alt_counts). Fixing this is going to be hard and make the already convoluted logic and code even more convoluted. You won't like this, but I will leave this bug unfixed and this ticket open for two reasons:

1. You are actually not meant to filter on base or mapping quality excessively in LoFreq! Filtering introduces calling biases (and -q has the biggest effect, apart from the bug) in LoFreq's core which was built to deal with and is robust against quality variations (even though you might lose your favourite variant every now and then :))
2. I have started LoFreq3 (even though as a hobby, because I left academia; a static binary for call will appear within weeks) already, which is not affected by design (because it separates pileup/filtering from calling by design)

Sorry. Hope this helps anyway
Andreas

It certainly gives a lot of insights, so thank you for taking time to inspect the issue.

@andreas-wilm I have the same issue

NC_045512.2     23709   .       C       T       49314   PASS    DP=9307;AF=0.751262;SB=23;DP4=45,18,4892,4304

When I run other variant callers (BCFTools) and look at IGV, they all report AF > 0.99. If I actually calculate from DP4 i.e. (4892+4304)/(45+18+4892+4304) allele frequency > 0.99. Therefore, I couldn't figure out where is 0.75 AF calculation is coming from
However, you mentioned that quality filters while running lofreq can cause biases in calculations. Please see my run below:

lofreq call --call-indels --min-cov 50 -q 30 -Q 30 -m 20 --verbose -f path_to_reference/genome_ref.fasta -o sample_bwa_mq20pp_mark_dup_lofreq_realign_indelq_alnq_call.vcf sample_bwa_mq20pp_mark_dup_lofreq_realign_indelq_alnq.bam

If we don't include mapping and base quality filters at the variant calling run, then how do we filter on base and mapping qualities after variants have been called using lofreq call?

Thanks, Rahil

Also, I really need to filter my variants on allele frequency to get high frequency variants (> 90%). Could you please suggest if I calculating from DP4 would be fine?

Thanks, Rahil