5.基于python的scRNA-seq细胞状态分析-差异表达

tzc_fly

于 2024-05-22 14:09:42 发布

阅读量654

点赞数 22

分类专栏：生物计算工具文章标签： python

本文链接：https://blog.csdn.net/qq_40943760/article/details/139105973

版权

生物计算工具专栏收录该内容

49 篇文章 62 订阅

订阅专栏

参考：
[1] https://github.com/Starlitnightly/single_cell_tutorial
[2] https://github.com/theislab/single-cell-best-practices

背景知识

在手动注释中，利用了一个先验知识：不同细胞类型特定的marker基因。现在我们反过来，当我们在面临新的分析场景时，假设已经得到了细胞类型，我们想利用单细胞的分辨率去分析更多的基因差异表达情况。这个过程也被称为差异分析。

差异分析方法最早是ttest，但是scRNA-seq和bulkRNA-seq相比，有更强的dropout，因此，一些样本会导致直接遍历每个细胞的差异分析不可靠，为了解决这个问题，我们利用pseudo-bulk方法来聚合单细胞数据，再进行差异分析（另一种是先插补，再ttest）。比如SEACell将单个细胞聚合成代表不同细胞状态的元细胞。

下面加载已经有标注的细胞数据，这里使用Kang，该数据集是来自8名红斑狼疮患者的10倍液滴式单细胞RNA测序技术测量的外周血单核细胞（PBMC）数据，在接受干扰素β（INF-β）治疗前后6小时进行采集（因此共16个样本）。这是大部分差异分析的应用场景：已知细胞类型，然后分析药物干预前后的基因差异表达。

包括预处理步骤和可视化：

import omicverse as ov
import scanpy as sc
import scvelo as scv

ov.utils.ov_plot_set()

adata = ov.read('./data/kang.h5ad')

#quantity control
adata=ov.pp.qc(adata,
              tresh={'mito_perc': 0.2, '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'])

print(adata)

ov.utils.embedding(adata,
                   basis='X_mde',
                    frameon='small',
                   color=['replicate','label','cell_type'])

"""
AnnData object with n_obs × n_vars = 24529 × 2000
    obs: 'nCount_RNA', 'nFeature_RNA', 'tsne1', 'tsne2', 'label', 'cluster', 'cell_type', 'replicate', 'nCount_SCT', 'nFeature_SCT', 'integrated_snn_res.0.4', 'seurat_clusters', 'nUMIs', 'mito_perc', 'detected_genes', 'cell_complexity', 'doublet_score', 'predicted_doublet', 'passing_mt', 'passing_nUMIs', 'passing_ngenes', 'n_genes'
    var: 'name', 'mt', 'n_cells', 'percent_cells', 'robust', 'mean', 'var', 'residual_variances', 'highly_variable_rank', 'highly_variable_features'
    uns: 'scrublet', 'log1p', 'hvg', 'scaled|original|pca_var_ratios', 'scaled|original|cum_sum_eigenvalues'
    obsm: 'X_pca', 'X_umap', 'scaled|original|X_pca', 'X_mde'
    varm: 'scaled|original|pca_loadings'
    layers: 'counts', 'scaled', 'lognorm'
"""

fig1

虽然看起来患者批次比较均衡，但在细胞类型的UMAP上，明显看到CD14+Monocytes具有严重批次效应，这肯定会干扰差异分析结果，因此，我们依然执行批次整合，方法使用scVI，批次标签为患者：

# 使用scVI整合
import scvi
adata1 = adata.copy()
scvi.model.SCVI.setup_anndata(adata1, layer="counts", batch_key="replicate")
vae = scvi.model.SCVI(adata1, n_layers=2, n_latent=30, gene_likelihood="nb")
vae.train()
adata.obsm["X_scVI"] = vae.get_latent_representation()

# 可视化
adata.obsm["X_mde_scVI"] = ov.utils.mde(adata.obsm["X_scVI"])
ov.utils.embedding(adata,
                   basis='X_mde_scVI',
                    frameon='small',
                   color=['replicate','label','cell_type'])

fig2
批次整合花费49m15.1s。

t-test差异分析

在下面的案例中，我们只分析CD4 T cells，即研究治疗前后CD4 T细胞的响应变化。

t-test是全细胞水平的差异分析方法，对于t-test而言，数据要求正态分布，所以我们使用对原始数据标准化后的数据：

test_adata=adata[adata.obs['cell_type'].isin(['CD4 T cells'])]
dds=ov.bulk.pyDEG(test_adata.to_df(layer='lognorm').T)
dds.drop_duplicates_index()
print('... drop_duplicates_index success')

选取label分别为Ctrl与Stim作为需要进行差异分析的两个组别：

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='ttest')

不指定差异表达倍数，由数据分布自动计算阈值，然后火山图可视化差异表达基因的分布：

# -1 means automatically calculates
dds.foldchange_set(fc_threshold=-1,
                   pval_threshold=0.05,
                   logp_max=10)

dds.plot_volcano(title='DEG Analysis',figsize=(4,4),
                 plot_genes_num=8,plot_genes_fontsize=12,)

fig3
或者使用箱线图来观察指定的差异表达基因在不同组的差异情况：

dds.plot_boxplot(genes=['NELL2','TNFSF13B'],treatment_groups=treatment_groups,
                control_groups=control_groups,figsize=(2,3),fontsize=12,
                 legend_bbox=(2,0.55))

fig4
也可以在UMAP中直接检验基因表达：

ov.utils.embedding(test_adata,
                   basis='X_mde_scVI',
                    frameon='small',
                   color=['label','NELL2','TNFSF13B'])

fig5

DEseq2差异分析

DEseq2广泛用于差异分析，不用考虑正态分布：

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()
print('... drop_duplicates_index success')

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')

# -1 means automatically calculates
dds.foldchange_set(fc_threshold=1.2,
                   pval_threshold=0.05,
                   logp_max=10)

dds.plot_volcano(title='DEG Analysis',figsize=(4,4),
                 plot_genes_num=8,plot_genes_fontsize=12,)

dds.plot_boxplot(genes=['IFI6','RGCC'],treatment_groups=treatment_groups,
                control_groups=control_groups,figsize=(2,3),fontsize=12,
                 legend_bbox=(2,0.55))

ov.utils.embedding(test_adata,
                   basis='X_mde_scVI',
                    frameon='small',
                   color=['label','IFI6','RGCC'])