Evaluation of stromal cell infiltration in the tumor microenvironment enable prediction of treatment sensitivity and prognosis in colon cancer

评价肿瘤微环境中的基质细胞浸润可预测结肠癌的治疗敏感性和预后

发表期刊：Comput Struct Biotechnol J

发表日期：2022 Apr 30

DOI:  10.1016/j.csbj.2022.04.037

期刊相关信息

一、背景

        结肠癌是一种以相当大的分子异质性为特征的疾病，与预后和治疗效果有关。大量文献强调了各种类型癌症的肿瘤微环境（TME）特征的临床重要性。对于结肠癌来说，一个共识的免疫分数，已经被证明是一个独立于TNM分期的强有力的预后指标，从而使结肠癌成为免疫分类的典型肿瘤。除了预后价值，免疫分数也被发现是ADJC性能的参考因素。

        共识免疫分数是通过免疫组化染色检测的，这种方法严重依赖有限的表型标志物库和足够大小的活检标本。换句话说，免疫组化染色只能反映有限的细胞类型，因此大部分关于TME全景特征的信息被忽略了。此外，免疫组化染色通常应用于分析特定的二维组织切片，所以它不能用来估计整个三维组织样本的TME含量。

二、材料与方法

1.数据来源

1） 结肠癌样本的公共转录组数据是从GEO和TCGA-COAD数据集中回顾性收集的：GSE39582、TCGA-COAD

2） 来自中山大学癌症中心队列的样本：30个组织学上诊断为非转移性结肠癌的样本的转录组数据

3） 从UCSC Xena Data Porta下载了TCGA数据集中的11种实体瘤的泛癌归一化RNA测序数据，以探索SIIS模型的广泛适用性

2.实验流程

1） 基于31种基因的定量指标的产生：设计了一个四阶段的筛选策略，以筛选出有信息量的基因来建立“基质浸润强度评分（SIIS）”；弹性网分析、Cox比例风险模型、LASSO；使用 "LIMMA "包进行集群之间的差异表达基因（DEG）分析

2） SubMap和最近的模板预测（NTP）分析：SubMap分析和NTP分析被用来探索SIIS和TMECS算法的内部关系，通过GenePattern网站的在线模块实现

3） 预测化疗反应：使用R软件包 "pRRophetic"，它分别基于CTRP数据库和GDSC数据库实现了内置的脊回归模型

4） 路径、免疫和代谢评分的估计：GSVA、ssGSEA；每个肿瘤样本中的新陈代谢激活状态和其他流行的基因特征水平是用IOBR软件包确定的

5） 多组学分析：突变分析、蛋白质组学分析

6） CRISPR协同激活介体库筛选鉴定化疗抗性基因、CRISPR/Cas9敲除库筛选

7） 连接图分析：根据Cmap网站提供的说明，使用折返变化最明显的150个基因(上调和下调)进行连接图(Cmap)分析。

8） 体外和体内实验

三、实验结果

01 - 识别与化疗获益相关的不同的TME群组

        利用无监督的分层聚类分析（图1A），根据meta-GEO队列中的990个样本确定了三个不同的TME集群，分别表示为肿瘤微环境细胞亚型1（TMECS1；n = 369）、肿瘤微环境细胞亚型2（TMECS2；n = 218）和肿瘤微环境细胞亚型3（TMECS3；n = 403）。如图1B和1C所示，TMECS1的特点是富含最高水平的基质细胞和具有免疫抑制特性的先天性免疫细胞（这些细胞类型被归入细胞群1），中等水平的细胞毒性淋巴细胞和活化树突状细胞（这些细胞类型被归入细胞群2），以及最低水平的Th2细胞、中央记忆T细胞和总T辅助细胞（这些细胞类型被归入细胞群3），这被认为是免疫排斥的表型。与TMECS1相比，TMECS2的细胞群2和3浸润的丰度明显较高，而细胞群1浸润明显减少，表明TMECS2的TME特征被认为是免疫炎症表型[13]。与TMECS1和TMECS2不同，TMECS3的细胞群1和2的丰度估计都是最低的，这与以前报道的免疫-缺失表型密切相关。而细胞群3的浸润与TMECS2的浸润相似。

        在TMECS与结肠癌分子特征的相关性方面，发现TMECS2在dMMR、CIMP（+）、BRAF突变、CMS1和C2（CIT）肿瘤中明显富集，而CIN（+）、CMS4、C4（CIT）和C6（CIT）肿瘤在TMECS1群中比其他两个群中更常见（图1D，E）。然而，在TME群中，KRAS突变或P53突变没有明显的差异。生存分析显示，在三个TME集群中，TMECS2患者的无复发生存期（RFS）和总生存期（OS）都表现出突出的生存优势，而TMECS1则与最差的预后有关（图1F）。

        接下来，评估了TME集群与化疗反应之间的关系。只有在TMECS3群中，ADJC的施用明显降低了患者的死亡风险；它没有给TMECS1或TMECS2群的患者带来生存优势（图1G）。考虑到TMECS1患者的高复发风险，推测ADJC未能给这些患者带来好处可能是由于存在原发性耐药。作者比较了TMECS1和TMECS2集群的氟尿嘧啶反应率和几种化疗药物的预测IC50值。结果显示，在所有集群中，TMECS1集群的患者表现出最低的氟尿嘧啶反应率（图1H），与TMECS2集群的患者相比，大多数化疗药物的估计IC50值明显更高（图1I）。最后，生物学分析显示，TMECS2集群有丰富的DNA损伤修复和免疫激活的途径，而TMECS1集群的特点是有多种与基质激活、化疗抵抗和癌症进展有关的途径，如刺猬信号、上皮-间质转化（EMT）、细胞色素P450代谢和转化生长因子（TGF）-β信号（图1J）。

图1 结肠癌中TME细胞的无监督聚类和TME亚型的临床特征

02 - 在meta-GEO队列中建立区分TMECS1患者的基因特征

        通过两步回归（图2A），产生了一个由31个基因组成的定量指标，称为SIIS（图2B）。箱线图显示SIIS的中位数在TMECS1群中最高，在TMECS2群中最低（图2C），随后的TME相关分析显示SIIS和细胞群1的浸润之间有很强的正相关（图2D，E），尤其是成纤维细胞和内皮细胞，表明SIIS特征可以作为测量结肠癌样本中基质细胞浸润强度的定量指标。一致的是，对公共结肠癌单细胞RNA测序数据集的分析也显示与其他细胞类型相比，成纤维细胞和内皮细胞的SIIS值分布更高（图S2A-B）。SIIS区分TMECS1和TMECS2以及所有其他非TMECS1（包括TMECS2和TMECS3）患者的曲线下面积分别为0.90和0.82（图2F）。为了进一步阐明SIIS和TMECS算法的内部关系，进行了NTP分析和SubMap分析。如图S3A和S3B所示，预测的TMECS1亚型样本的SIIS值明显高于预测的非TMECS1亚型（图S3A），高SIIS组的样本与TMECS1样本的标记基因表达有明显的相似性（图S3B）。

图S2 (A) 结肠癌单细胞 RNA 测序数据集GSE188711 细胞类型UMAP 图； (B) GSE188711 数据集中不同细胞类型中 SIIS 值分布的箱线图

图S3    GEO、TCGA 和 SYSUCC 数据集中预测的 TMECS1 和非 TMECS1 亚型中 SIIS 值分布的箱线图。 (B) 显示通过 SubMap 分析揭示的 SIIS 组和 TMECS 集群之间相似性比较的热图。

        作者试图确定SIIS模型与meta-GEO队列的临床相关性。根据单变量（图2G）和多变量分析，SIIS较高的患者的复发和死亡风险都明显增加，SIIS与OS之间的负相关关系在接受ADJC的患者中明显增加（图2G），表明这一标志物可能有能力表明化疗的好处。为了验证这一猜测，在GSE39582数据集中分析了SIIS、ADJC效益和化疗反应之间的关系。在氟尿嘧啶无应答组中SIIS的中位数明显较高，ADJC的效益随着SIIS的增加而明显下降，特别是在II期的患者中（图2H-J）。同样，在GDSC的8种化疗药物方面，高SIIS组表现出明显高于低SIIS组的估计IC50（图2K），并与晚期、pMMR、CIN（+）、C4（CIT）、CMS4以及基质和抗化疗途径的激活相关（图2L-M）。

图2 在meta-GEO队列中构建和验证SIIS模型

03 - 对SIIS特征的作用进行验证

        使用TCGA-COAD数据集（n = 382）进一步评估了SIIS特征的性能。与meta-GEO的结果一致，SIIS值在预测的TMECS1亚型中也明显升高，TCGA-COAD数据集中高SIIS组样本的转录组特征与meta-GEO中TMECS1亚型的样本密切相关。此外，在TCGA-COAD队列中，SIIS值也与基质成分浸润水平呈正相关（图3A），高SIIS组在TCGA-COAD队列中表现出与meta-GEO队列中相似的分子特征和代谢特征（图3B，E）。通过分析TCGA-COAD样本的蛋白质组数据，转录组水平上的路径富集结果在蛋白质水平上得到进一步证实（图3C）。

        利用从Thorsson等人的研究中获得的几个已知的特征，进一步测试了TCGA队列中SIIS的相关性。与研究结果一致，高SIIS组的患者表现出明显较高的EMT、TGF-β和瘤内异质性分数，而与肿瘤纯度相关的分数都下调（图3D）。然而，不同SIIS组之间在肿瘤突变负担或新抗原数量方面没有明显差异（图3D）。在临床上，较高的SIIS不仅表明死亡风险增加，而且化疗的生存效益也减少了（图3F，G）。最后，药物反应分析表明，在对TCGA-COAD数据集的分析中，较高的SIIS也明显与氟尿嘧啶无应答特征（图3H）和多种化疗药物的估计IC50升高（图3I）相关。这些结果强调了SIIS特征在预测结肠癌患者的不利预后和化疗抗性方面的稳健性。

        除了meta-GEO和TCGA-COAD队列外，还在SYSUCC队列中验证了SIIS的作用。与meta-GEO和TCGA-COAD数据库的结果一致，在预测的TMECS1中SIIS值分布较高（图S3A），高SIIS组与TMECS1的转录组相似（图S3B），在这个独立队列中也可以观察到细胞群1的浸润，基质相关和化疗抗性途径的激活，氟尿嘧啶不反应状态和CMS4亚型的丰富（图3J-M）。这些数据强调了SIIS反映结肠癌基质细胞浸润强度的性能，从而显示了强大的可重复性。

图3 在TCGA-COAD和SYSUCC队列中对SIIS模型的验证

04 - 探讨SIIS特征在泛癌症中的作用

        作者评估了SIIS在各种癌症中的作用，选择了11个TCGA队列（BLCA、BRCA、CESC、ESCA、HNSC、LIHC、LUAD、LUSC、READ、SKCM、STAD）的转录组数据进行泛癌症分析。在SIIS值的TME特征和转录组特征方面，发现SIIS值与基质细胞（成纤维细胞、内皮细胞）和肥大细胞的浸润丰度，以及TGF-β反应得分和EMT特征的激活呈正相关，但与分析的所有癌症队列中的T辅助细胞总数呈负相关（图4A）。这些结果表明，SIIS评分也可以作为泛癌基质成分的一个有效的定量工具。然而，与结肠癌队列分析结果不同的是，BRCA、CESC、HNSC、LUAD和SKCM的SIIS值与细胞毒性细胞浸润和IFNγ信号传递呈明显的负相关（图4A)。

        至于SIIS与化疗药物敏感性之间的关系，结果显示在BLCA、BRCA、CESC、HNSC、LUAD、LUSC、SKCM和STAD队列中，高SIIS组的6种以上药物的预测IC50值明显增加。值得注意的是，在这些肿瘤类型中，SIIS模型对化疗耐药性的指示作用在HNSC可能是最强的，因为所有化疗药物的t统计量都大于5（7种药物的t统计量都大于10）。有趣的是，多西他赛在高SIIS组的预测IC50的下调只能在BRCA、ESCA和LIHC队列中观察到，这表明多西他赛作为高SIIS组患者的潜在药物选择作用可能高度依赖于癌症类型。然后探讨了SIIS模型对这些癌症的生存影响。如图4B所示，10个队列（ESCA除外）当患者具有较高的SIIS时，死亡率风险增加；然而，其中只有三个队列（BLCA、SKCM和STAD）达到统计学意义（图4B）。最后，比较了SIIS在不同肿瘤类型中的分布（图4C），发现结肠癌和乳腺癌分别是SIIS分布最低和最高的肿瘤类型。进一步获得了TCGA-BRCA数据库中记载的PAM50亚型信息：腔镜A亚型患者的SIIS中值最高，而基底样疾病（也称为三阴性乳腺癌）患者的SIIS中值最低（图4D）。据报道，基底样亚型对化疗敏感；因此，这些结果进一步表明，SIIS可能能够预测乳腺癌患者的化疗疗效。

        最近有证据证明，基质细胞富集、TGF-β信号激活和EMT途径引发的免疫排斥表型与免疫治疗缺乏反应有关。考虑到泛癌症队列中SIIS和基质激活之间的紧密联系，作者接下来调查了SIIS特征是否能预测实体瘤的免疫治疗效果。分析了四个独立的免疫治疗队列，每个队列中有超过50个样本。在所有队列中，较高的SIIS与较低的治疗反应率明显相关（图4E-H）。此外，生存数据显示，低SIIS组接受免疫治疗的患者的无进展生存期（GSE176307）和总生存期（Imvigor210）比高SIIS组的患者明显延长。不同分子特征的转移性尿路癌患者的SIIS分布比较显示，免疫排斥表型、低TMB和低NEO组的SIIS明显高于其他组（图4I，J）。从TCGA亚型比较来看（图4J），SIIS在 "浸润 "亚型患者中最高，在基因组不稳定患者中最低。然而，尽管SIIS在FGFR突变体亚组（图4I）和免疫细胞上PD-L1高表达的亚组（图4J）中趋于下降，但没有达到统计学差异。最重要的是，这四个免疫治疗队列的结果表明，SIIS与免疫疗法的反应有关，可以进一步预测患者的预后。

图4 在TCGA泛癌症队列和免疫疗法治疗队列中验证SIIS模型

05 - 多组学分析

        为了揭示高SIIS组患者内在化疗抵抗的潜在机制，在TCGA-COAD队列中进行了基因组分析。有17个基因的突变频率在高、低SIIS组患者之间有明显差异，其中4个基因在高SIIS组的突变频率较高，其余13个基因的突变主要在低SIIS组富集（图5A）。这些基因的共同出现状况见图5B。通过t检验分析，发现在这13个低SIIS相关的差异突变基因中，有12个基因（除SMAD4外）的突变状态与氟尿嘧啶敏感性显著相关。进一步根据这12个基因的突变状态构建了一个综合突变评分：如果患者在这些基因中的任何一个有突变，将被计入1分。如图5D-F所示，随着突变分数的增加，SIIS值（图5D-E）和氟尿嘧啶无反应率（图5F）都有所下降。接下来的生存分析显示，这些基因没有突变的患者（突变分数为0）不能从ADJC中获益（图5G）。综上所述，这些结果表明，缺乏基因突变是造成SIIS高分患者化疗耐药的机制之一。

        除了基因组突变外，还探讨了基因启动子甲基化水平与SIIS值之间的相关性。通过对DNA启动子甲基化情况的全面调查，我们发现共有125个探针具有高变异性，其检测的β值与SIIS值显著相关。注意到有5个探针针对ADCY4基因的启动子区域，由cg14287235检测的甲基化水平与SIIS值显示出最强的相关性（图5H），在随后的分析中重点关注ADCY4。如图5I和5J所示，相关分析显示，ADCY4的甲基化与基质细胞浸润和EMT及TGFβ途径的激活水平显著负相关，而其转录表达与这些基质相关特征正相关。然而，ADCY4的转录表达和甲基化之间的关联很弱，这表明ADCY4甲基化对生物学行为的影响可能不是通过影响其转录来实现的。最后，化疗疗效分析结果也表明，对氟尿嘧啶无应答的患者，ADCY4的甲基化水平明显下降（图5K），只有ADCY4甲基化程度高的患者（β值>0.5）才倾向于从ADJC中获益（图5L）。

图5 SIIS相关的基因组改变和甲基化特征

06 - GPX3和PRICKLE2被确定为高SIIS组患者化疗耐药的驱动基因

        为了确定参与高SIIS患者化疗抗性的关键基因，作者利用人类基因组规模的CRISPR/CAS9 SAM2集合库对SW480细胞进行了功能增益筛选，40μM的奥沙利铂和2μM的5-氟尿嘧啶（5-FU）作为有效选择压力（图6A，B）。从这个筛选中，发现有488个sgRNAs子集在奥沙利铂处理的细胞中明显富集（图6C，上），而与对照相比，279个sgRNAs在5-FU处理的细胞中明显富集（图6C，下），表明这些sgRNAs靶向的基因可能是氟尿嘧啶或奥沙利铂抗性的潜在驱动因素。通过对SIIS和这些潜在驱动基因的表达进行相关性分析，在GSE39582和TCGA-COAD队列中共发现了31个与SIIS明显正相关的基因（图6C-E）。在这些基因中，GPX3和PRICKLE2引起了特别的注意，因为它们是sgRNA在奥沙利铂治疗和5-FU治疗的人群中同时显著增加的交叉基因，并且在GSE39582和TCGA-COAD队列中高SIIS组的样本中表达增加（图6F）。临床相关性分析表明，较高的GPX3和PRICKLE2表达反映了RFS和OS的不良结果，特别是对于接受ADJC的患者（图6G，左）。此外，ADJC增加了高GPX3和高PRICKLE2组患者的死亡风险（图6G，右）。Boxplots显示，GPX3和PRICKLE2的mRNA表达在氟尿嘧啶无反应组和CMS4亚型组明显升高（图6H，I）。通路和免疫分析证实，在GS39582、TCGA-COAD和SYSUCC队列中，随着GPX3和PRICKLE2表达的增加，细胞群1中成员的浸润丰度和基质通路的激活水平明显增加（图6J）。

        作者发现分析的所有三个数据库中，SIIS与WNT/PCP活性水平明显正相关（图6K），而SIIS与β-catenin依赖的典型WNT信号之间的相关性很弱（TCGA-COAD），甚至表现为明显负相关（GSE39582和SYSUCC）（图6K）。还观察到在GSE39582、TCGA-COAD和SYSUCC结肠癌队列的氟尿嘧啶无反应组和CMS4组中，WNT/PCP途径的GSVA得分持续增加（图6L，M）。总之，这些结果表明，WNT/PCP途径，而不是典型的WNT/β-catenin途径，主要促成了高SIIS患者的生物学特征。

图6 CRISPR-CAS9 SAM对高SIIS组的化疗抗性驱动基因的筛选

        由于GPX3是SIIS特征的成员之一，在接下来的研究中对GPX3进行了进一步研究。为了验证CRISPR文库筛选结果，进行了一项体外研究：将外源性GPX3质粒或针对GPX3的siRNA分别转染到SW480和CACO2结肠癌细胞系。一致的是，外源GPX3的过表达再现了SW480细胞中sgRNA介导的转录激活对5-FU和奥沙利铂处理的化疗抗性促进作用，这一点通过MTT和细胞集落形成确定（图7A）。而在体外（图7B，C）和体内（图7D，E）实验中，CACO2细胞中GPX3的下调导致对5-FU和奥沙利铂的反应急剧增加。使用石蜡包埋组织芯片检测了GPX3在结肠癌中的表达和临床意义（图7G）。基于蛋白质水平，还证实GPX3的高表达可以独立识别出无病生存率较低的患者（图7F）。在化疗效益方面，发现ADJC特别是高强度的化疗方案（奥沙利铂联合，6个月的化疗，或8个周期的Xeloda），往往会增加GPX3蛋白高表达患者的复发或死亡风险（图7H，I）。肿瘤分期的分层也观察到类似的趋势。

        为了确定可能是实现高SIIS患者化疗增敏的候选分子靶点和化合物，使用Cmap工具。在过滤结果后，在GSE39582和TCGA-COAD队列中进行的分析中确定了几个最丰富的候选分子靶点和化合物（图7J），化合物的作用机制见图7K。在交叉的分子靶点中，注意到SIAH2是一个上游调节器，需要在缺氧条件下维持HIFα的表达水平和激活其转录活性，在GSE39582和TCGA-COAD队列中都有最高的富集分数。巧合的是，通过GPX3启动子区域的HIF1α结合位点的存在，缺氧被确定为GPX3表达的强大转录调节器。因此，在转染了针对SIAH2的siRNA（siSIAH2）的CACO2细胞和转染了SIAH2过表达质粒（oxSIAH2）的SW480细胞中测试GPX3的表达，发现通过siSIAH2和oxSIAH2转染，GPX3的mRNA和蛋白水平明显下降或上升（图7L）。此外，siSIAH2转染明显增加了CACO2细胞中由5-FU或奥沙利铂诱导的凋亡细胞数量（图7M），进一步表明SIAH2调节GPX3介导的化疗抗性。此外，测试了富集度最高的化合物的化疗效果，并报告了几种抑制剂，如MLN8054、LY317615和SJ-72550，在CACO2细胞中拥有对5-FU和奥沙利铂的化疗敏感能力（图7M）。

图7 验证GPX3在推动结肠癌化疗抗性中的作用

四、结论

        总之，作者开发并验证了一个基于量化TME中肿瘤相关基质成分的新型个体评分系统，该系统可作为预测结肠癌生存率的实用和强大的生物标志物，并提供更精确的治疗方案。