gimbleprep filters out almost all SNPs

[PierreLesturgie]

Hi,

I think I may have encountered an issue with gimbleprep, as it appears to filter out almost all variants in my VCF:

• Initial number of variants = 6,384,223
• In gimble.vcf.gz = 64

This reduction from approximately 6.4 million variants to only 64 seems unexpectedly large. Do you have any idea which step of the preprocessing pipeline could be responsible for removing almost all the variants?

My pipeline:

(1) freebayes

I ran FreeBayes separately for each chromosome, splitting each chromosome into 100-kb regions for computational efficiency. For example, for chromosome 1:

GENOME=/reference_genome/merged_genomes.fasta
BAM_LIST=/bamfiles/data/bam.list
REGIONS=/freebayes/Chromosome1_100000_regions.txt
OUTDIR="Chromosome1"

mkdir -p "$OUTDIR"

parallel -j 128 \
    'freebayes \
        -f '$GENOME' \
        -L '$BAM_LIST' \
        -r {} \
        > '$OUTDIR'/{#}.vcf' \
    :::: $REGIONS

I then concatenated all regional and chromosome-level VCF files into a single compressed VCF using bcftools concat, producing: SNPs.vcf.gz

(2) gimbleprep

I ran:

GENOME=/reference_genome/merged_genomes.fasta
BAM_FOLDER=/bamfiles/data/
VCF=gimble/SNPs.vcf.gz

gimbleprep -f $GENOME -b $BAM_FOLDER -v $VCF --min_depth 10 --max_depth 6 -k

The VCF header of gimble.vcf.gz indicates that the steps were:

• bcftools norm
• vcfallelicprimitives
• filling AN and AC
• QUAL<1 filtering
• strand/read-position balance filtering
• --SnpGap 2
• removal of non-SNP records

I manually reproduced the normalization, decomposition, AN/AC, quality, and balance-filtering steps, and still retained several million variants.

Do you have any idea what could explain the discrepancy between the variants retained by gimbleprep and those retained when manually reproducing the filtering steps described in the intermediate VCF headers?

Thank you very much for your help.
Best,
Pierre