Integrated Multi-omics Profiling Yields a Clinically-Relevant Molecular Classification for Esophageal Squamous Cell Carcinoma

食管鳞状细胞癌（esophageal squamous cell carcinomas），WGS, whole-genome bisulfite sequencing (WGBS), RNA sequencing (RNA-seq), small RNA sequencing (sRNA-seq), and proteomics on 155 ESCC cases

前言结构：

流行病学，食管癌 5 年生存率 20%，大概 50 % 的病例发生在东亚和中亚，特别是中国，ESCC 是主要的组织学类型
治疗策略，目前主要是放疗和化疗，ESCC 没有特定的治疗靶点，VEGFR2 拮抗剂和 PDL1 抗体已被批准用作 ESCC 的治疗，但是反应率异质性非常大，目前没有有效的生物标志物
目前研究现状，病理学亚型而非分子亚型，大都是单组学研究，缺乏多组学表征
本文做的事

基于多组学的 ESCC 分型

155 个样本，进行了多组学测序，包括：WGS, WGBS, RNA-seq 和 sRNA-seq（小 RNA 测序，这个做的比较少），其中 73 个样本还进行了蛋白组测序。利用 WGS 分析突变和拷贝数，进一步利用突变分析突变特征，变异频率等，利用 RNA-seq 和蛋白组进行差异表达分析，甲基化数据分析低甲基化和高甲基化区域。

这部分的主要结论为：

对 ESCC 基因组贡献最大的 mutation signature 是 SBS1 （年龄相关）和 SBS2/13 （APOBEC 活性相关）。
发生频率最高的变异类型为 TP53, CDKN2A, and NFE2L2, amplifications in 11q13.3, and deletions in 9p21.3/CDKN2A and 4q35.2.

聚类方法：

数据预处理：对于 mRNA 数据，TMM（Trimmed Mean of M-values）表达在超过 70% 的样本中小于 1 的基因被剔除；对于 miRNA 数据，则是 CPM（Counts Per Million）在超过 70% 的样本中小于 1 的基因被剔除；前 25 % 的可变 mRNA 和 miRNA 被用来聚类。甲基化数据 beta 值的方差是通过在基因 promoter 区域（定义为 TSS 上下游 1000 bp）计算得到，有着最可变甲基化信号的前 1% 基因被用来聚类；对于 CNV，则使用 GISTIC2 得到的 q-value<0.25 作为筛选标准
对不同的组学分别用一致性聚类进行聚类，但是用的参数不同
多组学整合：使用的方法是 cluster of cluster assignments (COCA) integrative clustering analysis, 也就是整合不同组学的聚类信息，形成最终的样本类别

在聚类分析中，特别是在层次聚类（Hierarchical Clustering）的方法中，“linkage”指的是用于定义或度量不同簇之间相似性的方式，进而决定如何合并这些簇以形成更大的簇。Consensus Clustering（一致性聚类）是一种集成方法，它通过多次运行聚类算法并整合结果来提高聚类稳定性与可靠性，而在每次单独的层次聚类过程中也可能用到不同的linkage策略。

在层次聚类中，常见的linkage方法包括：

Single Linkage（单链接）：两个簇之间的距离定义为两个簇中最接近（最近邻）的两个对象之间的距离。

Complete Linkage（全链接）：两个簇之间的距离定义为两个簇中最远（最远邻）的两个对象之间的距离。

Average Linkage（平均链接）：又称为UPGMA（Unweighted Pair Group Method with Arithmetic mean），两个簇之间的距离定义为两个簇中所有配对对象距离的平均值。

Centroid Linkage（质心链接）：两个簇之间的距离是根据两个簇的质心（即所有对象坐标的均值向量）之间的距离计算的。

Ward's Linkage（沃德链接）：这是一种最小化簇内平方和变异的策略，合并过程选择使得簇间平方误差增量最小的簇进行合并。

在Consensus Clustering中，可能会结合不同的linkage方法来生成多个聚类解决方案，并通过共识过程确定最终的稳定聚类结构。

COCA 方法：文献 Multiple kernel learning for integrative consensus clustering of omic datasets

开始对于每个组学进行聚类，然后生成一个 Matrix Of Clusters (MOC)，这个矩阵的行数等于每种组学中的聚类个数，列数等于样本的数目，如果样本属于某个组学中的某类，那么矩阵的这个位置就是 1，否则就是 0，然后再将这个矩阵输入一致性聚类进行聚类

最终利用这种方法将样本分成四个亚型，并且和临床特征相关（TNM 分期，淋巴转移，和免疫评分），依据不同分型的特征将其命名：

The four subtypes consisted of the following: a cell cycle signaling activation (CCA) subtype (25.2%, n = 39) characterized with CNAs for CCND1 and CDKN2A/B (Figures 1A, 1D, and S4A); a NRF2 oncogenic activation (NRFA) subtype (24.5%, n = 38) clustering with upregulation of NFE2L2 and SOX2 genes (Figures 1A, 1E, and S4A); an immune suppression (IS) (19.3%, n = 30) subtype with higher levels of infiltrating B cells and NK CD56bright cells, and higher ERBB2 protein level; and an immune modulation (IM) subtype (31.0%, n = 48) with higher infiltrating of fibroblasts, CD8+ T cells and innate immune cells, and high immune score

CCA 亚型和激活的细胞周期信号相关

CCA 亚型的特征包括 11q13.3/CCND1 频繁扩增以及 9p21.3 CDKN2A/B 的纯合子缺失，这两个区域的扩展和缺失使得这些区域内的基因表达相对于正常样本发生显著的升高或降低，另外 CDKN2A/B 除了同源删除外，也通过其他的途径失活，如突变或启动子区域的超甲基化等，总的来说在 CCA 中，有 84.6 % 的样本含有直接参与细胞周期失调的基因组变异。

CDK4/6 抑制剂 palbociclib 是靶向细胞周期变异的药物，并且已批准作为乳腺癌的治疗药物，作者做了病人来源的类器官（分别来自 CDKN2A 删除，CDKN2A 删除+CCND1 扩增，以及两个基因都是 WT），比较在这些类器官中，palbociclib 药物的效果：

这个药物看起来浓度都是 uM 级别的，效果也不是太好。

总的来说：Thus, a significant proportion of CCA patients might benefit from CDK4/6 inhibitor therapy.

NRFA 亚型的特征是 NRF2 和 SOX2 信号的激活

NRFA 亚型的特点是 NRF2 信号通路基因的变异比较多，包括 NFE2L2，KEAP1 和 CUL3。并且 NFE2L2 的表达在 NRFA 样本中是显著升高的，即使在那些没有 NRF2 通路基因突变的样本中也是这样，同时 KEAP1 的蛋白水平在 NRFA 亚型中是下调的。并且在 NRFA 亚型中 SOX2 也是频繁扩增和高表达的，将 NRF2 通路和 SOX2 基因综合考虑，在 NRFA 亚型中有 73.7 % 的样本有着 NRF2 通路基因变异或/和 SOX2 基因扩增。假设：转录因子 SOX2 高表达，从而导致 NRF2 激活。

接下来作者验证了 SOX2 induced transcription of NFE2L2：

mRNA 水平两者的表达是显著正相关的，
在两个 ESCC 细胞系 ZEC145 and ZEC166 进行了 CHIP，发现 SOX2 是特异性结合到 NFE2L2 的启动子区域，并且也使用了公共数据的 Chip-seq 进行了验证
敲低 SOX2 显著降低 NRF2 和 NQO1 的 mRNA 和蛋白水平（NQO1 是被 NRF2 调控表达的）（NFE2L2 就是 NRF2）

作者进一步利用 IHC 来比较不同 NRFA 亚型中的不同组：

NFE2L2 变异的样本，包括 NFE2L2 amplifications and NFE2L2 mutations
调控 NFE2L2 的基因变异，包括 KEAP1/CUL3 mutations, KEAP1 hypermethylation, and SOX2 amplifications
在 NRF2 通路相关的基因中没有任何变异

总的来说，NRFA 亚型可能会受益于 NRF2 抑制剂（虽然目前没有这类药物）

CCA 和 NRFA 富集 CpG 超甲基化表型

进行常规的 WGBS 甲基化数据分析，发现和正常组织相比肿瘤样本表现出全局低甲基化而 CpG 岛区域（CGI）高甲基化的特征；分析差异甲基化区域（DMR），发现 HyperDMR 和 hypoDMR 的分布不同，HyperDMR主要是富集在启动子区域，特别是启动子的 CGI 区域，而 hypoDMR 则主要分布在基因间区，增强子，以及基因主体区域。作者定义了一个 esophageal CpG island methylator phenotype (E-CIMP) 表型，基于 208 个甲基化位点（这些位点在超过 50% 的样本中是超甲基化的），E-CIMP 阳性病人和临床表型显著相关，并且富集在 CCA 和 NRFA 亚型中。

IS 和 IM 亚型有着不同的肿瘤免疫浸润细胞分布

之前已经分析过 IS 和 IM 亚型有着较高的免疫评分：

为了系统研究不同亚型的 TME 差别，作者解析了主要的肿瘤浸润免疫细胞和基质细胞的分布（包括 25 种免疫细胞，内皮细胞和成纤维细胞）；基于不同细胞类型的 signature（marker），利用 ssgsea 计算 NES 值，作为细胞类型的丰度估计。发现 IS 和 IM 亚型浸润的免疫细胞较其他亚型要显著高，并通过多标记免疫荧光（multilabel immuno-fluorescence，MIF）进行验证，IS 亚型主要富集 B 细胞和 NK 细胞，而 IM 亚型富集 CD8 细胞和巨噬细胞：

为了研究不同浸润细胞的协同作用，作者开发了一种算法，TME module，该方法先计算了细胞类型之间的 NES 值相关性（Spearman’s rank correlation），保留相关系数大于 0.6 且 P 值小于 0.0001 的细胞对，构建一个细胞之间的相互作用网络。利用 cluster_edge_betweenness 算法来鉴别这个图中的细胞 communities/module（一个社区的细胞可能协同发挥作用），再对每个模块进行打分（下面式子中 A 表示细胞的丰度，n 表示细胞数量）：

通过这种方法发现，在 IM 中占主要的免疫细胞亚群模块是 “CD8+ T-Macrophage”，包含 CD8+ T cells,
myeloid-derived suppressor cells (MDSCs), and macrophages，而 IS 中主要的模块是 “CD4+ Tcm-Tem”，包含 CD4+ T central memory (Tcm) and T effector memory (Tem) cells。并且 dendritic cells (DCs), interstitial DCs and plasmacytoid DCs 富集在 IM 中，并且和巨噬细胞有着强相关性，但是巨噬细胞和 CD8 细胞没有出现在 IS 的主要模块中。

IM 亚型肿瘤表现出较高的肿瘤突变负担，而IS肿瘤表现出较高的 ERBB2 蛋白表达水平，这表明 IM 亚型的肿瘤更适用于免疫检查点抑制剂疗法，而 IS 亚型适用于 ERBB2 靶向疗法（ERBB2 又叫 HER2）。

IM 亚型分类器可以预测免疫治疗反应

43 个接受抗 PD-1 抗体治疗的 ESCC 病人（I 期临床试验），进行多组学测序，mRNA expression data for 38 samples, exome sequencing data for 23 samples, and DNA methylation data for 21 samples。

基于差异基因将这 38 个样本和前面的 155 个样本放在一起进行聚类，发现 PR 样本和 IM 亚型聚在一起，而 PD 则和 CCA/NRFA/IS 亚型聚在一起，这个结果显示 IM 亚型可能比其他亚型对 ICB 治疗反应更好。

基于这个结果，作者接下来开发了一个分类器来基于多组学特征对样本的 IM 亚型进行分类（是否为 IM 亚型）。设计流程是利用 155 个 ESCC 样本构建模型，在 TCGA 样本上进行外部验证，在免疫治疗数据集上进行预测效果。作者首先选择了 48 个多组学特征，通过特征选择步骤，最终得到 28 个特征用于构建随机森林模型，在验证集上的 AUC 达到 0.969，并且这 28 个特征都是重要的，因为删除任一基因都会显著降低模型性能。

分类模型训练步骤：

选择初始 48 个特征；前面的聚类分析和差异分析可以得到 12 个 CNV 特征，20 个基因的表达。两个免疫细胞毒性和耗竭相关的基因集也被纳入考虑，算这些基因集基因的平均表达作为特征。基于 IM-vs-NonIM 和 PR-vs-NonPR 的差异表达，对于 PR-vs-NonPR，选择显著上调和下调的两个基因集，对于 IM-vs-NonIM 选择显著上调和下调的基因分别进行富集分析，在显著富集的通路中出现频率最多的 10 个基因被选择（上调和下调分别 10 个），这样基于差异分析就可以得到 4 个基因集，对基因集中的基因进行 CPM 标准化后取平均作为特征值。差异甲基化也进行同样的分析，得到 10 个差异甲基化特征，最终得到 12 CNV + 26 mRNA + 10 Meth = 48 个特征

数据处理：首先基于上面选择的特征对样本进行聚类，从而鉴别 6 个离群样本（可以在聚类树上设定距离阈值），另外为了避免正负例不平衡（155 个样本中 IM 负例和正例样本的比例为 1:2.2），对 IM 正例样本进行过采样，方法是 SMOTE（Synthetic Minority Oversampling Technique），从而使得正负样本比例为 1:1

特征选择和模型训练：采用迭代特征选择的方法从上述 48 个特征中选择具有最佳分类能力的特征，利用随机森林（RF）构建模型；70% 训练集，30% 验证集；迭代特征选择的流程为：

首先使用所有的 48 个特征训练 RF

计算 48 个特征的模型重要性得分，并降序排列

移除排在后面的 4~5 个特征（大概占 10%）构成新的特征集合

使用新的特征集合训练新的 RF

计算 out of bag （OOB）误差，如果 OOB 上升，则将 c 步移除的特征中重要性最大的特征拿回来

重复 c, d, e 直到特征数量收敛，记录下所有可能的特征集合和相应的模型

选择具有最低 OOB 误差并且在验证集上表现最好的特征集合

随机森林是一种集成模型，使用的是 bagging（bootstrap aggregating 缩写）集成方法，也就是在训练集的随机选择的不同子集上进行训练，并且这种随机选择是有放回的抽样。在使用 bagging 时，有些实例可能会被抽到多次而另一些实例可能根本不会被抽到；在 BaggingClassifier 中默认是有放回地抽取 n 个训练实例（n 也是训练集的大小），因此平均来说，对于每个预测器大概只有 63% 的实例被抽到，剩下的 37% 的实例就叫做 out-of-bag (oob) 实例；由于这些实例没有被抽取作为训练样本，所以可以在这些实例上进行模型的评估

同时作者也使用单组学特征进行建模，经过上面同样的特征选择步骤，最终 mRNA-only 使用 20 个基因表达特征，CNV-only 模型使用 11 个 CNV 特征，methylation-only 模型使用 22 个差异甲基化特征，这些单组学模型的性能没有使用多组学特征的模型好：

作者将模型应用于免疫治疗数据集，其中有 17 个样本有着全部 28 个多组学特征，模型预测 7 个样本是 IM，10 个样本是 non-IM；IM 组有着更好的 OS 生存，同时通过 IHC 检测病人的 PDL1 表达，发现 PDL1 表达阳性和阴性样本间的生存没有显著差异：

在这 17 个病人中，6/7 个 PR 样本被模型识别为 IM 类，9/10 个 PD 或 SD 样本被模型识别成 non-IM，因此达到了 85.7% 的灵敏度和 90% 的特异性，AUC 为 0.857：

在 81 个 TCGA ESCC 样本中，模型鉴别出 17 个 IM 肿瘤，IM 样本显示出更高的免疫评分和 NK 细胞以及巨噬细胞的浸润，与本文中的 ESCC 样本类似：

Synthetic Minority Oversampling Technique

参考：Synthetic Minority Over-sampling TEchnique (SMOTE)

在训练分类模型时，经常会遇到数据集中样本类别不平衡的情况，可能会导致模型性能不佳。一种解决方法是对样本量较多的那一类进行降采样，但是这样会丢弃一些训练样本，这在样本量稀缺的领域可能不太适合（比如临床样本）。另一种方法则是对样本量较少的那一类进行过采样，即随机复制一些观测值加入训练数据，但是这种方法可能会导致模型由于过于关注某些样本从而过拟合，SMOTE 要解决的就是这类问题。SMOTE 一般分为四步：

对于一个样本，找到距离较近的邻近样本
计算样本和邻近样本之间的差值
将差值乘上一个 0-1 之间的随机数，将样本加上这个值就得到一个新的样本

#pip install imbalanced-learn
from imblearn.over_sampling import SMOTE
sm = SMOTE(random_state=42, k_neighbors=5) ##默认情况下，得到的正例和负例样本一样多
X_res, y_res = sm.fit_resample(X, y)

‍

目录CONTENT

ESCC 分型文献