VariantRecalibrator fails on 30x chr22 subset (GIAB NA12878)

splaisansplaisan Leuven (Belgium)Member ✭✭

I did a full workflow with the 300x depth data but now, retrying with a subset of 30x depth fails (required for a training session).

Is there a requirement for recalibration that is set too high by default and fails on my data?
Or is it because my data is present in the calibration collections.
This GIAB data was mapped then called using only chr22 reads (10% of the 300x) against HG38 and led to a gvcf from which I derived the vcf

-rw-r--r-- 1 root domain users 142M Aug  3 11:17 NA12878_0.1.g.vcf.gz
-rw-r--r-- 1 root domain users  60K Aug  3 11:17 NA12878_0.1.g.vcf.gz.tbi
-rw-r--r-- 1 root domain users 4.6M Aug  3 11:22 NA12878_0.1.vcf.gz
-rw-r--r-- 1 root domain users  24K Aug  3 11:22 NA12878_0.1.vcf.gz.tbi

Command and output are below. I can add any useful extract on your go.

Thanks a lot for your helping me find why this fails.

Stephane

# True sites training resource: HapMap
truetraining15=reference/hg38_v0_hapmap_3.3.hg38.vcf.gz

# True sites training resource: Omni
truetraining12=reference/hg38_v0_1000G_omni2.5.hg38.vcf.gz

# Non-true sites training resource: 1000G
nontruetraining10=reference/hg38_v0_Mills_and_1000G_gold_standard.indels.hg38.vcf.gz

# Known sites resource, not used in training: dbSNP
knowntraining2=reference/hg38_v0_Homo_sapiens_assembly38.dbsnp138.vcf.gz

# indels True sites training resource: Mills
truetrainingindel12=reference/hg38_v0_Mills_and_1000G_gold_standard.indels.hg38.vcf.gz

java -Xmx${maxmem} -jar $GATK/gatk.jar VariantRecalibrator \
-R $BWA_INDEXES/NCBI_GRCh38.fa \
-V ${mappings}/${samplename}_${p}.vcf.gz \
--resource hapmap,known=false,training=true,truth=true,prior=15.0:${truetraining15} \
--resource omni,known=false,training=true,truth=true,prior=12.0:${truetraining12} \
--resource 1000g,known=false,training=true,truth=false,prior=10.0:${nontruetraining10} \
--resource dbsnp,known=true,training=false,truth=false,prior=2.0:${knowntraining2} \
--resource Mills_and_1000G_gold,known=false,training=true,truth=true,prior=12.0:${truetrainingindel12} \
-an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -an DP \
--mode BOTH \
--output ${mappings}/output.recal_${p}.vcf \
--tranches-file ${mappings}/output.tranches_${p} \
--rscript-file ${mappings}/output.plots_${p}.R

stderr

java -Xmx${maxmem} -jar $GATK/gatk.jar VariantRecalibrator \
> -R $BWA_INDEXES/NCBI_GRCh38.fa \
> -V ${mappings}/${samplename}_${p}.vcf.gz \
> --resource hapmap,known=false,training=true,truth=true,prior=15.0:${truetraining15} \
> --resource omni,known=false,training=true,truth=true,prior=12.0:${truetraining12} \
> --resource 1000g,known=false,training=true,truth=false,prior=10.0:${nontruetraining10} \
> --resource dbsnp,known=true,training=false,truth=false,prior=2.0:${knowntraining2} \
> --resource Mills_and_1000G_gold,known=false,training=true,truth=true,prior=12.0:${truetrainingindel12} \
> -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -an DP \
> --mode BOTH \
> --output ${mappings}/output.recal_${p}.vcf \
> --tranches-file ${mappings}/output.tranches_${p} \
> --rscript-file ${mappings}/output.plots_${p}.R
11:33:58.090 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/opt/biotools/gatk-4.0.7.0/gatk-package-4.0.7.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
11:33:58.338 INFO  VariantRecalibrator - ------------------------------------------------------------
11:33:58.339 INFO  VariantRecalibrator - The Genome Analysis Toolkit (GATK) v4.0.7.0
11:33:58.339 INFO  VariantRecalibrator - For support and documentation go to https://software.broadinstitute.org/gatk/
11:33:58.339 INFO  VariantRecalibrator - Executing as [email protected] on Linux v4.4.0-131-generic amd64
11:33:58.339 INFO  VariantRecalibrator - Java runtime: OpenJDK 64-Bit Server VM v1.8.0_171-8u171-b11-0ubuntu0.16.04.1-b11
11:33:58.340 INFO  VariantRecalibrator - Start Date/Time: August 3, 2018 11:33:58 AM CEST
11:33:58.340 INFO  VariantRecalibrator - ------------------------------------------------------------
11:33:58.340 INFO  VariantRecalibrator - ------------------------------------------------------------
11:33:58.340 INFO  VariantRecalibrator - HTSJDK Version: 2.16.0
11:33:58.340 INFO  VariantRecalibrator - Picard Version: 2.18.7
11:33:58.340 INFO  VariantRecalibrator - HTSJDK Defaults.COMPRESSION_LEVEL : 2
11:33:58.340 INFO  VariantRecalibrator - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
11:33:58.340 INFO  VariantRecalibrator - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
11:33:58.340 INFO  VariantRecalibrator - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
11:33:58.340 INFO  VariantRecalibrator - Deflater: IntelDeflater
11:33:58.341 INFO  VariantRecalibrator - Inflater: IntelInflater
11:33:58.341 INFO  VariantRecalibrator - GCS max retries/reopens: 20
11:33:58.341 INFO  VariantRecalibrator - Using google-cloud-java patch 6d11bef1c81f885c26b2b56c8616b7a705171e4f from https://github.com/droazen/google-cloud-java/tree/dr_all_nio_fixes
11:33:58.341 INFO  VariantRecalibrator - Initializing engine
11:33:58.804 INFO  FeatureManager - Using codec VCFCodec to read file file:///mnt/freenas/NGS_Variant-Analysis-training2018/reference/hg38_v0_hapmap_3.3.hg38.vcf.gz
11:33:59.042 INFO  FeatureManager - Using codec VCFCodec to read file file:///mnt/freenas/NGS_Variant-Analysis-training2018/reference/hg38_v0_1000G_omni2.5.hg38.vcf.gz
11:33:59.177 INFO  FeatureManager - Using codec VCFCodec to read file file:///mnt/freenas/NGS_Variant-Analysis-training2018/reference/hg38_v0_Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
11:33:59.290 INFO  FeatureManager - Using codec VCFCodec to read file file:///mnt/freenas/NGS_Variant-Analysis-training2018/reference/hg38_v0_Homo_sapiens_assembly38.dbsnp138.vcf.gz
11:33:59.402 INFO  FeatureManager - Using codec VCFCodec to read file file:///mnt/freenas/NGS_Variant-Analysis-training2018/reference/hg38_v0_Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
11:33:59.513 INFO  FeatureManager - Using codec VCFCodec to read file file:///mnt/freenas/NGS_Variant-Analysis-training2018/bwa_mappings_10pc/NA12878_0.1.vcf.gz
11:33:59.633 INFO  VariantRecalibrator - Done initializing engine
11:33:59.645 INFO  TrainingSet - Found hapmap track:    Known = false   Training = true         Truth = true    Prior = Q15.0
11:33:59.646 INFO  TrainingSet - Found omni track:      Known = false   Training = true         Truth = true    Prior = Q12.0
11:33:59.646 INFO  TrainingSet - Found 1000g track:     Known = false   Training = true         Truth = false   Prior = Q10.0
11:33:59.646 INFO  TrainingSet - Found dbsnp track:     Known = true    Training = false        Truth = false   Prior = Q2.0
11:33:59.646 INFO  TrainingSet - Found Mills_and_1000G_gold track:      Known = false   Training = true         Truth = true    Prior = Q12.0
11:33:59.693 INFO  ProgressMeter - Starting traversal
11:33:59.693 INFO  ProgressMeter -        Current Locus  Elapsed Minutes    Variants Processed  Variants/Minute
11:34:08.098 INFO  ProgressMeter -       chr22:50362905              0.1                 81887         584628.7
11:34:08.098 INFO  ProgressMeter - Traversal complete. Processed 81887 total variants in 0.1 minutes.
11:34:08.113 INFO  VariantDataManager - QD:      mean = 19.68    standard deviation = 9.58
11:34:08.124 INFO  VariantDataManager - MQ:      mean = 59.83    standard deviation = 1.70
11:34:08.132 INFO  VariantDataManager - MQRankSum:       mean = -0.02    standard deviation = 0.29
11:34:08.145 INFO  VariantDataManager - ReadPosRankSum:          mean = 0.03     standard deviation = 1.00
11:34:08.157 INFO  VariantDataManager - FS:      mean = 1.97     standard deviation = 3.36
11:34:08.165 INFO  VariantDataManager - SOR:     mean = 1.02     standard deviation = 0.58
11:34:08.173 INFO  VariantDataManager - DP:      mean = 25.19    standard deviation = 7.16
11:34:08.276 INFO  VariantDataManager - Annotations are now ordered by their information content: [MQ, DP, QD, MQRankSum, FS, SOR, ReadPosRankSum]
11:34:08.284 INFO  VariantDataManager - Training with 27312 variants after standard deviation thresholding.
11:34:08.288 INFO  GaussianMixtureModel - Initializing model with 100 k-means iterations...
11:34:09.210 INFO  VariantRecalibratorEngine - Finished iteration 0.
11:34:09.860 INFO  VariantRecalibratorEngine - Finished iteration 5.    Current change in mixture coefficients = 2.94951
11:34:10.490 INFO  VariantRecalibratorEngine - Finished iteration 10.   Current change in mixture coefficients = 0.39301
11:34:11.065 INFO  VariantRecalibratorEngine - Finished iteration 15.   Current change in mixture coefficients = 0.00837
11:34:11.507 INFO  VariantRecalibratorEngine - Finished iteration 20.   Current change in mixture coefficients = 0.01513
11:34:12.037 INFO  VariantRecalibratorEngine - Finished iteration 25.   Current change in mixture coefficients = 0.01792
11:34:12.531 INFO  VariantRecalibratorEngine - Finished iteration 30.   Current change in mixture coefficients = 0.01851
11:34:12.998 INFO  VariantRecalibratorEngine - Finished iteration 35.   Current change in mixture coefficients = 0.02430
11:34:13.451 INFO  VariantRecalibratorEngine - Finished iteration 40.   Current change in mixture coefficients = 0.01579
11:34:13.891 INFO  VariantRecalibratorEngine - Finished iteration 45.   Current change in mixture coefficients = 0.00536
11:34:14.433 INFO  VariantRecalibratorEngine - Finished iteration 50.   Current change in mixture coefficients = 0.00169
11:34:14.433 INFO  VariantRecalibratorEngine - Convergence after 50 iterations!
11:34:14.549 WARN  VariantRecalibratorEngine - Model could not pre-compute denominators.
11:34:14.567 INFO  VariantDataManager - Selected worst 0 scoring variants --> variants with LOD <= -5.0000.
11:34:14.593 INFO  VariantRecalibrator - Shutting down engine
[August 3, 2018 11:34:14 AM CEST] org.broadinstitute.hellbender.tools.walkers.vqsr.VariantRecalibrator done. Elapsed time: 0.28 minutes.
Runtime.totalMemory()=5262802944
java.lang.IllegalArgumentException: No data found.
        at org.broadinstitute.hellbender.tools.walkers.vqsr.VariantRecalibratorEngine.generateModel(VariantRecalibratorEngine.java:34)
        at org.broadinstitute.hellbender.tools.walkers.vqsr.VariantRecalibrator.onTraversalSuccess(VariantRecalibrator.java:630)
        at org.broadinstitute.hellbender.engine.GATKTool.doWork(GATKTool.java:981)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:137)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:182)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:201)
        at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:160)
        at org.broadinstitute.hellbender.Main.mainEntry(Main.java:203)
        at org.broadinstitute.hellbender.Main.main(Main.java:289)

Best Answer

Answers

Sign In or Register to comment.