Combined HiFi and Illumina SNP calling not possible?

[astridboehne]

Hi all, we had raised this as an issue in June but did not get any reply:

Hello,

I am trying to use gatk/4.1.4.1 and picard/2.22.0 to do joint variant calling of PacBio HiFi reads and Illumina short-reads.
My pipeline is basically the recommended one (without base recalibration) where after HaplotypeCaller, I use GenomicsDBImport and GenotypeGVCFs and end up with the vcf file.

The HiFI reads have normal hifi phred scores up to 93 and the illumina are encoded with phred33 quality scores.

The output of the joint variant calling is strange and it causes issues when I try to use it with tools like plink and it makes me wonder if the joint variant calling went badly as well. This is even after variant filtration where I filter for QUAL<30 and QD<2

Here is an example of what the first few lines of my vcf calls look like. Notice the QUAL column that looks like Num,Num:

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Illumina_sample0 sample1 sample2 sample_3 sample4 sample5 sample6 hifi_sample
Chr1 10608 . GCA G 174,34 PASS AC=2;AF=0.143;AN=14;BaseQRankSum=0.00;DP=55;ExcessHet=3.3579;FS=7.404;MLEAC=2;MLEAF=0.143;MQ=57.97;MQRankSum=0.00;QD=15.85;ReadPosRankSum=1.00;SOR=3.611 GT:AD:DP:GQ:PL ./.:0,0:0:.:0,0,0 0/1:3,3:6:99:117,0,117 0/0:5,0:5:15:0,15,141 0/1:3,2:5:71:71,0,120 0/0:3,0:3:9:0,9,113 0/0:7,0:7:21:0,21,190 0/0:5,0:5:15:0,15,175 0/0:20,0:20:60:0,60,900
Chr1 10616 . G A 903,42 PASS AC=8;AF=0.571;AN=14;BaseQRankSum=0.00;DP=66;ExcessHet=0.7136;FS=9.277;MLEAC=9;MLEAF=0.643;MQ=59.80;MQRankSum=0.00;QD=20.08;ReadPosRankSum=1.00;SOR=0.251 GT:AD:DP:GQ:PL ./.:0,0:0:.:0,0,0 0/0:10,0:10:30:0,30,367 1/1:0,5:5:15:141,15,0 0/1:2,3:5:75:110,0,75 0/0:3,0:3:9:0,9,113 1/1:0,8:8:24:232,24,0 1/1:0,7:7:21:224,21,0 0/1:14,6:20:99:210,0,570

I believe this can be reproduced with any illumina+hifi samples. Are you aware of this problem? Could the snp calling be wrong because of the different phred score encoding? What do you suggest to do in these cases?

I have read your response in this issue:
https://gatk.broadinstitute.org/hc/en-us/community/posts/360072716972-Variant-calling-with-PacBio-HiFi-reads and have used minimap2 in a similar way. (minimap2 -acyYL --secondary=no --MD --eqx -x asm20 -k 20)

Any advice on this from the community?

[astridboehne]

here the GATK error log
Using GATK jar /share/scientific_bin/gatk/4.1.4.1/gatk-package-4.1.4.1-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /share/scientific_bin/gatk/4.1.4.1/gatk-package-4.1.4.1-local.jar IndexFeatureFile -I output/called/final/allsites.filtered.vcf
00:57:08.257 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/share/scientific_bin/gatk/4.1.4.1/gatk-package-4.1.4.1-local.jar!/com/intel/gkl/native/libgkl_compression.so
Sep 11, 2023 12:57:08 AM shaded.cloud_nio.com.google.auth.oauth2.ComputeEngineCredentials runningOnComputeEngine
INFO: Failed to detect whether we are running on Google Compute Engine.
00:57:08.789 INFO IndexFeatureFile - ------------------------------------------------------------
00:57:08.789 INFO IndexFeatureFile - The Genome Analysis Toolkit (GATK) v4.1.4.1
00:57:08.789 INFO IndexFeatureFile - For support and documentation go to https://software.broadinstitute.org/gatk/
00:57:08.790 INFO IndexFeatureFile - Executing as ychrysostomakis@compute-0-3.local on Linux v3.10.0-1160.53.1.el7.x86_64 amd64
00:57:08.790 INFO IndexFeatureFile - Java runtime: Java HotSpot(TM) 64-Bit Server VM v1.8.0_231-b11
00:57:08.790 INFO IndexFeatureFile - Start Date/Time: 11. September 2023 00:57:07 MESZ
00:57:08.790 INFO IndexFeatureFile - ------------------------------------------------------------
00:57:08.790 INFO IndexFeatureFile - ------------------------------------------------------------
00:57:08.790 INFO IndexFeatureFile - HTSJDK Version: 2.21.0
00:57:08.790 INFO IndexFeatureFile - Picard Version: 2.21.2
00:57:08.790 INFO IndexFeatureFile - HTSJDK Defaults.COMPRESSION_LEVEL : 2
00:57:08.790 INFO IndexFeatureFile - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
00:57:08.790 INFO IndexFeatureFile - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
00:57:08.790 INFO IndexFeatureFile - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
00:57:08.790 INFO IndexFeatureFile - Deflater: IntelDeflater
00:57:08.790 INFO IndexFeatureFile - Inflater: IntelInflater
00:57:08.790 INFO IndexFeatureFile - GCS max retries/reopens: 20
00:57:08.790 INFO IndexFeatureFile - Requester pays: disabled
00:57:08.790 INFO IndexFeatureFile - Initializing engine
00:57:08.790 INFO IndexFeatureFile - Done initializing engine
00:57:09.166 INFO FeatureManager - Using codec VCFCodec to read file file:///share/pool/CompGenomVert/phoxy_snp_calling/VC_HIFI/output/called/final/allsites.filtered.vcf
00:57:09.250 INFO ProgressMeter - Starting traversal
00:57:09.250 INFO ProgressMeter - Current Locus Elapsed Minutes Records Processed Records/Minute
00:57:09.251 INFO IndexFeatureFile - Shutting down engine
[11. September 2023 00:57:09 MESZ] org.broadinstitute.hellbender.tools.IndexFeatureFile done. Elapsed time: 0.02 minutes.
Runtime.totalMemory()=2114453504
java.lang.NumberFormatException: For input string: "175,24"
at sun.misc.FloatingDecimal.readJavaFormatString(FloatingDecimal.java:2043)
at sun.misc.FloatingDecimal.parseDouble(FloatingDecimal.java:110)
at java.lang.Double.parseDouble(Double.java:538)
at htsjdk.variant.vcf.VCFUtils.parseVcfDouble(VCFUtils.java:262)
at htsjdk.variant.vcf.AbstractVCFCodec.parseQual(AbstractVCFCodec.java:620)
at htsjdk.variant.vcf.AbstractVCFCodec.parseVCFLine(AbstractVCFCodec.java:422)
at htsjdk.variant.vcf.AbstractVCFCodec.decodeLine(AbstractVCFCodec.java:384)
at htsjdk.variant.vcf.AbstractVCFCodec.decode(AbstractVCFCodec.java:328)
at htsjdk.variant.vcf.AbstractVCFCodec.decode(AbstractVCFCodec.java:48)
at htsjdk.tribble.AsciiFeatureCodec.decode(AsciiFeatureCodec.java:70)
at htsjdk.tribble.AsciiFeatureCodec.decode(AsciiFeatureCodec.java:37)
at htsjdk.tribble.AbstractFeatureCodec.decodeLoc(AbstractFeatureCodec.java:43)
at org.broadinstitute.hellbender.utils.codecs.ProgressReportingDelegatingCodec.decodeLoc(ProgressReportingDelegatingCodec.java:46)
at htsjdk.tribble.index.IndexFactory$FeatureIterator.readNextFeature(IndexFactory.java:689)
at htsjdk.tribble.index.IndexFactory$FeatureIterator.(IndexFactory.java:606)
at htsjdk.tribble.index.IndexFactory.createDynamicIndex(IndexFactory.java:446)
at org.broadinstitute.hellbender.tools.IndexFeatureFile.createAppropriateIndexInMemory(IndexFeatureFile.java:118)
at org.broadinstitute.hellbender.tools.IndexFeatureFile.doWork(IndexFeatureFile.java:75)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:139)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:191)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:210)
at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:163)
at org.broadinstitute.hellbender.Main.mainEntry(Main.java:206)
at org.broadinstitute.hellbender.Main.main(Main.java:292)
Full Traceback (most recent call last):
File "/home/ychrysostomakis/.local/lib/python3.9/site-packages/snakemake/executors/init.py", line 2578, in run_wrapper
run(
File "/share/pool/CompGenomVert/phoxy_snp_calling/VC_HIFI/joint_snp_calling.smk", line 422, in __rule_index_feature_file
File "/home/ychrysostomakis/.local/lib/python3.9/site-packages/snakemake/shell.py", line 300, in new
raise sp.CalledProcessError(retcode, cmd)
subprocess.CalledProcessError: Command 'set -euo pipefail;
module load gatk/4.1.4.1
gatk IndexFeatureFile -I output/called/final/allsites.filtered.vcf' returned non-zero exit status 3.