简介
一般来说，肿瘤体细胞突变的分析都要求需要肿瘤与正常配对样本，采用如 MuTect2、MuSE、Varscan2、SomaticSniper、Strelka2 之类的工具来 call 体细胞突变。 对于得到的体细胞突变位点，以 vcf 文件的形式保存，需要进一步过滤，突变过滤主要有以下几种策略：
- 基于阈值：比如过滤掉 reads counts < 3，VAF < 0.05 等
- 基于数据库：比如过滤掉 1000G、dbSNP、ExAC、gnomAD 等数据库突变人群频率 > 0.001
- 基于功能：过滤掉同义突变位点，内含子位点等
上面的几种方法，可以在对 vcf 进行注释转 maf 之后，根据 maf 文件对应的注释列进行过滤，比较简单，这里介绍其他的方法。
多个体细胞突变检测工具取交集
很多文章会使用多个软件 call 突变，然后进行合并，常用见的方法就是多个软件得到的结果中，一个突变在两个或两个以上的软件检测到即保留，认为其可信度较高，这是在较多文章中出来的方法，也算是一种过滤策略。不过在这之后，还可能对保留的突变位点再进一步的过滤。
但是不同的突变检测工具，得到的 vcf 文件并不完全一致，因为不同工具定义的 format tag 不太一致。所以直接进行合并可能导致格式有问题。这里简单简介一种比较粗暴的 vcf 文件合并方法。主要是针对体细胞突变检测工具 Mutect2 和 Strelka2 的结果。
如果想对vcf文件有一个更加深入的了解，可以查看官方帮助文档：https://samtools.github.io/hts-specs/VCFv4.2.pdf
Mutect2 的 vcf 文件
Mutect2 的体细胞突变检测方法在GATK 的Somatic Mutation流程--肿瘤基因组测序数据分析专栏 已经有详细介绍了。如果不看 vcf 的 header 的话，主要突变信息的内容是：
代码语言：javascript复制$ less -S 6.mutect/case1_biorep_B_filter.vcf
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT case1_biorep_B case1_germline
chr1 1228610 . G T . PASS AS_FilterStatus=SITE;AS_SB_TABLE=49,22|1,1;DP=74;ECNT=2;GERMQ=93;MBQ=38,20;MFRL=157,100;MMQ=60,60;MPOS=26;NALOD=1.43;NLOD=7.48;POPAF=6.00;TLOD=5.80 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:43,2:0.049:45:26,2:17,0:27,16,1,1 0/0:28,0:0.036:28:15,0:12,0:22,6,0,0
chr1 3476004 . T C . PASS AS_FilterStatus=SITE;AS_SB_TABLE=58,28|2,1;DP=93;ECNT=2;GERMQ=93;MBQ=32,20;MFRL=145,91;MMQ=60,60;MPOS=15;NALOD=1.41;NLOD=7.52;POPAF=6.00;TLOD=4.49 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:57,3:0.056:60:28,3:29,0:39,18,2,1 0/0:29,0:0.037:29:17,0:12,0:19,10,0>
chr1 3635071 . G T . PASS AS_FilterStatus=SITE;AS_SB_TABLE=80,113|1,3;DP=206;ECNT=1;GERMQ=93;MBQ=32,33;MFRL=160,133;MMQ=60,60;MPOS=27;NALOD=1.78;NLOD=17.70;POPAF=6.00;TLOD=6.21 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:130,4:0.039:134:69,2:55,2:50,80,1,3 0/0:63,0:0.016:63:31,0:30,0:30,33,0>
chr1 3781221 . G T . PASS AS_FilterStatus=SITE;AS_SB_TABLE=107,58|4,0;DP=178;ECNT=1;GERMQ=93;MBQ=38,33;MFRL=142,154;MMQ=60,60;MPOS=27;NALOD=1.75;NLOD=16.55;POPAF=6.00;TLOD=5.96 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:104,4:0.048:108:46,2:57,2:74,30,4,0 0/0:61,0:0.017:61:31,0:30,0:33,28,0>
chr1 3848765 . G T . PASS AS_FilterStatus=SITE;AS_SB_TABLE=32,157|1,3;DP=205;ECNT=1;GERMQ=93;MBQ=37,28;MFRL=170,140;MMQ=60,60;MPOS=35;NALOD=1.80;NLOD=17.76;POPAF=6.00;TLOD=5.34 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:126,4:0.033:130:58,3:66,1:26,100,1,3 0/0:63,0:0.016:63:37,0:26,0>
chr1 6630273 . G T . PASS AS_FilterStatus=SITE;AS_SB_TABLE=131,44|2,2;DP=193;ECNT=1;GERMQ=93;MBQ=38,20;MFRL=144,108;MMQ=60,60;MPOS=17;NALOD=1.76;NLOD=16.55;POPAF=6.00;TLOD=3.81 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:115,4:0.028:119:68,0:46,4:83,32,2,2 0/0:60,0:0.017:60:25,0:35,0:48,12,0>
chr1 10324117 . G T . PASS AS_FilterStatus=SITE;AS_SB_TABLE=5,117|0,3;DP=131;ECNT=1;GERMQ=93;MBQ=39,33;MFRL=165,171;MMQ=60,60;MPOS=11;NALOD=1.59;NLOD=11.14;POPAF=6.00;TLOD=3.86 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:85,3:0.045:88:45,2:39,1:5,80,0,3 0/0:37,0:0.025:37:21,0:16,0>
chr1 13779345 . G T . PASS AS_FilterStatus=SITE;AS_SB_TABLE=67,9|2,0;DP=83;ECNT=1;GERMQ=93;MBQ=29,40;MFRL=173,164;MMQ=60,60;MPOS=28;NALOD=1.43;NLOD=7.77;POPAF=6.00;TLOD=4.01 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:50,2:0.056:52:26,1:24,1:44,6,2,0 0/0:26,0:0.036:26:11,0:14,0>
chr1 15487348 . G T . PASS AS_FilterStatus=SITE;AS_SB_TABLE=88,78|3,0;DP=181;ECNT=1;GERMQ=93;MBQ=30,39;MFRL=148,150;MMQ=60,60;MPOS=24;NALOD=1.87;NLOD=21.01;POPAF=6.00;TLOD=6.22 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:88,3:0.046:91:50,3:38,0:50,38,3,0 0/0:78,0:0.013:78:39,0:38,0>
......
Strelka2 的 vcf 文件
Strelka2 的体细胞突变检测方法在Strelka2 call Somatic 流程--肿瘤基因组测序数据分析专栏 已经有详细介绍了。其体细胞突变检测的结果分为 SNV 和 INDELs 两个 vcf 文件，需要先进行合并，这里为了后面方便操作，就不保留 vcf 文件的 header 信息：
代码语言：javascript复制$ zcat 6.strelka/case1_biorep_B/results/variants/*gz | grep '^chr' | grep 'PASS' |sort -V >6.strelka/case1_biorep_B/results/variants/case1_biorep_B_filter.vcf
$ zless -S 6.strelka/case1_biorep_B/results/variants/case1_biorep_B_filter.vcf
chr1 6264722 . G A . PASS SOMATIC;QSS=41;TQSS=1;NT=ref;QSS_NT=41;TQSS_NT=1;SGT=GG->AG;DP=107;MQ=60.00;MQ0=0;ReadPosRankSum=-0.76;SNVSB=0.00;SomaticEVS=8.56 DP:FDP:SDP:SUBDP:AU:CU:GU:TU 39:0:0:0:0,0:0,0:39,39:0,0 68:0:0:0:2,2:0,0:64,64:2,2
chr1 6620092 . C A . PASS SOMATIC;QSS=29;TQSS=1;NT=ref;QSS_NT=29;TQSS_NT=1;SGT=CC->AC;DP=38;MQ=60.00;MQ0=0;ReadPosRankSum=0.55;SNVSB=0.00;SomaticEVS=8.49 DP:FDP:SDP:SUBDP:AU:CU:GU:TU 17:0:0:0:0,0:17,17:0,0:0,0 21:0:0:0:2,2:18,18:0,0:1,1
chr1 19152288 . G A . PASS SOMATIC;QSS=120;TQSS=1;NT=ref;QSS_NT=120;TQSS_NT=1;SGT=GG->AG;DP=274;MQ=60.00;MQ0=0;ReadPosRankSum=-0.91;SNVSB=0.00;SomaticEVS=19.98 DP:FDP:SDP:SUBDP:AU:CU:GU:TU 57:0:0:0:0,0:0,0:57,57:0,0 216:0:0:0:57,58:0,0:159,159:0,0
chr1 20485269 . G A . PASS SOMATIC;QSS=46;TQSS=1;NT=ref;QSS_NT=46;TQSS_NT=1;SGT=GG->AG;DP=101;MQ=60.00;MQ0=0;ReadPosRankSum=1.72;SNVSB=0.00;SomaticEVS=12.70 DP:FDP:SDP:SUBDP:AU:CU:GU:TU 54:0:0:0:0,0:0,0:54,54:0,0 47:0:0:0:2,2:0,0:43,43:2,2
chr1 23050336 . G T . PASS SOMATIC;QSS=36;TQSS=1;NT=ref;QSS_NT=36;TQSS_NT=1;SGT=GG->GT;DP=228;MQ=60.00;MQ0=0;ReadPosRankSum=0.85;SNVSB=0.00;SomaticEVS=7.71 DP:FDP:SDP:SUBDP:AU:CU:GU:TU 81:0:0:0:0,0:0,0:81,81:0,0 147:0:0:0:0,0:0,0:143,143:4,4
chr1 31425849 . G T . PASS SOMATIC;QSS=35;TQSS=1;NT=ref;QSS_NT=35;TQSS_NT=1;SGT=GG->GT;DP=153;MQ=60.00;MQ0=0;ReadPosRankSum=1.20;SNVSB=0.00;SomaticEVS=7.15 DP:FDP:SDP:SUBDP:AU:CU:GU:TU 45:0:0:0:0,0:0,0:45,45:0,0 108:0:0:0:0,0:0,0:104,104:4,4
chr1 31586610 . C G . PASS SOMATIC;QSS=33;TQSS=1;NT=ref;QSS_NT=33;TQSS_NT=1;SGT=CC->CG;DP=70;MQ=60.00;MQ0=0;ReadPosRankSum=0.69;SNVSB=0.00;SomaticEVS=8.28 DP:FDP:SDP:SUBDP:AU:CU:GU:TU 27:0:0:0:0,0:27,27:0,0:0,0 43:0:0:0:0,0:40,40:2,2:1,1
chr1 43569766 . C G . PASS SOMATIC;QSS=84;TQSS=1;NT=ref;QSS_NT=84;TQSS_NT=1;SGT=CC->CG;DP=112;MQ=60.00;MQ0=0;ReadPosRankSum=-0.43;SNVSB=0.00;SomaticEVS=19.98 DP:FDP:SDP:SUBDP:AU:CU:GU:TU 31:0:0:0:0,0:31,31:0,0:0,0 81:0:0:0:0,0:62,62:19,19:0,0
chr1 56696122 . G T . PASS SOMATIC;QSS=49;TQSS=1;NT=ref;QSS_NT=49;TQSS_NT=1;SGT=GG->GT;DP=184;MQ=60.00;MQ0=0;ReadPosRankSum=-0.08;SNVSB=0.00;SomaticEVS=7.30 DP:FDP:SDP:SUBDP:AU:CU:GU:TU 47:0:0:0:0,0:0,0:46,46:1,1 137:0:0:0:0,0:1,1:131,131:5,5
......
合并两个 vcf 文件
如果突变检测只用了 Mutect2 和 Strelka2 ，可以用下面代码进行合并。首先是取出 Mutect2 的 vcf 文件的 header 信息，即最后合并后的 vcf 文件 header 来自于 Mutect2 。然后再用 awk 对两个vcf文件的突变坐标进行判断，如果坐标一致就保留，主体区的其他信息如 info format等，也是来自于 Mutect2：
代码语言：javascript复制cat 6.mutect/case1_biorep_B_filter.vcf | grep '^#' >6.somatic/case1_biorep_B_filter.vcf
awk 'NR==FNR {A[$1"t"$2]; next}(($1"t"$2) in A)' 6.strelka/case1_biorep_B/results/variants/case1_biorep_B_filter.vcf 6.mutect/case1_biorep_B_filter.vcf >>6.somatic/case1_biorep_B_filter.vcf
这样最后得到的 6.somatic/case1_biorep_B_filter.vcf就是 Mutect2 和 Strelka2 的 overlap 了。这种方法比较简单粗暴，实际上也有其他工具可以实现，欢迎大家在评论区提出自己的方法。
基于 igv 的人工筛查
对于体细胞突变的过滤，要在 igv 进行 “人工智能”过滤，已经有文章报道：Standard operating procedure for somatic variant refinement of sequencing data with paired tumor and normal samples https://www.nature.com/articles/s41436-018-0278-z
文章中作者开发了一个工具 IGVNav ，可以结合 igv 进行使用。在文章附件，作者给了许多例子，供读者进行参考。这里从附件中挑出来几个具有代表性的例子：
正常的体细胞突变示例
当变异在肿瘤样本中有足够的支持并且没有明显的测序伪影时，就会进行体细胞调用。在本例中，该突变被假定为真正的体细胞突变。在评估基因组特征部分中的参考序列时，参考等位基因是胞嘧啶 (C)。DNA 肿瘤样本中的比对和覆盖率显示，大约 20% 的reads支持变异腺嘌呤 (A) 等位基因（绿色）。重要的是，正常样本中没有 reads 支持该变异，表明该变异是体细胞变异而不是种系多态性。使用基因注释，我们可以预测这种 (C>A) 碱基变化将导致PANK4中的 ATG（M；蛋氨酸）到 ATT（I；异亮氨酸）错义突变。
注意：由于肿瘤样本不纯，体细胞变异通常具有低于 50% 的 VAF。然而，后者不是一个严格的规则，因为随机抽样、拷贝数改变、杂合性丢失和其他因素有时会产生 50% 或以上的体细胞 VAF。
正常的种系突变示例
当变异在正常样本中有足够的reads支持超出了被认为可归因于正常的肿瘤污染时，就会进行种系调用。在这个例子中，该变体被假定为种系多态性。参考等位基因是鸟嘌呤 (G)，但 DNA 正常/肿瘤轨迹中的读数支持胸腺嘧啶 (T) 等位基因。这表明该突变可能是纯合种系多态性。单核苷酸多态性 (SNP) 跟踪进一步支持所讨论的变体是常见的多态性。
单链reads假阳性
在这个例子中，该突变在肿瘤样本中得到了 14 个 reads 的支持，但大多数都在负链上 (93%)。此外，几个支持reads有多个不匹配，表明潜在的低质量reads。将需要更多信息来调用此突变体细胞或失败，否则认为其为假阳性。
突变支持的 reads 出现错配碱基
当包含变体的读数具有其他不匹配的碱基对时，使用 Multiple Mismatches 标签，这会降低读数质量的置信度。具体来说，考虑到高错误率和错误的随机分布，当错误在肿瘤样本中而不是在正常样本中时，可能会出现假变异。
低比对质量
在 igv 中，当reads按readstrand着色时，半透明/透明reads表示比对质量较低，不透明reads表示映射质量较高。低映射质量reads支持的突变则可能为假阳性。
单碱基重复
当在包含单个核苷酸重复序列（例如，AAAAAAA…）的参考序列区域附近发生突变时，在这种情况下，被调用的变异很可能是由reads与参考基因组的错配引起的。
