发布日期：2024-08-24 09:30 点击次数：59
图片希威社 姐妹花希威社 姐妹花
作家按
本章节主要老师了单细胞数据的各别抒发基因分析门径，详备对比了ttest与DEseq2在整个细胞进行各别抒发分析的异同，以及为什么要使用元细胞分析的原因。本教程首发于单细胞最好的华文教程[1]，未经授权许可，退却转载。
全翰墨数|瞻望阅读时间: 3000|5min
——Starlitnightly（星夜）
1. 配景咱们在前边留心的章节中，筹划了不同细胞的特异性marker（记号）基因，但好多时候，咱们更颐养在某一类细胞中，两种不同景色下的组别各别，举例药物调养与未经药物调养，肿瘤细胞与正常细胞（癌旁）细胞等。因此，咱们但愿能在单细胞水平上，进行各别抒发分析。
图片
各别抒发基因模子单细胞RNA-seq与批量RNA-seq数据比拟，有着更强的寥落性，举例dropout表象的宽绰出现，这使得咱们若是只是将每一个细胞视作样本进行各别抒发分析，由于权贵性熟习的P值关于样本量的敏锐度很高，因此若是使用沿途细胞进行各别抒发分析，那么分析末端中的权贵性各别可能是不成靠的。咱们鄙人面的对比实验会进一步评释这种表象。为了处分单细胞寥落性与样本过大的问题，咱们不错接受一种叫作念伪Bulk（pseudo-bulk）的门径来团聚单细胞数据，从而进行各别抒发分析。但这又引入另一个问题，最近的一项筹划强调了伪Bulk的问题，其中推行统计期骗于统计上不独处的生物复制。未能探讨类似（来自归并个体的细胞）的内在关系性会加多无剪发现率（FDR）。因此，应在各别抒发分析之前，应先期骗批量效应纠正或通过每个个体的总数、平均或立时效应（即伪Bulk生成）对个体内的细胞类型特异性抒发值进行团聚，以解说样本内关系性。
元细胞的见识（代表不同细胞景色的细胞组，其中元细胞内的变异是由于技能而非生物开端）被冷落手脚保握统计着力同期最大化灵验数据分辨率的一种方式。与伪Bulk不同，SEACells 寻求以与数据模态无关的方式将单个细胞团聚成代表不同细胞景色的元细胞。使用计数矩阵手脚输入，它提供每个元单位的每个单位权重、每个元单位的每个单位硬分拨以及每个元单位的团聚计数手脚输出，故在本教程中，咱们将展示怎样使用SEACells完成各别抒发分析。
import omicverse as ovimport scanpy as scimport scvelo as scvov.utils.ov_plot_set()
____ _ _ __ / __ \____ ___ (_)___| | / /__ _____________ / / / / __ `__ / / ___/ | / / _ / ___/ ___/ _ \ / /_/ / / / / / / / /__ | |/ / __/ / (__ ) __/ \____/_/ /_/ /_/_/\___/ |___/\___/_/ /____/\___/ Version: 1.5.5, Tutorials: https://omicverse.readthedocs.io/2. 加载数据
在这里，咱们使用pertpy的演示数据Kang 数据集，它是来自 8 名狼疮患者 INF-β 调养前后 6 小时前后的 10 倍基于 scRNA-seq 的外周血单核细胞 (PBMC) 数据（所有这个词 16 个样本）
下载地址：https://figshare.com/ndownloader/files/34464122
adata = ov.read('data/kang.h5ad')adata# AnnData object with n_obs × n_vars = 24673 × 15706 # obs: 'nCount_RNA', 'nFeature_RNA', 'tsne1', 'tsne2', 'label', 'cluster', 'cell_type', 'replicate', 'nCount_SCT', 'nFeature_SCT', 'integrated_snn_res.0.4', 'seurat_clusters' # var: 'name' # obsm: 'X_pca', 'X_umap'
adata.X.max()
与别的门径不雷同，咱们在这里使用ttest门径计较各别抒发基因，因此咱们需要使用log对数化处理后的数据，因此咱们对数据进行预处理
#quantity controladata=ov.pp.qc(adata, tresh={'mito_perc': 0.2, 'nUMIs': 500, 'detected_genes': 250})#normalize and high variable genes (HVGs) calculatedadata=ov.pp.preprocess(adata,mode='shiftlog|pearson',n_HVGs=2000,)#save the whole genes and filter the non-HVGsadata.raw = adataadata = adata[:, adata.var.highly_variable_features]#scale the adata.Xov.pp.scale(adata)#Dimensionality Reductionov.pp.pca(adata,layer='scaled',n_pcs=50)adata.obsm['X_mde']=ov.utils.mde(adata.obsm['scaled|original|X_pca'])ov.utils.embedding(adata, basis='X_mde', frameon='small', color=['replicate','label','cell_type'])
图片
降维可视化图3. 批次效应的纠正咱们发现数据存在着一定的批次效应，咱们这里使用scVU进行批次效应的纠正。虽然，你也不错使用别的门径，如harmony等，更多的批次效应纠正的门径见omicverse的教程：https://omicverse.readthedocs.io/en/latest/Tutorials-single/t_single_batch/
adata=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"])
# ...Begin using scVI to correct batch effect# GPU available: True (cuda), used: True# TPU available: False, using: 0 TPU cores# IPU available: False, using: 0 IPUs# HPU available: False, using: 0 HPUs# LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]# Epoch 326/326: 100%|██████████| 326/326 [13:22<00:00, 2.47s/it, loss=599, v_num=1] # `Trainer.fit` stopped: `max_epochs=326` reached.# Epoch 326/326: 100%|██████████| 326/326 [13:22<00:00, 2.46s/it, loss=599, v_num=1]# ```python# ov.utils.embedding(adata,# basis='X_mde_scVI',# frameon='small',# color=['replicate','label','cell_type'])
图片
批次效应纠正图4. 各别抒发分析（全细胞水平）在使用元细胞分析前，我想先使用全细胞水平进行分析，以展示全细胞水平下的各别抒发分析恶果，在教程中，咱们只分析CD4 T cells的各别抒发基因，但愿能通过CD4 T cells的各别来标明调养前后免疫应酬反映的变化。
咱们这里不错使用ttest或者是DESeq2两种模子，两种模子的恶果我将诀别演示：
bt核工厂发布器4.1 基于ttest进行各别抒发分析关于ttest统计门径而言，数据要求倨傲正态漫衍的基本形态。中心极纵容理，是指概率论中征询立时变量和漫衍渐近于正态漫衍的定理，是数理统计学和邪恶分析的表面基础，指出了多数立时变量近似战胜正态漫衍的要求。故土们使用原数据的对数圭臬化后的实践，来进行各别抒发分析。
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')# ... 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 calculatesdds.foldchange_set(fc_threshold=-1, pval_threshold=0.05, logp_max=10)
# ... Fold change threshold: 0.8351151943206787
咱们使用自带的火山图函数plot_volcano来不雅察各别抒发基因的漫衍
dds.plot_volcano(title='DEG Analysis',figsize=(4,4), plot_genes_num=8,plot_genes_fontsize=12,)
图片
火山图1咱们还不错使用自带的箱线图函数plot_boxplot来不雅察各别抒发基因在不同组的各别情况
dds.plot_boxplot(genes=['NELL2','TNFSF13B'],treatment_groups=treatment_groups, control_groups=control_groups,figsize=(2,3),fontsize=12, legend_bbox=(2,0.55))
图片
箱线图1除此除外，咱们还不错在低维流型图umap中，不雅察各别抒发基因在不同组的漫衍情况
ov.utils.embedding(test_adata, basis='X_mde_scVI', frameon='small', color=['label','NELL2','TNFSF13B'])
图片
降维图14.2 基于DEseq2进行各别抒发分析DEseq2被过去使用在各别抒发分析的各式时势，收货于其负二项漫衍模子的筹算。使得咱们不错使用原始计数Counts进活动直分析。DESeq2将对原始reads进行建模，使用圭臬化因子(scale factor)来解说库深度的各别。然后，DESeq2预计基因的闹翻度，并疲塌这些预计值以生成更准确的闹翻度预计，从而对reads count进行建模。终末，DESeq2拟合负二项漫衍的模子，并使用Wald熟习或似然比熟习进行假定熟习。
与ttest不同，使用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')# ... drop_duplicates_index success # Fitting size factors... # ... done in 0.20 seconds. # Fitting dispersions... # ... done in 3.69 seconds. # Fitting dispersion trend curve... # ... done in 0.60 seconds. # logres_prior=nan, sigma_prior=nan # Fitting MAP dispersions... # ... done in 73.18 seconds. # Fitting LFCs... # ... done in 299.53 seconds. # Refitting 5 outliers. # Fitting dispersions... # ... done in 0.09 seconds. # Fitting MAP dispersions... # ... done in 0.23 seconds. # Fitting LFCs... # ... done in 0.70 seconds. # Running Wald tests... # ... done in 71.47 seconds. # Log2 fold change & Wald test p-value: condition Treatment vs Control
# -1 means automatically calculatesdds.foldchange_set(fc_threshold=1.2, pval_threshold=0.05, logp_max=10)
#... Fold change threshold: 1.2
dds.plot_volcano(title='DEG Analysis',figsize=(4,4), plot_genes_num=8,plot_genes_fontsize=12,)
图片
火山图2dds.plot_boxplot(genes=['IFI6','RGCC'],treatment_groups=treatment_groups, control_groups=control_groups,figsize=(2,3),fontsize=12, legend_bbox=(2,0.55))
图片
箱线图2ov.utils.embedding(test_adata, basis='X_mde_scVI', frameon='small', color=['label','IFI6','RGCC'])
图片
降维图25. 基于元细胞的各别抒发分析咱们不错看到，不管是使用ttest与DEseq2，齐不错一定经由上计较出两个组别的细胞各别抒发基因，况兼也确乎在两个组别中竣事了分离。可是，这也带来了多数的生物学杂音。两种细胞并莫得竣事荒谬无缺的分离，各别抒发基因的抒发在两类细胞中也有着一定的丰采。
为了处分这一问题，扩大样本（细胞）的深度，咱们使用SEACells团聚细胞，然后在元细胞水平上，奉行各别抒发分析。值得一提的是，SEACells团聚的元细胞在信号团聚和细胞分辨率之间竣事了最好均衡，况兼它们拿获通盘表型谱中的细胞景色，包括陌生景色。元细胞保留了样本之间高明的生物学各别，这些各别通过替代门径手脚批次效应被摒除，因此，为数据集成提供了比寥落单个细胞更好的起初。
5.1 元细胞驱动化在这里，咱们使用scVI手脚元细胞特征向量的驱动化，你也不错使用X_pca或者是omicverse计较获取的scaled|original|X_pca来手脚元细胞的驱动输入。
meta_obj=ov.single.MetaCell(adata,use_rep='X_scVI',n_metacells=250, use_gpu=True)
# Welcome to SEACells GPU! # Computing kNN graph using scanpy NN ... # computing neighbors # finished: added to `.uns['neighbors']` # `.obsp['distances']`, distances for each pair of neighbors # `.obsp['connectivities']`, weighted adjacency matrix (0:00:03) # Computing radius for adaptive bandwidth kernel... # 0%| | 0/24528 [00:00<?, ?it/s] # Making graph symmetric... # Parameter graph_construction = union being used to build KNN graph... # Computing RBF kernel... # 0%| | 0/24528 [00:00<?, ?it/s] # Building similarity LIL matrix... # 0%| | 0/24528 [00:00<?, ?it/s] # Constructing CSR matrix...
咱们需要驱动化元细胞的原型，也即是其在特征空间scVI上的驱动位置，这里可能有一些难相识，若是咱们在umap图上来相识的话，不错示意为umap图上的元细胞的位置，不详会愈加平凡一些。
meta_obj.initialize_archetypes()
# Building kernel on X_scVI # Computing diffusion components from X_scVI for waypoint initialization ... # computing neighbors # finished: added to `.uns['neighbors']` # `.obsp['distances']`, distances for each pair of neighbors # `.obsp['connectivities']`, weighted adjacency matrix (0:00:02) # Done. # Sampling waypoints ... # Done. # Selecting 239 cells from waypoint initialization. # Initializing residual matrix using greedy column selection # Initializing f and g... # 100%|██████████| 21/21 [00:00<00:00, 42.44it/s] # Selecting 11 cells from greedy initialization.5.2 元细胞模子熟习
驱动完模子结构后，咱们就需要对元细胞进行团聚，在这里咱们通过熟习SEACells模子来完成这一主见
meta_obj.train(min_iter=10, max_iter=50)
# Randomly initialized A matrix. # Setting convergence threshold at 0.28753 # Starting iteration 1. # Completed iteration 1. # Converged after 8 iterations. # Converged after 9 iterations. # Starting iteration 10. # Completed iteration 10. # Converged after 10 iterations.
咱们提供了一个保存与读取函数，幸免类似计较
meta_obj.save('seacells/model.pkl')meta_obj.load('seacells/model.pkl')5.3 预测元细胞咱们不错使用 predicted 来预测原始 scRNA-seq 数据的元胞。有两种门径可供礼聘，一种是 "soft"门径，另一种是 "hard"门径。
在 "soft"门径中，汇总每个 SEACell 中的细胞，对整个原始数据乞降 x 为属于一个 SEACell 的整个细胞分拨权重。数据未圭臬化，伪原始团聚计数存储在 .layer['raw']中。与变量（.var）关系的属性会被复制过来，但每个 SEACell 的关系属性必须手动复制，因为某些属性可能需要乞降或取平均值等，具体取决于属性。
在 "hard"门径中，汇总每个 SEACell 中的单位格，对属于一个 SEACell 的整个单位格的整个原始数据乞降。数据未经圭臬化处理，原始汇共计数存储在 .layer['raw']中。与变量（.var）关系的属性会被复制过来，但每个 SEACell 的关系属性必须手动复制，因为某些属性可能需要乞降或取平均值等，具体取决于属性。
咱们最先需要串联细胞类型(cell_type)与组别(label)的标签用于预测。
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='lognorm')
# 100%|██████████| 250/250 [01:03<00:00, 3.94it/s]
ad
# AnnData object with n_obs × n_vars = 250 × 2000 # obs: 'Pseudo-sizes', 'celltype', 'celltype_purity'
import matplotlib.pyplot as pltfig, ax = plt.subplots(figsize=(4,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', )
图片
元细胞图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)
# If you pass `n_top_genes`, all cutoffs are ignored. # extracting highly variable genes # finished (0:00:00) # --> added # 'highly_variable', boolean vector (adata.var) # 'means', float vector (adata.var) # 'dispersions', float vector (adata.var) # 'dispersions_norm', float vector (adata.var) # computing PCA # on highly variable genes # with n_comps=50 # finished (0:00:00) # computing neighbors # finished: added to `.uns['neighbors']` # `.obsp['distances']`, distances for each pair of neighbors # `.obsp['connectivities']`, weighted adjacency matrix (0:00:00) # computing UMAP # finished: added # 'X_umap', UMAP coordinates (adata.obsm) (0:00:00)
import matplotlib.pyplot as pltfrom matplotlib import patheffectsfig, ax = 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=ax, 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=ax, adjust_kwargs=dict(arrowprops=dict(arrowstyle='-', color='black')), text_kwargs=dict(fontsize= 12 ,weight='bold', path_effects=[patheffects.withStroke(linewidth=2, foreground='w')] ),)图片
元细胞可视化test_adata=ad[ad.obs['celltype'].isin(['CD4 T cells-stim','CD4 T cells-ctrl'])]dds_meta=ov.bulk.pyDEG(test_adata.to_df().T)dds_meta.drop_duplicates_index()print('... drop_duplicates_index success')treatment_groups=test_adata.obs[test_adata.obs['celltype']=='CD4 T cells-stim'].index.tolist()control_groups=test_adata.obs[test_adata.obs['celltype']=='CD4 T cells-ctrl'].index.tolist()dds_meta.deg_analysis(treatment_groups,control_groups,method='ttest')# -1 means automatically calculatesdds_meta.foldchange_set(fc_threshold=-1, pval_threshold=0.05, logp_max=10)# ... drop_duplicates_index success # ... Fold change threshold: 1.4815978396550413
dds_meta.plot_volcano(title='DEG Analysis',figsize=(4,4), plot_genes_num=8,plot_genes_fontsize=12,)
图片
火山图3dds_meta.plot_boxplot(genes=['IFIT1','OSCP1'],treatment_groups=treatment_groups, control_groups=control_groups,figsize=(2,3),fontsize=12, legend_bbox=(2,0.55))
图片
箱线图3ov.utils.embedding(test_adata, basis='X_umap', frameon='small', color=['IFIT1','OSCP1','celltype'], cmap='RdBu_r')
图片
降维图3ov.utils.embedding(adata, basis='X_mde_scVI', frameon='small', color=['label','IFIT1','OSCP1'])
图片
降维图4咱们发现，IFIT1和OSCP1诀别在stim和ctrl组中特异性下调与上调，其中IFIT1编码一种含有四肽类似序列的卵白质，领先被浮滑为骚动素调养指导的。元细胞深切出了更优的各别基因的识别上风，况兼打消了dropout的骚动，咱们对比箱线图可获取沟通的不雅点。
test_adata=adata[adata.obs['cell_type'].isin(['CD4 T cells'])]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.plot_boxplot(genes=['IFIT1','OSCP1'],treatment_groups=treatment_groups, control_groups=control_groups,figsize=(2,3),fontsize=12, legend_bbox=(2,0.55))
图片
箱线图45. 回想咱们在对单细胞数据进行各别抒发分析的时候，不错从全细胞和元细胞两个角度去探讨，时常来说为了幸免生物学杂音和dropout的骚动，咱们会礼聘元细胞手脚分析战术来进行。但更多时候，咱们会同期使用两种门径，并进行互相印证，来标明分析的准确性与可靠性。
6. 想考咱们为什么要使用元细胞来进行各别抒发分析？ttest与DESeq2的各别抒发分析末端有什么不同？文中一语气[1]单细胞最好的华文教程: https://single-cell-tutorial.readthedocs.io/zh/latest/
图片
本站仅提供存储劳动，整个实践均由用户发布，如发现存害或侵权实践，请点击举报。