I am a technically-competent non-bioinformatician with a set of three WGS sequences for a grandmother, son, and daughter that were provided by Nebula Genomics (for each, I have two XXX.fq.gz files, along with XXX.cram, XXX.vcf, and XXX.tbi). While I have been able to inspect the genomes individually using tools like gene.iobio.io and IGV, and I have been able to get the GATK Docker container up and running, I have run into some errors while trying use the three genomes to generate phased haplotypes. In particular, when I run the following:
./gatk HaplotypeCaller -R /mnt/Backup/Medical/Genetic/hg38/resources-broad-hg38-v0-Homo_sapiens_assembly38.fasta -I /mnt/Backup/Medical/Genetic/person1.cram -I /mnt/Backup/Medical/Genetic/person2.cram -I /mnt/Backup/Medical/Genetic/person3.cram -O /mnt/Backup/Medical/Genetic/gatk_output.vcf.gz
I get the following output:
Using GATK jar /gatk/gatk-package-4.6.2.0-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 /gatk/gatk-package-4.6.2.0-local.jar HaplotypeCaller -R /mnt/Backup/Medical/Genetic/hg38/resources-broad-hg38-v0-Homo_sapiens_assembly38.fasta -I /mnt/Backup/Medical/Genetic/.person1.cram -I /mnt/Backup/Medical/Genetic/person2.cram -I /mnt/Backup/Medical/Genetic/person3.cram -O /mnt/Backup/Medical/Genetic/gatk_output.vcf.gz
18:45:26.985 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/gatk/gatk-package-4.6.2.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
18:45:27.109 INFO HaplotypeCaller - ------------------------------------------------------------
18:45:27.112 INFO HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.6.2.0
18:45:27.112 INFO HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/
18:45:27.112 INFO HaplotypeCaller - Executing as root@739efec904b5 on Linux v6.10.14-linuxkit amd64
18:45:27.112 INFO HaplotypeCaller - Java runtime: OpenJDK 64-Bit Server VM v17.0.12+7-Ubuntu-1ubuntu222.04
18:45:27.112 INFO HaplotypeCaller - Start Date/Time: May 28, 2025 at 6:45:26 PM GMT
18:45:27.113 INFO HaplotypeCaller - ------------------------------------------------------------
18:45:27.113 INFO HaplotypeCaller - ------------------------------------------------------------
18:45:27.114 INFO HaplotypeCaller - HTSJDK Version: 4.2.0
18:45:27.114 INFO HaplotypeCaller - Picard Version: 3.4.0
18:45:27.114 INFO HaplotypeCaller - Built for Spark Version: 3.5.0
18:45:27.117 INFO HaplotypeCaller - HTSJDK Defaults.COMPRESSION_LEVEL : 2
18:45:27.117 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
18:45:27.118 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
18:45:27.118 INFO HaplotypeCaller - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
18:45:27.118 INFO HaplotypeCaller - Deflater: IntelDeflater
18:45:27.118 INFO HaplotypeCaller - Inflater: IntelInflater
18:45:27.119 INFO HaplotypeCaller - GCS max retries/reopens: 20
18:45:27.119 INFO HaplotypeCaller - Requester pays: disabled
18:45:27.120 INFO HaplotypeCaller - Initializing engine
WARNING 2025-05-28 18:45:27 CRAMFileReader CRAM index file /mnt/Backup/Medical/Genetic/person1.cram.crai is older than CRAM /mnt/Backup/Medical/Genetic/person1.cram
WARNING 2025-05-28 18:45:28 CRAMFileReader CRAM index file /mnt/Backup/Medical/Genetic/person2.cram.crai is older than CRAM /mnt/Backup/Medical/Genetic/person2.cram
WARNING 2025-05-28 18:45:28 CRAMFileReader CRAM index file /mnt/Backup/Medical/Genetic/person3.cram.crai is older than CRAM /mnt/Backup/Medical/Genetic/person3.cram
18:45:31.908 INFO HaplotypeCaller - Done initializing engine
18:45:32.023 INFO NativeLibraryLoader - Loading libgkl_utils.so from jar:file:/gatk/gatk-package-4.6.2.0-local.jar!/com/intel/gkl/native/libgkl_utils.so
18:45:32.029 INFO NativeLibraryLoader - Loading libgkl_smithwaterman.so from jar:file:/gatk/gatk-package-4.6.2.0-local.jar!/com/intel/gkl/native/libgkl_smithwaterman.so
18:45:32.030 INFO SmithWatermanAligner - Using AVX accelerated SmithWaterman implementation
18:45:32.034 INFO HaplotypeCallerEngine - Disabling physical phasing, which is supported only for reference-model confidence output
18:45:32.044 INFO NativeLibraryLoader - Loading libgkl_pairhmm_omp.so from jar:file:/gatk/gatk-package-4.6.2.0-local.jar!/com/intel/gkl/native/libgkl_pairhmm_omp.so
18:45:32.059 INFO IntelPairHmm - Flush-to-zero (FTZ) is enabled when running PairHMM
18:45:32.060 INFO IntelPairHmm - Available threads: 4
18:45:32.061 INFO IntelPairHmm - Requested threads: 4
18:45:32.062 INFO PairHMM - Using the OpenMP multi-threaded AVX-accelerated native PairHMM implementation
18:45:32.230 INFO ProgressMeter - Starting traversal
18:45:32.231 INFO ProgressMeter - Current Locus Elapsed Minutes Regions Processed Regions/Minute
18:45:39.171 INFO VectorLoglessPairHMM - Time spent in setup for JNI call : 0.0
18:45:39.171 INFO PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 0.0
18:45:39.172 INFO SmithWatermanAligner - Total compute time in native Smith-Waterman : 0.00 sec
18:45:39.174 INFO HaplotypeCaller - Shutting down engine
[May 28, 2025 at 6:45:39 PM GMT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.20 minutes.
Runtime.totalMemory()=700448768
htsjdk.samtools.SAMException: No sequence dictionary mapping available for header: SAMFileHeader{VN=1.4, SO=coordinate}
at htsjdk.samtools.SamFileHeaderMerger.getMergedSequenceIndex(SamFileHeaderMerger.java:741)
at htsjdk.samtools.MergingSamRecordIterator$MergedSequenceDictionaryCoordinateOrderComparator.getReferenceIndex(MergingSamRecordIterator.java:234)
at htsjdk.samtools.MergingSamRecordIterator$MergedSequenceDictionaryCoordinateOrderComparator.fileOrderCompare(MergingSamRecordIterator.java:214)
at htsjdk.samtools.SAMRecordCoordinateComparator.compare(SAMRecordCoordinateComparator.java:48)
at htsjdk.samtools.SAMRecordCoordinateComparator.compare(SAMRecordCoordinateComparator.java:43)
at htsjdk.samtools.ComparableSamRecordIterator.compareTo(ComparableSamRecordIterator.java:75)
at htsjdk.samtools.ComparableSamRecordIterator.compareTo(ComparableSamRecordIterator.java:36)
at java.base/java.util.PriorityQueue.siftUpComparable(PriorityQueue.java:647)
at java.base/java.util.PriorityQueue.siftUp(PriorityQueue.java:639)
at java.base/java.util.PriorityQueue.offer(PriorityQueue.java:330)
at htsjdk.samtools.MergingSamRecordIterator.addIfNotEmpty(MergingSamRecordIterator.java:161)
at htsjdk.samtools.MergingSamRecordIterator.<init>(MergingSamRecordIterator.java:94)
at org.broadinstitute.hellbender.engine.ReadsPathDataSource.prepareIteratorsForTraversal(ReadsPathDataSource.java:429)
at org.broadinstitute.hellbender.engine.ReadsPathDataSource.iterator(ReadsPathDataSource.java:336)
at org.broadinstitute.hellbender.engine.MultiIntervalLocalReadShard.iterator(MultiIntervalLocalReadShard.java:134)
at org.broadinstitute.hellbender.engine.AssemblyRegionIterator.<init>(AssemblyRegionIterator.java:88)
at org.broadinstitute.hellbender.engine.AssemblyRegionWalker.processReadShard(AssemblyRegionWalker.java:188)
at org.broadinstitute.hellbender.engine.AssemblyRegionWalker.traverse(AssemblyRegionWalker.java:173)
at org.broadinstitute.hellbender.engine.GATKTool.doWork(GATKTool.java:1119)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:150)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:203)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:222)
at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:166)
at org.broadinstitute.hellbender.Main.mainEntry(Main.java:209)
at org.broadinstitute.hellbender.Main.main(Main.java:306)
So it seems that (maybe) I need to regenerate new .crai files? Using the suggestion in this post : How to produce a single joint-called VCF file using as input three WGS samples (VCF,CRAM,FASTQ) that the reference files used to generate the .cram files are not the same, I found that the file used for person1 was different from the other two:
more person1.vcf | grep "##ref"
##reference=file:///mnt/ssd/MegaBOLT_scheduler/reference/hg38.fa
more person2.vcf | grep "##ref"
##reference=file:///mnt/ssd/MegaBOLT_scheduler/reference/G42_refdata/GCA_000001405.15_GRCh38_no_alt_plus_hs38d1_analysis_set.fna
more person3.vcf | grep "##ref"
##reference=file:///mnt/ssd/MegaBOLT_scheduler/reference/G42_refdata/GCA_000001405.15_GRCh38_no_alt_plus_hs38d1_analysis_set.fna
I'm looking for suggestions on what the most reasonable straight line path would be. Should I try to generate new .cram/.crai/.vcf files from the source FASTA data or is there an easier way to rebase these downstream analyses on a new reference genome? If so, what is the most widely adopted/robust human reference genome? Any pointers that might help me accomplish this?
Thanks!
Matthias
•
written
by
biostars •
I don't have too much experience running tools like this but my 2c
regenerating the CRAM index (.crai) files is easy and relatively fast (minutes), you can run "samtools index -@8 file.cram" (to use 8 threads) for each of your cram files. unless that error was related to something really weird, I wouldn't suspect this to be the issue necessarily though.
it is not generally good to mix reference genomes, but both "hg38.fa" and "GCA_000001405.15_GRCh38_no_alt_plus_hs38d1_analysis_set.fna" are both "versions" of GRCh38 but subtle differences can really confuse things, particularly with cram (because CRAM files are reference compressed, so you need the exact reference genome to properly decompress/interpret them). you may want to ensure all your data is on the same reference, by potentially performing realignment of e.g. the one that is aligned to "hg38.fa" to "GCA_000001405.15_GRCh38_no_alt_plus_hs38d1_analysis_set.fna" (aligning can take hours to maybe a day). this may be a good thing to do because the error from GATK is during "SAM header merging", while I don't know the exact thing going on in their code, it triggers my spidey sense because info about the reference sequences is in the SAM header (BAM and CRAM files, while being compressed versions of SAM, have "SAM headers") so if the info in these headers are different, potentially it could create that error
GRCh38 is 'good enough' for most purposes. The latest genome is the "T2T" genome aka hs1 aka CHM13-T2T-v2.0, but it is not yet commonly adopted for most analysis pipelines, and the lessons that were learned while sequencing it have already helped the hg38 annotations. It's possible that they are similar enough that it is not a worry. if it does worry you, then you can redo alignment from the fastq to produce your own cram files.
as long as you are doing stuff with trio's I might as well plug this little tutorial I made. I don't walk through the creation of phased variant calls, but once you have them, you can visualize them in jbrowse https://jbrowse.org/jb2/docs/user_guides/analyze_trio/
you could also try bcftools mpileup+whatshap as an alternative pipeline, though again, i don't have much practical experience running these tools so I could be off base
I would ask Nebula to reanalyze person1.vcf for you, since you show it was aligned against a different reference so is not compatible.
If you have a beefy local workstation/cluster you might be able to do the reanalysis yourself but gatk is known to be slow and resource hungry.
I'm not sure if your error relates to a lacking sequence dict. eg Building Dict File for GATK
Also - check to see if the SAMPLE eg person1 2 3 is present in the last header line of the VCF file (header starting with #, not ##. This is important for some later eg trio operations.
I believe I found the correct reference genome for the odd man out. It is hg38 which seems like one of the original older reference genomes. when I have some time, I am going to try redoing the alignment using GATK. I do have a reasonably fast workstation and a functioning Docker image. Any suggestions on effectively accomplishing? This would be most welcome thanks.
You'll need to align the reads using a tool like bwa mem first. Try doing a tutorial to get the concepts or read through at least some of these - https://training.galaxyproject.org/training-material/topics/variant-analysis/ . GATK has good tutorials too.
Log in to answer this question.
