假阳性突变的出现居然是因为duplicates mark的不够？

生信菜鸟团

发布于 2022-04-08 17:28:50

1K0

文章被收录于专栏：生信菜鸟团生信菜鸟团

介绍

大家好，我是sssimon yang，菜鸟团的新成员。我目前的研究方向是基因组学，欢迎跟我一起，挖掘基因组里的宝藏。

做基因组分析，每天就是和突变打交道。拿到突变文件后，我习惯在bam里实际的看一看突变，好有个印象，对突变的质量有一个整体的把握。但是不看不知道，一看，就总是吓一跳。

问题的出现

样本走gatk best practice流程拿到somatic vcf后，查看关注区域，显示tumor样本中，7号染色体138235828处有一个C到A的SNV，AD为11。

看一下bam。

$ hg19=~/project/0-reference/2-GATK/hg19/human_g1k_v37_decoy.fasta
$ samtools tview -p 7:138235828 tumor.bam ${hg19}

往下翻，找支持突变的reads。

总共看下来，只有这6条reads支持该突变，vcf里写着11，现在就只有6个了，真行。而且这个右端，上面三个一致，下面三个一致，有点像duplicates。

左端也是一致的，reads上也没有其他突变，那就是duplicates了，duplicates被去除后这不就只剩2，更离谱了。

把duplicates从bam当中去除，再看一下。

bam的处理我更喜欢用sambamba，速度快。

$ sambamba view -f bam -F 'not duplicate' tumor.bam 7:138235828 > view.bam
$ sambamba index view.bam
$ samtools tview -p 7:138235828 view.bam ${hg19}

没区别啊，显然，这6个reads并没有被标记成duplicates，怪不得这个突变被call出来了。

gatk best practice是包含mark duplicatse步骤的，这一步我应该做了.

$ sambamba view -f bam -F 'not duplicate' -c tumor.bam 7:138235828-138235828
126
$ sambamba view -f bam -F '' -c tumor.bam 7:138235828-138235828
272

-c输出reads数目。前后的数量对比能看出，确实有一半的duplicates已被去除，是做了啊。

那为什么这些reads没有被标记成duplicates呢？

问题的解决

找到GATK MarkDuplicates (Picard)[1]的文档，扫了一下，发现了重点。

“The program can take either coordinate-sorted or query-sorted inputs, however the behavior is slightly different. When the input is coordinate-sorted, unmapped mates of mapped records and supplementary/secondary alignments are not marked as duplicates. However, when the input is query-sorted (actually query-grouped), then unmapped mates and secondary/supplementary reads are not excluded from the duplication test and can be marked as duplicate reads.

在coordinate-sorted的时候，unmapped mate read和supplementary read是不会被标记为duplicates的。emm？这也行？

考虑到在我的处理中bam一直都是坐标排序的，也就是coordinate-sorted。那极有可能这些支持的reads都是supplementary read，没有被去除。

检查一下。

$ sambamba view -F 'supplementary' -f bam view.bam > supplementary.bam
$ sambamba index supplementary.bam
$ samtools tview -p 7:138235828 supplementary.bam ${hg19}

果然，6条支持突变的reads都是supplemtary read。

那有没有办法解决这个问题呢，有，gatk说，query-sorted可以，好，试一下。

为了方便查看结果，就直接将检测到的duplicates去除了。

$ sambamba sort -t 24 -n tumor.bam -o query_sorted.bam
$ gatk MarkDuplicates -I query_sorted.bam -O gatk_markduplicates.bam -M temp.metrics --REMOVE_DUPLICATES true
$ sambamba sort -t 24 gatk_markduplicates.bam -o gatk_markduplicates_sorted.bam
$ sambamba index gatk_markduplicates_sorted.bam
$ samtools tview -p 7:138235828 gatk_markduplicates_sorted.bam ${hg19}

成功。gatk能够处理supplemtary的duplicates问题，但是需要query-sorted，是有些麻烦了。

其他方案

samtools[2]也提供了remove duplicates的方法，下面进行一下测试。

阅读文档后，我发现samtools的markdup提供了一个-S选项，可以在primary read被标记为duplicates时，其supplementary read也被标记为duplicates。

那么问题是，我所遇到的情况中，primary read有没有被标记为duplicates。

挑一个出来验证一下。

$ sambamba view supplementary.bam 7:138235828-138235828 | tail -1
A00583:225:H2FV5DSXY:1:2528:10926:19727 2131    7       138235810       60      71H79M  =       138235610       -279    TACATACCGGTTACGGCAACTCCTTCGTGCAAGGTGTGATGCATCCCCAGTGACCAACAACACCATCCAATTTCACTGT 99<99:>7;;;996;;9:9=;==::5:9:;99:9898:789:68:::978796:965865859857986577679>:;9 SA:Z:7,138235650,-,79M71S,60,0; MD:Z:18C60      PG:Z:MarkDuplicates     RG:Z:T191162_PT_S2012-18117_R19058217   NM:i:1  AS:i:74XS:i:20

MD信息显示18C，即138235810+18有个C，有突变，是我们要的read。找到该read name的所有reads。

$ sambamba view ~/project/2-comb/3-combine/1-bam/2-tumor/T191162_PT_S2012-18117_R19058217.bam 7 | grep A00583:225:H2FV5DSXY:1:2528:10926:19727
A00583:225:H2FV5DSXY:1:2528:10926:19727 1187    7       138235610       60      119M31S =       138235650       119     GGTTAAATAGTGCAGCTGAGTTACGCCTCAAAAGTTTTTAACTTGTTTTGGAAAGAAAATTGTGTTGTATGTCCAACAAACCTATATATGATTGTTTATTGTGTATTGATGTACATACCGGTTACGGCAACTCCTTCGTGCAAGGTGTGA 476745645948769768795658388797888:6778879988:7888:;799<::::89;8;89;888;8;:9::9:9998777777:97796776679696667875855765577295657297677687679369878::7:797   SA:Z:7,138235810,+,111S39M,60,1;        MD:Z:119        PG:Z:MarkDuplicates     RG:Z:T191162_PT_S2012-18117_R19058217  NM:i:0   AS:i:119        XS:i:20
A00583:225:H2FV5DSXY:1:2528:10926:19727 1107    7       138235650       60      79M71S  =       138235610       -119    ACTTGTTTTGGAAAGAAAATTGTGTTGTATGTCCAACAAACCTATATATGATTGTTTATTGTGTATTGATGTACATACCGGTTACGGCAACTCCTTCGTGCAAGGTGTGATGCATCCCCAGTGACCAACAACACCATCCAATTTCACTGT 6;989959899777977765787887866789;:87:887;;7777778;8:9:::88::;:;99;:<9:;99<99:>7;;;996;;9:9=;==::5:9:;99:9898:789:68:::978796:965865859857986577679>:;9   SA:Z:7,138235810,-,71S79M,60,1; MD:Z:79 PG:Z:MarkDuplicates     RG:Z:T191162_PT_S2012-18117_R19058217   NM:i:0  AS:i:79XS:i:0
A00583:225:H2FV5DSXY:1:2528:10926:19727 2211    7       138235810       60      111H39M =       138235650       -83     TACATACCGGTTACGGCAACTCCTTCGTGCAAGGTGTGA 45776366294657297677687779369878::696:9 SA:Z:7,138235610,+,119M31S,60,0;        MD:Z:18C20      PG:Z:MarkDuplicates     RG:Z:T191162_PT_S2012-18117_R19058217   NM:i:1  AS:i:34 XS:i:0
A00583:225:H2FV5DSXY:1:2528:10926:19727 2131    7       138235810       60      71H79M  =       138235610       -279    TACATACCGGTTACGGCAACTCCTTCGTGCAAGGTGTGATGCATCCCCAGTGACCAACAACACCATCCAATTTCACTGT 99<99:>7;;;996;;9:9=;==::5:9:;99:9898:789:68:::978796:965865859857986577679>:;9 SA:Z:7,138235650,-,79M71S,60,0; MD:Z:18C60      PG:Z:MarkDuplicates     RG:Z:T191162_PT_S2012-18117_R19058217   NM:i:1  AS:i:74XS:i:20

有四条reads。查询数字对应SAM Flag[3]的含义，1187=1+2+32+128+1024，是duplicates，1107=1+2+16+64+1024，是duplicates。

所以duplicated supplementary read是会对应到duplicated primary read的，从supplementary read的产生机制上也可以理解。

-S选项对于解决问题应该是有效的，麻烦之处就是markdup要求处理的bam必须先行经过samtools fixmate -m处理，照着他的example试一下。

$ samtools fixmate -@ 24 -m query_sorted.bam fixmate.bam
$ samtools sort -@ 24 -o positionsort.bam fixmate.bam
$ samtools markdup -@ 24 -r positionsort.bam samtools_markdup.bam
$ sambamba index -t 24 samtools_markdup.bam
$ samtools markdup -@ 24 -r -S positionsort.bam samtools_markdup_S.bam
$ sambamba index -t 24 samtools_markdup_S.bam

在不加-S选项时和原来一样，使用-S选项后，成功。

fixmate步骤是必须的，直接运行会出现以下错误。

$ samtools markdup tumor.bam markdup.bam
[markdup] error: no MC tag. Please run samtools fixmate on file first.
[markdup] error: unable to assign pair hash key.

作者的建议是，应当在比对完成后立即进行fixmate步骤，然后sort，这样可以最大限度的减少额外处理步骤。

以下方案均已尝试：

sambamba markdup 作者说标准与picard一致失败。
samtools rmdup 该子命令已废弃失败。
samtools rmdup -s 在单端数据上去除duplicates 失败。
samtools rmdup -S 将双端数据视作单端数据处理成功。

samtools rmdup在处理时检测条件较为简单，只考虑其位置，这也是被废弃的原因。在使用-S时，supplementary read被视作单端read处理，因此成功。

总结

总的来说，我推荐使用samtools markdup配合-S选项进行处理。值得一提的是，并不只有一个突变出现这种情况，因此，对supplementary read的额外考虑是十分必要的。

参考资料

[1]GATK MarkDuplicates (Picard): https://gatk.broadinstitute.org/hc/en-us/articles/360057439771-MarkDuplicates-Picard-

[2]samtools: http://www.htslib.org/doc/samtools-markdup.html

[3]SAM Flag: https://www.samformat.info/sam-format-flag

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2022-03-03，如有侵权请联系 cloudcommunity@tencent.com 删除

python

本文分享自生信菜鸟团微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

python

登录后参与评论

0 条评论

热度