I've been running snippy with M. tuberculosis samples, and I've come across the following issue.
I run snippy like this:
snippy --outdir 'QC_2025_nr14_filt' --R1 'QC_2025_nr14_filt_1.fastp.fastq.gz' --R2 'QC_2025_nr14_filt_2.fastp.fastq.gz' --mincov 10 --mapqual 10 --basequal 5 --minqual 30 --ref AncMTB.fasta --cpus 5 &
This process generates first a sorted .bam file that is then used by freebayes for variant calling, which is run like this:
freebayes-parallel reference/ref.txt 5 -p 2 -P 0 -C 2 -F 0.05 --min-coverage 0 --min-repeat-entropy 1.0 -q 5 -m 10 --strict-vcf -f reference/ref.fa snps.bam > snps.raw.vcf
This process generates a snps.raw.vcf file that contains all variants detected for such sample, before performing any successive filtering procedures. bcftools view is run afterwards, which is when these filters are applied (in this case, I'm aware that they can be modified with the corresponding parameters), generating the final snps.vcf:
bcftools view --include 'FMT/GT="1/1" && QUAL>=30 && FMT/DP>=0 && (FMT/AO)/(FMT/DP)>=0' snps.raw.vcf`
The issue here is that this command discards all variants that have been detected in heterozygosity, no matter the proportion existing between the reference and alternative alleles. Check this example:
SAMPLE EQA7-25:
$ cat snps.raw.vcf | grep 104712
MTB_anc 104712 . T C 3578.25 . AB=0;ABP=0;AC=2;AF=1;AN=2;AO=123;CIGAR=1X;DP=132;DPB=132;DPRA=0;EPP=15.8802;EPPR=20.3821 ;GTI=0;LEN=1;MEANALT=2;MQM=54.4715;MQMR=38.125;NS=1;NUMALT=1;ODDS=118.12;PAIRED=0.910569;PAIREDR=1;PAO=0;PQA=0;PQR=0;PRO=0;QA=4397;QR=296;RO=8;R PL=35;RPP=52.6011;RPPR=20.3821;RPR=88;RUN=1;SAF=69;SAP=6.98251;SAR=54;SRF=8;SRP=20.3821;SRR=0;TYPE=snp GT:DP:AD:RO:QR:AO:QA:GL 1/1:132:8,123:8: 296:123:4397:-360.795,-14.8772,0
$ cat snps.filt.vcf | grep 104712
MTB_anc 104712 . T C 3578.25 . AB=0;AO=123;DP=132;QA=4397;QR=296;RO=8;TYPE=snp GT:DP:RO:QR:AO:QA:GL 1/1:132:8:296:12 3:4397:-360.795,-14.8772,0
$ cat snps.vcf | grep 104712
MTB_anc 104712 . T C 3578.25 . AB=0;AO=123;DP=132;QA=4397;QR=296;RO=8;TYPE=snp GT:DP:RO:QR:AO:QA:GL 1/1:132:8:296:12 3:4397:-360.795,-14.8772,0
SAMPLE EQA5-25:
$ cat snps.raw.vcf | grep 104712
MTB_anc 104712 . T C 4549.05 . AB=0.879121;ABP=230.227;AC=1;AF=0.5;AN=2;AO=160;CIGAR=1X;DP=182;DPB=182;DPRA=0;EPP=16.90 77;EPPR=48.6112;GTI=0;LEN=1;MEANALT=2;MQM=54.8125;MQMR=37.4286;NS=1;NUMALT=1;ODDS=97.0257;PAIRED=0.925;PAIREDR=1;PAO=0;PQA=0;PQR=0;PRO=0;QA=5745 ;QR=739;RO=21;RPL=53;RPP=42.5854;RPPR=48.6112;RPR=107;RUN=1;SAF=91;SAP=9.579;SAR=69;SRF=21;SRP=48.6112;SRR=0;TYPE=snp GT:DP:AD:RO:QR:AO:QA:GL0 /1:182:21,160:21:739:160:5745:-451.732,0,-6.13472
$ cat snps.filt.vcf | grep 104712
$ cat snps.vcf | grep 104712
$
As you can tell, in this position the reference allele is T, while the alternative allele is C.
- In the sample called EQA7-25, freebayes estimates that this sample is homozygous in this position for the alternative allele. Therefore, this position passes all the filters that are applied when running
bcftools view and appears in the final snps.vcf file, which is then used by snippy-core.
- In EQA5-25, freebayes estimates that this sample is heterozygous in this position.
- However, in this sample, the number of reads that have the reference allele is 21 (RO=21), and the number of reads that have the alternative allele is 160 (AO=160) so, even if there is heterozygosity in this position, the alternative allele is majority in this position (since it represents the 88% of the reads).
- The alternative allele should be called in this position for this sample since it's majority, but this does not happen.
- Since there is heterozygosity in this position, the
bcftools view --include 'FMT/GT="1/1" parameter discards this position regardless of the allele balance, no variant is reported in this site for this sample, even though there is an actual variant in that position, which can be checked in IGV when viewing the corresponding .bam file. This happens no matter what the allele ratio is: heterozygous sites are discarded even if there is a 50/50, 90/10 or 10/90 ratio, and this is the cause for false negatives existing for a certain allele in a certain position.
This is problematic to us since positions in the snps.raw.vcf that do not comply with the --include 'FMT/GT="1/1" criterion are discarded in the final snps.vcf file plus these positions are marked with "N" in the snps.aligned.fa file that is later on used by snippy-core. These positions do not appear in the core files and cannot be used for the construction of neither an SNP matrix nor a phylogenetic tree.
This should be taken into account for evaluation since this causes false negatives in our analyses, which means that not all true variants present in a sample are being reported by snippy, and therefore SNP matrices are not as reliable as they should be.
Even though the bcftools view --include 'FMT/GT="1/1" command is not wrong, the way snippy considers and discards variants should be done in another way, so that balanced heterozygous sites are discarded but unbalanced heterozygous sites are not discarded, so that snippy chooses the majority allele in order to avoid that this site is not present in the core files. For example, a certain threshold for the allele frequencies might be established for this purpose, so that an "N" is added in that position in the snps.aligned.fa file in cases of true heterozygosity, or the majority allele is added when the heterozygosity in unbalanced.
Please let me know if you need any more info from my side. Thanks a lot in advance!
I've been running snippy with M. tuberculosis samples, and I've come across the following issue.
I run
snippylike this:This process generates first a sorted .bam file that is then used by freebayes for variant calling, which is run like this:
This process generates a
snps.raw.vcffile that contains all variants detected for such sample, before performing any successive filtering procedures.bcftools viewis run afterwards, which is when these filters are applied (in this case, I'm aware that they can be modified with the corresponding parameters), generating the finalsnps.vcf:The issue here is that this command discards all variants that have been detected in heterozygosity, no matter the proportion existing between the reference and alternative alleles. Check this example:
SAMPLE EQA7-25:
SAMPLE EQA5-25:
As you can tell, in this position the reference allele is T, while the alternative allele is C.
bcftools viewand appears in the finalsnps.vcffile, which is then used bysnippy-core.bcftools view --include 'FMT/GT="1/1"parameter discards this position regardless of the allele balance, no variant is reported in this site for this sample, even though there is an actual variant in that position, which can be checked in IGV when viewing the corresponding .bam file. This happens no matter what the allele ratio is: heterozygous sites are discarded even if there is a 50/50, 90/10 or 10/90 ratio, and this is the cause for false negatives existing for a certain allele in a certain position.This is problematic to us since positions in the
snps.raw.vcfthat do not comply with the--include 'FMT/GT="1/1"criterion are discarded in the finalsnps.vcffile plus these positions are marked with "N" in thesnps.aligned.fafile that is later on used bysnippy-core. These positions do not appear in the core files and cannot be used for the construction of neither an SNP matrix nor a phylogenetic tree.This should be taken into account for evaluation since this causes false negatives in our analyses, which means that not all true variants present in a sample are being reported by snippy, and therefore SNP matrices are not as reliable as they should be.
Even though the
bcftools view --include 'FMT/GT="1/1"command is not wrong, the way snippy considers and discards variants should be done in another way, so that balanced heterozygous sites are discarded but unbalanced heterozygous sites are not discarded, so that snippy chooses the majority allele in order to avoid that this site is not present in the core files. For example, a certain threshold for the allele frequencies might be established for this purpose, so that an "N" is added in that position in thesnps.aligned.fafile in cases of true heterozygosity, or the majority allele is added when the heterozygosity in unbalanced.Please let me know if you need any more info from my side. Thanks a lot in advance!