发表日期：2020.02

影响因子：42.778

癌症是全球第二大常见死因，每年超过800万人因癌症丧命。预计在未来十年，癌症发生率将增加50%以上。癌症是体细胞亚克隆自主发展和扩散类疾病的总称。癌症克隆控制多个细胞通路，打破正常细胞的生长和调控等限制，获取自主发展和扩散的特征。单个细胞通路改变不足以引发癌症。每个癌症由潜在的致病异常“池”中的多个异常通路组合而引发。

肿瘤异质性来自于达尔文进化的随机性。达尔文进化的三个先决条件：（1）群体中的特征是变化的；（2）变异从亲本遗传到子代；（3）群体为了生存进行竞争。一部分突变改变细胞表型，一部分突变使克隆获取逃逸正常生理控制的优势。提供选择优势的突变称为驱动突变，反之称为乘客突变。

选用2834个患者人全基因组测序数据（WGS），去除176个患者低质量数据，共计2658个患者的WGS数据，其中有2583个患者高质量数据。2658个患者共取2605个原发肿瘤和173个转移或复发肿瘤，正常样本平均测序深度为39×，肿瘤测序深度分别为38×和60×。研究群体包括1469男性（55%）和1189女性（45%），平均年龄56岁，覆盖38种肿瘤类型。其中，1222个患者具有RNA-seq数据。

利用以上数据分析somatic SNVs, somatic Indels, somatic CNVs, somatic SVs,体细胞逆转录事件，线粒体DNA突变、端粒长度以及germline SNV, Indel, SVs等事件。

利用3个核心变异检测流程和额外10个变异检测流程，对63对tumor-normal变异检测，估测3个核心流程的敏感度和精确度。并对其中50对进行高深度靶向测序验证。3个核心流程检测到真实变异的敏感度为80~90%，每个流程检测的95%以上变异是真实的somatic mutations。针对Indel检测，3个核心流程的敏感度是40~50%，精确度是70~95%。SV检测算法的精确度在80~95%。

对3个核心流程的变异结果合并，评估合并集合中突变的属性：Somatic SNVs敏感度为95%（90%置信区间，88~98%），精确度为95%（90%置信区间，71~99%）。Somatic Indels 检测敏感度为60%（34~72%）和精确度91%（73~96%）。合并的Somatic SVs 敏感度为90%，精确度为97.5%。多种方法检测变异提高了低频突变检出的准确性。

分析2583个患者数据，共检测到43,778,859个somatic SNVs，410,123个somatic 多核酸突变，2,418,247个somatic Indels，288,416个somatic SVs，19,166 体细胞逆转录事件，8,185个新线粒体突变。通过相关性分析，发现诊断年龄和体细胞突变数量相关：年龄每增长一年，增加约190个SNVs，约22个Indels。

3.1癌症驱动突变全景图

根据突变的性质和已知癌症相关基因，预测肿瘤的驱动基因；利用已知的启动子和增强子分析非编码驱动突变。结果发现，91%的肿瘤至少有1个驱动突变，每个肿瘤平均有4.6个驱动突变（癌种之间变化较大）。对于编码区点突变，每个肿瘤平均有2.6个驱动突变。除此之外， 13%（785/5913）的驱动点突变是非编码突变，而且1/3（237/785）突变发生在 TERT 启动子上；25%肿瘤具有非编码驱动突变。说明：非编码区驱动点突变频率较编码区低；与 TERT 启动子相比，其他启动子和增强子并不常发生驱动突变。

根据肿瘤类型，SVs和点突变致力于不同的癌症发生机制。驱动SVs常发生在乳腺癌和卵巢腺癌；驱动点突变常出现在在结肠腺癌和成熟B细胞淋巴瘤。

文章发现抑癌基因的驱动突变多为二次打击事件。例如，954个肿瘤具有 TP53 突变，736（77%）个肿瘤样本的两个等位基因均发生突变，其中96%（707/736）是等位基因突变和等位基因缺失同时发生。17%的病人在癌症相关基因上具有稀少的胚系蛋白截断体突变，5.4%病人由于somatic mutations导致以上基因次等位基因失活。

3.2没有驱动突变的PCAWG肿瘤数据分析

90%以上的PCAWG样本鉴定到驱动突变，仍有181个样本未检测到驱动突变。分析肿瘤样本未找到驱动突变的原因，有以下几点：（1）样本质量低：4/181个样本的正常对照被肿瘤DNA污染，每个对照含有超过5%的肿瘤DNA；同理，肿瘤样本中肿瘤细胞含量较低也会影响突变检出；（2）驱动突变位点覆盖度较低无法满足突变检出：6个肝细胞癌和2个胆管癌在高深度靶向测序后检测到 TERT 突变；（3）生信分析方法：35个骨髓增生性肿瘤未检测到 JAK2 V617F 突变，由于利用Panels of normals作为对照去除测序影响导致。2~5%的健康人群具有造血克隆，可能涵盖了驱动突变；（4）驱动基因检测力不足，说明某些肿瘤中存在未被发现的基因富集；（5）染色体变异：19/43肾细胞癌和18/81前列腺癌缺少驱动突变，但发生染色体异常，有可能单凭染色体扩增或缺失足以引发癌症。

3.3成簇突变和SVs模式

癌症中，单个灾难性事件可产生多个聚集性突变，导致基因组大量重组。主要包含：（1）染色体重排：不同染色体的DNA双链断裂修复导致重排发生；（2）Kataegis（雷雨）：单链DNA局部超突变，导致聚集性核苷酸替换；（3）染色体碎裂：数十数百个DNA断裂同时发生在一个或者几个染色体，产生的碎片随机组合在一起。

467个样本（17.8%）发生染色体重排和平衡易位，主要发生在前列腺癌、淋巴系统恶性肿瘤和甲状腺癌。重排事件导致甲状腺癌的部分融合基因的产生，例如 RET 、 NTRK3 和 IGF2BP3 等等。

60.5%癌症中发生Kataegis事件，例如肺鳞癌、膀胱癌、肢端黑色素瘤和肉瘤等。Kataegis主要包含（1）由APOBEC活性导致TpC的C>N 突变；（2）聚合酶导致 T pT或Cp T 的T >N突变。81.7%的Kataegis事件与 APOBEC3B 表达水平相关，5.7%与易错聚合酶相关，以及2.3%事件是GpC 或 CpC的胞嘧啶脱氨导致的。Kataegis事件与SV断点相关，尤其是缺失和复杂重排事件，包括在缺失附近10-25kb内Cp T pT的T>N 突变。

Kataegis事件包含4种局部超突变类型：（1）脱靶体细胞超突变和局部Cp T pT的T>N 突变；（2）与复杂重排相关的APOBEC；（3）后随链和早期复制区域的APOBEC；（4）后两种类型混合。

587（22.3%）个染色体碎裂样本，主要为肉瘤、脑胶质瘤、肺鳞癌、黑色素瘤和乳腺癌样本。染色体碎裂伴随全基因组重复，相关的驱动基因为 TP53 。肉瘤和B细胞淋巴瘤患者中，女性发生染色体碎裂的频率高于男性；前列腺患者中，晚期患者具有更高频率的染色体碎裂。染色体碎裂区域包含3.6%驱动基因和7%拷贝数驱动。

3.4进化中时间聚集性突变

根据分子时钟分析每个肿瘤的进化史：主克隆发生在早期，亚克隆突变发生在后期；拷贝数扩增区域，分子时间根据突变发生在拷贝之前或者之后进行划分。染色体碎裂通常发生在主克隆，特别是在脂肪肉瘤、前列腺癌和肺鳞癌说明是癌症进化早期事件。在黑色素瘤中，染色体碎裂扩增涉及到较多的癌症相关基因，例如 CCND1 ,  TERT ,  CDKN2A ,  TP53 和 MYC 。

在扩增的染色体碎裂事件中，利用SNV的拷贝数目计算扩增发生的时间，SNV发生在扩增之前，将会有很高比例的reads携带SNVs。相反，SNV发生在拷贝数变异之后，将只有一条染色体携带SNV，具有较低的变异频率。肢端黑色素瘤的 CCND1 扩增区域具有较少的高频突变，而皮肤黑色素瘤更多突变发生在扩增之前。

3.5胚系突变对somatic mutations的影响

根据检测到的胚系突变分析胚系突变对体细胞突变率和模式的影响作用。利用欧洲群体中MAF>5%的胚系突变位点进行GWAS分析，发现 APOBEC3B 突变机制可以利用22q13.1预测，信号最强位点是rs12628403。该位点标记了常见的30kb胚系 APOBE3B 编码序列缺失和 APOBEC3B 的3’非翻译区域 APOBE3A 编码序列融合。除此，文章在22q13.1位置发现一个新的突变位点rs2142833，并验证其与 APOBEC3B 突变相关性。rs12628403和 rs2142833在欧洲群体中是独立遗传的，rs2142833是 APOBEC3B 的eQTL。

利用稀有突变（MAF<0.5%）分析欧洲群体中胚系蛋白截短体（PTVs）和体细胞DNA重排相关性。胚系BRCA2和BRCA1蛋白截短体和小于10kb的体细胞缺失和串联重复负荷相关。BRCA1蛋白截短体和模板插入具有显著相关。20/21个BRCA1相关肿瘤出现模板插入表型，且胚系突变和体细胞突变均发生在该基因上。说明 BRCA1 基因的次等位基因失活驱动模板插入SV表型。

稀有突变关联分析发现胚系MBD4蛋白截短体突变增加CpG位置的体细胞C>T突变。 MBD4 编码DNA修复基因，移除甲基化CpG上的T:G错配的胸腺嘧啶。

评估LINE调控体细胞反转座子事件，验证114个胚系LINE对体细胞反转座激活能力，包含70个人类基因组相关插入和53个连锁不平衡SNP。16个L1元件介导67%（2440/3669）的转座事件，以两种形式进行体细胞激活，称为Strombolian和Plinian；Strombolian在人群中分布频率较高，引发中小规模的体细胞L1激活；Plinian在群体中频率很低，引发严重的体细胞L1激活。

3.6复制的永生

癌症特征之一是逃避细胞衰老，保持端粒长度是癌症永久复制的因素之一。16%的肿瘤在 ATRX ,  DAXX 和 TERT 基因上发生突变。聚类端粒序列的12个特征得到4个肿瘤亚型，说明 ALT 和 TERT 介导的端粒变异的不同。

体细胞驱动突变在四个亚型中分布不同。C1主要富集 RB1 突变和影响 ATRX 的SV，C2主要富集 ATRX 和 DAXX 的体细胞点突变，C3样本主要发生 TERT 启动子突变。 RB 基因缺失与端粒延长相关。高频发生端粒异常机制的肿瘤主要由于组织中低复制活性。

总结

利用泛癌全基因组测序数据对驱动突变、结构变异、克隆进化以及转座子事件和端粒模式进行详细分析，绘制泛癌基因组特征和阐明引发癌症的多样性因素。

参考文献

ICGC/TCGA Pan-Cancer Analysis of Whole Genomes Consortium. Pan-cancer analysis of whole genomes. Nature. 2020, 578(7793): 82-93.

原文链接：https://www.nature.com/articles/s41586-020-1969-6

随着癌症基因组学的进步，突变注释格式（MAF）被广泛接受并用于存储检测到的体细胞变体。 癌症基因组图谱项目对30多种不同的癌症进行了测序，每种癌症类型的样本量超过200种。由体细胞变体组成的结果数据以MAF格式形式存储。 只要数据采用MAF格式，该软件包就会尝试从TCGA源或任何内部研究中有效地汇总，分析，注释和可视化MAF文件.

使用前要先将文件转换为maf格式,对于VCF格式文件，可以使用 vcf2maf 进行格式转换.

maf文件包含的内容：

注： 安装过程特别麻烦，按了好几天，R版本要求3.3以上，也不要使用最新版本，可能有的包新版本还没同步 。我使用的是：

annovarToMaf函数说明

Converts variant annotations from Annovar into a basic MAF.将annovar格式转换为maf格式

| 参数 |详细解释 |

| annovar | input annovar annotation file.|

| Center | Center field in MAF file will be filled with this value. Default NA.（MAF文件中的中心字段将填充此值。 默认NA）|

| refBuild | NCBI_Build field in MAF file will be filled with this value. Default hg19.（MAF文件中的NCBI_Build字段将填充此值。 默认hg19）|

| tsbCol | column name containing Tumor_Sample_Barcode or sample names in input file.（列名包含Tumor_Sample_Barcode或输入文件中的示例名称） |

| table | reference table used for gene-based annotations. Can be 'ensGene' or 'refGene'. Default 'refGene'（用于基于基因的注释的参考表。 可以是'ensGene'或'refGene'。 默认'refGene)|

| basename | If provided writes resulting MAF file to an output file. （将结果MAF文件写入输出文件）|

| sep | field seperator for input file. Default tab seperated.|

| MAFobj | If TRUE, returns results as an [MAF](http://127.0.0.1:37698/help/library/maftools/help/MAF object.|

| sampleAnno | annotations associated with each sample/Tumor_Sample_Barcode in input annovar file. If provided it will be included in MAF object. Could be a text file or a data.frame. Ideally annotation would contain clinical data, survival information and other necessary features associated with samples. Default NULL.（与输入annovar文件中的每个样本/ Tumor_Sample_Barcode相关联的注释。 如果提供，它将包含在MAF对象中。 可以是文本文件或data.frame。 理想情况下，注释将包含临床数据，生存信息和与样本相关的其他必要特征。 默认为NULL）|

然后用linux处理掉那些无义突变，也可以在后续设置参数去掉无义突变

Takes tab delimited MAF (can be plain text or gz compressed) file as an input and summarizes it in various ways. Also creates oncomatrix - helpful for visualization.

该文件将每个样本中的变体数显示为堆积条形图，将变体类型显示为Variant_Classification汇总的箱形图。 我们可以在堆积的条形图中添加平均线或中线，以显示整个群组中变体的平均值/中值数

Plots maf summary.

Oncoplot函数使用“ComplexHeatmap”来绘制oncoplots2。 具体来说，oncoplot是ComplexHeatmap的OncoPrint功能的包装器，几乎没有任何修改和自动化，使绘图更容易。 侧面条形图和顶部条形图可分别由drawRowBar和drawColBar参数控制。

top的值需要视情况而定

takes output generated by read.maf and draws an oncoplot

通过包括与样本相关的注释（临床特征），改变变体分类的颜色并包括显着性的q值（从MutSig或类似程序生成），可以进一步改善Oncoplots。

[图片上传失败...(image-fc6334-1536734778754)]

使用oncostrip函数可视化任何一组基因，它们在每个样本中绘制类似于cBioPortal上的OncoPrinter工具的突变。 oncostrip可用于使用top或gene参数绘制任意数量的基因

titv函数将SNP分类为 Transitions_vs_Transversions ，并以各种方式返回汇总表的列表。 汇总数据也可以显示为一个箱线图，显示六种不同转换的总体分布，并作为堆积条形图显示每个样本中的转换比例。

takes output generated by read.maf and classifies Single Nucleotide Variants into Transitions and Transversions.

棒棒糖图是简单且最有效的方式，显示蛋白质结构上的突变点。许多致癌基因具有比任何其他基因座更频繁突变的优先位点。这些斑点被认为是突变热点，棒棒糖图可以用于显示它们以及其他突变。我们可以使用函数lollipopPlot绘制这样的图。这个功能要求我们在maf文件中有氨基酸改变信息。然而，MAF文件没有关于命名氨基酸变化字段的明确指南，不同的研究具有不同的氨基酸变化的字段（或列）名称。默认情况下，lollipopPlot查找列AAChange，如果在MAF文件中找不到它，则会打印所有可用字段并显示警告消息。对于以下示例，MAF文件包含字段/列名称“Protein_Change”下的氨基酸变化。我们将使用参数AACol手动指定它。此函数还将绘图作为ggplot对象返回，如果需要，用户稍后可以修改该对象。

maftools还可以制作很多图，比如

还可以用函数geneCloud绘制突变基因的词云图。 每个基因的大小与其突变/改变的样品总数成比例。

癌症中的许多引起疾病的基因共同发生或在其突变模式中显示出强烈的排他性。 可以使用somaticInteractions函数检测这种相互排斥或共同发生的基因组，其执行成对的Fisher's Exact检验以检测这种显着的基因对。 somaticInteractions函数还使用cometExactTest来识别涉及>2个基因的潜在改变的基因集

maftools包 功能很强大，具体可参考：

http://bioconductor.org/packages/release/bioc/vignettes/maftools/inst/doc/maftools.html

《深度学习精要（基于R语言）》学习笔记

机器学习主要用于开发和使用那些从原始数据中学习、总结出来的用于进行预测的算法。

深度学习是一种强大的多层架构，可以用于模式识别、信号检测以及分类或预测等多个领域。

神经网络包括一系列的神经元，或者叫作节点，它们彼此连结并处理输入。神经元之间的连结经过加权处理，权重取决于从数据中学习、总结出的使用函数。一组神经元的激活和权重（从数据中自适应地学习）可以提供给其他的神经元，其中一些最终神经元的激活就是预测。

经常选择的激活函数是sigmoid函数以及双曲正切函数tanh，因为径向基函数是有效的函数逼近，所以有时也会用到它们。

权重是从每个隐藏单元到每个输出的路径，对第i个的输出通过(w_i)表示。如创建隐藏层的权重，这些权重也是从数据中学习得到的。分类会经常使用一种最终变换，softmax函数。线性回归经常使用恒等（identity）函数，它返回输入值。权重必须从数据中学习得到，权重为零或接近零基本上等同于放弃不必要的关系。

R中神经网络相关包：

一旦集群完成初始化，可以使用R或本地主机（127.0.0.1:54321）提供的Web接口与它连接。

如果数据集已经加载到R，使用as.h2o()函数：

如果数据没有载入R，可以直接导入到h2o中：

也可以直接导入网络上的文件：

导入基于图片识别手写体数字，数据集的每一列（即特征），表示图像的一个像素。每张图像都经过标准化处理，转化成同样的大小，所以所有图像的像素个数都相同。第一列包含真实的数据标签，其余各列是黑暗像素的值，它用于分类。

使用caret包训练模型：

生成数据的一组预测，查看柱状图：

跟训练集数据柱状图对比，很明显模型不是最优的。

通过混淆矩阵检查模型性能：

No Information Rate（无信息率）指不考虑任何信息而仅仅通过猜测来决定最频繁的类的准确度期望。在情形“1”中，它在11.16%的时间中发生。P值（P-Value [Acc >NIR]）检验了观测准确度（Accuracy : 0.3674）是否显著不同于无信息率（11.16%）。

Class: 0的灵敏度（Sensitivity）可以解释为：89.07%的数字0被正确地预测为0。特异度（Specificity）可以解释为：95.14%的预测为非数字0被预测为不是数字0。

检出率（Detection Rate）是真阳性的百分比，而最后的检出预防度（detection prevalence）是预测为阳性的实例比例，不管它们是否真的为阳性。

平衡准确度（balanced accuracy）是灵敏度和特异度的平均值。

接下来我们通过增加神经元的个数来提升模型的性能，其代价是模型的复杂性会显著增加：

隐藏神经元的数量从5个增加到10个，样本内性能的总准确度从36.74% 提升到了 65.4%。我们继续增加隐藏神经元的数量：

增加到40个神经元后准确度跟10个神经元的一样，还是65.4%。如果是商业问题，还需要继续调节神经元的数量和衰变率。但是作为学习，模型对数字9的表现比较差，对其他数字都还行。

RSNNS包提供了使用斯图加特神经网络仿真器（Stuttgart Neural Network Simulator ， SNNS）模型的接口，但是，对基本的、单隐藏层的、前馈的神经网络，我们可以使用mlp()这个更为方便的封装函数，它的名称表示多层感知器（multi-layer perceptron）。

RSNNS包要求输入为矩阵、响应变量为一个哑变量的 矩阵 ，因此每个可能的类表示成矩阵列中的 0/1 编码。

通过decodeClassLabels()函数可以很方便的将数据转换为哑变量矩阵。

预测结果的值为1-10，但是实际值为0-9，所以在生成混淆矩阵时，需要先减去1：

RSNNS包的学习算法使用了相同数目的隐藏神经元，计算结果的性能却有极大提高。

函数I()有两个作用：

1.在对data.frame的调用中将对象包含在I()中来保护它，防止字符向量到factor的转换和名称的删除，并确保矩阵作为单列插入。

2.在formula函数中，它被用来禁止将“+”、“-”、“*”和“^”等运算符解释为公式运算符，因此它们被用作算术运算符。

从RSNNS包返回的预测值（pred.ml4）中可以看到，一个观测可能有40%的概率成为“5”，20%的概率成为“6”，等等。最简单的方法就是基于高预测概率来对观测进行分类。RSNNS包有一种称为赢者通吃（winner takes all，WTA）的方法，只要没有关系就选择概率最高的类，最高的概率高于用户定义的阈值（这个阈值可以是0），而其他类的预测概率都低于最大值减去另一个用户定义的阈值，否则观测的分类就不明了。如果这两个阈值都是0（缺省），那么最大值必然存在并且唯一。这种方法的优点是它提供了某种质量控制。

但是在实际应用中，比如一个医学背景下，我们收集了病人的多种生物指标和基因信息，用来分类确定他们是否健康，是否有患癌症的风险，是否有患心脏病的风险，即使有40%的患癌概率也需要病人进一步做检查，即便他健康的概率是60%。RSNNS包中还提供一种分类方法称为“402040”，如果一个值高于用户定义的阈值，而所有的其他值低于用户定义的另一个阈值。如果多个值都高于第一个阈值，或者任何值都不低于第二个阈值，我们就把观测定性为未知的。这样做的目的是再次给出了某种质量控制。

“0”分类表示未知的预测。

通常来说，过拟合指模型在训练集上的性能优于测试集。过拟合发生在模型正好拟合了训练数据的噪声部分的时候。因为考虑了噪声，它似乎更准确，但一个数据集和下一个数据集的噪声不同，这种准确度不能运用于除了训练数据之外的任何数据 — 它没有一般化。

使用RSNNS模型对样本外数据预测：

模型在第一个5000行上的准确度为85.1%，在第二个5000行上的准确度减少为80%，损失超过5%，换句话说，使用训练数据来评价模型性能导致了过度乐观的准确度估计，过度估计是5%。

这个问题我们后面再处理。