satijalabsctransform包的github網(wǎng)址：https://github.com/ChristophH/sctransform
sctransform作者發(fā)表的論文：https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1874-1
sctransform中的vst方法：https://rawgit.com/ChristophH/sctransform/supp_html/supplement/variance_stabilizing_transformation.html
在seurat如何使用sctransform：https://satijalab.org/seurat/articles/sctransform_vignette.html
注：vst方法(Variance stabilizing transformation)，可通過'?sctransform::vst'詳細(xì)了解

一唆樊、sctransform簡介

seurat包的sctransform函數(shù)調(diào)用sctransform::vst仿便。sctransform包是由紐約基因組中心Rahul Satija實驗室的Christoph Hafemeister開發(fā)(也是satijalab實驗室出品)绊诲，使用正則化負(fù)二項式回歸(regularized negative binomial regression)對單細(xì)胞UMI表達(dá)數(shù)據(jù)進(jìn)進(jìn)行建模硫豆，以消除由于測序深度引起的變化袭艺。該方法首先使用測序深度作為自變量和UMI計數(shù)作為響應(yīng)或因變量為每個基因構(gòu)建廣義線性模型 (GLM)蹬癌。 然后根據(jù)基因表達(dá)對參數(shù)估計進(jìn)行正則化（或調(diào)整）。 使用正則化參數(shù)應(yīng)用第二輪負(fù)二項式回歸粗卜。 該模型的輸出（殘差）是每個基因的標(biāo)準(zhǔn)化表達(dá)水平刀森。 SCTransform方法學(xué)習(xí)基因集特定的特征( gene-group specific factors)杠览，而不是使用常數(shù)因子來標(biāo)準(zhǔn)化所有基因。 這些基因集特征分別針對低当宴、中和高表達(dá)基因畜吊，消除了技術(shù)變異的影響，同時保留了真正的生物異質(zhì)性户矢。

該軟件包的核心功能已集成到Seurat包中玲献。與傳統(tǒng)的lognormalize歸一化方法相比，它集成了NormalizeData(), FindVariableFeatures()，ScaleData()三個函數(shù)的功能青自，是單細(xì)胞基因count矩陣預(yù)處理的替代方法株依。 結(jié)果保存在一個新的assay中（默認(rèn)命名為 SCT），包含：

• seurat_obj[['SCT']]@counts：矯正后的counts延窜；
• seurat_obj[['SCT']]@data：矯正后的counts通過log1p(counts)處理后的data恋腕；

• seurat_obj[['SCT']]@scale.data：pearson殘差；
  sctransform::vst 中間結(jié)果保存在seurat_obj[['SCT']]@misc中逆瑞。

  
  
  image.png

我們先從比較淺的教材出發(fā)荠藤，一層層去了解和熟悉sctransform算法背后的邏輯，來回答以下幾個問題

1. 什么情況下SCTransform是優(yōu)于全局scale歸一化方法（global scaling normalization）
2. SCTransform的工作原理
3. 在執(zhí)行SCTransform分析中要注意哪些

二获高、seurat分析使用sctransform

教程鏈接：https://satijalab.org/seurat/articles/sctransform_vignette.html
單細(xì)胞RNA-seq數(shù)據(jù)中的生物異質(zhì)性通常會受到包括測序深度在內(nèi)的技術(shù)因素的影響哈肖。在每個細(xì)胞中檢測到的分子數(shù)量在細(xì)胞之間可能存在顯著的差異，即使在同一細(xì)胞類型中也是如此念秧。scRNA-seq數(shù)據(jù)需要進(jìn)行有效的預(yù)處理和標(biāo)準(zhǔn)化以消除這種實驗技術(shù)來源的變異性淤井。在Hafemeister and Satija, 2019文章中引入了一個建模框架摊趾，用scRNA-seq 實驗的UMI計數(shù)數(shù)據(jù)進(jìn)行歸一化和方差穩(wěn)定化轉(zhuǎn)化币狠。同時也改進(jìn)了常見的下游分析，如高可變基因選擇砾层、降維和差異基因分析漩绵。

library(Seurat)
library(ggplot2)
library(sctransform)

加載數(shù)據(jù)并創(chuàng)建 Seurat 對象

pbmc_data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
pbmc <- CreateSeuratObject(counts = pbmc_data)

應(yīng)用sctransform歸一化

• 這一命令替換NormalizeData()，ScaleData()和FindVariableFeatures();
• 轉(zhuǎn)換后的數(shù)據(jù)將存儲在SCT assay中肛炮，運(yùn)行sctransform后將'SCT'設(shè)置為默認(rèn)assay;
• 在標(biāo)準(zhǔn)化處理中止吐，我們還可以去除混雜的變異來源，例如線粒體百分比;

# store mitochondrial percentage in object meta data
pbmc <- PercentageFeatureSet(pbmc, pattern = "^MT-", col.name = "percent.mt")

# run sctransform
pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt", verbose = FALSE)

最新版本sctransform還支持使用glmGamPoi包侨糟，這大大提高了算法學(xué)習(xí)過程的速度碍扔。可以指定調(diào)用方法method="glmGamPoi"粟害。

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
  BiocManager::install("glmGamPoi")

pbmc <- SCTransform(pbmc, method = "glmGamPoi", vars.to.regress = "percent.mt", verbose = FALSE)

怎么理解glmGamPoi蕴忆？
gmlGamPoi的設(shè)計目的：在任意大小的數(shù)據(jù)集上擬合Gamma-Poisson廣義線性模型颤芬；比替代方法更快悲幅，例如DESeq2或edgeR。

執(zhí)行PCA和UMAP降維

# These are now standard steps in the Seurat workflow for visualization and clustering
pbmc <- RunPCA(pbmc, verbose = FALSE)
pbmc <- RunUMAP(pbmc, dims = 1:30, verbose = FALSE)
pbmc <- FindNeighbors(pbmc, dims = 1:30, verbose = FALSE)
pbmc <- FindClusters(pbmc, verbose = FALSE)
DimPlot(pbmc, label = TRUE) + NoLegend()

image.png

問題1：為什么我們在使用sctransform時可以選擇更多的PC杰妓？

標(biāo)準(zhǔn)Seurat分析流程，我們對此數(shù)據(jù)集只關(guān)注前10個PC碘勉，但我們強(qiáng)調(diào)nPC=10巷挥，與PC數(shù)目設(shè)置得更高，下游結(jié)果是相似的验靡。但是倍宾，我們發(fā)現(xiàn)在使用 sctransform時，通常會將此參數(shù)設(shè)置得更高而越有利胜嗓。這是因為sctransform分析流程執(zhí)行了更有效的標(biāo)準(zhǔn)化高职，從數(shù)據(jù)中強(qiáng)烈地去除技術(shù)帶來的影響。

使用log-normalization歸一化后辞州，測序深度的變化仍然是一個混雜因素怔锌，這種影響可以微妙地影響更高的PC。在 sctransform 中变过，這種影響大大減輕产禾，這意味著更高的PC更有可能代表微妙但與生物學(xué)相關(guān)的異質(zhì)性來源——因此sctransform可能會改善下游分析。

此外牵啦，sctransform默認(rèn)返回 3,000個高可變基因亚情，而不是 2,000個」基本原理是相似的楞件，額外的高可變基因不太可能由跨細(xì)胞的技術(shù)差異所引起的，而可能是代表更微妙的生物波動裳瘪。一般而言土浸，我們發(fā)現(xiàn)使用sctransform產(chǎn)生的結(jié)果對這些參數(shù)的依賴性較小（實際上彭羹，當(dāng)使用轉(zhuǎn)錄組中的所有基因時黄伊，我們獲得了幾乎相同的結(jié)果，盡管使用所有基因這確實降低了計算效率）派殷。這可以幫助用戶生成更可靠的結(jié)果还最，此外，還可以搭建具有相同參數(shù)的標(biāo)準(zhǔn)分析管道毡惜，快速應(yīng)用于新數(shù)據(jù)集：

pbmc <- CreateSeuratObject(pbmc_data) %>%
    PercentageFeatureSet(pattern = "^MT-", col.name = "percent.mt") %>%
    SCTransform(vars.to.regress = "percent.mt") %>%
    RunPCA() %>%
    FindNeighbors(dims = 1:30) %>%
    RunUMAP(dims = 1:30) %>%
    FindClusters()

我的理解是拓轻，sctransform分析流程不需要關(guān)注PC數(shù)目，和高可變基因數(shù)目這些參數(shù)经伙，而反復(fù)調(diào)參扶叉，使用默認(rèn)的參數(shù)即可，便于搭建標(biāo)準(zhǔn)分析流程。

問題2：sctransform處理后的標(biāo)準(zhǔn)化值存儲在哪里枣氧？

如論文中所述溢十，sctransform使用“正則化負(fù)二項式回歸”計算scRNA-seq數(shù)據(jù)中的技術(shù)噪聲模型。該模型的殘差是歸一化數(shù)值达吞，可以是正值也可以是負(fù)值茶宵。給定細(xì)胞中給定基因的正殘差表明，鑒于該基因在細(xì)胞群中的平均表達(dá)和細(xì)胞測序深度宗挥，我們觀察到的UMI比預(yù)期的多乌庶，而負(fù)殘差則表明相反。

sctransfrom 的結(jié)果存儲在'SCT' assay中契耿。

• pbmc[["SCT"]]@scale.data包含殘差（歸一化值）瞒大，并直接用作PCA 的輸入。請注意搪桂，這個矩陣是非稀疏矩陣透敌，因此如果為所有基因存儲，可能會占用大量內(nèi)存踢械。為了節(jié)省內(nèi)存酗电，我們只為高可變基因存儲這些值，方法是在SCTransform()函數(shù)調(diào)用中默認(rèn)設(shè)置 return.only.var.genes = TRUE 内列。

• 協(xié)助基因可視化和解釋撵术。我們還將Pearson殘差轉(zhuǎn)換回“校正后的”UMI 計數(shù)counts』扒疲“校正后的”UMI 計數(shù)counts可以這樣理解嫩与，假設(shè)所有細(xì)胞的測序深度一致，我們希望觀察到的 UMI計數(shù)（原始counts）交排，通過sctransfrom算法剔除測序?qū)嶒炞儺惡蠡蹋C正后的counts。

• “矯正后”的UMI計數(shù)存儲在pbmc[["SCT"]]@counts中埃篓。我們將這些校正計數(shù)進(jìn)行l(wèi)og-normalized处坪，歸一化版本存儲在pbmc[["SCT"]]@data，這對基因表達(dá)量可視化非常有幫助架专。

• 可以使用校正的log-normalized計數(shù)counts（也就是pbmc[["SCT"]]@data）進(jìn)行差異表達(dá)分析和整合分析同窘。但是，原則上胶征，最好直接對殘差（存儲在scale.data中）執(zhí)行這些計算塞椎。目前，Seurat v3不支持此功能睛低，但很快就會支持應(yīng)用scale.data做差異基因等分析。

三、sctransform vs lognormalize

我們通過閱讀作者的論文钱雷，了解sctransform和lognormalize方法上的差異骂铁。

3.1 為什么要Normalization?

在單細(xì)胞RNA-seq (scRNA-seq) 數(shù)據(jù)的分析中，進(jìn)行有效的預(yù)處理和標(biāo)準(zhǔn)化非常關(guān)鍵罩抗，但是也很有挑戰(zhàn)性拉庵。 原始UMI計數(shù)不能直接用于比較細(xì)胞之間的基因表達(dá)，因為它們會被技術(shù)和“無趣”的生物變異所混淆套蒂。 特別是钞支，觀察到的測序深度（每個細(xì)胞檢測到的基因或分子的數(shù)量）在細(xì)胞之間可能存在顯著差異，即使在同一細(xì)胞類型內(nèi)操刀，分子計數(shù)的變化可能跨越一個數(shù)量級烁挟。 雖然現(xiàn)在在scRNA-seq中廣泛使用唯一分子標(biāo)識符 (UMI) 消除了 PCR 擴(kuò)增偏差，但仍需要通過標(biāo)準(zhǔn)化以消除其他技術(shù)變化的影響骨坑，如測序深度撼嗓、細(xì)胞裂解和逆轉(zhuǎn)錄效率等帶來的變異。其實在bulk RNA-seq分析中同樣存在欢唾，但由于scRNA-seq數(shù)據(jù)的極度稀疏性且警，這類問題會更嚴(yán)重。

單細(xì)胞標(biāo)準(zhǔn)化的主要目標(biāo)是消除技術(shù)效應(yīng)對潛在分子計數(shù)的影響礁遣，同時保留真正的生物學(xué)變異斑芜。具體來說，作者建議經(jīng)過有效歸一化分析流程處理的數(shù)據(jù)集應(yīng)具有以下特征：

1. 一般來說祟霍，基因的標(biāo)準(zhǔn)化表達(dá)水平不應(yīng)與細(xì)胞的總測序深度相關(guān)押搪。 下游分析流程（降維、差異表達(dá)）也不應(yīng)受測序深度變化的影響浅碾。
2. 歸一化基因（跨細(xì)胞）的方差應(yīng)主要反映生物異質(zhì)性大州，與基因豐度或測序深度無關(guān)。
   例如垂谢，標(biāo)準(zhǔn)化后具有高方差的基因應(yīng)在不同細(xì)胞類型之間差異表達(dá)厦画，而看家基因應(yīng)表現(xiàn)出低方差。 此外滥朱，在考慮深度測序的細(xì)胞或淺測序的細(xì)胞時根暑，基因的方差應(yīng)該是相似的。

3.2 兩種主流的標(biāo)準(zhǔn)化方法

方法1：全局scale歸一化方法（global scaling normalization）：即基于縮放因子的歸一化方法徙邻。確定單個細(xì)胞的“大小因素”（size factors）排嫌，通過對每個細(xì)胞進(jìn)行統(tǒng)一縮放，這些方法假設(shè)數(shù)據(jù)集中所有細(xì)胞的潛在 RNA 含量是恒定的缰犁，并且可以對所有基因應(yīng)用單個縮放因子淳地，以便將技術(shù)噪音與生物細(xì)胞間的變異性區(qū)分開來怖糊。lognormalize標(biāo)準(zhǔn)化屬于此類。

方法2：使用概率方法對分子計數(shù)進(jìn)行建模颇象，sctransform屬于此類型伍伤。他們認(rèn)為不同的基因組別(gene group)不能被相同的常數(shù)因子歸一化，質(zhì)疑基于縮放因子的歸一化方法遣钳，縮放因子是如何計算的扰魂。sctransform的“正則化負(fù)二項式回歸”的殘差代表了有效歸一化的數(shù)據(jù)值，不再受技術(shù)特征的影響蕴茴，但保留了由不同生物狀態(tài)驅(qū)動的異質(zhì)性劝评。最后，證明這些歸一化值能夠進(jìn)行下游分析倦淀，例如降維和差異表達(dá)測試蒋畜，結(jié)果不會被細(xì)胞測序深度混淆。

這兩種方法也是目前seurat主要的標(biāo)準(zhǔn)化處理方法晃听。

3.3 單一縮放因子不能有效地標(biāo)準(zhǔn)化低表達(dá)和高表達(dá)的基因

為了探究單細(xì)胞測序深度的影響程度百侧，作者從10x Genomics官網(wǎng)下載了不同組織的5個數(shù)據(jù)集，共33,148 個人類外周血單核細(xì)胞 (PBMC) 能扒，查看log-normalization方法標(biāo)準(zhǔn)化處理的有效性佣渴。

圖A：細(xì)胞總UMI計數(shù)（“測序深度”）的分布圖；每個細(xì)胞平均 1891個UMI（中位數(shù)）初斑。
圖B：根據(jù)基因在數(shù)據(jù)集中的平均表達(dá)量將基因分為6組辛润，即Gene group1-group6高中低表達(dá)組。
圖C：x軸是每個細(xì)胞的總UMI數(shù)（“測序深度”）见秤，y軸是觀察到的基因UMI數(shù)目砂竖，為每個基因擬合一條平滑線，黑色線表示均值擬合線鹃答，彩色區(qū)域表示四分位距乎澄。能很明顯的看出，這6個Gene group均表現(xiàn)出非標(biāo)準(zhǔn)化表達(dá)量（基因UMI計數(shù)测摔，y軸）和細(xì)胞測序深度（x軸）之間存在很強(qiáng)的線性關(guān)系置济。
圖D：我們應(yīng)用scaled log-normalized處理（后面簡稱log-normalization）后的標(biāo)準(zhǔn)化表達(dá)值，和圖C一樣繪圖锋八，y軸是縮放的log歸一化值浙于，而非原始的基因UMI計數(shù)。6個Gene group處理后的均值擬合線趨勢不一致挟纱。
圖E：觀察基因變異方差來源與細(xì)胞測序深度的關(guān)系羞酗；作者的做法是將6個Gene group，又按照細(xì)胞的總UMI計數(shù)將細(xì)胞群分成5個大小相同的組(Cell Group), Cell Group1-5紊服，其中Cell Group1的測序深度最大檀轨，Group2-5其他依次降低胸竞。評價一個有效標(biāo)準(zhǔn)化的數(shù)據(jù)是，每個細(xì)胞對Gene group的方差貢獻(xiàn)率應(yīng)該是一樣的裤园，20%撤师。

結(jié)論：log-normalization這種方法雖然減輕了測序深度和基因表達(dá)量之間的關(guān)系剂府，但我們發(fā)現(xiàn)具有不同總體豐度的基因在log-normalization處理后表現(xiàn)出不同的模式拧揽，并且只有底部三組（Gene group4-group6）的低/中豐度基因被有效地歸一化，高豐度基因歸一化無效腺占。作者在人腦RNA數(shù)據(jù)集也發(fā)現(xiàn)了相同的模式淤袜。默認(rèn)的lognormalize歸一化僅適用于中低表達(dá)基因。

image.png

3.4 將正則化負(fù)二項式回歸應(yīng)用于單細(xì)胞歸一化

作者發(fā)現(xiàn)使用負(fù)二項式分布對單細(xì)胞數(shù)據(jù)進(jìn)行建模會導(dǎo)致過度擬合衰伯，再此不細(xì)說铡羡。
此后，作者引入正則化負(fù)二項式回歸模型應(yīng)用于單細(xì)胞歸一化意鲸，pearson殘差有效地標(biāo)準(zhǔn)化了技術(shù)差異烦周，同時保留了生物變異。說明一下怎顾，pearson殘差可以理解為sctransform標(biāo)準(zhǔn)化读慎。
圖A和B：同3.3上圖，但是y軸是Pearson殘差槐雾，對scRNA-seq 數(shù)據(jù)的歸一化處理明顯優(yōu)于log-normalization夭委。特別是圖B，每個cell group對Gene group的方差貢獻(xiàn)率幾乎一致募强。
圖C：比較log-normalization和Pearson殘差兩種標(biāo)準(zhǔn)化處理方法株灸，計算6個Gene group的Pearson殘差和細(xì)胞總UMI計數(shù)之間Pearson相關(guān)性，繪制箱線圖擎值。Pearson殘差這個方法的中位數(shù)幾乎為0慌烧，也就是說Pearson殘差和總細(xì)胞 UMI計數(shù)（測序深度）幾乎無關(guān)。
結(jié)論：圖A-C都表明鸠儿，與log-normalization相比屹蚊，Pearson殘差的表達(dá)量和方差與測序深度無關(guān)，運(yùn)用正則化NB回歸的Pearson殘差捆交，對scRNA-seq數(shù)據(jù)進(jìn)行了有效的歸一化淑翼。這應(yīng)該消除技術(shù)變異但保留生物異質(zhì)性，同時避免了模型過度擬合數(shù)據(jù)品追。

image.png

四玄括、小結(jié)

4.1 為什么我們要標(biāo)準(zhǔn)化scRNA-seq基因表達(dá)值？

• 我們根據(jù)基因表達(dá)譜的差異對細(xì)胞進(jìn)行聚類
• 基因表達(dá)值的差異應(yīng)反映細(xì)胞間的生物學(xué)變異
• 單細(xì)胞基因表達(dá)值是噪音數(shù)據(jù)
• 測序深度（每個細(xì)胞的UMI數(shù)量）在細(xì)胞之間存在顯著差異
• 標(biāo)準(zhǔn)化的基因表達(dá)值應(yīng)與測序深度無關(guān)

4.2 全局scale標(biāo)準(zhǔn)化(Global scaling normalization)

計算方法：

• 將細(xì)胞中基因的UMI計數(shù)除以該細(xì)胞中的UMI總數(shù)
• 將上面的比值乘以scale 因子（默認(rèn)為10,000）
• 對結(jié)果取自然對數(shù)進(jìn)行轉(zhuǎn)換肉瓦，這將每個細(xì)胞縮放到10,000這個轉(zhuǎn)錄本總數(shù)

問題：全局scale標(biāo)準(zhǔn)化的問題是對高表達(dá)基因不能有效的歸一化

• 默認(rèn)的lognormalize歸一化僅適用于中低表達(dá)基因
• 高表達(dá)基因的表達(dá)值與測序深度相關(guān)

4.3 Sctransform

方法：使用正則化負(fù)二項式回歸對單細(xì)胞 RNA-seq數(shù)據(jù)進(jìn)行歸一化和方差穩(wěn)定化處理
優(yōu)勢：對高表達(dá)基因遭京，Sctransform可以更好地處理這個問題

下面胃惜，我們回答第一個問題：
問：什么情況下SCTransform是優(yōu)于全局scale歸一化方法（global scaling normalization）
答： 對高表達(dá)基因，全局scale標(biāo)準(zhǔn)化不能有效的歸一化哪雕，Sctransform可以更好地處理這個問題船殉；