Nature methods|scPerturb:单细胞数据扰动怎么办？看看这篇文章或许你就有思路

忙活了一周，终于有一个下午的时间抽空出来看看文献，写写总结了，今天分享的是单细胞数据扰动的处理，由于数据互操作性差，对越来越多的单细胞扰动数据集的分析受到阻碍。为了促进计算方法的开发和基准测试，该篇文献收集了一组 44 个公开可用的单细胞扰动响应数据集，其中包含分子读数，包括转录组学、蛋白质组学和表观基因组学。文献中应用统一的质量控制流程并协调特征注释。由此产生的信息资源 scPerturb 可以开发和测试计算方法，并促进跨数据集的比较和集成。该文献描述了用于量化扰动效应和显著性检验的能量统计（E 统计），并展示了 E 距离作为单细胞表达谱组之间的一般距离测量。最后也说明了 E 统计在量化扰动的相似性和功效中的应用。
文章地址：https://www.nature.com/articles/s41592-023-02144-y
代码地址：https://github.com/sanderlab/scPerturb
扰动实验

扰动实验探测细胞或细胞系统对条件变化的反应。传统的方法是通过改变温度或添加药物等方式平等地作用于实验室实验中的所有细胞。如今，随着功能基因组学技术的进步，可以对单个细胞成分进行单细胞遗传扰动。使用不同技术的扰动针对蛋白质生产层次的不同层。

在最低层，CRISPR-cas9 直接作用于基因组，利用插入-缺失多态性诱导移码突变，从而有效敲除一个或多个指定基因。比较新的CRISPRi和 CRISPRa 技术分别抑制或激活转录。CRISPR-cas13 作用于蛋白质生产层次的下一层，促进 RNA 降解。小分子药物通常直接作用于酶和受体等蛋白质产品。当这些技术应用于大规模筛选时，它们会在基因型、转录本、蛋白质、染色质可及性以及在某些情况下的表型之间创建一张图谱。

图中展示单细胞扰动实验分析

扰动实验数据分析

对越来越大的扰动数据集进行可靠的分析需要高效的统计工具来处理大量细胞和扰动。扰动-响应数据固有的高维性使扰动之间距离的计算变得复杂，细胞间差异和数据稀疏性也使计算扰动之间距离变得复杂。理想情况下，扰动之间的统计比较和扰动强度的量化应基于细胞组之间的多元距离测量。这种距离测量描述了用不同扰动处理的细胞表达谱之间的差异或相似性，从而推断独特或共享的机制或识别扰动目标，这些目标往往会产生分子谱的类似变化。

（1）一些研究通过组合给定扰动设置中的所有细胞来进行伪批量计算，意味着会丢失有关细胞间响应差异的信息。

（2）对多种细胞类型的研究已经开发出复杂的方法来量化异质细胞群之间的相似性。

（3）基于细胞的统计测量可用于识别成功扰动的细胞，但目前不能量化扰动效应的相似性。

（4）单细胞社区已经探索了 scRNA-seq 的多种距离测量，包括在最佳传输框架中计算的 Wasserstein 距离、最大平均差异、基于邻域的测量和能量距离 (E 距离)。E距离可用于统计测试以识别强或弱扰动以及区分影响不同细胞子过程的扰动。

扰动实验数据库

数据库地址：http://projects.sanderlab.org/scperturb/。下图展示了数据库中扰动数据的详细信息。

E距离计算

E 距离是高维分布之间的统计距离，已用于定义多变量双样本检验，称为能量检验 (E 检验)。它更常见的名称是能量距离，源于物理学中使用引力能量的原始解释。从形式上讲，它阐述了这样一种概念：如果高维空间中的两个点分布与两个分布的宽度相比相距很远，则它们是可区分的。计算如图所示，计算出E值。E 距离的定义，将高维分子谱的细胞分布宽度与它们之间的距离联系起来。受扰动细胞与未受扰动细胞之间的 E 距离较大，表明扰动引起的分子谱发生了强烈变化。下图a展示了E距离之间的计算。

代码实现及分析结果展示

#####python包安装
pip install scperturb
#####代码实现
import scanpy as sc
import matplotlib.pyplot as pl
import anndata as ad
import pandas as pd
import numpy as np
import seaborn as sns
from src.scperturb import *
from tqdm.notebook import tqdm
from IPython.display import clear_output
from IPython.core.display import display, HTML
##############         1.数据载入         #######
file = '/scratch/peidli/scPerturb/DatlingerBock2021.h5ad'
adata = sc.read_h5ad(file)
if 'processed' in adata.uns.keys():
    print('The dataset is already processed. Skipping processing...')
else:
    adata.layers['counts'] = adata.X.copy()
    # basic qc and pp
    sc.pp.filter_cells(adata, min_counts=1000)
    sc.pp.normalize_per_cell(adata)
    sc.pp.filter_genes(adata, min_cells=50)
    sc.pp.log1p(adata)
    # high class imbalance
    adata = equal_subsampling(adata, 'perturbation', N_min=200)
    sc.pp.filter_genes(adata, min_cells=3)  # sanity cleaning
    # select HVGs
    n_var_max = 2000  # max total features to select
    sc.pp.highly_variable_genes(adata, n_top_genes=n_var_max, subset=False, flavor='seurat_v3', layer='counts')
    sc.pp.pca(adata, use_highly_variable=True)
    sc.pp.neighbors(adata)
    adata.uns['processed'] = True
##############         2.E-distance计算       #######
estats = edist(adata, obs_key='perturbation', obsm_key='X_pca', dist='sqeuclidean')
with sns.axes_style('whitegrid'):
    sns.violinplot(data=estats.loc['control'], inner=None)
    sns.swarmplot(data=estats.loc['control'], size=10, color='white')
##############       2.1 小提琴图可视化距离分布图    #######
pl.xticks([0], [''])
pl.ylabel('E-distance to control')
pl.title('Distribution of E-distances to control')
pl.show()

##############       2.2 热图可视化距离大小    #######
fig, ax = pl.subplots(1,1, figsize=[10,8])
order = estats.sort_index().index
sns.heatmap(estats.loc[order, order])
pl.show()

##############         3.E-testing计算         #######
df = etest(adata, obs_key='perturbation', obsm_key='X_pca', dist='sqeuclidean', control='control', alpha=0.05, runs=1000, n_jobs=-1)
df.loc[df.index=='control', 'significant_adj'] = 'control'
df['neglog10_pvalue_adj'] = -np.log10(df['pvalue_adj'])
with sns.axes_style('whitegrid'):
    sns.scatterplot(data=df, y='neglog10_pvalue_adj', x='edist', hue='significant_adj', palette={True: 'tab:green', False: 'tab:red', 'control': 'tab:orange'}, s=30)
pl.title('E-test results')
pl.xlabel('E-distance from control')
pl.ylabel('E-test neg log10 of adjusted p-value')
pl.show()

with sns.axes_style('whitegrid'):
    sns.violinplot(data=df, y='edist', inner=None)
    sns.swarmplot(data=df, x=[0]*len(df), y='edist', size=10, hue='significant_adj', dodge=False, palette={True: 'tab:green', False: 'tab:red', 'control': 'tab:orange'})
pl.xticks([0], [''])
pl.ylabel('E-distance to control')
pl.title('Distribution of E-distances to control')
pl.show()

欢迎大家留言讨论，如对结果有疑问或者方法实现有疑问的可以联系管理员进行解决：kriswcyYQ。