最近安排团队一起学习肿瘤外显子流程，处理的数据来源于2017年CELL杂志的文章《Tumor and Microenvironment Evolution during Immunotherapy with Nivolumab》，跑完上游流程拿到了maf格式的肿瘤特异性的somatic突变位点后，第一张图就应该是突变全景图，所以有如下代码：

```R
library(maftools)
library(data.table)

vep_maf = fread("7.annotation/vep/vep_merge.maf")
oncoplot( read.maf(vep_maf) )
```

起初看到下面的图表，团队成员是一脸懵逼的，因为跟原文差异太大了，这些top基因，比如TTN和MUC16其实是因为它本身基因长度以及复杂度等特性导致它容易被算作是突变位点。


如果看文章仔细一点，就会发现，首先它并没有显示全部的病人队列，仅仅是，免疫治疗前的基线数据哦。而且呢，其实作者是挑选了指定的基因进行展示突变情况，并不是简单的看突变的频率来挑选排名靠前的基因。

原文主图如下所示：


所以，我们需要有两个操作，首先，对病人进行筛选，其次，对基因进行筛选。自己复制粘贴上面图片里面的基因名字，然后写入到 gene.list 文本文件里面，代码如下所示：

```r
maf <- read.maf(vep_maf)
gene = read.table("gene.list")[,1]
barcode = unique(vep_maf$Tumor_Sample_Barcode)
pre.barcode = barcode[str_detect(barcode,'pre')]
pre.maf =subsetMaf(maf,tsb = pre.barcode)
oncoplot(
pre.maf,
genes = gene,
showTumorSampleBarcodes = T,
GeneOrderSort = gene,
keepGeneOrder = T,
removeNonMutated = F
)
```

再次出图如下：


可以看到， BRAF，NRAS, NF1基因的突变情况基本上跟原位是一模一样的啦。可以从侧面说明我们的肿瘤外显子流程是OK的。

### 本文的研究方法解读

- 病人和样本：85 名（实际上符合要求的只有 68 位）黑色素瘤患者，每位患者在治疗前取了基线样本，治疗后重复取样活检，其中 35 名是 IPi-P，即接受 IPi 治疗，33名是 Ipi-N

- 测序：

wes：试剂盒是 Agilent SureSelect All Exon V2，测序仪是 HiSeq 2000 V3 或 HiSeq 2500，双端 76 bp，测序深度 150x，

- wes 数据处理：
- 采用 4 个软件 call somatic SNV：MuTect 1.1.4、SomaticSniper 1.0.4、VarScan 2.3.7 和 Strelka 1.0.13，并且过滤掉等位基因 reads 小于 5 或相应正常覆盖率 reads 小于 7 。用 VarScan 2.3.7 和 Strelka 1.0.13 call INDEL。TMB 高低的阈值为 100 个突变。
- 使用 deconstructSigs 包进行mutation signature 分析。
- 新抗原分析是先用 SOAP-HLA 预测 HLA 表型，然后用 MAC 进行校正（因为是黑色素瘤，存在紫外线照射引起的 SNV），用 NetMHC v4.0 版本预测突变肽与 HLA 的结合力，保留 level < 2 的突变肽，如果一个 突变位点产生了多个合格的 9肽，这认为太突变产生一种新抗原。
- 拷贝数变异数据分析用 FACETS 。
- 克隆分析用 Pyclone ，用到了 facet 评估的肿瘤纯度，得到每个 SNV 的 CCF ，并且根据给出的 CI 置信区间划分克隆和亚克隆，阈值为 95%。（clone 分析的方法值得参考，不过我在 pyclone 的结果中，没有发现这个置信区间 CI ）

本文是MSKCC那边的团队的成果，非常值得学习哦！