1. 背景
在前面的教程中恼五，我們從數(shù)據(jù)集中刪除了低質(zhì)量的細(xì)胞磨淌，包括計數(shù)較差以及雙細(xì)胞给涕，并將數(shù)據(jù)存放在anndata文件中鬼佣。由于單細(xì)胞測序技術(shù)的限制驶拱，我們在樣本中獲得RNA的時候，經(jīng)過了分子捕獲晶衷，逆轉(zhuǎn)錄還有測序蓝纲。這些步驟會影響同一種細(xì)胞的細(xì)胞間的測序計數(shù)深度的變異性，故單細(xì)胞測序數(shù)據(jù)中的細(xì)胞間差異可能會包含了這部分測序誤差晌纫，等價于計數(shù)矩陣中包含了變化很大的方差項税迷。但在目前的統(tǒng)計方法中，絕大部分模型都預(yù)先假定了數(shù)據(jù)具有相同的方差結(jié)構(gòu)锹漱。
伽馬-泊松分布
從理論上和經(jīng)驗上建立的 UMI 數(shù)據(jù)模型是 Gamma-Poisson 分布箭养，即，其中
代表UMI平均值哥牍，
代表細(xì)胞UMI的過度離散值毕泌。若
時，意味著此時UMI的分布為泊松分布嗅辣。
“歸一化”的預(yù)處理步驟旨在通過將“UMI的方差”縮放到指定范圍撼泛，來調(diào)整數(shù)據(jù)集中的原始UMI計數(shù)以實現(xiàn)模型建模。而在真實的單細(xì)胞分析中澡谭，有不同的歸一化方法以解決不同的分析問題愿题。但經(jīng)驗發(fā)現(xiàn)，移位對數(shù)在大部分?jǐn)?shù)據(jù)中的表現(xiàn)良好蛙奖，這在2023年4月的Nature Method上的基準(zhǔn)測試中有提到潘酗。
本章將向讀者介紹兩種不同的歸一化技術(shù)：移位對數(shù)變換和皮爾遜殘差的解析近似。移位對數(shù)有利于穩(wěn)定方差雁仲，以利于后續(xù)降維和差異表達(dá)基因的識別仔夺。皮爾森近似殘差可以保留生物學(xué)差異，并鑒定稀有細(xì)胞類型攒砖。
我們首先導(dǎo)入我們所需要的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
我們首先檢查原始計數(shù)UMI的分布祭衩，一般在后續(xù)的分析中我們會忽略這一步，但對該分布的認(rèn)識有利于我們理解歸一化的意義阅签。
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)槊恳粋€細(xì)胞確定自己的尺寸因子精算，以滿足同時考慮采樣效果和不同細(xì)胞尺寸的變換瓢宦。細(xì)胞的尺寸因子可以計算為：
其中 代表不同的基因，
代表基因的計數(shù)總和灰羽。確定尺寸因子的方法有很多驮履，在
scanpy中，我們默認(rèn)使用原始計數(shù)深度的中位數(shù)來計算廉嚼，而在seruat中使用固定值玫镐，而在
omicverse的預(yù)處理中，我們將設(shè)定為
怠噪。不同的值會使得過度離散值
的不同恐似。
過度離散值
過度離散值
移位對數(shù)是一種快速歸一化技術(shù)傍念，優(yōu)于其他揭示數(shù)據(jù)集潛在結(jié)構(gòu)的方法（特別是在進(jìn)行主成分分析時）矫夷，并且有利于方差的穩(wěn)定性，以進(jìn)行后續(xù)的降維和差異表達(dá)基因的識別憋槐。我們現(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ù)中的細(xì)胞間變異包括了生物異質(zhì)性以及技術(shù)效應(yīng)。而移位計數(shù)并不能很好的排除兩種不同變異來源的混淆誤差顺又。皮爾森近似殘差利用了“正則化負(fù)二項式回歸”的皮爾森殘差來計算數(shù)據(jù)中潛在的技術(shù)噪音更卒，將計數(shù)深度添加為廣義線性模型中的協(xié)變量，而在不同的歸一化方法的測試中稚照，皮爾森殘差法可以消除計數(shù)效應(yīng)帶來的誤差蹂空，并且保留了數(shù)據(jù)集中的細(xì)胞異質(zhì)性。
此外果录，皮爾森殘差法不需要進(jìn)行啟發(fā)式步驟上枕，如偽計數(shù)加法/對數(shù)變化，該方法的輸出就是歸一化后的值弱恒，包括了正值和負(fù)值辨萍。細(xì)胞和基因的負(fù)殘差表明與基因的平均表達(dá)和細(xì)胞測序深度相比，觀察到的計數(shù)少于預(yù)期返弹。正殘差分別表示計數(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中提供了預(yù)處理函數(shù)pp.preprocess椅棺，該方法可直接計算移位對數(shù)或皮爾森殘差犁罩，方法內(nèi)同時包括了基于移位對數(shù)/皮爾森殘差的高可變基因的選擇方法，高可變基因會在下一節(jié)的教程中進(jìn)行講解土陪。
此外昼汗，我們omicverse的歸一化方法還提供了穩(wěn)定基因的識別與過濾。
4.1 移位對數(shù)
在omicverse中鬼雀，我們設(shè)置mode='shiftlog|pearson'即可完成移位對數(shù)的計算顷窒，一般來說，默認(rèn)的target_sum=50*1e4鞋吉，同時高可變基因定義為前2000個,需要注意的是，當(dāng)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.
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 皮爾森近似殘差
我們也可以設(shè)置mode='pearson|pearson'來完成皮爾森近似殘差的計算，此時我們不需要輸入target_sum坛掠，需要注意的是赊锚，當(dāng)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.
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()
以上屉栓，就是本章所要介紹的歸一化內(nèi)容舷蒲，通過benchmark的測試，我們發(fā)現(xiàn)移位對數(shù)適用于大多數(shù)任務(wù)友多。但是如果我們的分析目標(biāo)是尋找稀有細(xì)胞的時候牲平，可以考慮采用皮爾森殘差法來進(jìn)行歸一化。
