文献复现：单细胞转录组中的免疫原性细胞死亡特征与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回归分析的结果以及这些基因之间的相互关系。

最近成功复现了一篇PMID为37275552的文献，复现率高达95%，感觉收获满满，迫不及待来和大家分享一下其中的关键内容。我还准备了实例数据视频，大家可以更直观地感受整个过程。

代码主要内容介绍

1. 计算ICD活性得分

在本研究里，使用ssGSEA算法为每个TCGA-KIRC样本获取一个ICD活性得分，这个得分将作为后续WGCNA分析的表型数据。这里假设我们用Python的相关库来实现，示例代码可能长这样：

import pandas as pd
from ssGSEA_package import ssGSEA  # 假设存在这样一个包

# 读取样本数据
samples_data = pd.read_csv('TCGA-KIRC_samples.csv')
# 假设我们有一个预先定义好的基因集用于ssGSEA计算
icd_gene_set = ['gene1', 'gene2', 'gene3']  # 实际应该是完整的与ICD相关基因集

icd_scores = ssGSEA(samples_data, icd_gene_set)
icd_scores_df = pd.DataFrame(icd_scores, columns=['ICD_activity_score'])

这里通过ssGSEA算法，基于样本数据和定义好的基因集计算出了ICD活性得分，并保存到了DataFrame里。

2. 构建共表达网络

使用在单细胞测序水平上识别出的836个与ICD相关的DEGs（差异表达基因），在移除异常样本后构建共表达网络。

# 假设我们已经读取了836个与ICD相关的DEGs数据到deg_df
deg_df = pd.read_csv('ICD_related_DEGs.csv')
# 移除异常样本，这里假设异常样本的判断标准是某一列数据大于某个阈值
abnormal_threshold = 100
filtered_deg_df = deg_df[deg_df['specific_column'] <= abnormal_threshold]

3. 选择最佳软阈值

选择最佳的软阈值power = 7（R2 = 0.874）以确保一个无尺度的拓扑网络。在实际操作中，可能会通过循环不同的power值，并计算对应的R2来找到这个最佳值。

import networkx as nx
import numpy as np

# 构建网络相关操作，假设已经有构建好的表达矩阵expression_matrix
possible_powers = range(1, 10)
r2_values = []

for power in possible_powers:
    adjacency_matrix = calculate_adjacency_matrix(expression_matrix, power)  # 假设存在这样一个函数
    G = nx.from_numpy_array(adjacency_matrix)
    # 计算R2相关操作，这里简化表示
    r2 = calculate_R2(G)  # 假设存在这样一个函数
    r2_values.append(r2)

best_power_index = np.argmax(r2_values)
best_power = possible_powers[best_power_index]

通过这样的操作，我们就找到了能确保无尺度拓扑网络的最佳软阈值。

4. 模块获取

通过将最小模块基因计数设置为60和MEDissThres设置为0.25，共获得了四个模块。

# 假设已经有基于上述步骤构建好的网络G
from networkx.algorithms import community

communities_generator = community.greedy_modularity_communities(G, weight='weight')
communities = list(communities_generator)
filtered_communities = []

for community in communities:
    if len(community) >= 60:
        filtered_communities.append(community)

# 根据MEDissThres进一步筛选，这里简化处理
final_communities = []
for community in filtered_communities:
    # 假设存在计算MEDissThres相关函数
    medissthres_value = calculate_MEDissThres(community)  
    if medissthres_value <= 0.25:
        final_communities.append(community)

5. 模块与ICD得分相关性分析

发现表明，MEblue模块与整体RNA测序中的ICD得分强烈相关（cor = 0.7）。而且蓝色模块的GS与MM散点图显示了显著的相关性（cor = 0.89，p = 3.9e57），这暗示蓝色模块内的基因可能与免疫原性细胞死亡有功能上的关联。

import seaborn as sns
import matplotlib.pyplot as plt

# 假设已经有蓝色模块基因数据blue_module_genes，ICD得分数据icd_scores_df
blue_module_icd_correlation = blue_module_genes['expression'].corr(icd_scores_df['ICD_activity_score'])
print(f"MEblue模块与ICD得分相关性: {blue_module_icd_correlation}")

# 绘制GS与MM散点图
sns.scatterplot(x='GS', y='MM', data=blue_module_genes)
plt.title('蓝色模块GS与MM散点图')
plt.show()

6. 差异表达基因分析

火山图展示了在TCGA-KIRC整体RNA-seq中肿瘤和正常肾组织间的差异表达基因。将蓝色模块中的164个基因与整体RNA-seq的DEGs取交集，最终鉴定出131个基因，这些基因被命名为免疫原性细胞死亡相关基因（ICDRgenes）。

# 假设已经有蓝色模块基因列表blue_module_gene_list，整体RNA-seq的DEGs列表rna_seq_deg_list
icdr_genes = list(set(blue_module_gene_list).intersection(set(rna_seq_deg_list)))
print(f"鉴定出的ICDRgenes数量: {len(icdr_genes)}")

7. GO富集分析

对ICDRgenes的GO富集分析显示在BP（生物学过程）中包括T细胞激活、白细胞介导的免疫以及抗原处理和呈递；在CC（细胞组成）中如内吞囊泡和MHC蛋白复合体；以及在MF（分子功能）中如酰胺结合、免疫受体活性和MHC蛋白复合体的显著富集。

文献复现：单细胞转录组中的免疫原性细胞死亡特征与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回归分析的结果以及这些基因之间的相互关系。

这里可能会用到一些专门的富集分析工具包，比如clusterProfiler：

from clusterProfiler import enrichGO

# 假设icdr_genes是前面鉴定出的基因列表，并且已经转换为合适的格式
ego = enrichGO(gene=icdr_genes, OrgDb='org.Hs.eg.db', ont='BP', pAdjustMethod='BH')
print(ego)

8. 单因素Cox回归分析及模型构建

随后，对131个ICDR基因进行了单因素Cox回归分析，识别出39个P值小于0.05的显著基因。为了进一步构建和验证模型，将TCGA基因列表与外部数据集E - MTAB - 1980进行交集，并提取了37个共有基因。

import lifelines
from lifelines.statistics import logrank_test

# 假设已经有ICDR基因表达数据icdr_gene_expression，生存数据survival_data
cox_results = []
for gene in icdr_genes:
    gene_expression = icdr_gene_expression[gene]
    results = lifelines.CoxPHFitter().fit(pd.DataFrame({'gene_expression': gene_expression,'survival_time': survival_data['time'],'survival_status': survival_data['status']}), 'gene_expression')
    if results.p_value < 0.05:
        cox_results.append(gene)

print(f"P值小于0.05的显著基因数量: {len(cox_results)}")

# 假设已经有TCGA基因列表tcga_gene_list，外部数据集E - MTAB - 1980基因列表external_gene_list
common_genes = list(set(tcga_gene_list).intersection(set(external_gene_list)))
print(f"共有的基因数量: {len(common_genes)}")

通过这些步骤，我们基本完成了文献中的关键分析过程，整个复现过程虽然复杂，但每一步都充满了挑战与乐趣，希望对大家有所启发。