A novel prognostic signature of immune-related genes for patients with colorectal cancer
01
02
—
目前，结直肠癌（CRC）是全世界第三常见的诊断癌症，死亡率方面排名第二。尽管，新的治疗方法使晚期大肠癌患者的总生存率（OS）增加一倍，然而，大肠癌的诊断率却仍然很低，患者往往到了晚期才能被诊断。免疫疗法现在是许多肿瘤标准方案的一部分，抗肿瘤免疫治疗方法分为六类：溶瘤病毒治疗、辅助免疫治疗和细胞因子、过继细胞治疗、治疗性疫苗、检查点抑制剂和单克隆抗体。然而，免疫反应不足一直是一个问题，特别是针对CRC中PD-1和PD-L1s的检查点抑制剂，免疫检查点抑制剂（ICIs）的使用在大多数转移性CRC患者中很少或没有起作用。
近年来，基因芯片和RNA测序已经鉴定出来与几种人类癌症的预后相关的生物标志物，而且针对多种癌症类型建立了几种与免疫相关的预后标志物。正是基于这些研究，本研究的目的是在CRC患者中建立与免疫相关的预后标志物。
03
—
1.数据来源
-TCGA数据：568例CRC组织和44例正常组织的RNA-seq数据以及相应的临床随访数据。
-GEO数据：GSE38832包含115例CRC样本表达谱芯片数据与临床随访信息，注释平台为GPL570.
-临床样本数据：来源于浙江大学医学院附属第二医院25对CRC组织和匹配的正常组织
-特定数据集：2498个与人体免疫相关的基因（IRGs）来源于ImmPort数据库
2.数据分析
-差异分析：R程序包limma筛选肿瘤组织与正常组织之间的差异基因，阈值：FDR<0.05和log2|fold change|>1。
-富集分析：R程序包clusterProfiler进行GO与KEGG富集分析，阈值：FDR<0.01
-筛选与预后相关IRG:将435例癌症患者按照8：2分为训练集与验证集，使用R程序包survival进行单因素cox回归分析，阈值p<0.05.
-构建与验证基因风险模型：将单因素cox回归分析与预后显著的IRGs进行多因素回归分析，建立预后基因风险模型。根据风险得分中值将训练组的患者分为低风险组和高风险组，使用Kaplan‐Meier高低风险组患者之间的生存差异。使用R软件包survival ROC计算模型1年、3年和5年的AUC。最后使用RTqPCR使用进行风险基因临床验证
-风险得分与患者临床因素关联分析：使用R程序包rms构建列线图，R程序包regplot计算碎石图。
-生物信息学分析：使用R程序包GSVA进行GSEA分析，筛选与风险得分相关的通路，阈值：FDR < 0.001和相关系数>0.3；Spearman相关系数分析免疫相关基因与高低风险患者之间的相关性；STRING数据库用于PPI分析，阈值：confidence>0.4,nodes>10；使用Cistrome数据库分析差异转录因子（TFs），阈值：FDR<0.05和log2|fold change|>1；TarBase和LncBase数据库构建了风险基因与ceRNA调控网络，并使用Cytoscape软件进行了可视化。
统计分析：Wilcoxon检验，Kruskal-Wallis检验，log‐rank检验，阈值：p<0.05
04
—
1.筛选与免疫相关的差异基因
——–
在TCGA正常样本中和肿瘤样本筛选出476个IRG的差异基因，其中后上调基因177个，下调基因299个（图2A-B）。然后，作者对IRG的差异基因进行功能分析，GO功能富集分析结果这些基因参与了补体激活，体液免疫反应，蛋白激活，免疫球蛋白介导的免疫反应，抗原结合，吞噬作用等功能中显著富集（图2C）；KEGG通路富集分析结果这些基因参与了细胞因子-细胞因子受体相互作用中，IL-17信号通路，自然杀伤细胞介导的细胞毒性，趋化因子信号通路，Rap1信号和MAPK信号等通路中显著富集（图2D）。
图2. 免疫相关基因差异分析与差异基因功能分析
2.构建基因风险模型
——–
作者首先使用训练集中携带的生存数据进行单因素cox回归分析,共得到18个与患者预后显著相关的基因（图3A）。然后作者将单因素显著的基因进行多因素cox回归分析，得到8个基因与患者的预后显著相关，8个基因多因素cox回归分析结果见表2。
表2. 8个分析基因多因素cox回归结果
接下来，作者基于每个基因多因素cox回归系数和每个的表达值构建了一个由8个基因组成的风险评分公式：risk score =(0.639 * SLC10A2 ) + (0.387 * FGF2 ) + (−0.094 * CCL28 ) + (0.012 * NDRG1 ) + (0.124 * ESM1 ) + (0.378 * UCN ) + (0.254 * UTS2 ) + (0.129 * TRDC )，然后基于这个风险评分公式对训练集的每一位患者进行风险评分，并根据所有样本的风险得分的中值（0.95）将患者分为高低风险两组（图3B），从图中可以看出，高风险组患者的生存时间明显低于低风险组患者（log‐rank P<0.0001）。此外，作者使用ROC曲线来评估模型的预测性能，结果显示，3年和5年风险得分的AUC为0.79与0.76（图3C），图3D展示了训练集中患者的生存状态，风险评分和风险基因表达量之间的关系。
图3 构建与免疫相关的基因风险模型
3.验证集验证风险模型的鲁棒性
——–
为了验证免疫基因风险模型的鲁棒性，作者使用内部验证集、整个TCGA数据集对基因模型进行验证。结果显示，在内部验证集中低风险组患者的生存时间显著高于高风险组(n = 91, P =0.006, 图4A)，风险模型3年和5年的AUC为0.75和0.82（图4B）。在整个TCGA数据集中低风险组患者的生存时间显著高于高风险组(n = 435, P <0.0001, 图4C)，风险模型3年和5年的AUC为0.77和0.78（图4B）。因此，8个IRG风险模型能准确的预测CRC患者的生存状况。另外，在外部数据集GSE38832中低风险组患者的生存时间显著高于高风险组(n = 115, P =0.016, 图4F)，并且高低两组患者的风险评分存在显著差异(秩和检验, P<0.0001,图4F)。最后，作者使用qRT-PCR技术验证了25对肿瘤样本和正常样本中基因ESM1和基因SLC10A2表达情况，得到与生信分析类似的结果，在肿瘤组织中ESM1表达量显著升高，SLC10A2表达量显著下降（P<0.05）。
图4.验证集验证基因风险模型的鲁棒性
4.IRG风险得分与对不同临床类型CRC患者的预测能力
——–
作者首先将临床因素包括年龄、pStage、pT、pN、pM和免疫风险评分进行单因素cox回归分析，结果如图5A所示。进一步使用多因素cox回归模型验证了风险模型能够独立与其他临床变量，对患者的预后进行独立的预测（P<0.001，图5B和表3）。8个IRG风险基因的表达谱如5C所示，从图中可以得出结论随着风险得分的增加肿瘤越来越靠近肿瘤晚期。
表3. 基因风险模型与临床病理因素单因素与多因素cox回归分析
图5.基因风险模型与临床因素的相关性
此外，作者分析CCL28的表达水平与pT（P=0.048），ESM1的pT（P=0.008）和pStage（P=0.01），具有pM的UTS2（P=0.039）和具有pN的TRDC（P=0.012），pM（P<0.0001）和pStage（P<0.0001，图6A）这些临床指标显著相关。作者结合临床病理特征和IRG基因风险模型构建列线图和碎石图综合评估了IRG基因风险模型对CRC患者的预测性能（图6B-C）.
图6 风险模型中的基因与临床因素的相关性，列线图，碎石图
5.IRG基因风险模型与生物学通路及调控网络的相关性
——–
作者使用TCGA表达谱数据，使用GSEA评估了基因风险模型与生物学通路之间得到联系，结果显示有38条KEGG富集通路与高低风险患者预后情况显著相关。调节癌症代谢相关功能的途径，如性粘连、肌动蛋白细胞骨架、ECM受体相互作用和内吞作用与风险得分呈正相关，而胞浆DNA传感、亚油酸代谢、花生四烯酸代谢和醚脂代谢途径与风险得分呈负相关。
最后，作者分析了在低风险组和高风险组中几种CRC生物标记物的在肿瘤组织中表达情况。BRAF、NRAS和PIK3CA在高低风险组中表达存在显著差异，这表明IRG风险评分与CRC进展密切相关（图7A-C）。在此基础上，利用TarBase和LncBase数据库构建了一个复杂的ceRNA网络来调控IRGs（图7D）。
图7. 免疫相关基因与生物标志物的相关性以及构建结直肠癌竞争性内源性RNA（ceRNA）网络
