报错代码

05:28:18.297 INFO  BaseRecalibrator - Shutting down engine
[June 16, 2022 5:28:18 AM CST] org.broadinstitute.hellbender.tools.walkers.bqsr.BaseRecalibrator done. Elapsed time: 711.85 minutes.
Runtime.totalMemory()=7918321664
htsjdk.samtools.SAMException: Query asks for data past end of contig. Query contig chrUn_KI270509v1 start:2248 stop:2396 contigLength:2318
        at org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile.getSubsequenceAt(CachingIndexedFastaSequenceFile.java:330)
        at org.broadinstitute.hellbender.engine.ReferenceFileSource.queryAndPrefetch(ReferenceFileSource.java:78)
        at org.broadinstitute.hellbender.utils.recalibration.BaseRecalibrationEngine.calculateIsSNPOrIndel(BaseRecalibrationEngine.java:372)
        at org.broadinstitute.hellbender.utils.recalibration.BaseRecalibrationEngine.processRead(BaseRecalibrationEngine.java:131)
        at org.broadinstitute.hellbender.tools.walkers.bqsr.BaseRecalibrator.apply(BaseRecalibrator.java:191)
        at org.broadinstitute.hellbender.engine.ReadWalker.lambda$traverse$0(ReadWalker.java:96)
        at java.util.stream.ForEachOps$ForEachOp$OfRef.accept(ForEachOps.java:184)
        at java.util.stream.ReferencePipeline$3$1.accept(ReferencePipeline.java:193)
        at java.util.stream.ReferencePipeline$2$1.accept(ReferencePipeline.java:175)
        at java.util.stream.ReferencePipeline$3$1.accept(ReferencePipeline.java:193)
        at java.util.Iterator.forEachRemaining(Iterator.java:116)
        at java.util.Spliterators$IteratorSpliterator.forEachRemaining(Spliterators.java:1801)
        at java.util.stream.AbstractPipeline.copyInto(AbstractPipeline.java:482)
        at java.util.stream.AbstractPipeline.wrapAndCopyInto(AbstractPipeline.java:472)
        at java.util.stream.ForEachOps$ForEachOp.evaluateSequential(ForEachOps.java:151)
        at java.util.stream.ForEachOps$ForEachOp$OfRef.evaluateSequential(ForEachOps.java:174)
        at java.util.stream.AbstractPipeline.evaluate(AbstractPipeline.java:234)
        at java.util.stream.ReferencePipeline.forEach(ReferencePipeline.java:418)
        at org.broadinstitute.hellbender.engine.ReadWalker.traverse(ReadWalker.java:94)
        at org.broadinstitute.hellbender.engine.GATKTool.doWork(GATKTool.java:1049)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:140)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:192)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:211)
        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)

real    711m57.046s
user    539m17.576s
sys     21m29.580s

这其实就是说有一个read比对到contig上的时候末尾超出去了，那只要简单粗暴删掉这对reads就不会再报错啦~

寻找报错行

将排序重比对后的bam文件转sam，根据第3列染色体名称和第4列起始位置找到导致报错的行

samtools view ￥{file}.bam | \
awk 'BEGIN {FS="\t"} $3=="chrUn_KI270509v1" && $4=="2248" {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8}'

sam分隔符是tab而非awk默认的空格，所以如果不设置BEGIN {FS="\t"}的话，第一行将不能正确识别

执行结果

ST-E00200:183:HLJJHCCXX:7:1111:21816:17799      147     chrUn_KI270509v1        2248    19      29M1I120M       =       1970

第一列是read名称，第二列是flag（147=128【read2】+16【比对到负链】+2【正常比对】+1【双端测序】），第三列是染色体名称，第四列是起始位置，第五列是MAPQ值，第六列是CIGAR字符串，第七列是双端测序的对应read比对到的染色体名称（在这里=意思是比对到同一条染色体），第八列是双端测序的对应read比对的起始位置

我的59个样本里有7个出现这个报错，并且报错的read全都是比对到负链的，也不知道为什么会这样。。。

但是总之我们可以通过第一列read名称找到这一行并删除

删除报错行

#bam转sam先，一定要记得加-h，不然会报错没有SQ信息，默认是不显示header的！！！
samtools view -@ 6 -h ${file} > to_del_error_line.sam
#删除含有这个read名称的行
sed -i "/ST-E00200:183:HLJJHCCXX:7:1111:21816:17799/d" to_del_error_line.sam
#sam转回bam，就可以重新跑BaseRecalibrator啦
samtools view -@ 6 -bS -o del_error_line.bam to_del_error_line.sam

也可能后面还有一行有这个问题，重复上述操作就可以