癌症基因的somatic mutation calling 流程的评价体系
文章是：A comprehensive assessment of somatic mutation detection in cancer using whole-genome sequencing
癌症研究通常都是获取tumor和normal组织一起进行WGS或者WES测序，这样能过滤掉germline mutation，达到获取真正属于癌症的somatic mutation部分。当然，仅仅是这样也是不够的，癌症是一个持续进化的过程，从而会累计各种各样的突变，而只有那些真正导致癌症的突变才是科学家们想关心的，那部分突变称为driver mutation。
- 从call mutation的角度来说，早在七八年前就有个各种各样的benchmark测试啦，这么多年的科研路上，留下来的就是GATK,varscan，freebayes啦。
- 从call somatic mutation的角度来说，也是三五年前就出现了不少benchmark测试，经过了时间的考验的就是TCGA项目采用的mutect,muse,strella,varscan,somaticSnipper.
- 只剩下call driver mutation这一块处女地了，应该可以在这一块做一些软件测评，发一些文章的。
本文选取了chronic lymphocytic leukaemia和medulloblastoma这两种癌症的WGS数据来分析，CLL和MB的基因组特性很大区别，具有代表性。CLL 有着一些已知基因组倒位以及其它大的结构变异，而MB的基因组有着很高的四倍体基因组特性，却并没有多少结构变异。对CLL我们主要比较MiSeq和IonTorrent测序仪的表现，主要看看它们在有限的测序深度条件下假阴性率的高低。而对MB，我们产生了近300X的高测序深度，再加上它的1, 4, 8, 15 ,17号染色体的四倍体特性，我们可以探究其低频率变异位点。
虽然理论上癌症基因组研究是需要对癌症样品以及其配对的正常组织一起测序，测序深度最好是超过30X，而且是150个碱基长度的双端测序，才能更好的找到肿瘤特异性的somatic mutation。但在实际操作过程中，总会有这样那样的条件无法满足。导致后续的生物信息学分析困难重重，包括somatic single-base mutations (SSM), somatic insertion/deletion mutations (SIM)，larger structural changes 的确定，简称SNV,INDEL,CNV,SV 这4种常见的变异类型。为了更好的解释数据分析的结果，有必要进行系统性的评价测序数据的产出对分析结果的影响。本文就首先对下面四点进行分步探讨。
- 文库构建
- 肿瘤测序深度的选择
- 肿瘤纯度的选择
- 配对的normal样本的测序深度
文库构建
这一块知识是大部分生物信息学分析着的短板，因为大家都是直接拿到了NGS测序数据，最多就是对其做一个简单的QC及过滤而已。很少有人会愿意追溯其具体产出的细节。本文对MB样品构建了8种不同的文库，在样品起始量，打断片段大小, 片段选择，文库试剂盒，是否PCR，选择的测序仪，测序深度这些条件都做了调整。具体指标如下所示：
| Library | Starting DNA (μg) | Fragment Size (bp) | Size selection | Library protocol | PCR cycles | Sequencing machine | Chemistry (Illumina) | Depth (×) control:tumour |
|---|---|---|---|---|---|---|---|---|
| L.A | 4 | ∼400 | 2% Agarose gel | KapaBio | 0 | HiSeq 2500 HiSeq 2000MiSeq | V1 (RR)V3V2 | 29.6 : 40.5 |
| L.B | 1 | ∼400 | 2% Agarose gel, Invitrogen E-gel | TrueSeq DNA | 10 | HiSeq 2000 | V3 | 44.9 : 62.8 |
| L.C | 2.5 | ∼500 | 2% Agarose gel | NEBNext | 12 | HiSeq 2500HiSeq 2000 | V1 (RR)V3 | 58.9 : 66.8 |
| L.D | 1 | ∼550 | Agarose gel | TrueSeq DNA | 10 | HiSeq 2000 | V3 | 35.3 : 39.1 |
| L.E | 2.8 | ∼620 | 1.5% Agarose gel pippin | NEBNext | 0 | HiSeq 2000 | V3 | 40.5 : 60.4 |
| L.F | 1 | ∼400 | AMPureXP beads | NEBDNA | 10 | HiSeq 2000 | V3 | 38.7 : 37.9 |
| L.G | 1 | ∼350 | AMPureXP beads | TrueSeq DNA PCR-Free | 0 | HiSeq 2000 | V3 | 19.4 : 19.3 |
| L.H | 0.5 | ∼175 | AMPureXP beads | SureSelect WGS | 10 | HiSeq 2500 | V3 | 28.7 : 26.5 |
因为L.G和L.H文库测序深度低于28，所以在分析测序深度的评价过程中剔除了这两个文库的数据，把其它文库的数据都用picard的downsample功能挑选了28X的数据进行比较。
可以看出，文库A, E，G 这3个PCR-free的都受GC%含量的影响较小，是比较好的选择。
肿瘤测序深度的选择
把上面的8个文库数据合并，得到了314X的肿瘤测序数据和272X的正常组织测序数据。再随机把它们抽样得到250, 200,150, 100, 50, 30，20X的测序深度数据。再对这些数据分别进行snp-calling ，得到的结果如下：
可以很清楚的看到，大部分的突变(SSM)在30X的时候就已经被找出来了。
但是随着测序深度的提高，能找到的SSM也逐步增多。但是100X之后增加幅度很低，近乎达到饱和，能检测到95%的变异了。所以推荐100X为比较理想的测序深度，可以平衡测序花费，和得到有意义的生物学结果。
肿瘤纯度的选择
主流的测序仍然是bulk测序，没办法达到单细胞水平的测序，必然是百万级别的细胞一起处理，这样就几乎不可能得到100%纯度的肿瘤样品，而且肿瘤细胞其实也处于数量不少的免疫细胞的包围之中。这些正常细胞的含量会对肿瘤测序找somatic mutation产生干扰。
所以本文特意模拟了各种纯度的肿瘤样品，在肿瘤测序数据里面掺入 17，33，50%的正常组织测序数据并重新找变异。发现如果正常组织含量过高，哪怕是把测序深度提高到250X也会漏掉很多真正的变异位点。但总的来说，提高测序深度是可以增加变异的检出率的，具体指标如下：
- At 100 ×, the detected proportions of mutation calls from the pure sample were 95%, 90% and 85%, respectively, for 17%, 33% and 50% ‘contamination’
- At 30 ×, only 92%, 83% or 68% of the calls from the 30 × pure sample were called when adding 17%, 33% or 50% normal reads, respectively
配对的normal样本的测序深度
做肿瘤测序研究很容易陷入一个误区，就是以为自己的重点是肿瘤样品，之所以测正常组织只是为了过滤germline mutation，所以不需要给它太高的测序深度，标准的30X就已经足够了。所以作者在这里选取了250X的肿瘤测序数据，加上不同测序深度的正常样品测序结果，来分别进行somatic mutation的calling。正常样品的测序深度包括，200, 150, 100, 50,30。
结果发现，在高深度的正常组织测序数据的分析过程中得到的SSM要少于低深度的，分析那些少掉的变异位点，发现有显著的GpTpG情况，可能是测序过程中的人工错误，说明提高正常样本的测序深度，可以减少假阳性。
金标准
经过这一系列的测评，作者给出了一个经过专家认证的somatic mutation集合的金标准数据集放在文章附件供下载。这里选取的是314X的tumor和272X的normal测序数据，经过DKFZ calling pipeline得到的somatic mutation数据集，并且把这些位点分成了5个级别，如下：
- Tier 1 mutations have a mutant AF≥10%
- Tier 2 mutations have AF≥5%
- Tier 3 includes all verified mutations supported by unambiguous alignments
- Tier 4 includes additional mutations with more complicated or ambiguous local alignments
- Tier 5 includes those with unusually high or low read depth.
这个金标准里面的不同级别的SSM和SIM的数量简单的统计如下：
| Definition | MB SSM | MB SIM | |
|---|---|---|---|
| Class 1 | Mutant AF≥0.10 | 962 | 337 |
| Class 2 | 0.05≤Mutant AF<0.10 | 139 | |
| Class 3 | Mutant AF<0.05 | 154 | |
| Class 4 | Ambiguous alignment | 8 | 10 |
| Class 5 | High or low depth | 29 | |
| Tier 1 | Class 1 | 962 | 337 |
| Tier 2 | Classes 1 and 2 | 1,101 | |
| Tier 3 | Classes 1, 2 and 3 | 1,255 | |
| Tier 4 | Classes 1, 2, 3 and 4 | 1,263 | 347 |
| Tier 5 | Classes 1, 2, 3, 4 and 5 | 1,292 |
值得注意的是，这套数据集里面的1620个比较好的位点中32%的位点在重复区域，9%位于串联重复区域，还有4.4%位于串联重复区域的附件。接近四分之一的位点的AF低于10%，这些位点是在数据分析过程中很容易被漏检。而SIMs里面，83%都是在串联重复区域。
流程测评
上面的测评用的都是DKFZ calling pipeline，但是其中数据分析的流程也需要进行测评！所以作者针对CLL和MB两套数据集都用了16个不同的找SSM和SIM的流程，这些流程来自于世界各地的科研实验室自己用的最佳流程。每个流程的各个步骤所用的软件，参数都在本文的超百页的PDF附件里面列出。比如：
MB.O
FASTQ PROCESSING
Lane Quality control:Each lane is checked with our quality control pipeline. All lanes were kept to create the sample bamfile.
MAPPING AND BAM PROCESSING
Alignment:
Core mapping with BWA v0.7.5a
bwa aln -l 32 -t 6 -q 20 (6 CPUs used)
bwa sampe -P
Presently 'bwa sampe' is called within a Perl wrapper and pairs with both reads unmapped arewritten to a separate BAM file.
Bam processing:
Sorting + indexing lane BAM file (Picard tools v1.73) Indel realignment at known sites (GATK2 v2.0-23) Duplicates Marking (Picard) Recalibration (GATK2): Calculating table and generating recalibrated BAM file Picard MergeSamFile Picard Markduplicates (with remove duplicates option) Picard and GATK are used with standard parameters.In step picard MergeSamFile, we merged all the lanes that passed the quality controls. We then removed the duplicates to obtain a BAM file at the sample level.The final estimated coverage were : 28x for the normal sample and 38x for the tumor sample
MUTATION CALLING
SSM
We called somatic mutations with MuTect (version 1.1.4, Broad tool). This tool compares normal and tumor samplesto retrieve somatic mutations.Default values for all parameters were used.
SIM
We used SomaticIndelDetector, a tool included in GATK2 to call small indels.
MUTATION FILTRATION
SSM
Among MuTect calls, 40064 candidate calls were retained. These candidate mutations were then filtered according tothese filters: TAC : minimum number of read carrying mutant allele DP : read depth + base quality filter
MQ : mapping quality of reads carrying mutant allele PB : read position filter SB : strand bias filter GL : germline filter
RE : mutations found nearby repeat regions BL : mutant allele fount in satellite regions CU2 : the presence of the mutant allele in a panel of 18 other germline samples sequenced with Illumina protocol
SIM
Calls were further filtered by our pipeline: min number of reads carrying indel mapping quality of reads carrying indel germline filtering
strand bias filtering mismatches in mapping of reads carrying indel (number and windows around indel filtering) base quality around indel the position of indels within boundaries of known repeated, centromeric and telomeric regions mutation filtering (indel starts at a SNV position) position of indels in the carrying reads the presence of the indel in a panel of 20 other germline samples sequenced with Illumina protocol
COMPUTATIONAL RESOURCES
Almost the same as described in CLL.O1.
把流程拆解开来就是测序数据fastq的预处理，与参考基因组比对，比对文件的处理，变异的检出，somatic mutation的检出，变异的过滤。
既然做了这么复杂的评价，那么结果的解读也必然非常难以描述，作者别出心裁的创造了这个rainfall图，如下：
作者发现，只有205个SSMs和仅1个SIM被所有的流程都发现了，所有的数据如下：
| SSM calls | Aligner | SSM Detection Software | TP | FP | FN | P | R | F1 |
|---|---|---|---|---|---|---|---|---|
| MB.GOLD | BWA, GEM | Curated | 1,255 (8) | 0 | 0 | 1.00 | 1.00 | 1.00 |
| MB.A | BWA | In-house | 775 (0) | 147 | 480 | 0.84 | 0.62 | 0.71 |
| MB.B | BWA | samtools, Varscan | 788 (1) | 12 | 467 | 0.99 | 0.63 | 0.77 |
| MB.C | GEM | samtools, bcftools | 766 (3) | 1,025 | 489 | 0.43 | 0.61 | 0.50 |
| MB.D | n.a. | SMuFin | 737 (4) | 1,086 | 518 | 0.41 | 0.59 | 0.48 |
| MB.E | BWA | SomaticSniper | 750 (4) | 229 | 505 | 0.77 | 0.60 | 0.67 |
| MB.F | BWA | Strelka | 884 (2) | 165 | 371 | 0.84 | 0.70 | 0.77 |
| MB.G | BWA | Caveman, Picnic | 899 (3) | 140 | 356 | 0.87 | 0.72 | 0.78 |
| MB.H | Novoalign | MuTect | 947 (3) | 6,296 | 308 | 0.13 | 0.76 | 0.22 |
| MB.I | BWA | samtools | 879 (7) | 129 | 376 | 0.87 | 0.70 | 0.78 |
| MB.J | None, BWA | SGA+freebayes | 856 (1) | 62 | 399 | 0.93 | 0.68 | 0.79 |
| MB.K | BWA | Atlas2-snp | 945 (8) | 7,923 | 310 | 0.11 | 0.75 | 0.19 |
| MB.L1 | BWA | MuTect, Strelka | 385 (0) | 3 | 870 | 0.99 | 0.31 | 0.47 |
| MB.L2 | BWA | MuTect, Strelka | 900 (1) | 253 | 355 | 0.78 | 0.72 | 0.75 |
| MB.M | BWA mem | samtools, GATK+MuTect | 937 (4) | 1,695 | 318 | 0.36 | 0.75 | 0.48 |
| MB.N | BWA | Strelka | 847 (1) | 289 | 408 | 0.75 | 0.68 | 0.71 |
| MB.O | BWA | MuTect | 944 (3) | 272 | 311 | 0.78 | 0.75 | 0.76 |
| MB.P | BWA | Sidron | 833 (3) | 256 | 422 | 0.77 | 0.66 | 0.71 |
| MB.Q | BWA | qSNP+GATK | 842 (2) | 25 | 413 | 0.97 | 0.67 | 0.79 |
上面的每一个流程所找到的somatic mutation都跟作者得到的近标准进行比较，计算出了各种指标如上表，需要理解一些基础知识。
- TP代表真实值是阳性且预测出来也是阳性的数量。
- FN代表真实值是阳性但是预测出来是阴性的数量。
- FP代表真实值是阴性但是预测出来是阳性的数量。
- P这里是precision，就是准确率，公式是(TP)/(TP+FP)
- R这里是recall，就是召回率，公式是(TP)/(TP+FN)
- F这个才是衡量每个流程的标准，公式是 2 × (Precision × Recall)/(Precision+Recall).
做这样的 统计也非常简单，就是把每个流程最后得到的somatic mutation list的位点拿去跟金标准比较，看看它们的overlap情况即可。
可以看到表现最好的是MB.Q和MB.J这两套流程，它们有着最高的F1值，用得是SGA+freebayes和qSNP+GATK这两套流程。得到的突变数量越多，并不意味着召回率就很高。因为大部分的流程的召回率都低于80%，因为在40X的测序深度的限制下，必然有很多AF低的位点是无法检出的。
同理，对CLL数据集也进行了同样的分析。
因为这些流程用到了近百款软件，有些流程之间的差异仅仅是比对工具的不一样，比如4个主流的比对工具，Novoalign2, BWA, BWA-mem 和 GEM，控制这一个变量的不同，就可以对之进行探究，看看它们到底有哪些不同的表现。还有3种不同的参考基因组b37, b37+decoy and ‘hg19r’以及基因组不同区域的特性，这里作者并没有对最新版参考基因组进行比较。还有TCGA项目纳入的4个somatic mutation caller，包括MuTect，muse，varscan，somaticsnipper，还有同样比较流行的Strelka。
拿到高质量的somatic mutation之后做什么呢？
跟TCGA的样本一起比较TNB，我博客有教程，下载TCGA所有癌症的maf文件计算TMB分析mutation signature，我博客有教程，下载TCGA所有癌症的maf文件做signature分析
最后的结论是：
PCR-free建库方法，测序深度100X左右，N-T测序深度一致，是最佳选择。