肺癌是21世纪全球发生率和死亡率最高的恶性肿瘤，每年约有209万新发病例和176万死亡病例，严重威胁着人类生命健康[1-2]。其中，肺腺癌（lung adenocarcinoma，LUAD）是最普遍的组织病理学类型，发病率逐年上升[3]。目前，手术切除仍是治疗LUAD最有效的方法。虽然随着治疗方法的多元化发展，晚期肺癌患者的生存状况显著改善，然而，其5年总生存率依然低于20%[4]。肿瘤微环境是由肿瘤、基质和内皮细胞组成的异质、复杂的组织，其特点是肿瘤与先天性和适应性免疫细胞之间的交互作用[5]。同时，研究[6]发现免疫微环境可通过促进肿瘤血管生成或调节肿瘤干细胞活性等方式直接或间接影响肿瘤的发生、发展。近年来，免疫治疗作为一种新兴治疗模式，在LUAD治疗中的潜力较大。目前针对免疫检查点PD1/PD-L1、CTLA4的抗体已经广泛应用于LUAD的治疗[7]。同时，随着免疫治疗的进展，已有多项临床试验CheckMate-032 （NCT01928394）[8]、KEYNOTE-010（NCT01905657）[9]证实了Nivolumab、Pembrolizumab、Atezolizumab 较化疗均显著延长患者总生存期。然而，LUAD的免疫治疗仍存在许多困境，如耐药性增多、不良反应频发、特异性低、疗效出现瓶颈等问题，阻碍了其进一步发展。因此，挖掘安全有效的免疫相关的潜在作用靶点对于改善LUAD免疫治疗困境来说具有重要意义。
近年来，利用生物信息学分析肿瘤免疫微环境，挖掘肿瘤生物标志物是一种流行而有效的方法。例如，利用WGCNA对乳腺癌进展和预后的研究[10]表明，有一个基因模块与患者生存时间相关，并且该研究通过共表达网络分析确定了8个候选生物标志物，为进一步了解乳腺癌的分子发病机制提供了依据。也有研究[11]通过ssGSEA将The Cancer Genome Atlas（TCGA）数据库中舌鳞状细胞癌患者分类为不同免疫细胞侵蚀的集群，并利用ESTIMATE和CIBERSORT分析对集群之间的免疫状况进行评估，最终确定了15个免疫相关基因可以作为舌鳞状细胞癌患者预后的潜在生物标志物。由此可见，通过生物信息学的分析方法挖掘癌症治疗的潜在生物标志物已受到众多研究者的关注。
本研究通过WGCNA和ESTIMATE等生物信息学分析方法对TCGA数据库来源的LUAD表达数据进行分析，最终获得与免疫细胞相关的hub基因。同时，基于hub基因构建了预后风险模型，并评估了RiskScore对免疫治疗益处的预测能力。总的来说，本研究的结果有利于寻找LUAD潜在的诊断或治疗生物标志物，这对于临床为患者提供个性化的治疗策略具有重要意义。
1 资料与方法
1.1 数据来源
LUAD mRNA表达数据集（FPKM）TCGA-LUAD和对应的临床信息从TCGA数据库（
1.2 肺腺癌样本免疫微环境分析
利用R包“estimate”[12]对TCGA-LUAD数据集中LUAD肿瘤样本肿瘤微环境进行评估。
1.3 加权基因共表达网络分析
首先通过R包“WGCNA”来建立基因共表达网络[13]。然后通过使用 Pearson 相关函数测量 TCGA-LUAD 数据中基因之间的共表达。通过将共表达相似度提高到β=5来构建邻接矩阵。从邻接矩阵出发，构建拓扑重叠矩阵（TOM）来考虑拓扑相似性，并建立相应的不相似矩阵（1-TOM）来形成聚类。使用hclust函数利用不同矩阵进行分层聚类。最后基于拓扑重叠矩阵进行层次聚类，按照动态混合切割方法，每个模块最少基因数目设置为50，相异性阈值设置为0.25。计算由第一主成分（PC1）表示的模块特征基因表达式（ME），确定每个模块的表达式。临床资料与模块特征基因表达相关，并确定基因显著性（gene significance，GS）。通过将模块的ME与模块的基因表达相关联，确定模块成员的定量方法，颜色被随机分配给模块，最终选择与免疫评分、基质评分和ESTIMATE评分显著相关的模块作为后续的研究对象。为了进一步探究特征模块内基因的生物学功能，利用R包“clusterProfiler”[14]对特征模块内基因进行富集分析，特征模块基因最显著富集terms的筛选标准为P<0.05和Q<0.05。
1.4 获取hub基因
本研究以cor.geneModuleMembership>0.8 和cor.geneTraitSignificance>0.6 为阈值筛选模块内候选关键基因。同时，利用特征模块内的所有基因在STRING （http://string-db.org/cgi/input.pl）网站构建 蛋白质-蛋白质相互作用（PPI）网络，以相互作用评分>0.4 作为阈值，将获得的候选hub基因与PPI 网络中degree前30的基因取交集，最终获得hub 基因。
1.5 hub基因的进一步验证
首先比较hub基因在正常组中的表达差异。然后通过Pearson相关性分析hub 基因与各Stage分期之间的关系。Cancer Cell Line Encyclopedia（CCLE）数据库（
1.6 评估RiskScore对免疫治疗益处的预测能力
为了进一步分析RiskScore在预测治疗益处中的作用。我们在获得LUAD患者IPS数据的基础上，将患者按风险评分的中位值分为高风险组和低风险组，利用t检验比较anti-PD1和anti-CTLA4 的治疗效果在高低风险组中的差异（P<0.05），并利用R包“ggpubr”（
2 结果
2.1 WGCNA和特征模块的识别
利用R包“WGCNA”对预处理后的511例肿瘤样本MAD值最大的前5 000个基因构建加权基因共表达网络。本研究中是以软阈值为β=5来建立无尺度化网络（图1a～d）。接下来，基于平均层次聚类和动态树切割，最终共获得13个基因模块（图1e，表1）。同时，13个基因模块和临床性状的相关性分析结果表明，brown模块与StromalScore（r=0.67，P=6e-68）、ImmuneScore（r=0.91，P=1e-198）、ESTIMATEScore（r=0.86，P =3e-153）呈显著正相关，与TumorPurity呈负相关（r=–0.88，P=3e-169）（图1f），因此我们假设用于构建brown模块的基因可能与免疫力有关。我们选择brown模块作为后续分析的研究对象。
图1
肺腺癌的加权基因共表达网络
a：不同软阈值功率下的无标度拟合指数分析（β）；X轴表示软阈值，Y轴表示不同软阈值下的无标度拟合指数；红线是对应于WGCNA分析的无标度拟合指数，最优软阈值为5；b：不同软阈值功率下的平均连通性；X轴表示软阈值功率，Y轴反映平均连通性，表示不同软阈值下的网络连接性；c：β=5时连通性分布直方图；d：检查β=5时的无标度拓扑；e：基于不同度量（1-TOM）聚类的所有差异表达基因的树状图，共13个不同的模块，不同的颜色作为标识符标识不同的模块；f：热图表示模块hub基因与临床性状之间的相关性。ME表示模块本征基因，它是在对模块中的所有基因进行PCA分析后获得的主成分1（PC1）的值；红色表示与表型呈正相关，蓝色表示与表型负相关
表1
13个基因模块及其相对基因数
2.2 特征模块内基因的生物功能和通路富集分析
对brown模块内的521个基因进行富集分析，结果表明该模块内基因大多与中性粒细胞活化参与免疫反应、先天免疫反应的调节、细胞激活的正向调节、中性粒细胞活化和白细胞增殖等免疫功能相关；与lumenal side of endoplasmic reticulum membrane、integral component of lumenal side of endoplasmic reticulum membrane等内质网细胞组分相关；与MHC class Ⅱ protein complex binding和MHC class Ⅱ receptor activity等MHC受体活性有关（图2a）。KEGG富集分析结果显示该模块内基因主要富集到NF-kappa B signaling pathway、NOD-like receptor signaling pathway、primary immunodeficiency、complement and coagulation cascade和human T-cell leukemia virus 1 infection等免疫和癌症相关通路中（图2b）。上述结果表明brown模块中的基因主要与免疫功能相关，并在免疫相关途径中富集。
图2
Brown模块基因的富集分析
a：brown模块中hub基因的GO功能富集分析，纵坐标表示不同的GO项，横坐标表示基因富集数；b：brown模块特征基因KEGG通路富集分析，纵坐标为KEGG通路注释，横坐标为基因富集数。 图中的点的大小表示富集基因的数量；
2.3 与免疫相关的hub基因的鉴定
在brown基因模块中筛选出与ImmuneScore具有较高相关性的27个基因作为候选hub 基因（图3a）。同时，将brown模块中的所有基因输入STRING网站并构建PPI网络，根据degree值筛选出前30基因（图3b）。以相互作用评分>0.4作为阈值，将27个候选基因与前30基因取交集，最终获得5个hub 基因，分别为CD53、PLEK、SPI1、IL10RA和C3AR1（图3c）。
图3
与免疫相关的hub基因的鉴定
a：散点图表示brown模块中与组织学免疫评分相关的hub基因，X轴值表示基因和模块之间的相关性；如果绝对值接近0，则该基因不是该模块的一部分；如果这个值的绝对值接近1，那么该基因与该模块高度相关。Y轴值表示基因和免疫评分之间的相关性；0表示该基因与免疫力不相关，1表示该基因与免疫力高度相关，图中每个点代表一个基因；b：brown模块中PPI网络节点的前30个基因；c：韦恩图筛选最终与免疫相关的关键基因；27个候选hub基因与PPI中degree前30个基因取交集
2.4 正常组与肿瘤组hub基因的表达差异及其与肿瘤分期的相关性及hub基因的Kaplan-Meier生存分析
首先利用t检验分析了5个hub基因在正常组和肿瘤组中的表达差异，结果显示，与正常组相比，5个hub 基因在肿瘤组中显著低表达（图4a～e）。同时，分析了5个hub基因的表达与Stage分期之间的相关性，结果显示， CD53和IL10RA的表达水平与Stage分期有关联（图4f）。CCLE数据库有大量的人类癌症细胞系和独特的数据集，能够对大量人类癌症细胞系进行详细的遗传表征[18]。因此，为了探索关键基因在大量LUAD细胞系中的表达差异，我们利用CCLE数据库发现SPI1、PLEK、IL10RA、CD53和C3AR1在LUAD细胞株中的表达存在差异（图5）。同时对患者的高低表达组进行生存分析，结果发现5个hub 基因的高表达预示着患者较好的预后（图6）。
图4
分析5个关键基因在正常组和肿瘤组的表达差异及在不同肿瘤分期中的mRNA表达变化
a～e：正常组和肿瘤组中C3AR1、CD53、IL10RA、PLEK和SPI1表达的差异；f：箱线图表示5个hub基因在不同肿瘤分期中的mRNA表达情况。
图5
5个hub基因在大量肺腺癌细胞系中的表达
通过CCLE分析肺腺癌细胞系中SPI1（a）、PLEK（b）、IL10RA（c）、CD53（d）和C3AR1（e）的表达水平；红色表示高表达，蓝色表示低表达
图6
5个hub基因的Kaplan-Meier生存分析
a～e：C3AR1，CD53，IL10RA，PLEK 以及 SPI1C3AR1的Kaplan-Meier生存分析
2.5 分析hub 基因与免疫细胞和免疫检查点的相关性
为了进一步分析5个hub基因与免疫细胞和免疫检查点的相关性，通过TIMER网站对5个hub 基因进行评估。结果表明，这5个hub 基因的表达水平与6种免疫细胞的浸润丰度均呈现显著正相关，与肿瘤纯度呈现显著负相关（图7）。同时，与免疫检查点的相关性分析表明这5个hub 基因的表达水平与免疫检查点CTLA4和PDCD1的表达量同样呈现显著正相关（图8）。
图7
5个hub基因与免疫细胞及肿瘤纯度的相关性分析
a～e：C3AR1，CD53，IL10RA，PLEK 以及 SPI1C3AR1与免疫细胞和肿瘤纯度的相关性；X轴表示免疫细胞的浸润水平和肿瘤纯度，Y轴表示基因表达水平
图8
5个hub基因与免疫检查点（CTLA4和PDCD1）的相关性分析
a：5个hub基因与CTLA4的相关性研究；b：5个hub基因与PDCD1的相关性研究
2.6 Hub 基因的GSEA富集分析
对5个hub 基因的高低表达组进行GSEA富集分析，结果表明，5个hub 基因均富集到Primary immunodeficiency、PD-L1 expression and PD-1 checkpoint pathway in cancer、NF-kappa B signaling pathway和Toll-like receptor signaling pathway等免疫相关的通路（图9）。
图9
5个关键基因的GSEA富集分析
a～e：C3AR1，CD53，IL10RA，PLEK 以及 SPI1C3AR1的GSEA富集分析。上图为通路的富集分数折线，下图中的线对应各通路的基因
2.7 多因素cox回归分析和生存分析
基于5个免疫相关的hub 基因进行多因素cox回归分析并构建了预后风险模型，RiskScore=（−0.04×ExpressionSPI1）+（−0.23×ExpressionPLEK）+（−0.13×ExpressionIL10RA）+（−0.06×ExpressionCD53）+（0.3×ExpressionC3AR1）（图10a）。同时，对患者的高低风险组进行了Kaplan-Meier生存分析。结果显示，高风险组的患者预后较差（图10b）。随后，使用外部数据集GSE72094验证模型，ROC曲线结果显示，1年、3年和5年生存期的AUC值分别为0.63、0.64和0.65（图10d）。生存曲线显示，高风险组和低风险组的生存状态存在显著差异，高风险组的预后较差（图10c）。以上结果证明，由5个特征基因构建的模型具有合理性。
图10
5个hub基因的多因素Cox回归分析和 Kaplan-Meier生存分析
a：多因素Cox回归分析森林图；b：高风险和低风险组的Kaplan-Meier生存分析；c：外部验证集中高风险组和低风险组的Kaplan-Meier生存分析；d： 外部验证集的ROC曲线
2.8 评估RiskScore对免疫治疗益处的预测能力
为了进一步分析RiskScore在预测治疗益处中的作用，我们分析了高低风险组LUAD患者的IPS数据。IPS代表患者免疫原性，高IPS意味着患者有更高的免疫性能[19]。我们发现，IPS、IPS-CTLA4、IPS-PD1和IPS-CTLA4 and PD1在RiskScore组内差异具有统计学意义（P均<0.05）（图11a～d）。先前的研究发现，TIDE可以准确预测患者的预后[20]。因此，我们计算了高低风险组患者的TIDE评分，结果显示，低风险组的TIDE得分显著低于高风险组（图11e），这表明低风险组更能从免疫检查点阻断治疗中受益。基于上述结果，我们认为我们建立的RiskScore在预测免疫治疗方面有很大的潜力。
图11
肺腺癌患者的IPS、TIDE和RiskScore之间的关系
a～d：RiskScore组的IPS、IPS-CTLA4、IPS-PD1、IPS-CTL A4和PD1有显著差异（均
3 讨论
近年来，随着诊断和治疗方法的改进，LUAD的死亡率和发病率有了一定改善，但由于LUAD易转移、易复发等原因导致该病患者的长期预后仍不令人满意。目前，多项研究[21-22]发现肿瘤微环境的免疫细胞浸润会影响肿瘤免疫治疗的预后和临床结果。因此，寻找免疫细胞浸润相关关键基因为寻找有效诊断或治疗LUAD的潜在靶点提供了一定的理论基础。
本研究中，我们通过生存分析，发现5个hub基因的表达水平和LUAD患者的预后密切相关。同样地，有研究[23]发现CD53的表达与三阴性乳腺癌患者的 无病生存期正相关，并且CD53在免疫细胞中表达，是淋巴细胞活化所必需的。相似地，有研究[24]发现CD53与骨肉瘤患者的生存时间显著相关，并且在骨肉瘤转移中发挥了重要作用，可以作为骨肉瘤转移的候选生物标志物。Bai等[25]的研究发现，在LUAD中，PLEK是与肿瘤纯度负相关的肿瘤纯度共表达基因，并且可以作为预后风险因素预测LUAD患者的临床结局。也有研究[26]表明SPI1诱导的lncRNA SNHG6上调可以通过miR-485-3p/VPS45轴促进前列腺癌进展。同时，Cheng等[27]的研究表明IL10RA的表达与黑色素瘤患者的预后显著相关，并且体外实验发现IL10RA能显著抑制黑色素瘤细胞的增殖、迁移和侵袭，在黑色素瘤的进展中发挥了重要作用。Qu等[28]的研究发现，C3AR1与食管鳞状细胞癌患者的生存时间相关，并且C3AR1可能通过影响巨噬细胞向M2表型的极化而导致免疫抑制微环境，从而进一步影响食管鳞状细胞癌的进展，提示C3AR1可能是免疫治疗的潜在靶点。综上所述，这些研究表明我们获得的hub 基因对肿瘤的发生发展和患者预后中都产生了重要影响。
多项研究[29-30]已证实，肿瘤的发生和进展会受到TME中免疫细胞和因子的调节，这些都可能会成为癌症治疗有希望的潜在靶点。在本研究中，我们发现5个hub 基因的表达水平与6种免疫细胞的浸润丰度是显著正相关的。已有研究[31-32]表明CD53是一种four-transmembrane protein，主要表达于骨髓淋巴系统，对B细胞功能是至关重要的，因为CD53促进BCR依赖性蛋白激酶C信号传导，使其能够磷酸化其底物，从而进一步促进B细胞活化。PLEK已被认为与免疫逃逸有关，并且是多种癌症的潜在治疗靶点[33]。SPI1原癌基因对于骨髓细胞、淋巴细胞和巨噬细胞的发育、成熟、分化和调控至关重要[34]。同时，研究[35]表明，SPI1的高表达与胃癌的免疫细胞浸润相关，并参与细胞周期调控，导致预后不良。在膀胱癌中，SPI1可靶向调控CXCL12的表达，并募集肿瘤相关巨噬细胞[36]。IL10RA缺陷是原发性免疫缺陷病的病因[37]。IL10RA和IL10结合，出发JAK1-STAT3信号通路，抑制促炎细胞因子，还降低巨噬细胞活性，抑制单核细胞分化以及抗原呈递细胞上共刺激分子的下调[38]。在癌症中，IL10RA的抗炎和抗免疫反应也已经被报道多次[39-40]。C3AR1也在卵巢癌[41]，肾透明细胞癌[42]，骨肉瘤[43]等多种癌症中发现与免疫细胞浸润、介导肿瘤免疫微环境有关。由此可见，我们获得的5个hub基因与免疫细胞的功能及激活密切相关，可以作为免疫细胞相关的潜在生物标志物进一步应用于LUAD的临床免疫治疗。
值得注意的是，我们对5个hub 基因进行了GSEA富集分析，结果发现这5个hub基因主要富集到Intestinal immune network for IgA production、Th1 and Th2 cell differentiation、PD-L1 expression and PD-1 checkpoint pathway in cancer、Th17 cell differentiation、NF-kappa B signaling pathway、T cell receptor signaling pathway和Toll-like receptor signaling pathway等免疫相关的通路。Yang等[44]的研究发现Bufalin可以通过APOBEC3F诱导的intestinal immune network for IgA production 信号通路抑制HCC细胞的增殖和迁移。Liu等[45]的研究发现，用从草珊瑚中纯化的酸性多糖处理树突状细胞后，会使Th1 and Th2 cell differentiation通路中的DLL4基因显著上调，从而进一步增强抗肿瘤免疫。Sun等[46]的研究表明，肿瘤外泌体传递的CRNDE-h通过抑制结直肠癌中Itch介导的泛素化和 RORγt降解来促进Th17 cell differentiation，从而进一步影响结直肠癌的发生和发展。Jin等[47]的研究则发现，TRIM7可以通过负调控NF-kappa B 信号通路显著抑制肺癌细胞的增殖和迁移，为肺癌的治疗提供了潜在靶点。同时，有研究[48]报道Toll-like receptor signaling pathway参与激活先天性和适应性免疫反应，在炎症诱导的疾病中起着关键作用，该信号通路的失调可导致上皮层止血紊乱、慢性炎症、过度修复反应和结直肠癌的发生。这些结果提示本研究所获得的这5个hub基因可能会通过其相关的免疫信号通路在LUAD的发生和发展中发挥重要作用。
多项研究[49-50]证实，免疫治疗为癌症患者带来了新的希望，尤其是免疫检查点抑制剂的发现和应用在癌症免疫治疗中发挥了不可或缺的作用。如PD1 和 CTLA4等抑制剂有广泛而多样的机会通过调节免疫反应来增强抗肿瘤免疫[51]。在本研究中，我们发现基于5个hub 基因所构建的预后风险模型在预测免疫治疗方面有很大的潜力，这为临床实践中治疗不同风险患者的免疫治疗策略的发展铺垫了道路。在本研究中发现，低风险组的患者更可能从免疫检查点治疗当中获益，这是基于理论层次上的研究。在实际情况当中，由于高风险患者的病情更为紧迫，病情更为复杂，医生更倾向于尽快解决高风险患者的危急情况，稳定患者病情[52]。而对低风险的患者，尽管理论上他们更可能从治疗当中获益，在实际治疗当中也应该考虑到药物带来的副作用和成本[53]。
综上所述，我们的研究从理论上揭示了与免疫细胞密切相关的5个枢纽基因的巨大潜力，其不仅可以被选为LUAD免疫治疗的预后生物标志物，而且可以被选作LUAD免疫疗法的治疗靶点。然而，本研究也存在一定的不足之处，本研究的结果均属于理论层次上的研究，缺少临床样本和相关实验的验证。因此，对于使用外部验证集验证出hub基因的诊断效果不是十分理想，后期通过相关实验对本研究的结果进行进一步验证是必要的。
利益冲突：无。
作者贡献：何东元负债设计并实施研究，分析数据，撰写文章；陈波提供研究思路，分析研究方案可行性，并参与文章的修订与审核；梁靖瑶参与提出文章的撰写思路，收集并分析数据；叶海波：收集并分析数据，设计论文框架，并参与文章的修订与终版文章的审核；易小杏、梁广妮收集并分析数据，并参与文章的修订；
肺癌是21世纪全球发生率和死亡率最高的恶性肿瘤，每年约有209万新发病例和176万死亡病例，严重威胁着人类生命健康[1-2]。其中，肺腺癌（lung adenocarcinoma，LUAD）是最普遍的组织病理学类型，发病率逐年上升[3]。目前，手术切除仍是治疗LUAD最有效的方法。虽然随着治疗方法的多元化发展，晚期肺癌患者的生存状况显著改善，然而，其5年总生存率依然低于20%[4]。肿瘤微环境是由肿瘤、基质和内皮细胞组成的异质、复杂的组织，其特点是肿瘤与先天性和适应性免疫细胞之间的交互作用[5]。同时，研究[6]发现免疫微环境可通过促进肿瘤血管生成或调节肿瘤干细胞活性等方式直接或间接影响肿瘤的发生、发展。近年来，免疫治疗作为一种新兴治疗模式，在LUAD治疗中的潜力较大。目前针对免疫检查点PD1/PD-L1、CTLA4的抗体已经广泛应用于LUAD的治疗[7]。同时，随着免疫治疗的进展，已有多项临床试验CheckMate-032 （NCT01928394）[8]、KEYNOTE-010（NCT01905657）[9]证实了Nivolumab、Pembrolizumab、Atezolizumab 较化疗均显著延长患者总生存期。然而，LUAD的免疫治疗仍存在许多困境，如耐药性增多、不良反应频发、特异性低、疗效出现瓶颈等问题，阻碍了其进一步发展。因此，挖掘安全有效的免疫相关的潜在作用靶点对于改善LUAD免疫治疗困境来说具有重要意义。
近年来，利用生物信息学分析肿瘤免疫微环境，挖掘肿瘤生物标志物是一种流行而有效的方法。例如，利用WGCNA对乳腺癌进展和预后的研究[10]表明，有一个基因模块与患者生存时间相关，并且该研究通过共表达网络分析确定了8个候选生物标志物，为进一步了解乳腺癌的分子发病机制提供了依据。也有研究[11]通过ssGSEA将The Cancer Genome Atlas（TCGA）数据库中舌鳞状细胞癌患者分类为不同免疫细胞侵蚀的集群，并利用ESTIMATE和CIBERSORT分析对集群之间的免疫状况进行评估，最终确定了15个免疫相关基因可以作为舌鳞状细胞癌患者预后的潜在生物标志物。由此可见，通过生物信息学的分析方法挖掘癌症治疗的潜在生物标志物已受到众多研究者的关注。
本研究通过WGCNA和ESTIMATE等生物信息学分析方法对TCGA数据库来源的LUAD表达数据进行分析，最终获得与免疫细胞相关的hub基因。同时，基于hub基因构建了预后风险模型，并评估了RiskScore对免疫治疗益处的预测能力。总的来说，本研究的结果有利于寻找LUAD潜在的诊断或治疗生物标志物，这对于临床为患者提供个性化的治疗策略具有重要意义。
1 资料与方法
1.1 数据来源
LUAD mRNA表达数据集（FPKM）TCGA-LUAD和对应的临床信息从TCGA数据库（
1.2 肺腺癌样本免疫微环境分析
利用R包“estimate”[12]对TCGA-LUAD数据集中LUAD肿瘤样本肿瘤微环境进行评估。
1.3 加权基因共表达网络分析
首先通过R包“WGCNA”来建立基因共表达网络[13]。然后通过使用 Pearson 相关函数测量 TCGA-LUAD 数据中基因之间的共表达。通过将共表达相似度提高到β=5来构建邻接矩阵。从邻接矩阵出发，构建拓扑重叠矩阵（TOM）来考虑拓扑相似性，并建立相应的不相似矩阵（1-TOM）来形成聚类。使用hclust函数利用不同矩阵进行分层聚类。最后基于拓扑重叠矩阵进行层次聚类，按照动态混合切割方法，每个模块最少基因数目设置为50，相异性阈值设置为0.25。计算由第一主成分（PC1）表示的模块特征基因表达式（ME），确定每个模块的表达式。临床资料与模块特征基因表达相关，并确定基因显著性（gene significance，GS）。通过将模块的ME与模块的基因表达相关联，确定模块成员的定量方法，颜色被随机分配给模块，最终选择与免疫评分、基质评分和ESTIMATE评分显著相关的模块作为后续的研究对象。为了进一步探究特征模块内基因的生物学功能，利用R包“clusterProfiler”[14]对特征模块内基因进行富集分析，特征模块基因最显著富集terms的筛选标准为P<0.05和Q<0.05。
1.4 获取hub基因
本研究以cor.geneModuleMembership>0.8 和cor.geneTraitSignificance>0.6 为阈值筛选模块内候选关键基因。同时，利用特征模块内的所有基因在STRING （http://string-db.org/cgi/input.pl）网站构建 蛋白质-蛋白质相互作用（PPI）网络，以相互作用评分>0.4 作为阈值，将获得的候选hub基因与PPI 网络中degree前30的基因取交集，最终获得hub 基因。
1.5 hub基因的进一步验证
首先比较hub基因在正常组中的表达差异。然后通过Pearson相关性分析hub 基因与各Stage分期之间的关系。Cancer Cell Line Encyclopedia（CCLE）数据库（
1.6 评估RiskScore对免疫治疗益处的预测能力
为了进一步分析RiskScore在预测治疗益处中的作用。我们在获得LUAD患者IPS数据的基础上，将患者按风险评分的中位值分为高风险组和低风险组，利用t检验比较anti-PD1和anti-CTLA4 的治疗效果在高低风险组中的差异（P<0.05），并利用R包“ggpubr”（
2 结果
2.1 WGCNA和特征模块的识别
利用R包“WGCNA”对预处理后的511例肿瘤样本MAD值最大的前5 000个基因构建加权基因共表达网络。本研究中是以软阈值为β=5来建立无尺度化网络（图1a～d）。接下来，基于平均层次聚类和动态树切割，最终共获得13个基因模块（图1e，表1）。同时，13个基因模块和临床性状的相关性分析结果表明，brown模块与StromalScore（r=0.67，P=6e-68）、ImmuneScore（r=0.91，P=1e-198）、ESTIMATEScore（r=0.86，P =3e-153）呈显著正相关，与TumorPurity呈负相关（r=–0.88，P=3e-169）（图1f），因此我们假设用于构建brown模块的基因可能与免疫力有关。我们选择brown模块作为后续分析的研究对象。
图1
肺腺癌的加权基因共表达网络
a：不同软阈值功率下的无标度拟合指数分析（β）；X轴表示软阈值，Y轴表示不同软阈值下的无标度拟合指数；红线是对应于WGCNA分析的无标度拟合指数，最优软阈值为5；b：不同软阈值功率下的平均连通性；X轴表示软阈值功率，Y轴反映平均连通性，表示不同软阈值下的网络连接性；c：β=5时连通性分布直方图；d：检查β=5时的无标度拓扑；e：基于不同度量（1-TOM）聚类的所有差异表达基因的树状图，共13个不同的模块，不同的颜色作为标识符标识不同的模块；f：热图表示模块hub基因与临床性状之间的相关性。ME表示模块本征基因，它是在对模块中的所有基因进行PCA分析后获得的主成分1（PC1）的值；红色表示与表型呈正相关，蓝色表示与表型负相关
表1
13个基因模块及其相对基因数
2.2 特征模块内基因的生物功能和通路富集分析
对brown模块内的521个基因进行富集分析，结果表明该模块内基因大多与中性粒细胞活化参与免疫反应、先天免疫反应的调节、细胞激活的正向调节、中性粒细胞活化和白细胞增殖等免疫功能相关；与lumenal side of endoplasmic reticulum membrane、integral component of lumenal side of endoplasmic reticulum membrane等内质网细胞组分相关；与MHC class Ⅱ protein complex binding和MHC class Ⅱ receptor activity等MHC受体活性有关（图2a）。KEGG富集分析结果显示该模块内基因主要富集到NF-kappa B signaling pathway、NOD-like receptor signaling pathway、primary immunodeficiency、complement and coagulation cascade和human T-cell leukemia virus 1 infection等免疫和癌症相关通路中（图2b）。上述结果表明brown模块中的基因主要与免疫功能相关，并在免疫相关途径中富集。
图2
Brown模块基因的富集分析
a：brown模块中hub基因的GO功能富集分析，纵坐标表示不同的GO项，横坐标表示基因富集数；b：brown模块特征基因KEGG通路富集分析，纵坐标为KEGG通路注释，横坐标为基因富集数。 图中的点的大小表示富集基因的数量；
2.3 与免疫相关的hub基因的鉴定
在brown基因模块中筛选出与ImmuneScore具有较高相关性的27个基因作为候选hub 基因（图3a）。同时，将brown模块中的所有基因输入STRING网站并构建PPI网络，根据degree值筛选出前30基因（图3b）。以相互作用评分>0.4作为阈值，将27个候选基因与前30基因取交集，最终获得5个hub 基因，分别为CD53、PLEK、SPI1、IL10RA和C3AR1（图3c）。
图3
与免疫相关的hub基因的鉴定
a：散点图表示brown模块中与组织学免疫评分相关的hub基因，X轴值表示基因和模块之间的相关性；如果绝对值接近0，则该基因不是该模块的一部分；如果这个值的绝对值接近1，那么该基因与该模块高度相关。Y轴值表示基因和免疫评分之间的相关性；0表示该基因与免疫力不相关，1表示该基因与免疫力高度相关，图中每个点代表一个基因；b：brown模块中PPI网络节点的前30个基因；c：韦恩图筛选最终与免疫相关的关键基因；27个候选hub基因与PPI中degree前30个基因取交集
2.4 正常组与肿瘤组hub基因的表达差异及其与肿瘤分期的相关性及hub基因的Kaplan-Meier生存分析
首先利用t检验分析了5个hub基因在正常组和肿瘤组中的表达差异，结果显示，与正常组相比，5个hub 基因在肿瘤组中显著低表达（图4a～e）。同时，分析了5个hub基因的表达与Stage分期之间的相关性，结果显示， CD53和IL10RA的表达水平与Stage分期有关联（图4f）。CCLE数据库有大量的人类癌症细胞系和独特的数据集，能够对大量人类癌症细胞系进行详细的遗传表征[18]。因此，为了探索关键基因在大量LUAD细胞系中的表达差异，我们利用CCLE数据库发现SPI1、PLEK、IL10RA、CD53和C3AR1在LUAD细胞株中的表达存在差异（图5）。同时对患者的高低表达组进行生存分析，结果发现5个hub 基因的高表达预示着患者较好的预后（图6）。
图4
分析5个关键基因在正常组和肿瘤组的表达差异及在不同肿瘤分期中的mRNA表达变化
a～e：正常组和肿瘤组中C3AR1、CD53、IL10RA、PLEK和SPI1表达的差异；f：箱线图表示5个hub基因在不同肿瘤分期中的mRNA表达情况。
图5
5个hub基因在大量肺腺癌细胞系中的表达
通过CCLE分析肺腺癌细胞系中SPI1（a）、PLEK（b）、IL10RA（c）、CD53（d）和C3AR1（e）的表达水平；红色表示高表达，蓝色表示低表达
图6
5个hub基因的Kaplan-Meier生存分析
a～e：C3AR1，CD53，IL10RA，PLEK 以及 SPI1C3AR1的Kaplan-Meier生存分析
2.5 分析hub 基因与免疫细胞和免疫检查点的相关性
为了进一步分析5个hub基因与免疫细胞和免疫检查点的相关性，通过TIMER网站对5个hub 基因进行评估。结果表明，这5个hub 基因的表达水平与6种免疫细胞的浸润丰度均呈现显著正相关，与肿瘤纯度呈现显著负相关（图7）。同时，与免疫检查点的相关性分析表明这5个hub 基因的表达水平与免疫检查点CTLA4和PDCD1的表达量同样呈现显著正相关（图8）。
图7
5个hub基因与免疫细胞及肿瘤纯度的相关性分析
a～e：C3AR1，CD53，IL10RA，PLEK 以及 SPI1C3AR1与免疫细胞和肿瘤纯度的相关性；X轴表示免疫细胞的浸润水平和肿瘤纯度，Y轴表示基因表达水平
图8
5个hub基因与免疫检查点（CTLA4和PDCD1）的相关性分析
a：5个hub基因与CTLA4的相关性研究；b：5个hub基因与PDCD1的相关性研究
2.6 Hub 基因的GSEA富集分析
对5个hub 基因的高低表达组进行GSEA富集分析，结果表明，5个hub 基因均富集到Primary immunodeficiency、PD-L1 expression and PD-1 checkpoint pathway in cancer、NF-kappa B signaling pathway和Toll-like receptor signaling pathway等免疫相关的通路（图9）。
图9
5个关键基因的GSEA富集分析
a～e：C3AR1，CD53，IL10RA，PLEK 以及 SPI1C3AR1的GSEA富集分析。上图为通路的富集分数折线，下图中的线对应各通路的基因
2.7 多因素cox回归分析和生存分析
基于5个免疫相关的hub 基因进行多因素cox回归分析并构建了预后风险模型，RiskScore=（−0.04×ExpressionSPI1）+（−0.23×ExpressionPLEK）+（−0.13×ExpressionIL10RA）+（−0.06×ExpressionCD53）+（0.3×ExpressionC3AR1）（图10a）。同时，对患者的高低风险组进行了Kaplan-Meier生存分析。结果显示，高风险组的患者预后较差（图10b）。随后，使用外部数据集GSE72094验证模型，ROC曲线结果显示，1年、3年和5年生存期的AUC值分别为0.63、0.64和0.65（图10d）。生存曲线显示，高风险组和低风险组的生存状态存在显著差异，高风险组的预后较差（图10c）。以上结果证明，由5个特征基因构建的模型具有合理性。
图10
5个hub基因的多因素Cox回归分析和 Kaplan-Meier生存分析
a：多因素Cox回归分析森林图；b：高风险和低风险组的Kaplan-Meier生存分析；c：外部验证集中高风险组和低风险组的Kaplan-Meier生存分析；d： 外部验证集的ROC曲线
2.8 评估RiskScore对免疫治疗益处的预测能力
为了进一步分析RiskScore在预测治疗益处中的作用，我们分析了高低风险组LUAD患者的IPS数据。IPS代表患者免疫原性，高IPS意味着患者有更高的免疫性能[19]。我们发现，IPS、IPS-CTLA4、IPS-PD1和IPS-CTLA4 and PD1在RiskScore组内差异具有统计学意义（P均<0.05）（图11a～d）。先前的研究发现，TIDE可以准确预测患者的预后[20]。因此，我们计算了高低风险组患者的TIDE评分，结果显示，低风险组的TIDE得分显著低于高风险组（图11e），这表明低风险组更能从免疫检查点阻断治疗中受益。基于上述结果，我们认为我们建立的RiskScore在预测免疫治疗方面有很大的潜力。
图11
肺腺癌患者的IPS、TIDE和RiskScore之间的关系
a～d：RiskScore组的IPS、IPS-CTLA4、IPS-PD1、IPS-CTL A4和PD1有显著差异（均
3 讨论
近年来，随着诊断和治疗方法的改进，LUAD的死亡率和发病率有了一定改善，但由于LUAD易转移、易复发等原因导致该病患者的长期预后仍不令人满意。目前，多项研究[21-22]发现肿瘤微环境的免疫细胞浸润会影响肿瘤免疫治疗的预后和临床结果。因此，寻找免疫细胞浸润相关关键基因为寻找有效诊断或治疗LUAD的潜在靶点提供了一定的理论基础。
本研究中，我们通过生存分析，发现5个hub基因的表达水平和LUAD患者的预后密切相关。同样地，有研究[23]发现CD53的表达与三阴性乳腺癌患者的 无病生存期正相关，并且CD53在免疫细胞中表达，是淋巴细胞活化所必需的。相似地，有研究[24]发现CD53与骨肉瘤患者的生存时间显著相关，并且在骨肉瘤转移中发挥了重要作用，可以作为骨肉瘤转移的候选生物标志物。Bai等[25]的研究发现，在LUAD中，PLEK是与肿瘤纯度负相关的肿瘤纯度共表达基因，并且可以作为预后风险因素预测LUAD患者的临床结局。也有研究[26]表明SPI1诱导的lncRNA SNHG6上调可以通过miR-485-3p/VPS45轴促进前列腺癌进展。同时，Cheng等[27]的研究表明IL10RA的表达与黑色素瘤患者的预后显著相关，并且体外实验发现IL10RA能显著抑制黑色素瘤细胞的增殖、迁移和侵袭，在黑色素瘤的进展中发挥了重要作用。Qu等[28]的研究发现，C3AR1与食管鳞状细胞癌患者的生存时间相关，并且C3AR1可能通过影响巨噬细胞向M2表型的极化而导致免疫抑制微环境，从而进一步影响食管鳞状细胞癌的进展，提示C3AR1可能是免疫治疗的潜在靶点。综上所述，这些研究表明我们获得的hub 基因对肿瘤的发生发展和患者预后中都产生了重要影响。
多项研究[29-30]已证实，肿瘤的发生和进展会受到TME中免疫细胞和因子的调节，这些都可能会成为癌症治疗有希望的潜在靶点。在本研究中，我们发现5个hub 基因的表达水平与6种免疫细胞的浸润丰度是显著正相关的。已有研究[31-32]表明CD53是一种four-transmembrane protein，主要表达于骨髓淋巴系统，对B细胞功能是至关重要的，因为CD53促进BCR依赖性蛋白激酶C信号传导，使其能够磷酸化其底物，从而进一步促进B细胞活化。PLEK已被认为与免疫逃逸有关，并且是多种癌症的潜在治疗靶点[33]。SPI1原癌基因对于骨髓细胞、淋巴细胞和巨噬细胞的发育、成熟、分化和调控至关重要[34]。同时，研究[35]表明，SPI1的高表达与胃癌的免疫细胞浸润相关，并参与细胞周期调控，导致预后不良。在膀胱癌中，SPI1可靶向调控CXCL12的表达，并募集肿瘤相关巨噬细胞[36]。IL10RA缺陷是原发性免疫缺陷病的病因[37]。IL10RA和IL10结合，出发JAK1-STAT3信号通路，抑制促炎细胞因子，还降低巨噬细胞活性，抑制单核细胞分化以及抗原呈递细胞上共刺激分子的下调[38]。在癌症中，IL10RA的抗炎和抗免疫反应也已经被报道多次[39-40]。C3AR1也在卵巢癌[41]，肾透明细胞癌[42]，骨肉瘤[43]等多种癌症中发现与免疫细胞浸润、介导肿瘤免疫微环境有关。由此可见，我们获得的5个hub基因与免疫细胞的功能及激活密切相关，可以作为免疫细胞相关的潜在生物标志物进一步应用于LUAD的临床免疫治疗。
值得注意的是，我们对5个hub 基因进行了GSEA富集分析，结果发现这5个hub基因主要富集到Intestinal immune network for IgA production、Th1 and Th2 cell differentiation、PD-L1 expression and PD-1 checkpoint pathway in cancer、Th17 cell differentiation、NF-kappa B signaling pathway、T cell receptor signaling pathway和Toll-like receptor signaling pathway等免疫相关的通路。Yang等[44]的研究发现Bufalin可以通过APOBEC3F诱导的intestinal immune network for IgA production 信号通路抑制HCC细胞的增殖和迁移。Liu等[45]的研究发现，用从草珊瑚中纯化的酸性多糖处理树突状细胞后，会使Th1 and Th2 cell differentiation通路中的DLL4基因显著上调，从而进一步增强抗肿瘤免疫。Sun等[46]的研究表明，肿瘤外泌体传递的CRNDE-h通过抑制结直肠癌中Itch介导的泛素化和 RORγt降解来促进Th17 cell differentiation，从而进一步影响结直肠癌的发生和发展。Jin等[47]的研究则发现，TRIM7可以通过负调控NF-kappa B 信号通路显著抑制肺癌细胞的增殖和迁移，为肺癌的治疗提供了潜在靶点。同时，有研究[48]报道Toll-like receptor signaling pathway参与激活先天性和适应性免疫反应，在炎症诱导的疾病中起着关键作用，该信号通路的失调可导致上皮层止血紊乱、慢性炎症、过度修复反应和结直肠癌的发生。这些结果提示本研究所获得的这5个hub基因可能会通过其相关的免疫信号通路在LUAD的发生和发展中发挥重要作用。
多项研究[49-50]证实，免疫治疗为癌症患者带来了新的希望，尤其是免疫检查点抑制剂的发现和应用在癌症免疫治疗中发挥了不可或缺的作用。如PD1 和 CTLA4等抑制剂有广泛而多样的机会通过调节免疫反应来增强抗肿瘤免疫[51]。在本研究中，我们发现基于5个hub 基因所构建的预后风险模型在预测免疫治疗方面有很大的潜力，这为临床实践中治疗不同风险患者的免疫治疗策略的发展铺垫了道路。在本研究中发现，低风险组的患者更可能从免疫检查点治疗当中获益，这是基于理论层次上的研究。在实际情况当中，由于高风险患者的病情更为紧迫，病情更为复杂，医生更倾向于尽快解决高风险患者的危急情况，稳定患者病情[52]。而对低风险的患者，尽管理论上他们更可能从治疗当中获益，在实际治疗当中也应该考虑到药物带来的副作用和成本[53]。
综上所述，我们的研究从理论上揭示了与免疫细胞密切相关的5个枢纽基因的巨大潜力，其不仅可以被选为LUAD免疫治疗的预后生物标志物，而且可以被选作LUAD免疫疗法的治疗靶点。然而，本研究也存在一定的不足之处，本研究的结果均属于理论层次上的研究，缺少临床样本和相关实验的验证。因此，对于使用外部验证集验证出hub基因的诊断效果不是十分理想，后期通过相关实验对本研究的结果进行进一步验证是必要的。
利益冲突：无。
作者贡献：何东元负债设计并实施研究，分析数据，撰写文章；陈波提供研究思路，分析研究方案可行性，并参与文章的修订与审核；梁靖瑶参与提出文章的撰写思路，收集并分析数据；叶海波：收集并分析数据，设计论文框架，并参与文章的修订与终版文章的审核；易小杏、梁广妮收集并分析数据，并参与文章的修订；
