GATK CollectAlignmentSummaryMetrics等工具无输出报错求助

问题描述 #

在变异检测流程中，尝试运行GATK的CollectAlignmentSummaryMetrics和CollectInsertSizeMetrics工具以获取比对及插入片段长度指标，但未生成输出文件且出现报错。

使用命令 #

gatk CollectAlignmentSummaryMetrics -I /home/osh/Documents/aligned_reads/M2/S2_sorted_dedup_bqsr.bam -R /home/osh/Documents/data/ref/resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta -O /home/osh/Documents/aligned_reads/M2/Alignment_metrics_S2.txt

gatk CollectInsertSizeMetrics -I /home/osh/Documents/aligned_reads/M2/S2_sorted_dedup_bqsr.bam -O /home/osh/Documents/aligned_reads/M2/S2_insert_size_metrics.txt -H /home/osh/Documents/bam_reads/insert_size_histogram.pdf

报错信息 #

使用GATK jar包路径：/home/osh/Documents/gatk-4.4.0.0/gatk-package-4.4.0.0-local.jar
执行命令：
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 /home/osh/Documents/gatk-4.4.0.0/gatk-package-4.4.0.0-local.jar CollectAlignmentSummaryMetrics -I /home/osh/Documents/aligned_reads/M2/S2_sorted_dedup_bqsr.bam -R /home/osh/Documents/data/ref/resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta -O /home/osh/Documents/aligned_reads/M2/Alignment_metrics_S2.txt
22:58:27.695 INFO NativeLibraryLoader - 从jar包加载libgkl_compression.so：file:/home/osh/Documents/gatk-4.4.0.0/gatk-package-4.4.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
[2024年1月16日 周二 22:58:27 IST] CollectAlignmentSummaryMetrics 参数：--INPUT /home/osh/Documents/aligned_reads/M2/S2_sorted_dedup_bqsr.bam --OUTPUT /home/osh/Documents/aligned_reads/M2/Alignment_metrics_S2.txt --REFERENCE_SEQUENCE /home/oshi/Documents/data/ref/resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta --MAX_INSERT_SIZE 100000 --EXPECTED_PAIR_ORIENTATIONS FR --ADAPTER_SEQUENCE AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT --ADAPTER_SEQUENCE AGATCGGAAGAGCTCGTATGCCGTCTTCTGCTTG --ADAPTER_SEQUENCE AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT --ADAPTER_SEQUENCE AGATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCTCGTATGCCGTCTTCTGCTTG --ADAPTER_SEQUENCE AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT --ADAPTER_SEQUENCE AGATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNNNATCTCGTATGCCGTCTTCTGCTTG --METRIC_ACCUMULATION_LEVEL ALL_READS --IS_BISULFITE_SEQUENCED false --COLLECT_ALIGNMENT_INFORMATION true --ASSUME_SORTED true --STOP_AFTER 0 --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false
[2024年1月16日 周二 22:58:27 IST] 执行用户：oshi@oshi-VirtualBox，系统：Linux 5.15.0-91-generic amd64；OpenJDK 64-Bit Server VM 17.0.9+9-Ubuntu-120.04；压缩器：Intel；解压器：Intel；GCS可用；Picard版本：4.4.0.0
[2024年1月16日 周二 22:58:29 IST] picard.analysis.CollectAlignmentSummaryMetrics 执行完成。耗时：0.03分钟。
Runtime.totalMemory()=70488064
错误核心：
htsjdk.samtools.util.SequenceUtil$SequenceListsDifferException: 序列字典大小不一致（26, 3366）
at htsjdk.samtools.util.SequenceUtil.assertSequenceListsEqual(SequenceUtil.java:259)
at htsjdk.samtools.util.SequenceUtil.assertSequenceDictionariesEqual(SequenceUtil.java:342)
at htsjdk.samtools.util.SequenceUtil.assertSequenceDictionariesEqual(SequenceUtil.java:328)
at picard.analysis.SinglePassSamProgram.makeItSo(SinglePassSamProgram.java:117)
at picard.analysis.SinglePassSamProgram.doWork(SinglePassSamProgram.java:94)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:289)
at org.broadinstitute.hellbender.cmdline.PicardCommandLineProgramExecutor.instanceMain(PicardCommandLineProgramExecutor.java:37)
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)

问题原因 #

报错核心是序列字典大小不匹配：你的BAM文件的序列字典仅包含26条记录（对应标准的22条常染色体+X/Y/MT等主要染色体），但指定的参考基因组fasta的序列字典包含3366条记录（包含大量contig或备选contig）。GATK在默认的严格验证模式下，要求BAM文件与参考基因组的序列字典完全一致，否则会直接终止运行。

解决步骤 #

• 使用与BAM匹配的参考基因组：找到生成该BAM文件时使用的同版本参考基因组fasta，替换命令中的-R参数路径，确保两者序列字典完全一致。这是最稳妥的解决方案。
• 临时降低验证严格性：如果确认BAM与参考基因组的核心染色体序列一致，只是contig数量不同，可以在命令中添加--VALIDATION_STRINGENCY LENIENT参数跳过严格校验。例如：

  gatk CollectAlignmentSummaryMetrics -I /home/osh/Documents/aligned_reads/M2/S2_sorted_dedup_bqsr.bam -R /home/osh/Documents/data/ref/resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta -O /home/osh/Documents/aligned_reads/M2/Alignment_metrics_S2.txt --VALIDATION_STRINGENCY LENIENT
  

  注意：该方法可能忽略其他潜在的不匹配问题，仅作为临时应急方案使用。
• 重新生成匹配的BAM文件：如果必须使用当前的完整参考基因组，需要重新用该参考基因组比对原始reads，生成新的BAM文件后再执行指标收集操作。

内容的提问来源于stack exchange，提问作者user22752471