| 发布日期：2024-08-09 03:47 点击次数：104 |
图片软件开发价格
作家按
本章节主要培植了单细胞数据的互异抒发基因分析关节，详备对比了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两种模子，两种模子的扫尾我将差别演示：
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/
图片
本站仅提供存储劳动，系数践诺均由用户发布，如发现存害或侵权践诺，请点击举报。