scRNA——差异基因表达分析

最新推荐文章于 2025-03-04 17:17:18 发布

以待成追忆

最新推荐文章于 2025-03-04 17:17:18 发布

阅读量1.3k

点赞数 6

分类专栏： scRNA 文章标签： python 数据分析

本文链接：https://blog.youkuaiyun.com/qq_24680545/article/details/145204420

版权

scRNA 专栏收录该内容

8 篇文章

订阅专栏

单细胞RNA-seq由于自身测序特点的缘故，与RNA-seq数据相比，有明显的稀疏性（漏检），既是优势也是缺点，在做差异基因表达分析这方面就是缺点。为了尽量避免这个问题出现，引入了伪Bulk（pseudo-bulk）的方法来聚合单细胞数据，从而进行差异表达分析，差异表达分析之前，通过批量效应校正或通过每个个体的总和、平均或随机效应（即伪Bulk生成）对个体内的细胞类型特异性表达值进行聚合，以解释样本内相关性。

预处理

这里使用的是Kang 数据集，来自 8 名狼疮患者 INF-β 治疗前后 6 小时前后的 10 倍基于 scRNA-seq 的外周血单核细胞 (PBMC) 数据（总共 16 个样本），
下载地址：Kang.h5ad

adata = sc.read(input_path)
adata = ov.pp.qc(adata, tresh={'mito_perc': 0.05, 'nUMIs': 500, 'detected_genes': 250})
# normalize and high variable genes (HVGs) calculated
adata = ov.pp.preprocess(adata, mode='shiftlog|pearson', n_HVGs=2000, )

# save the whole genes and filter the non-HVGs
adata.raw = adata
adata = adata[:, adata.var.highly_variable_features]

# scale the adata.X
ov.pp.scale(adata)
# Dimensionality Reduction
ov.pp.pca(adata, layer='scaled', n_pcs=50)
adata.obsm['X_mde'] = ov.utils.mde(adata.obsm['scaled|original|X_pca'])
model = ov.single.batch_correction(adata, batch_key='replicate',
                                           methods='scVI', n_layers=2, n_latent=30, gene_likelihood="nb")
adata.obsm["X_mde_scVI"] = ov.utils.mde(adata.obsm["X_scVI"])

差异表达分析（全细胞水平）

暂时分析CD4 T cells的差异表达基因，当然也可以是别的类型的细胞。差异表达分析有对照组和实验组，那么这个在该数据集的label有记录，['ctrl', 'stim']

test_adata = adata[adata.obs['cell_type'].isin(['CD4 T cells'])]
dds = ov.bulk.pyDEG(test_adata.to_df(layer='counts').T)
dds.drop_duplicates_index()
treatment_groups = test_adata.obs[test_adata.obs['label'] == 'stim'].index.tolist()
control_groups = test_adata.obs[test_adata.obs['label'] == 'ctrl'].index.tolist()
dds.deg_analysis(treatment_groups, control_groups, method='DEseq2')
dds.foldchange_set(fc_threshold=1.2,pval_threshold=0.05,logp_max=10)

计算结果如下：

Log2 fold change & Wald test p-value: condition Treatment vs Control
                baseMean  log2FoldChange     lfcSE          stat  pvalue  padj
index                                                                         
MALAT1         99.966697   -6.098812e-01  0.012706 -4.800122e+01     0.0   0.0
B2M            38.339137    2.885976e-11  0.004293  6.721816e-09     1.0   1.0
TMSB4X         41.790238   -5.382882e-01  0.012457 -4.321212e+01     0.0   0.0
FTH1           31.313126   -1.255170e+00  0.029690 -4.227620e+01     0.0   0.0
RPL13          15.451456   -7.726909e-01  0.015888 -4.863298e+01     0.0   0.0
...                  ...             ...       ...           ...     ...   ...
RP11-102N12.3   0.000000             NaN       NaN           NaN     NaN   NaN
RP11-297B17.3   0.000000             NaN       NaN           NaN     NaN   NaN
CCNJL           0.000000             NaN       NaN           NaN     NaN   NaN
RP5-887A10.1    0.000000             NaN       NaN           NaN     NaN   NaN
RP11-407H12.8   0.000000             NaN       NaN           NaN     NaN   NaN

火山图

DEG_Analysis = dds.plot_volcano(title='DEG Analysis', figsize=(6, 8),
                     plot_genes_num=8, plot_genes_fontsize=12, )
DEG_Analysis.figure.savefig('DEG_Analysis.png', dpi=300, bbox_inches='tight')

挤在一起不太好分辨！！！，可以控制figsize的大小来调整图片大小，本来是figsize=(4, 4)
在这里插入图片描述

箱线图

dds.plot_boxplot(genes=['IFI6','RGCC'],treatment_groups=treatment_groups,
                control_groups=control_groups,figsize=(2,3),fontsize=12,
                 legend_bbox=(2,0.55))
#要保存就同上抄一下

在这里插入图片描述
这里可以对指定基因绘制embedding，这里指定IFI6和RGCC，当然也可以是别的

dds_embed = ov.utils.embedding(test_adata,
                       basis='X_mde_scVI',
                       frameon='small',
                       color=['label', 'IFI6', 'RGCC'],
                       show=False)
dds_embed[0].figure.savefig('dds_embed.png', dpi=300, bbox_inches='tight')

在这里插入图片描述

通路图

在omicverse中，使用的是gseapy来计算通路，包括wiki，go，kegg等等通路，可以从Enrichr知晓适配哪些通路。
以wiki为例

ov.utils.download_pathway_database()
pathway_dict = ov.utils.geneset_prepare('genesets/WikiPathways_2024_Human.txt', organism='Human')
# 替换文件就行
rnk = dds.ranking2gsea()
gsea_obj = ov.bulk.pyGSEA(rnk, pathway_dict)
enrich_res = gsea_obj.enrichment()
gsea_fig = gsea_obj.plot_enrichment(num=10, node_size=[10, 20, 30],
                             cax_fontsize=10,
                             fig_title='Wiki Pathway Enrichment', 
                             fig_xlabel='Fractions of genes',
                             figsize=(2, 4), cmap='YlGnBu',
                             text_knock=2, text_maxsize=30,
                             cax_loc=[2.0, 0.45, 0.5, 0.02],node_diameter=6)
gsea_fig.figure.savefig('gsea_fig.png', dpi=300, bbox_inches='tight')

在这里插入图片描述

差异表达分析（元细胞水平）

基于细胞水平分析的带有大量的生物学噪音，两种细胞并没有实现非常完美的分离，差异表达基因的表达在两类细胞中也有着一定的丰度。

为了解决这一问题，扩大样本（细胞）的深度，SEACells聚合细胞，然后在元细胞水平上，执行差异表达分析。

SEACells聚合的元细胞在信号聚合和细胞分辨率之间实现了最佳平衡，并且它们捕获整个表型谱中的细胞状态，包括罕见状态。元细胞保留了样本之间微妙的生物学差异，这些差异通过替代方法作为批次效应被消除，因此，为数据集成提供了比稀疏单个细胞更好的起点。

在这里，我们使用scVI作为元细胞特征向量的初始化，你也可以使用X_pca或者是omicverse计算得到的scaled|original|X_pca来作为元细胞的初始输入，并初始化训练。

import omicverse as ov
import matplotlib.pyplot as plt
import pandas as pd
meta_obj = ov.single.MetaCell(adata, use_rep='X_scVI', n_metacells=250,
                                  use_gpu=True)
meta_obj.initialize_archetypes()
meta_obj.train(min_iter=10, max_iter=50)

# 有的时候训练时间比较长，可以将他们保存到本地，然后再读取。
meta_obj.save('model.pkl')
meta_obj.load('model.pkl')

训练完成之后就进行预测。

meta_obj.adata.obs['celltype-label'] = [i + '-' + j for i, j in zip(meta_obj.adata.obs['cell_type'],
                                                                        meta_obj.adata.obs['label'])]
ad = meta_obj.predicted(method='soft', celltype_label='celltype-label',
                            summarize_layer='counts')
fig, ax = plt.subplots(figsize=(8, 4))
ov.utils.embedding(
        meta_obj.adata,
        basis="X_mde_scVI",
        color=['cell_type'],
        frameon='small',
        title="Meta cells",
        # legend_loc='on data',
        legend_fontsize=14,
        legend_fontoutline=2,
        size=10,
        ax=ax,
        alpha=0.2,
        # legend_loc='',
        add_outline=False,
        # add_outline=True,
        outline_color='black',
        outline_width=1,
        show=False,
        # palette=ov.utils.blue_color[:],
        # legend_fontweight='normal'
    )
ov.single._metacell.plot_metacells(ax, meta_obj.adata,use_rep='X_mde_scVI', color='#CB3E35',)
fig.savefig("output_image.png", dpi=300, bbox_inches="tight")
plt.tight_layout()
plt.show()

在这里插入图片描述

sc.pp.log1p(ad)
sc.pp.highly_variable_genes(ad, n_top_genes=1500, inplace=True)
sc.tl.pca(ad, use_highly_variable=True)
sc.pp.neighbors(ad, use_rep='X_pca')
sc.tl.umap(ad)
fig2, ax2 = plt.subplots(figsize=(4, 4))
ov.utils.embedding(
        ad,
        basis="X_umap",
        color=['celltype'],
        frameon='small',
        title="metacells celltypes",
        # legend_loc='on data',
        legend_fontsize=14,
        legend_fontoutline=2,
        # size=10,
        ax=ax2,
        legend_loc=None, add_outline=False,
        # add_outline=True,
        outline_color='black',
        outline_width=1,
        show=False,
        # palette=ov.palette()[12:],
        # legend_fontweight='normal'
    )
ov.utils.gen_mpl_labels(
        ad,
        'celltype',
        exclude=("None",),
        basis='X_umap',
        ax=ax2,
        adjust_kwargs=dict(arrowprops=dict(arrowstyle='-', color='black')),
        text_kwargs=dict(fontsize=12, weight='bold',
                         path_effects=[patheffects.withStroke(linewidth=2, foreground='w')]),
    )

fig2.savefig("output_image2.png", dpi=300, bbox_inches="tight")
plt.tight_layout()
plt.show()