1. 背景

在前面的教程中，我們從數(shù)據(jù)集中刪除了低質(zhì)量的細胞，包括計數(shù)較差以及雙細胞，并將數(shù)據(jù)存放在 anndata文件中。由于單細胞測序技術(shù)的限制，我們在樣本中獲得RNA的時候，經(jīng)過了分子捕獲，逆轉(zhuǎn)錄還有測序。這些步驟會影響同一種細胞的細胞間的測序計數(shù)深度的變異性，故單細胞測序數(shù)據(jù)中的細胞間差異可能會包含了這部分測序誤差，等價于計數(shù)矩陣中包含了變化很大的方差項。但在目前的統(tǒng)計方法中，絕大部分模型都預先假定了數(shù)據(jù)具有相同的方差結(jié)構(gòu)。

伽瑪-柏松分布

從理論上和經(jīng)驗上建立的 UMI 數(shù)據(jù)模型是 Gamma-Poisson 分布，即，其中代表UMI平均值，代表細胞UMI的過度離散值。若 時，意味著此時UMI的分布為泊松分布。

“歸一化”的預處理步驟旨在通過將“UMI的方差”縮放到指定范圍，來調(diào)整數(shù)據(jù)集中的原始UMI計數(shù)以實現(xiàn)模型建模。而在真實的單細胞分析中，有不同的歸一化方法以解決不同的分析問題。但經(jīng)驗發(fā)現(xiàn)，移位對數(shù)在大部分數(shù)據(jù)中的表現(xiàn)良好，這在2023年4月的Nature Method上的基準測試中有提到。

本章將向讀者介紹兩種不同的歸一化技術(shù)：移位對數(shù)變換和皮爾遜殘差的解析近似。移位對數(shù)有利于穩(wěn)定方差，以利于后續(xù)降維和差異表達基因的識別。皮爾森近似殘差可以保留生物學差異，并鑒定稀有細胞類型。

我們首先導入我們所需要的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:///ndownloader/files/40014331",
)

   # try downloading from url
   # https:///ndownloader/files/40014331
   # ... this may take a while but only happens once

我們首先檢查原始計數(shù)UMI的分布，一般在后續(xù)的分析中我們會忽略這一步，但對該分布的認識有利于我們理解歸一化的意義。

import seaborn as sns
p1 = sns.histplot(adata.obs["total_counts"], bins=100, kde=False)

2. 移位對數(shù)

我們將介紹的第一個歸一化方法是基于delta方法的移位對數(shù)。Delta方法應(yīng)用非線性函數(shù)，使得原始計數(shù) 中的差異更加相似。

我們定義非線性函數(shù)的變換如下：

其中是原始的計數(shù)，是尺寸因子，是偽計數(shù)。我們?yōu)槊恳粋€細胞確定自己的尺寸因子，以滿足同時考慮采樣效果和不同細胞尺寸的變換。細胞的尺寸因子可以計算為：

其中 代表不同的基因，代表基因的計數(shù)總和。確定尺寸因子的方法有很多，在scanpy中，我們默認使用原始計數(shù)深度的中位數(shù)來計算，而在seruat中使用固定值，而在omicverse的預處理中，我們將設(shè)定為。不同的值會使得過度離散值 的不同。

過度離散值

過度離散值 描述了數(shù)據(jù)集中存在著比期望更大的變異性。

移位對數(shù)是一種快速歸一化技術(shù)，優(yōu)于其他揭示數(shù)據(jù)集潛在結(jié)構(gòu)的方法（特別是在進行主成分分析時），并且有利于方差的穩(wěn)定性，以進行后續(xù)的降維和差異表達基因的識別。我們現(xiàn)在將檢查如何將此歸一化方法應(yīng)用于我們的數(shù)據(jù)集。我們可以使用pp.normalized_total來使用 scanpy 調(diào)用移位對數(shù)。并且我們設(shè)置target_sum=None,inplace=False來探索兩種不同的歸一化技術(shù)。

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圖來對比歸一化前后的計數(shù)變化

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

我們發(fā)現(xiàn)nUMI的最大值在1000左右，經(jīng)過移位對數(shù)化后，并且移位對數(shù)化后的nUMI的分布近似正態(tài)分布。

3. 皮爾森近似殘差

我們觀察到，scRNA-seq數(shù)據(jù)中的細胞間變異包括了生物異質(zhì)性以及技術(shù)效應(yīng)。而移位計數(shù)并不能很好的排除兩種不同變異來源的混淆誤差。皮爾森近似殘差利用了“正則化負二項式回歸”的皮爾森殘差來計算數(shù)據(jù)中潛在的技術(shù)噪音，將計數(shù)深度添加為廣義線性模型中的協(xié)變量，而在不同的歸一化方法的測試中，皮爾森殘差法可以消除計數(shù)效應(yīng)帶來的誤差，并且保留了數(shù)據(jù)集中的細胞異質(zhì)性。

此外，皮爾森殘差法不需要進行啟發(fā)式步驟，如偽計數(shù)加法/對數(shù)變化，該方法的輸出就是歸一化后的值，包括了正值和負值。細胞和基因的負殘差表明與基因的平均表達和細胞測序深度相比，觀察到的計數(shù)少于預期。正殘差分別表示計數(shù)越多。解析 Pearon 殘差在 scanpy 中實現(xiàn)，可以直接在原始計數(shù)矩陣上計算。

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

如果我們設(shè)置inplace=False時，我們歸一化的計數(shù)矩陣會取代原anndata文件中的計數(shù)矩陣，即更改adata.X的內(nèi)容。

4. 一鍵式歸一化

我們在omicverse中提供了預處理函數(shù)pp.preprocess，該方法可直接計算移位對數(shù)或皮爾森殘差，方法內(nèi)同時包括了基于移位對數(shù)/皮爾森殘差的高可變基因的選擇方法，高可變基因會在下一節(jié)的教程中進行講解。

此外，我們omicverse的歸一化方法還提供了穩(wěn)定基因的識別與過濾。

4.1 移位對數(shù)

在omicverse中，我們設(shè)置mode='shiftlog|pearson'即可完成移位對數(shù)的計算，一般來說，默認的target_sum=50*1e4，同時高可變基因定義為前2000個,需要注意的是，當omicverse的版本小于1.4.13時，mode的參數(shù)只能設(shè)置為scanpy或pearson，scanpy與shiftlog的意義是相同的

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.

4.2 皮爾森近似殘差

我們也可以設(shè)置mode='pearson|pearson'來完成皮爾森近似殘差的計算，此時我們不需要輸入target_sum，需要注意的是，當omicverse的版本小于1.4.13時，mode的參數(shù)只能設(shè)置為scanpy或pearson

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.

我們對比一下兩種分析的結(jié)果

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

以上，就是本章所要介紹的歸一化內(nèi)容，通過benchmark的測試，我們發(fā)現(xiàn)移位對數(shù)適用于大多數(shù)任務(wù)。但是如果我們的分析目標是尋找稀有細胞的時候，可以考慮采用皮爾森殘差法來進行歸一化。

5. 思考

為了加深你對本章的理解，我們提出了以下思考題，如有興趣作答者，可將答案發(fā)送至郵箱starlitnightly@163.com，郵件標題為姓名/昵稱-單細胞最好教程（二）思考題

• 我們在進行移位對數(shù)分析的時候，sc.pp.normalize_total(adata, target_sum=None, inplace=False) 中的target_sum 使用了默認值，在 seurat 中默認值是10,000，在一些教程中設(shè)定為1,000,000，我們雖然對這個值的意義進行了簡單介紹，但你認為不同的值背后的含義是什么？
• 我們?yōu)槭裁磿褂闷柹瓪埐顏碛嬎銡w一化值，相對于移位對數(shù)而言有什么更好的地方？
• 你可以找出別的歸一化方法，并比較其與移位對數(shù)，皮爾森殘差的好壞嗎？

如果你對omicverse的使用有什么疑問，或者也想加入一起開發(fā)擅長，想在生物信息學工具開發(fā)路上找到同道中人，那么這個Python的轉(zhuǎn)錄組學分析框架與生態(tài)交流群請不要錯過。感謝Jimmy老師對omicverse的大力支持。