文献复现：单细胞转录组中的免疫原性细胞死亡特征与101种机器学习组合 PMID:37275552 复现率达95成。 实例数据视频 代码主要内容介绍 在本研究中，使用ssGSEA算法为每个TCGA-KIRC样本获得一个ICD活性得分，作为后续WGCNA分析的表型数据。 使用在单细胞测序水平上识别的836个与ICD相关的DEGs，在移除异常样本后构建了一个共表达网络）。 选择最佳的软阈值power=7（R2=0.874）以确保一个无尺度的拓扑网络。 通过将最小模块基因计数设置为60和MEDissThres设置为0.25，共获得了四个模块。 发现表明，MEblue模块与整体RNA测序中的ICD得分强烈相关（cor=0.7）。 此外，蓝色模块的GS与MM散点图显示了显著的相关性（cor=0.89，p=3.9e57，这表明蓝色模块内的基因可能与免疫原性细胞死亡有功能上的关联。 火山图展示了在TCGA-KIRC整体RNA-seq中肿瘤和正常肾组织间的差异表达基因。 将蓝色模块中的164个基因与整体RNA-seq的DEGs取交集，最终鉴定出131个基因。 这些基因被命名为免疫原性细胞死亡相关基因（ICDRgenes）。 对ICDRgenes的GO富集分析显示在BP中包括T细胞激活、白细胞介导的免疫以及抗原处理和呈递；在CC中如内吞囊泡和MHC蛋白复合体；以及在MF中如酰胺结合、免疫受体活性和MHC蛋白复合体的显著富集。 随后，对131个ICDR基因进行了单因素Cox回归分析，识别出39个P值小于0.05的显著基因。 为了进一步构建和验证模型，将TCGA基因列表与外部数据集E-MTAB-1980进行交集，并提取了37个共有基因。 单因素Cox回归分析的结果以及这些基因之间的相互关系。

最近复现了一篇挺有意思的文献，核心是把免疫原性细胞死亡(ICD)特征从单细胞数据搬到bulk RNA测序场景，再套上101种机器学习组合做预后模型。整个过程堪称数据处理的极限拉扯，咱们直接上干货。

文献复现：单细胞转录组中的免疫原性细胞死亡特征与101种机器学习组合 PMID:37275552 复现率达95成。 实例数据视频 代码主要内容介绍 在本研究中，使用ssGSEA算法为每个TCGA-KIRC样本获得一个ICD活性得分，作为后续WGCNA分析的表型数据。 使用在单细胞测序水平上识别的836个与ICD相关的DEGs，在移除异常样本后构建了一个共表达网络）。 选择最佳的软阈值power=7（R2=0.874）以确保一个无尺度的拓扑网络。 通过将最小模块基因计数设置为60和MEDissThres设置为0.25，共获得了四个模块。 发现表明，MEblue模块与整体RNA测序中的ICD得分强烈相关（cor=0.7）。 此外，蓝色模块的GS与MM散点图显示了显著的相关性（cor=0.89，p=3.9e57，这表明蓝色模块内的基因可能与免疫原性细胞死亡有功能上的关联。 火山图展示了在TCGA-KIRC整体RNA-seq中肿瘤和正常肾组织间的差异表达基因。 将蓝色模块中的164个基因与整体RNA-seq的DEGs取交集，最终鉴定出131个基因。 这些基因被命名为免疫原性细胞死亡相关基因（ICDRgenes）。 对ICDRgenes的GO富集分析显示在BP中包括T细胞激活、白细胞介导的免疫以及抗原处理和呈递；在CC中如内吞囊泡和MHC蛋白复合体；以及在MF中如酰胺结合、免疫受体活性和MHC蛋白复合体的显著富集。 随后，对131个ICDR基因进行了单因素Cox回归分析，识别出39个P值小于0.05的显著基因。 为了进一步构建和验证模型，将TCGA基因列表与外部数据集E-MTAB-1980进行交集，并提取了37个共有基因。 单因素Cox回归分析的结果以及这些基因之间的相互关系。

数据炼金术第一步：用ssGSEA算法给TCGA-KIRC样本打ICD活性分。这里有个坑要注意——不同样本的基因长度校正。实际操作时用GSVA包的gsva()函数时记得开启kcdf="Poisson"参数：

library(GSVA)
icd_signature <- read.csv("ICD_genes.csv")$Gene
expr_matrix <- as.matrix(readRDS("TCGA_KIRC_expr.rds"))
ssgsea_score <- gsva(expr_matrix, list(ICD=icd_signature), 
                    method="ssgsea", kcdf="Poisson", parallel.sz=4)
pheno_data <- data.frame(Sample=colnames(expr_matrix), ICD_score=ssgsea_score[1,])

WGCNA模块挖掘这块最刺激。选power值的时候别迷信自动推荐，得手动验证无标度拓扑。用pickSoftThreshold函数跑完记得画个拟合曲线图，作者用的power=7确实在R²=0.87的位置最稳（图1）。构建网络时参数设置要够暴力：

library(WGCNA)
enableWGCNAThreads()
net <- blockwiseModules(expr_matrix, power=7, 
                       minModuleSize=60, 
                       mergeCutHeight=0.25,
                       numericLabels=TRUE)
module_colors <- labels2colors(net$colors)

模块基因关联分析才是重头戏。MEblue模块和ICD得分的皮尔逊相关系数飙到0.7，这可比普通转录组分析高出一大截。画基因显著性(GS)与模块成员(MM)的散点图时，用verboseScatterplot函数直接出带统计量的图：

blue_module <- module_colors == "blue"
gene_signif <- cor(t(expr_matrix[blue_module,]), pheno_data$ICD_score)
module_membership <- signedKME(expr_matrix[blue_module,], MEs[, "MEblue"])
verboseScatterplot(module_membership, gene_signif, 
                  col="steelblue", abline=TRUE)

基因交集操作需要点技巧。作者在TCGA的bulk RNA数据找到164个差异基因，和单细胞的836个DEGs取交集居然还剩131个。用VennDiagram包画韦恩图时记得调透明度：

library(VennDiagram)
venn.plot <- venn.diagram(
  x = list(TCGA_DEGs=deg_names, scDEGs=sc_deg_list),
  filename = NULL, fill=c("#1f78b4","#33a02c"),
  alpha=0.5, cat.cex=1.2)
grid.draw(venn.plot)

当看到GO富集结果里蹦出"抗原呈递"、"T细胞激活"这些词，就知道这波稳了。用clusterProfiler做富集时注意设置universe参数避免假阳性：

enrich_res <- enrichGO(gene = icdr_genes, 
                      OrgDb = org.Hs.eg.db,
                      universe = rownames(expr_matrix),
                      ont = "BP", 
                      pAdjustMethod = "BH")
dotplot(enrich_res, showCategory=15) + 
  theme(axis.text.y=element_text(size=8))

最后的Cox模型搭建才是大戏。39个显著基因经过数据集交集剩下37个，用glmnet做LASSO回归时要玩转交叉验证：

library(glmnet)
cvfit <- cv.glmnet(x=train_matrix, y=Surv(train_time,train_status),
                  family="cox", nfolds=10, alpha=1)
plot(cvfit) # 选lambda时看1SE规则
risk_score <- predict(cvfit, newx=test_matrix, s="lambda.1se")

整个流程跑下来，发现最大的挑战其实是数据清洗——原始TCGA数据里有5%的样本因批次效应被踢出局。不过当看到验证集上C-index冲到0.81时，只能说这波机器学习组合拳打得漂亮。