核心内容摘要
探寻日本ZZJJ的魅力：一场穿越时空的文化盛宴
转换数据
1 数据转换概述在生物信息学数据分析中我们需要将AnnData数据转换为DataFrame数据DataFrame数据的行为位点细胞列为基因。
这种转换使得数据更容易处理和分析。
2 数据预处理流程数据转换包含以下几个关键步骤数据读取从H5AD格式文件中读取单细胞RNA测序数据数据标准化对数据进行归一化处理对数变换对数据进行对数变换减少数据偏态高变基因筛选筛选出变化较大的基因减少维度数据格式转换将稀疏矩阵转换为密集矩阵并转换为DataFrame格式
3 数据转换代码实现importscanpyasscimportnumpyasnpimportpandasaspdfromscipyimportstatsdeftransform_dataframe(h5ad_file): 将AnnData数据转换为DataFrame数据 参数: h5ad_file: H5AD格式文件路径 返回: adata: AnnData对象 result_dataframe: 转换后的DataFrame数据 # 读取H5AD文件adatasc.read_h5ad(h5ad_file)# 标准化数据使每个细胞的总reads数为10,000sc.pp.normalize_total(adata,target_sum1e
# 对数变换使用log1p函数对数据进行对数变换sc.pp.log1p(adata)# 识别高变基因sc.pp.highly_variable_genes(adata,n_top_genes
# 只保留高变基因adataadata[:,adata.var[highly_variable]]# 将稀疏矩阵转换为密集矩阵Xadata.X.toarray()# 获取索引和列名index_listadata.obs[ground_truth].tolist()columnadata.var.index.tolist()# 创建DataFrameresult_dataframepd.DataFrame(dataX,indexindex_list,columnscolumn)returnadata,result_dataframe# 使用示例adata,datatransform_dataframe(
h5ad)# 查看数据前10行print(data.head(
)# 查看特定细胞类型的数据print(data.loc[Layer1,:])
4 数据转换结果分析数据转换后的DataFrame具有以下特征行索引表示不同的细胞类型或样本如’Layer1’等列名表示基因名称如’LINC00115’、‘FAM41C’、SAMD11’等数据值表示基因在细胞中的表达量经过标准化和对数变换后数值范围相对较小数据维度817行×3000列表示817个细胞和3000个高变基因数据示例LINC00115 FAM41C SAMD11 PERM1 ISG15 AL
3
2 B3GALT6 Layer1
0
0
0
0
446558
0
000000 Layer1
0
0
0
0
000000
0
000000 ...
寻找差异表达基因
1 差异表达基因分析概述差异表达基因分析是生物信息学中的重要分析任务用于识别在不同条件下表达水平显著不同的基因。
通过分析差异表达基因我们可以理解细胞类型、疾病状态、处理条件等对基因表达的影响。
2 统计检验方法
2.
1 Wilcoxon秩和检验Wilcoxon秩和检验Mann-Whitney U检验是一种非参数检验方法用于比较两组独立样本的分布差异。
在生物信息学中Wilcoxon检验被广泛应用于差异表达基因分析原因包括不需要数据满足正态分布假设对异常值不敏感适用于样本量较小的情况适合处理基因表达数据这种偏态分布数据
3 底层计算函数下面设计的函数为嵌入最深的函数或者说是含有底层数据处理逻辑的函数。
该函数用于计算单个基因在目标组和背景组之间的差异表达情况。
defcompute(target_gene,background_gene,gene_name,cluster): 计算单个基因的差异表达统计量 参数: target_gene: 目标组的基因表达值 background_gene: 背景组的基因表达值 gene_name: 基因名称 cluster: 存储结果的字典 返回: None结果直接存储在cluster字典中 # 初始化基因信息cluster[gene_name]{}# 计算Wilcoxon秩和检验的p值stat,p_valuestats.mannwhitneyu(target_gene,background_gene,alternativetwo-sided)# 计算log2 fold changepseudo_count1# 避免除以零mean1np.mean(target_gene)pseudo_count mean2np.mean(background_gene)pseudo_count log_fcnp.log2(mean1/mean
# 存储结果cluster[gene_name][p_value]p_value cluster[gene_name][log_fc]log_fc# 使用示例cluster{}gene_namerand# 创建模拟数据target_genenp.random.normal(loc5,scale1,size(10,))background_genenp.random.randn(
# 计算差异表达compute(target_gene,background_gene,gene_name,cluster)# 查看结果print(cluster)输出示例{ rand: { p_value:
0631739856825891e-07, log_fc:
6310103643667038 } }
4 统计量解释
2.
1 p值p_value定义在零假设成立的情况下观察到当前统计量或更极端情况的概率意义p值越小说明两组数据差异越显著常用阈值p
05认为差异显著
2.
2 log2 fold changelog_fc定义目标组平均表达量与背景组平均表达量的比值的对数意义表示基因表达水平的倍数变化解释log_fc 0基因在目标组中上调log_fc 0基因在目标组中下调|log_fc| 1表达水平变化超过2倍
5 批量差异表达分析在底层处理函数的基础上往外走一步。
这一步要处理的是两个矩阵目标细胞类型的表达矩阵和背景细胞类型的表达矩阵对所有基因进行批量差异表达分析。
defbatch_de_analysis(data,target_cluster): 批量差异表达分析 参数: data: DataFrame格式的表达数据 target_cluster: 目标细胞类型名称 返回: results: 所有基因的差异表达结果字典 # 提取目标细胞类型的数据target_datadata.loc[data.indextarget_cluster,:]# 提取背景细胞类型的数据除目标细胞类型外的所有数据background_datadata.loc[data.index!target_cluster,:]# 初始化结果字典results{}# 对每个基因进行差异表达分析forgeneindata.columns:target_genetarget_data[gene].values background_genebackground_data[gene].values compute(target_gene,background_gene,gene,results)returnresults# 使用示例target_clusterLayer1resultsbatch_de_analysis(data,target_cluster)# 查看部分结果fori,(gene,stats)inenumerate(list(results.items())[:10]):print(f{gene}: p_value{stats[p_value]:.2e}, log_fc{stats[log_fc]:.4f})输出示例LINC00115: p_value
88e-01, log_fc
0042 FAM41C: p_value
43e-01, log_fc
0031 SAMD11: p_value
98e-01, log_fc
0039 PERM1: p_value
23e-01, log_fc-
0021 ISG15: p_value
22e-09, log_fc-
1414 AL
3
2: p_value
48e-01, log_fc
0002 B3GALT6: p_value
49e-02, log_fc-
0289 MXRA8: p_value
17e-01, log_fc
0238 VWA1: p_value
97e-02, log_fc-
0256 CDK11B: p_value
51e-04, log_fc-
0.
0
6 多重检验校正在进行大量基因的差异表达分析时由于进行了多次假设检验需要进行多重检验校正以控制假阳性率。
常用的校正方法包括
2.
1 Bonferroni校正fromstatsmodels.stats.multitestimportmultipletestsdefbonferroni_correction(results): 使用Bonferroni方法进行多重检验校正 参数: results: 差异表达结果字典 返回: results: 添加了校正后p值的结果字典 # 提取所有p值p_values[results[gene][p_value]forgeneinresults]# 进行Bonferroni校正reject,p_corrected,_,_multipletests(p_values,alpha
05,methodbonferroni)# 将校正后的p值添加到结果中fori,geneinenumerate(results):results[gene][p_corrected]p_corrected[i]results[gene][significant]reject[i]returnresults
2.
2 Benjamini-Hochberg校正FDRdeffdr_correction(results): 使用Benjamini-Hochberg方法进行FDR校正 参数: results: 差异表达结果字典 返回: results: 添加了FDR校正后p值的结果字典 # 提取所有p值p_values[results[gene][p_value]forgeneinresults]# 进行FDR校正reject,p_corrected,_,_multipletests(p_values,alpha
05,methodfdr_bh)# 将校正后的p值添加到结果中fori,geneinenumerate(results):results[gene][fdr]p_corrected[i]results[gene][significant]reject[i]returnresults
7 差异表达基因筛选根据统计检验结果和设定的阈值筛选出显著的差异表达基因。
deffilter_de_genes(results,p_threshold
05,logfc_threshold
1.
: 筛选差异表达基因 参数: results: 差异表达结果字典 p_threshold: p值阈值 logfc_threshold: log2 fold change阈值 返回: up_genes: 上调基因列表 down_genes: 下调基因列表 up_genes[]down_genes[]forgene,statsinresults.items():# 判断基因是否显著is_significantstats.get(significant,stats[p_value]p_threshold)ifis_significant:ifstats[log_fc]logfc_threshold:up_genes.append({gene:gene,p_value:stats[p_value],log_fc:stats[log_fc]})elifstats[log_fc]-logfc_threshold:down_genes.append({gene:gene,p_value:stats[p_value],log_fc:stats[log_fc]})# 按log_fc排序up_genes.sort(keylambdax:x[log_fc],reverseTrue)down_genes.sort(keylambdax:x[log_fc])returnup_genes,down_genes# 使用示例# 首先进行多重检验校正resultsfdr_correction(results)# 筛选差异表达基因up_genes,down_genesfilter_de_genes(results,p_threshold
05,logfc_threshold
0.
print(f上调基因数量:{len(up_genes)})print(f下调基因数量:{len(down_genes)})# 显示前10个上调基因print(\n前10个上调基因:)fori,gene_infoinenumerate(up_genes[:10],
:print(f{i}.{gene_info[gene]}: log_fc{gene_info[log_fc]:.4f}, p{gene_info[p_value]:.2e})
8 结果可视化
2.
1 火山图Volcano Plotimportmatplotlib.pyplotaspltimportseabornassnsdefplot_volcano(results,p_threshold
05,logfc_threshold
0.
: 绘制火山图 参数: results: 差异表达结果字典 p_threshold: p值阈值 logfc_threshold: log2 fold change阈值 # 准备数据geneslist(results.keys())log_fcs[results[gene][log_fc]forgeneingenes]p_values[-np.log10(results[gene][p_value])forgeneingenes]# 判断基因是否显著colors[]forgeneingenes:ifresults[gene][log_fc]logfc_thresholdandresults[gene][p_value]p_threshold:colors.append(red)# 上调elifresults[gene][log_fc]-logfc_thresholdandresults[gene][p_value]p_threshold:colors.append(blue)# 下调else:colors.append(gray)# 不显著# 绘制火山图plt.figure(figsize(10,
)plt.scatter(log_fcs,p_values,ccolors,alpha
0.
# 添加阈值线plt.axvline(xlogfc_threshold,colorred,linestyle--,alpha
0.
plt.axvline(x-logfc_threshold,colorblue,linestyle--,alpha
0.
plt.axhline(y-np.log10(p_threshold),colorgray,linestyle--,alpha
0.
# 添加标签plt.xlabel(log2 Fold Change)plt.ylabel(-log10(p_value))plt.title(Volcano Plot)# 添加图例frommatplotlib.linesimportLine2D legend_elements[Line2D([0],[0],markero,colorw,markerfacecolorred,labelUp-regulated),Line2D([0],[0],markero,colorw,markerfacecolorblue,labelDown-regulated),Line2D([0],[0],markero,colorw,markerfacecolorgray,labelNot significant)]plt.legend(handleslegend_elements)plt.tight_layout()plt.show()# 使用示例plot_volcano(results,p_threshold
05,logfc_threshold
0.
5)
2.
2 热图Heatmapdefplot_heatmap(data,top_genes,n_samples
: 绘制差异表达基因热图 参数: data: DataFrame格式的表达数据 top_genes: 差异表达基因列表 n_samples: 显示的样本数量 # 提取top基因的表达数据top_gene_names[gene[gene]forgeneintop_genes[:20]]heatmap_datadata[top_gene_names].head(n_samples)# 绘制热图plt.figure(figsize(12,
)sns.heatmap(heatmap_data,cmapviridis,cbar_kws{label:Expression Level})plt.title(Differentially Expressed Genes Heatmap)plt.xlabel(Genes)plt.ylabel(Samples)plt.tight_layout()plt.show()# 使用示例all_top_genesup_genes[:10]down_genes[:10]plot_heatmap(data,all_top_genes,n_samples
20)
完整分析流程
1 端到端分析流程将上述步骤整合为一个完整的分析流程defcomplete_de_analysis_pipeline(h5ad_file,target_cluster): 完整的差异表达分析流程 参数: h5ad_file: H5AD格式文件路径 target_cluster: 目标细胞类型名称 返回: results: 差异表达分析结果 up_genes: 上调基因列表 down_genes: 下调基因列表 #
数据转换print(Step 1: 数据转换...)adata,datatransform_dataframe(h5ad_file)print(f数据维度:{data.shape})#
批量差异表达分析print(\nStep 2: 批量差异表达分析...)resultsbatch_de_analysis(data,target_cluster)print(f分析基因数量:{len(results)})#
多重检验校正print(\nStep 3: 多重检验校正...)resultsfdr_correction(results)significant_countsum([1forrinresults.values()ifr.get(significant,False)])print(f显著基因数量:{significant_count})#
筛选差异表达基因print(\nStep 4: 筛选差异表达基因...)up_genes,down_genesfilter_de_genes(results,p_threshold
05,logfc_threshold
0.
print(f上调基因数量:{len(up_genes)})print(f下调基因数量:{len(down_genes)})#
结果可视化print(\nStep 5: 结果可视化...)plot_volcano(results,p_threshold
05,logfc_threshold
0.
all_top_genesup_genes[:10]down_genes[:10]plot_heatmap(data,all_top_genes,n_samples
returnresults,up_genes,down_genes# 使用示例results,up_genes,down_genescomplete_de_analysis_pipeline(
h5ad,Layer
1)
2 结果解释通过完整分析流程我们可以获得以下信息数据概况数据集的大小、细胞类型数量、基因数量等基本信息差异表达基因在不同细胞类型之间表达水平显著不同的基因统计显著性通过p值和FDR校正评估基因差异的统计学意义生物学意义通过log2 fold change评估基因表达变化的幅度可视化结果通过火山图和热图直观展示差异表达基因
4、
注意事项与最佳实践
1 数据质量控制在进行差异表达分析前需要进行严格的质量控制检查细胞质量、基因质量、批次效应等问题对低质量细胞和基因进行过滤
2 参数选择p值阈值通常设置为
05但可根据研究需求调整log2 fold change阈值通常设置为
5-
0表示
1.