Identification of prognostic signature of non-small cell lung cancer based on TCGA methylation data
基于TCGA甲基化数据构建预后标志物来预测非小细胞肺癌患者的生存情况
发表期刊：Sci Rep
发表日期：2020 May 22
影响因子：4.01
DOI: 10.1038/s41598-020-65479-y
01
02
普通人群中与癌症相关的死亡中，肺癌是主要原因之一，非小细胞肺癌（NSCLC）占将近83％，但是其预后要好于小细胞肺癌。尽管接受手术治疗的早期NSCLC患者预后良好，但这些患者中仍然有40％会因发生转移容易复发，即使通过包括手术，化学疗法，靶向治疗，免疫治疗等治疗但是这些患者的预后仍然很差。因此，迫切需要找到一个预后标志物对NSCLC进行分类。
DNA甲基化可以表观遗传修饰基因组表达，并长期与NSCLC以及其他几种癌症的发展和进展相关。且大量的文献中的证明使一系列甲基化基因座的进行预后组合，可以作为NSCLC的预后标志物。这项研究基于NSCLC的TCGA甲基化谱，试图构建一个甲基化位点对患者预后有预测的预后模型。
03
1.数据下载
–
TCGA数据集：使用TCGA GDC API从TCGA数据库下载504名NSCLC患者临床随访信息，415名NSCLC患者450k甲基化数据以及551名NSCLC患者RNA-seq数据
2.数据预处理
–
甲基化数据：去除GpG缺失值大于70%的甲基化位点，排除含有交叉反应性和靶向多态性CpG 的探针，使用R软件包KNN补全甲基化表达谱中的缺失值最终得到208022个甲基化位点。
临床数据：去除生存时间少于30天的样本，最终得到349名患者用于后续分析
3.样本分组
–
作者将349随机分为训练集（n = 174）和验证集（n = 175），训练集与验证集样本的年龄分布，临床分期、随访时间和生存情况见表1。
表1 训练集与验证集临床信息统计
4.识别与患者预后相关的甲基化位点
–
作者在训练集中使用单因素Cox比例回归分析识别与患者预后相关的甲基化位点，然后，对单因素显著的甲基化位点和临床变量结合进行多因素cox比例回归分析，认为P<0.05具有统计学意义。
5.层次聚类
作者对甲基化位点进行了无监督的层次聚类，这些甲基化位点与患者的预后显著相关，使用R软件包Consensus Cluster Plus根据欧氏距离计算样本之间的距离，最佳聚类数由累积分布函数确定（CDF；图1A）。本次抽样采用双抽样方案，每次抽样80%，重复100次。图1B显示了当簇数（K）达到6时，结果的稳定性可以达到，并且这选为进一步分析的亚型数目。
图1. 一致性聚类参数选择
5.KEGG富集分析
–
R软件包clusterProfiler进行KEGG富集分析，阈值P <0.05
6.WGCNA共表达分析
–
使用R软件包WGCNA进行了WGCNA共表达分析，使用训练集中，WGCNA共表达算法分析与预后显著的甲基化位点中CpG位点的共表达。使用Pearson相关系数计算各CpG位点之间的距离，并构建加权共表达网络，作者选择的软阈值为7（图2A-B），筛选共表达模块，共表达网络符合无尺度网络，即出现连接度为k的节点的对数log(k)与该节点出现的概率的对数log(P(k))要负相关，且相关系数要大于0.8，然后再将邻接矩阵转换成拓扑矩阵，基于TOM，使用average-linkage层次聚类法对基因进行聚类，按照混合动态剪切树的标准，并设置每个lncRNA网络模块最少的基因数目30，最后对模块进行聚类分析，将具有相似性的模块合并为新模块（height = 0.25, deepSplit = 2, and minModuleSize = 30）。
图2. WGCNA网络共表达分析
04
—
1.识别对预后有影响的甲基化位点
——–
在训练集中对所有的甲基化位点进行单因素回归比例分析，获得了9201个对预后有显著影响的甲基化位点（P <0.05），对患者临床信息（年龄、性别、T分期、N分期、stage）进行单因素回归分析，只有T分期和stage与预后显著相关（P分别为0.0002197和0.005511）。进一步，使用T分期、stage分期和甲基化位点进行多因素回归分析，发现有4000个甲基化位点与NSCLC患者的OS显著相关。
2.甲基化位点层次聚类分析
———-
作者对多因素回归分析后4000个与预后显著相关甲基化位点使用一致性聚类将训练集样本聚成六个亚型（图3A），通过热图可以看出，6个亚型样本中甲基化程度有不同的差异，且大多数甲基化位点属于低甲基化位点。
图3. 训练集样本聚类分析
然后，作者对这六个亚型患者进行预后分析，结果显示6亚型患者之间的生存状态存在显著差异（p = 4e-5，图4A），其中第四亚型预后最好，第六亚型预后最差，最后，作者比较了六个分子亚型在T，N，stage等临床信息上存在的差异（图4B），结果显示第六亚型患者倾向于癌症晚期，肿瘤较大，淋巴结转移较多。
图4. 六个亚型患者之间的生存分析与临床信息差异分析
3.功能富集分析
———-
作者将4000个甲基化位点定位到3482个基因上。KEGG功能富集分析表明，3482个基因主要富集在与肿瘤发生发展相关的通路上，如MAPK信号途径、VEGF信号途径、肿瘤中心碳代谢、肿瘤转录失调等（图5A），进一步，作者在TCGA数据库中下载RNA-seq表达谱数据，从中提取这些基因对应的表达谱，并对训练集中174个样本进行热图分析。结果显示，在训练集中共提取到2747个基因表达谱数据，热图如图5B所示，从热图中可以看出mRNA表达水平的高低与甲基化位点的高低相一致，从而表明这些基因的DNA甲基化水平和基因表达是一致的。
图5 甲基化位点对应基因的功能富集分析与表达量分析
4.甲基化位点的WGCNA分析
———-
经过WGCNA分析以后，得到13个模块，灰色模块表示无法聚合到其他模块中的甲基化位点（表3），1833个CpG分成了12个模块。由于第1亚型、5亚型、6亚型样本较少，因此作者选取第2亚型、3亚型、4亚型并通过Pearson相关分析对三个亚型与各模块的相关性。结果显示：Cluster3与大多数模块正相关，而Cluster4与大多数模块呈负相关（图2D）。
表3 模块中包含甲基化位点的个数
5.选择Hub-CpG位点
———-
蓝色、棕色、绿色、黑色和绿松石色模块与Cluster3和Cluster4均表现出显著的相关性，我们选择了这些模块中的甲基化位点，并计算了它们与相应模块（MM）和Cluster3表型（GS）的相关性。通过MM > 0.9和GS > 0.2鉴定了hub-CpG位点（图6A）。结果显示共有17个hub-CpG位点，大部分位于黑色模块，这些CpG位点定位到了16个基因，14个CpG位点位于基因启动子CpG岛上。此外，还通过层次聚类分析分析了17个CpG基因座之间的甲基化相关性（图6B）。
图6 识别hub-CpG位点
6．构建和验证预后模型
———-
通过上述WGCNA分析得到了17个Hub-CpG位点，接下来选择了这17个CpG位点的甲基化谱，并测定了每个样品中每个CpG位点的甲基化强度，连接节点数的对数（log（I））与该节点出现的概率的对数（p（I））呈负相关，确定加权系数的值。然后使用多元回归分析的建立风险评分模型，公式为：
根据上述公式对训练集中每个样本进行分险得分，并且观察不同RiskScore患者对应的CpG高低与患者生存状况的关系，结果显示，随着RiskScore的升高，17个CpG位点在样本甲基化水平逐渐升高，而患者的生存时间呈下降趋势（图7A），最后，根据RiskScore的中值将样本分为分为高低分险两组，高风险组患者的生存时间显著低于低风险组（P = 0.00178，图7B）。
图7 训练集中构建预后模型
在验证集中，同样根据RiskScore模型对每个样本进行风险评分，验证集中RiskScore与甲基化的强弱和患者的生存情况的相关性与训练集中的相关性相似（图8A）。此外根据RiskScore的中值将样本分为分为高低风险两组，高风险组患者的生存时间显著低于低风险组（P <0.001，图8B）。由这17个CpG位点的甲基化谱构建的预后模型可以可靠地预测NSCLC患者的预后。
图8 验证集验证预后模型
