【摘要】 目的 分析肺腺癌单细胞转录组测序数据与TCGA数据，通过单细胞转录组测序数据的分组解析与TCGA肺腺癌样本的WGCNA,筛选与预后相关的核心基因，构建肺腺癌预后分子预测模型并验证其预测效能，通过肺腺癌细胞实验验证其临床转化潜力。方法 (1)自NCBI基因表达数据库获取GSE149655单细胞转录组测序数据，采用FindVariableFeatures函数筛选高变异基因后进行主成分分析，采用SingleR软件进行亚群标注，采用FindAllMarkers函数进行亚群分析筛选标记基因。自TCGA数据库下载肺腺癌样本的FPKM基因表达数据和临床信息，筛选与单细胞分析筛选的标记基因方差>0.5的FPKM基因，采用cibersort函数预测TCGA数据集中每个样本的亚群评分，比较肿瘤组织与正常组织样本各亚群评分。(2)采用Pearson相关性系数评估基因两两间的共表达关系，生成基因相似性矩阵；采用WGCNA构建加权共表达模块，筛选关键基因。(3)采用coxph函数对关键基因进行单因素Cox分析，筛选与预后相关的基因。采用nmf函数对相关基因进行聚类，采用非负矩阵分解聚类法划分肿瘤样本亚型，采用ESTIMATE算法对亚型进行免疫评分，采用Pearson相关法计算免疫评分与亚群评分的相关性。比较不同亚型患者的临床特征，绘制Kaplan-Meier生存曲线，比较不同亚型患者的3年总生存率；富集分析各样本中差异最显著的前20条通路。(4)采用limma包筛选不同亚型肿瘤组织与正常组织间差异表达基因，采用单因素Cox分析差异表达基因与患者预后的关系，采用LASSO回归筛选基因变量，采用stepAIC函数进一步筛选核心基因。根据回归系数构建风险模型，计算风险评分=∑(回归系数×核心基因表达量),对风险评分进行Z-score标准化后，将风险评分>0分定义为高风险，风险评分≤0分定义为低风险。绘制ROC曲线，分析风险评分对数据集肺腺癌患者1、2、3年预后的预测效能。绘制Kaplan-Meier生存曲线，比较高风险与低风险者3年总生存率。(5)以GSE31210为外部验证数据集，绘制风险评分分布图和生存状态散点图，采用Z-score标准化方法评估数据集中样本的风险评分分布。绘制ROC曲线，评估分子模型预测数据集样本1、3、5年预后的效能。绘制Kaplan-Meier生存曲线，比较高风险与低风险者5年总生存率。(6)取第3代对数生长期肺成纤维细胞系WI-38及肺腺癌细胞系NCI-H1395、A549、NCI-H1975,采用Western blot法检测4种细胞NID2、ANGPTL4、AKAP12、GJB3、CCL20、TPSB2蛋白相对表达量。(7)取对数生长期NCI-H1395、A549细胞分别分为NC-NCI-H1395组(转染shRNA-NC慢病毒)和shRNA-NCI-H1395组(转染shRNA-TPSB2慢病毒),NC-A549组(转染shRNA-NC慢病毒)和sh-A549组(转染shRNA-TPSB2慢病毒)。转染48 h,采用Transwell小室实验检测4组细胞侵袭能力，计算平均侵袭细胞数；采用EdU实验检测4组细胞转染48 h时EdU阳性细胞数，以EdU阳性细胞数表示细胞增殖能力。结果 (1)GSE149655单细胞数据中包含2份肺腺癌原发肿瘤组织和2份癌旁组织，4份样本中共包含22 357个基因和12 554个细胞。SingleR共注释出11个亚群，筛选出每个亚群中特异性表达最高的前5个基因作为Top5标记基因。TCGA的肺腺癌FPKM基因表达数据及临床信息包含513份肿瘤样本和59份正常样本，共筛选出19 611个基因，其中472份肺腺癌样本包含生存时间和生存状态数据。肿瘤组织与正常组织样本在C0和C2亚群评分比较差异均无统计学意义(P均>0.05),在其余9个亚群中评分比较差异均有统计学意义(t=2.356～8.741,P均<0.05)。(2)WGCNA分析获得6个共表达模块，从与C5亚群评分相关性最强的模块中筛选出关键基因333个。(3)coxph函数对333个关键基因进行单因素回归Cox分析，共筛选出25个与预后显著相关的基因。nmf函数对25个基因进行非负矩阵分解聚类，将472份肿瘤样本划分为L1和L2亚型。ESTIMATE算法预测2种亚型的免疫评分，结果显示L2亚型免疫评分[(2 318.75±25.63)分]高于L1亚型[(1 125.36±18.42)分](t=10.234,P<0.001),免疫评分与C5亚群评分呈正相关(r=0.613,P=0.001),L2亚型患者T3-4期、N2-3期、临床分级Ⅲ～Ⅳ级占比及C5亚群评分均高于L1亚型患者(χ~2=6.325～15.782,P均<0.05)。L2亚型患者3年总生存率低于L1亚型患者(χ~2=12.453,P=0.001)。功能富集分析筛选出不同临床特征患者差异最显著的前20条通路。(4)L1、L2亚型中共鉴定出337个基因表达上调及401个基因表达下调，单因素Cox分析筛选出90个与预后有关的基因，LASSO回归分析10倍交叉验证确定最优λ值为0.031 579 93时18个基因最具泛化能力的特征变量，采用stepAIC进一步回归分析筛选出6个核心基因(NID2、ANGPTL4、AKAP12、GJB3、CCL20、TPSB2)。风险评分=0.156×NID2+0.089×ANGPTL4+0.091×AKAP12+0.121×GJB3+0.084×CCL20-0.212×TPSB2。ROC曲线分析结果显示，分子模型风险评分预测肺腺癌患者1、3、5年预后的AUC分别为0.71(95%CI:0.721～0.849,P=0.001)、0.71(95%CI:0.753～0.871,P=0.001)、0.65(95%CI:0.730～0.862,P=0.001)。Kaplan-Meier生存分析结果显示，高风险者3年总生存率低于低风险者(χ~2=18.642,P<0.001）。TPSB2基因高表达者3年总生存率高于低表达者（χ~2=8.765,P=0.002）。（5）风险评分分布分析结果显示，高风险组与低风险组样本在验证集中分布不同，死亡样本多集中于高风险组。ROC曲线分析结果显示，分子模型预测外部验证数据集样本1、3、5年预后的AUC分别为0.79(95%CI:0.678～0.828,P=0.001)、0.62(95%CI:0.705～0.847,P=0.001)、0.73(95%CI:0.689～0.835,P=0.001)。Kaplan-Meier生存分析结果显示，高风险者5年总生存率低于低风险者（χ~2=10.325,P=0.003）。（6）NCI-H1395、A549、NCI-H1975细胞NID2、ANGPTL4、AKAP12、GJB3、CCL20蛋白相对表达量均高于WI-38细胞（t=4.559～25.770,P均<0.05）,TPSB2蛋白相对表达量均低于WI-38细胞（t=-17.668～-5.516,P均<0.05）。（7）转染48h,NC-NCI-H1395、NC-A549组侵袭细胞数及EdU阳性细胞数分别少于shRNA-NCI-H1395组、shRNA-A549组（t=8.240～16.310,P均<0.05）。结论 肺腺癌细胞NID2、ANGPTL4、AKAP12、GJB3、CCL20基因表达上调，TPSB2基因表达下调，基于单细胞转录组测序与TCGA多组学数据分析建立的包含6个核心基因的分子模型对肺腺癌患者的预后具有较好的预测效能，TPSB2可能是治疗肺腺癌的新靶点。更多还原

【Abstract】 Objective To analyze single-cell transcriptome sequencing data and TCGA data, to screen core prognosis-related genes of lung adenocarcinoma(LUAD) patients through subgroup analysis of single-cell sequencing data and WGCNA of the TCGA-LUAD cohort, to construct a molecular prognostic model for LUAD to validate its predictive performance, and to verify their clinical transformation potential via LUAD cell experiments. Methods Single-cell RNA sequencing data of GSE149655 were obtained from the NCBI Gene Expression Omnibus database. Highly variable genes were screened using the FindVariableFeatures function, followed by principal component analysis. Subgroup annotation was conducted using the SingleR software. Marker genes were screened for subgroup analysis using the FindAllMarkers software. FPKM gene expression data and clinical information of LUAD samples were downloaded from the TCGA database. FPKM genes with a variance > 0.5 among the marker genes identified from the single-cell analysis were selected. The cibersort function was used to predict the subgroup scores for each sample in the TCGA dataset, and the subgroup scores were compared between tumor tissue and normal tissue.(2) Pearson correlation coefficients were used to assess pairwise co-expression relationships between genes, generating a gene similarity matrix. WGCNA was employed to construct weighted co-expression modules and identify key genes.(3) Univariate Cox analysis was performed on the key genes using the coxph function to screen for prognosis-related genes. Non-negative matrix factorization(NMF) clustering of these related genes was performed using the nmf function to categorize tumor samples into subtypes. The ESTIMATE algorithm was used to calculate immune scores for each subtype, and the correlation between immune scores and subtype scores was assessed using Pearson correlation. Clinical characteristics were compared between different tumor subtypes. Kaplan-Meier survival curves were plotted to compare the 3-year overall survival rates among patients of different tumor subtypes. The top 20 most significantly enriched pathways were identified through enrichment analysis.(4) The limma package was used to identify differentially expressed genes between tumor and normal tissues. Univariate Cox regression analysis was performed to assess the relationship between these differentially expressed genes and prognosis. LASSO regression was applied to select gene variables, and the stepAIC function was used for further screening core genes. A risk model was constructed based on the regression coefficients: Risk Score = Σ(regression coefficient × core gene expression). The risk score was standardized using Z score. Samples with a risk score > 0 were defined as high-risk, and those with a risk score ≤ 0 were defined as low-risk. ROC curves were plotted to analyze the predictive performance of the risk score for 1-, 3-, and 5-year prognosis in the LUAD patient dataset. Kaplan-Meier survival curves were plotted to compare the 3-year overall survival rates between high-risk and low-risk groups.(5) The GSE31210 dataset was used as an external validation set. The distribution plots of risk scores and scatter plots of survival status were generated. The risk scores of samples in the validation dataset were assessed using Z-score standardization. ROC curves were plotted to evaluate the model’s predictive performance for 1-, 3-, and 5-year prognosis in the validation dataset samples. Kaplan-Meier survival curves were plotted to compare the 5-year overall survival rates between high-risk and low-risk groups.(6) The lung fibroblast cell line WI-38 and LUAD cell lines NCI-H1395, A549 and NCI-H1975 in the logarithmic growth phase were selected. Western blot was used to detect the relative protein expressions of NID2, ANGPTL4, AKAP12, GJB3, CCL20, and TPSB2 in these four cell lines.(7) NCI-H1395 and A549 cells in logarithmic growth phase were divided into NC-NCI-H1395 group(transfected with shRNA-NC lentivirus), shRNA-NCI-H1395 group(transfected with shRNA-TPSB2 lentivirus), NC-A549 group(transfected with shRNA-NC lentivirus), and sh-A549 group(transfected with shRNA-TPSB2 lentivirus), respectively. At 48 h post-transfection, cell invasion ability was detected using the Transwell chamber assay, and the average count of invading cells was calculated. The EdU assay was performed to detect the EdU-positive cell count at 48 h post-transfection in the four groups, with the EdU-positive cell count representing cell proliferation capacity. Results(1) The GSE149655 single-cell data included 2 copies of primary LUAD tumor tissue and 2 of adjacent normal tissue, encompassing 22 357 genes and 12 554 cells across the four samples. SingleR annotation identified 11 cell subgroups. The top 5 specifically expressed genes in each subgroup were selected as Top5 marker genes. The TCGA LUAD FPKM gene expression data and clinical information included 513 tumor samples and 59 normal tissue samples, from which 19 611 genes were screened. Among them, 472 LUAD samples had complete survival time and status data. No statistically significant differences were found in the C0 and C2 subgroup scores between tumor and normal tissues(all P>0.05), whereas significant differences were observed in the other 9 subgroups(t=2.356-8.741, all P<0.05).(2) WGCNA analysis identified 6 co-expression modules, screening out 333 key genes from the most significantly C5 subgroup score-related modules.(3) Univariate Cox analysis of the 333 key genes using coxph function identified 25 prognosis-related genes. NMF clustering based on these 25 genes categorized the 472 tumor samples into L1 and L2 subtypes. Immune scoring using the ESTIMATE algorithm revealed a significantly higher immune score in the L2 subtype(2 318.75±25.63) compared to the L1 subtype(1 125.36±18.42)(t=10.234, P<0.001). Immune score was positively correlated with C5 subgroup score(r=0.613, P=0.001). L2 subtype patients had significantly higher proportions of T3-4 stage, N2-3 stage and clinical grade Ⅲ-Ⅳ, and a higher C5 subgroup score than L1 subtype patients(χ~2=6.325-15.782, all P<0.05). The 3-year overall survival rate was significantly lower in L2 subtype patients compared to L1 subtype patients(χ~2=12.453, P=0.001). Functional enrichment analysis identified the top 20 most significantly enriched pathways.(4) A total of 337 up-regulated and 401 down-regulated genes were identified in the L1 and L2 subtypes. Univariate Cox analysis identified 90 prognosis-related genes. LASSO regression with 10-fold cross-validation determined that 18 genes served as the most generalizable feature variables when the optimal λ value was 0.031 579 93. Subsequently, stepAIC regression analysis screened 6 core genes(NID2, ANGPTL4, AKAP12, GJB3, CCL20, and TPSB2). The risk score was calculated as follows: Risk Score=0.156×NID2 +0.089×ANGPTL4+0.091×AKAP12+0.121×GJB3+0.084×CCL20-0.212×TPSB2. ROC curve analysis showed that the AUCs of the molecular model for predicting 1-, 3-, and 5-year prognosis were 0.71(95%CI: 0.721-0.849, P=0.001), 0.71(95%CI: 0.753-0.871, P=0.001) and 0.65(95%CI: 0.730-0.862, P=0.001), respectively. Kaplan-Meier survival analysis revealed the 3-year overall survival rate was significantly lower in the high-risk group compared to the low-risk group(χ~2=18.642, P<0.001). Patients with high TPSB2 expression had a significantly higher 3-year overall survival rate than those with low expression(χ~2=8.765, P=0.002).(5) The risk score distribution analysis results showed that the samples in the high-risk group and low-risk group were distributed differently in the validation set, with death samples predominantly concentrated in the high-risk group. The ROC curve analysis results indicated that the molecular model’s AUC values for predicting 1-, 3-, and 5-year prognosis in the external validation set were 0.79(95%CI: 0.678-0.828, P=0.001), 0.62(95%CI: 0.705-0.847, P=0.001), and 0.73(95%CI: 0.689-0.835, P=0.001), respectively. The Kaplan-Meier survival analysis results demonstrated that the 5-year overall survival rate was lower in the high-risk group than that in the low-risk group(χ~2=10.325, P=0.003).(6) The relative protein expression levels of NID2, ANGPTL4, AKAP12, GJB3, and CCL20 in NCI-H1395, A549, and NCI-H1975 cells were all higher than those in WI-38 cells(t=4.559-25.770, all P<0.05), and the relative protein expression of TPSB2 was lower in NCI-H1395, A549, and NCI-H1975 cells than that in WI-38 cells(t=-17.668 to-5.516, all P<0.05).(7) At 48 h post-transfection, the counts of invading cells and EdU-positive cells were significantly lower in the NC-NCI-H1395 and NC-A549 groups than those in the sh-NCI-H1395 and sh-A549 groups(t=8.240-16.310, all P<0.05). Conclusions The expressions of NID2, ANGPTL4, AKAP12, GJB3, and CCL20 genes are up-regulated, while the expression of TPSB2 gene is down-regulated in LUAD cells. The molecular prognostic model including these 6 genes based on the single-cell transcriptome sequencing data and TCGA database has a good predictive performance for LUAD patients, and TPSB2 may be a potential target.更多还原