单细胞测序最好的教程(二):归一化

1. 背景

在前面的教程中,我们从数据集中删除了低质量的细胞,包括计数较差以及双细胞,并将数据存放在anndata文件中。由于单细胞测序技术的限制,我们在样本中获得RNA的时候,经过了分子捕获,逆转录还有测序。这些步骤会影响同一种细胞的细胞间的测序计数深度的变异性,故单细胞测序数据中的细胞间差异可能会包含了这部分测序误差,等价于计数矩阵中包含了变化很大的方差项。但在目前的统计方法中,绝大部分模型都预先假定了数据具有相同的方差结构。

伽马-泊松分布
从理论上和经验上建立的 UMI 数据模型是 Gamma-Poisson 分布,即Var[Y] = \mu + \alpha \mu^2,其中\mu代表UMI平均值,\alpha代表细胞UMI的过度离散值。若 \alpha=0 时,意味着此时UMI的分布为泊松分布。

“归一化”的预处理步骤旨在通过将“UMI的方差”缩放到指定范围,来调整数据集中的原始UMI计数以实现模型建模。而在真实的单细胞分析中,有不同的归一化方法以解决不同的分析问题。但经验发现,移位对数在大部分数据中的表现良好,这在2023年4月的Nature Method上的基准测试中有提到。

本章将向读者介绍两种不同的归一化技术:移位对数变换和皮尔逊残差的解析近似。移位对数有利于稳定方差,以利于后续降维和差异表达基因的识别。皮尔森近似残差可以保留生物学差异,并鉴定稀有细胞类型。

我们首先导入我们所需要的Python包,以及上一个教程分析所得到的anndata文件。

import omicverse as ov
import scanpy as sc

ov.utils.ov_plot_set()
adata = sc.read(
    filename="s4d8_quality_control.h5ad",
    backup_url="https://figshare.com/ndownloader/files/40014331",
)
try downloading from url
https://figshare.com/ndownloader/files/40014331
... this may take a while but only happens once

我们首先检查原始计数UMI的分布,一般在后续的分析中我们会忽略这一步,但对该分布的认识有利于我们理解归一化的意义。

import seaborn as sns
p1 = sns.histplot(adata.obs["total_counts"], bins=100, kde=False)
原始计数分布

2. 移位对数

我们将介绍的第一个归一化方法是基于delta方法的移位对数。Delta方法应用非线性函数f(Y),使得原始计数 Y 中的差异更加相似。

我们定义非线性函数f(Y)的变换如下:

f(y) = \log(\frac{y}{s}+y_0)

其中y是原始的计数,s是尺寸因子,y_0是伪计数。我们为每一个细胞确定自己的尺寸因子,以满足同时考虑采样效果和不同细胞尺寸的变换。细胞的尺寸因子可以计算为:

s_c = \frac{\sum_g y_{gc}}{L}

其中 g代表不同的基因,L代表基因的计数总和。确定尺寸因子的方法有很多,在scanpy中,我们默认使用原始计数深度的中位数来计算,而在seruat中使用固定值L=10^4,而在omicverse的预处理中,我们将L设定为50*10^4。不同的值会使得过度离散值 \alpha的不同。

过度离散值 \alpha
过度离散值 \alpha描述了数据集中存在着比期望更大的变异性。

移位对数是一种快速归一化技术,优于其他揭示数据集潜在结构的方法(特别是在进行主成分分析时),并且有利于方差的稳定性,以进行后续的降维和差异表达基因的识别。我们现在将检查如何将此归一化方法应用于我们的数据集。我们可以使用pp.normalized_total来使用 scanpy 调用移位对数。并且我们设置target_sum=None,inplace=False来探索两种不同的归一化技术。

scales_counts = sc.pp.normalize_total(adata, target_sum=None, inplace=False)
# log1p transform
adata.layers["log1p_norm"] = sc.pp.log1p(scales_counts["X"], copy=True)
normalizing counts per cell
    finished (0:00:00)

我们使用hist图来对比归一化前后的计数变化

import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(8, 4))
p1 = sns.histplot(adata.obs["total_counts"], bins=100, kde=False, ax=axes[0])
axes[0].set_title("Total counts")
p2 = sns.histplot(adata.layers["log1p_norm"].sum(1), bins=100, kde=False, ax=axes[1])
axes[1].set_title("Shifted logarithm")
plt.show()
归一化前后数据分布对比

我们发现nUMI的最大值在1000左右,经过移位对数化后,并且移位对数化后的nUMI的分布近似正态分布。

3. 皮尔森近似残差

我们观察到,scRNA-seq数据中的细胞间变异包括了生物异质性以及技术效应。而移位计数并不能很好的排除两种不同变异来源的混淆误差。皮尔森近似残差利用了“正则化负二项式回归”的皮尔森残差来计算数据中潜在的技术噪音,将计数深度添加为广义线性模型中的协变量,而在不同的归一化方法的测试中,皮尔森残差法可以消除计数效应带来的误差,并且保留了数据集中的细胞异质性。

此外,皮尔森残差法不需要进行启发式步骤,如伪计数加法/对数变化,该方法的输出就是归一化后的值,包括了正值和负值。细胞和基因的负残差表明与基因的平均表达和细胞测序深度相比,观察到的计数少于预期。正残差分别表示计数越多。解析 Pearon 残差在 scanpy 中实现,可以直接在原始计数矩阵上计算。

from scipy.sparse import csr_matrix
analytic_pearson = sc.experimental.pp.normalize_pearson_residuals(adata, inplace=False)
adata.layers["analytic_pearson_residuals"] = csr_matrix(analytic_pearson["X"])
computing analytic Pearson residuals on adata.X
    finished (0:00:15)
fig, axes = plt.subplots(1, 2, figsize=(8, 4))
p1 = sns.histplot(adata.obs["total_counts"], bins=100, kde=False, ax=axes[0])
axes[0].set_title("Total counts")
p2 = sns.histplot(
    adata.layers["analytic_pearson_residuals"].sum(1), bins=100, kde=False, ax=axes[1]
)
axes[1].set_title("Analytic Pearson residuals")
plt.show()
皮尔森残差归一化数据分布

如果我们设置inplace=False时,我们归一化的计数矩阵会取代原anndata文件中的计数矩阵,即更改adata.X的内容。

4. 一键式归一化

我们在omicverse中提供了预处理函数pp.preprocess,该方法可直接计算移位对数或皮尔森残差,方法内同时包括了基于移位对数/皮尔森残差的高可变基因的选择方法,高可变基因会在下一节的教程中进行讲解。

此外,我们omicverse的归一化方法还提供了稳定基因的识别与过滤。

4.1 移位对数

omicverse中,我们设置mode='shiftlog|pearson'即可完成移位对数的计算,一般来说,默认的target_sum=50*1e4,同时高可变基因定义为前2000个,需要注意的是,当omicverse的版本小于1.4.13时,mode的参数只能设置为scanpypearsonscanpyshiftlog的意义是相同的

adata_log=ov.pp.preprocess(adata,mode='shiftlog|pearson',n_HVGs=2000,)
adata_log
Begin log-normalization, HVGs identification
After filtration, 20171/20171 genes are kept. Among 20171 genes, 20171 genes are robust.
End of log-normalization, HVGs identification.
Begin size normalization and pegasus batch aware HVGs selection or Perason residuals workflow
normalizing counts per cell The following highly-expressed genes are not considered during normalization factor computation:
[]
    finished (0:00:00)
2000 highly variable features have been selected.
End of size normalization and pegasus batch aware HVGs selection or Perason residuals workflow.





AnnData object with n_obs × n_vars = 14814 × 20171
    obs: 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_20_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'total_counts_hb', 'log1p_total_counts_hb', 'pct_counts_hb', 'outlier', 'mt_outlier', 'scDblFinder_score', 'scDblFinder_class'
    var: 'gene_ids', 'feature_types', 'genome', 'mt', 'ribo', 'hb', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells', 'percent_cells', 'robust', 'highly_variable_features', 'mean', 'var', 'hvf_loess', 'hvf_rank'
    uns: 'log1p'
    layers: 'counts', 'soupX_counts'

4.2 皮尔森近似残差

我们也可以设置mode='pearson|pearson'来完成皮尔森近似残差的计算,此时我们不需要输入target_sum,需要注意的是,当omicverse的版本小于1.4.13时,mode的参数只能设置为scanpypearson

adata_pearson=ov.pp.preprocess(adata,mode='pearson|pearson',n_HVGs=2000,)
adata_pearson
Begin log-normalization, HVGs identification
After filtration, 20171/20171 genes are kept. Among 20171 genes, 20171 genes are robust.
End of log-normalization, HVGs identification.
Begin size normalization and pegasus batch aware HVGs selection or Perason residuals workflow
extracting highly variable genes
--> added
    'highly_variable', boolean vector (adata.var)
    'highly_variable_rank', float vector (adata.var)
    'highly_variable_nbatches', int vector (adata.var)
    'highly_variable_intersection', boolean vector (adata.var)
    'means', float vector (adata.var)
    'variances', float vector (adata.var)
    'residual_variances', float vector (adata.var)
computing analytic Pearson residuals on adata.X
    finished (0:00:04)
End of size normalization and pegasus batch aware HVGs selection or Perason residuals workflow.





AnnData object with n_obs × n_vars = 14814 × 20171
    obs: 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_20_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'total_counts_hb', 'log1p_total_counts_hb', 'pct_counts_hb', 'outlier', 'mt_outlier', 'scDblFinder_score', 'scDblFinder_class'
    var: 'gene_ids', 'feature_types', 'genome', 'mt', 'ribo', 'hb', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells', 'percent_cells', 'robust', 'mean', 'var', 'residual_variances', 'highly_variable_rank', 'highly_variable_features'
    uns: 'hvg', 'pearson_residuals_normalization'
    layers: 'counts', 'soupX_counts'
import matplotlib.pyplot as plt
import seaborn as sns
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
p1 = sns.histplot(adata.obs["total_counts"], bins=100, kde=False, ax=axes[0])
axes[0].set_title("Total counts")

p2 = sns.histplot(
    adata_log.X.sum(1), bins=100, kde=False, ax=axes[1]
)
axes[1].set_title("Shifted logarithm")

p3 = sns.histplot(
    adata_pearson.X.sum(1), bins=100, kde=False, ax=axes[2]
)
axes[2].set_title("Analytic Pearson residuals")
plt.show()
两种不同数据分布的对比

以上,就是本章所要介绍的归一化内容,通过benchmark的测试,我们发现移位对数适用于大多数任务。但是如果我们的分析目标是寻找稀有细胞的时候,可以考虑采用皮尔森残差法来进行归一化。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,179评论 5 476
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,229评论 2 380
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 150,032评论 0 336
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,533评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,531评论 5 365
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,539评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,916评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,574评论 0 256
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,813评论 1 296
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,568评论 2 320
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,654评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,354评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,937评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,918评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,152评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,852评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,378评论 2 342

推荐阅读更多精彩内容