歸一化和回歸不必要的變異
現(xiàn)在我們有了高質(zhì)量的細(xì)胞旭愧，我們首先需要探索我們的數(shù)據(jù)并識(shí)別任何不需要的變異來(lái)源。然后我們需要對(duì)數(shù)據(jù)進(jìn)行歸一化置逻，執(zhí)行方差穩(wěn)定并回歸對(duì)我們的數(shù)據(jù)有影響的任何協(xié)變量的影響推沸。
目標(biāo)：
- 準(zhǔn)確標(biāo)準(zhǔn)化和縮放基因表達(dá)值，以解決測(cè)序深度和過(guò)度分散計(jì)數(shù)值的差異券坞。
- 找出最有可能指示不同細(xì)胞類型的變異基因鬓催。
挑戰(zhàn)：
- 檢查并消除不必要的變異，避免下游人為原因?qū)е碌募?xì)胞聚集
建議：
- 在執(zhí)行聚類之前恨锚，要很好地了解您對(duì)要呈現(xiàn)的細(xì)胞類型的期望宇驾。了解你所期望的是低復(fù)雜性的細(xì)胞類型還是高線粒體含量的細(xì)胞類型，以及這些細(xì)胞是否正在分化猴伶。
- 如果需要并且適合于實(shí)驗(yàn)课舍，回歸出UMIs的數(shù)量（默認(rèn)使用 sctransform）、線粒體含量和細(xì)胞周期他挎，因此不會(huì)用于下游聚類
開(kāi)始
讓我們首先創(chuàng)建一個(gè)新腳本進(jìn)行規(guī)范化和集成步驟筝尾。創(chuàng)建一個(gè)新腳本（文件 -> 新文件 -> R 腳本），并將其保存為SCT_integration_analysis.R.
對(duì)于工作流程的其余部分雇盖，我們將主要使用 Seurat 包中提供的函數(shù)忿等。因此，除了 tidyverse 庫(kù)和下面列出的其他一些庫(kù)之外崔挖，我們還需要加載 Seurat 庫(kù)贸街。
# Single-cell RNA-seq - normalization
# Load libraries
library(Seurat)
library(tidyverse)
library(RCurl)
library(cowplot)
此分析的輸入是一個(gè)seurat對(duì)象。我們將使用我們?cè)?QC 課程中創(chuàng)建的名為filtered_seurat.
探索不必要的變異來(lái)源
對(duì)生物協(xié)變量的校正目的是為了挑選出感興趣的特定生物信號(hào)狸相，而技術(shù)協(xié)變量的校正對(duì)于揭示潛在的生物信號(hào)可能至關(guān)重要薛匪。最常見(jiàn)的生物學(xué)數(shù)據(jù)校正是去除細(xì)胞周期對(duì)轉(zhuǎn)錄組的影響。這種數(shù)據(jù)校正可以通過(guò)對(duì)細(xì)胞周期評(píng)分的簡(jiǎn)單線性回歸來(lái)執(zhí)行脓鹃，我們將在下面演示逸尖。
第一步是查看數(shù)據(jù)，看看我們是否觀察到數(shù)據(jù)中的是否有任何影響。細(xì)胞之間的原始計(jì)數(shù)不具有可比性娇跟，我們不能將它們用于我們的探索性分析岩齿。因此，我們將通過(guò)除以每個(gè)細(xì)胞的總計(jì)數(shù)并取自然對(duì)數(shù)來(lái)執(zhí)行粗略的歸一化苞俘。這種標(biāo)準(zhǔn)化僅用于探索我們數(shù)據(jù)中的變化來(lái)源盹沈。
注意：Seurat 最近引入了一種稱為sctransform的新歸一化方法 ，它同時(shí)執(zhí)行方差穩(wěn)定并回歸出不必要的變化吃谣。這是我們?cè)诠ぷ髁鞒讨惺褂玫囊?guī)范化方法乞封。
# Normalize the counts
seurat_phase <- NormalizeData(filtered_seurat)
接下來(lái)，我們獲取這些歸一化數(shù)據(jù)并檢查是否需要數(shù)據(jù)校正方法岗憋。
評(píng)估細(xì)胞周期的影響
為了根據(jù)每個(gè)細(xì)胞的 G2/M 和 S 期標(biāo)記的表達(dá)為其分配一個(gè)分?jǐn)?shù)肃晚，我們可以使用 Seuart 的CellCycleScoring()函數(shù)。此函數(shù)需要輸入規(guī)范標(biāo)記計(jì)算細(xì)胞周期階段分?jǐn)?shù)仔戈。
我們?cè)?code>data文件夾中為您提供了人類細(xì)胞周期標(biāo)記的列表关串，作為一個(gè)名為cycle.rda. 但是，如果您不使用人類數(shù)據(jù)杂穷，我們還有其他物種詳細(xì)介紹了如何獲取其他感興趣生物的細(xì)胞周期標(biāo)記悍缠。
# Load cell cycle markers
load("data/cycle.rda")
# Score cells for cell cycle
seurat_phase <- CellCycleScoring(seurat_phase,
g2m.features = g2m_genes,
s.features = s_genes)
# View cell cycle scores and phases assigned to cells
View(seurat_phase@meta.data)
在對(duì)細(xì)胞的細(xì)胞周期進(jìn)行評(píng)分后，我們想使用 PCA 確定細(xì)胞周期是否是我們數(shù)據(jù)集中變異的主要來(lái)源耐量。要執(zhí)行 PCA飞蚓，我們需要首先選擇高度可變基因，然后對(duì)數(shù)據(jù)進(jìn)行縮放廊蜒。由于高表達(dá)的基因表現(xiàn)出最高的變異量趴拧，我們不希望我們的“高變異基因”只反映高表達(dá)，我們需要縮放數(shù)據(jù)反映出表達(dá)水平的變化山叮。Seurat的ScaleData()函數(shù)進(jìn)行數(shù)據(jù)縮放：
- 調(diào)整每個(gè)基因的表達(dá)著榴，使跨細(xì)胞的平均表達(dá)為 0
- 縮放每個(gè)基因的表達(dá)以使細(xì)胞間的方差為 1
# Identify the most variable genes
seurat_phase <- FindVariableFeatures(seurat_phase,
selection.method = "vst",
nfeatures = 2000,
verbose = FALSE)
# Scale the counts
seurat_phase <- ScaleData(seurat_phase)
注意：對(duì)于
selection.method和nfeatures參數(shù)，指定的值是默認(rèn)設(shè)置屁倔。因此脑又，代碼中可以不包含這些。我們將其包含在此處以提高透明度并告知您正在使用的內(nèi)容锐借。
現(xiàn)在问麸，我們可以執(zhí)行 PCA 分析并將前兩個(gè)主成分繪制出來(lái)。我們還按細(xì)胞周期階段拆分圖钞翔，以評(píng)估相似性或差異性严卖。由于細(xì)胞周期階段，我們沒(méi)有看到大的差異布轿∠剩基于這個(gè)圖来颤，我們將不進(jìn)行回歸，未發(fā)現(xiàn)細(xì)胞周期引起的變化稠肘。
# Perform PCA
seurat_phase <- RunPCA(seurat_phase)
# Plot the PCA colored by cell cycle phase
DimPlot(seurat_phase,
reduction = "pca",
group.by= "Phase",
split.by = "Phase")
細(xì)胞周期階段應(yīng)該何時(shí)回歸福铅？
下面是從處理“細(xì)胞周期評(píng)分和回歸”的 Seurat 案例中截取的兩個(gè) PCA 圖。
第一個(gè)圖類似于我們上面繪制的圖启具，它是回歸之前的 PCA本讥，用于評(píng)估細(xì)胞周期是否在驅(qū)動(dòng) PC1 和 PC2 中發(fā)揮重要作用。
image很明顯鲁冯，在這種情況下，細(xì)胞是按細(xì)胞類型分開(kāi)的色查，所以該案例建議回歸細(xì)胞周期產(chǎn)生的影響薯演。
第二個(gè) PCA 圖是post-regression，顯示了消除我們觀察到的影響方面回歸后的效果秧了。
image
線粒體表達(dá)是另一個(gè)極大影響聚類的因素跨扮。通常，對(duì)線粒體表達(dá)引起的變化進(jìn)行回歸是很有用的验毡。然而衡创，如果線粒體基因表達(dá)的差異代表了一種可能有助于區(qū)分細(xì)胞簇的生物學(xué)現(xiàn)象，那么我們建議不要將其回歸晶通。在本練習(xí)中璃氢，我們可以執(zhí)行類似于查看細(xì)胞周期的快速檢查，并決定是否要對(duì)其進(jìn)行回歸狮辽。
- 首先一也，將線粒體比率變量轉(zhuǎn)換為基于四分位數(shù)的新分類變量（使用以下代碼）：
# Check quartile values
summary(seurat_phase@meta.data$mitoRatio)
# Turn mitoRatio into categorical factor vector based on quartile values
seurat_phase@meta.data$mitoFr <- cut(seurat_phase@meta.data$mitoRatio,
breaks=c(-Inf, 0.0144, 0.0199, 0.0267, Inf),
labels=c("Low","Medium","Medium high", "High"))
接下來(lái)，繪制類似于我們?nèi)绾翁幚砑?xì)胞周期回歸的 PCA喉脖。提示：使用新
mitoFr變量拆分單元格并相應(yīng)地為它們著色椰苟。評(píng)估 #2 中生成的 PCA 圖。
1. 確定您是否觀察到效果树叽。
1\. Describe what you see.
1\. Would you regress out mitochndrial fraction as a source of unwanted variation?
使用 SCTransform 歸一化和回歸出不必要的變化源
現(xiàn)在我們可以使用 sctransform 方法作為更準(zhǔn)確的歸一化方法舆蝴，估計(jì)原始過(guò)濾數(shù)據(jù)的方差，并識(shí)別最易變的基因题诵。sctransform 方法使用正則化負(fù)二項(xiàng)式模型對(duì)UMI 計(jì)數(shù)進(jìn)行建模洁仗，以消除由于測(cè)序深度（每個(gè)細(xì)胞的總 nUMI）引起的變化，同時(shí)根據(jù)具有相似豐度的基因之間的匯集信息進(jìn)而調(diào)整方差（類似于一些bulk RNA-seq 方法）仇轻。
圖片來(lái)源： Hafemeister C 和 Satija R京痢。使用正則化負(fù)二項(xiàng)式回歸對(duì)單細(xì)胞 RNA-seq 數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化和方差穩(wěn)定化，bioRxiv 2019 (https://doi.org/10.1101/576827)
模型的輸出（殘差）是每個(gè)測(cè)試轉(zhuǎn)錄本的標(biāo)準(zhǔn)化表達(dá)水平篷店。
Sctransform 通過(guò)回歸測(cè)序深度 (nUMI) 自動(dòng)考慮細(xì)胞測(cè)序深度祭椰。但是臭家，如果在探索步驟期間在數(shù)據(jù)中發(fā)現(xiàn)了其他無(wú)用的變化來(lái)源，我們也可以包括這些方淤。由于在細(xì)胞周期階段钉赁，我們幾乎沒(méi)有觀察到影響，因此我們選擇不從我們的數(shù)據(jù)中回歸携茂。我們觀察到線粒體表達(dá)的一些影響你踩，因此我們選擇從數(shù)據(jù)中回歸。
為了運(yùn)行 SCTransform讳苦，我們以下面的代碼為例带膜。不要運(yùn)行此代碼，因?yàn)槲覀兏矚g在下面的下一部分中為每個(gè)示例分別運(yùn)行此代碼鸳谜。
## DO NOT RUN CODE ##
# SCTranform
seurat_phase <- SCTransform(seurat_phase, vars.to.regress = c("mitoRatio"))
迭代數(shù)據(jù)集中的樣本
由于我們的數(shù)據(jù)集中有兩個(gè)樣本（來(lái)自兩個(gè)條件）膝藕，我們希望將它們保留為單獨(dú)的對(duì)象并按照集成所需的方式轉(zhuǎn)換它們。我們首先將seurat_phase對(duì)象中的細(xì)胞拆分為“ctrl”和“stim”組：
# Split seurat object by condition to perform cell cycle scoring and SCT on all samples
split_seurat <- SplitObject(seurat_phase, split.by = "sample")
split_seurat <- split_seurat[c("ctrl", "stim")]
通過(guò)使用“for 循環(huán)”在每個(gè)樣本上運(yùn)行SCTransform()咐扭，并在函數(shù)SCTransform()的vars.to.regress參數(shù)中指定來(lái)回歸線粒體表達(dá)芭挽。
在我們運(yùn)行for循環(huán)之前，我們需要知道生成的R對(duì)象需要占用大量的內(nèi)存蝗肪。如果我們有一個(gè)大型數(shù)據(jù)集袜爪，那么我們需要使用以下代碼調(diào)整 R 內(nèi)允許的對(duì)象大小的限制（默認(rèn)為 500 * 1024 ^ 2 = 500 Mb）：
options(future.globals.maxSize = 4000 * 1024^2)
現(xiàn)在，我們運(yùn)行以下循環(huán)對(duì)所有樣本執(zhí)行 sctransform薛闪。這可能需要一些時(shí)間（約 10 分鐘）：
for (i in 1:length(split_seurat)) {
split_seurat[[i]] <- SCTransform(split_seurat[[i]], vars.to.regress = c("mitoRatio"))
}
注意：默認(rèn)情況下辛馆，在歸一化、調(diào)整方差和回歸不感興趣的變異源后逛绵，SCTransform 將按殘差對(duì)基因進(jìn)行排序怀各，并輸出 3000 個(gè)變異最大的基因。如果數(shù)據(jù)集具有更大的細(xì)胞數(shù)术浪，則使用該
variable.features.n參數(shù)將該參數(shù)調(diào)整得更高效果將會(huì)更好瓢对。
請(qǐng)注意，輸出的最后一行指定“將默認(rèn)檢測(cè)設(shè)置為 SCT”胰苏。我們可以查看存儲(chǔ)在 seurat 對(duì)象中的不同分析硕蛹。
# Check which assays are stored in objects
split_seurat$ctrl@assays
現(xiàn)在我們可以看到，除了原始 RNA 計(jì)數(shù)之外硕并，我們現(xiàn)在在我們的assays插槽中有一個(gè) SCT 組件法焰。最可變的特征將是存儲(chǔ)在 SCT 分析中的唯一基因。當(dāng)我們進(jìn)行 scRNA-seq 分析時(shí)倔毙，我們將選擇最合適的檢測(cè)方法用于分析中的不同步驟埃仪。
訓(xùn)練
-
split_seurat對(duì)象內(nèi)的“stim”樣本是否可以進(jìn)行相同的分析？你用來(lái)運(yùn)行的代碼是什么陕赃？ - 對(duì)于“ctrl”與“stim”的“前 10個(gè)feature：”和“最大的10 個(gè)variable feature：”中列出的基因或特征是否有任何觀察卵蛉？
保存對(duì)象颁股！
在完成之前，讓我們將此對(duì)象保存到data/文件夾中傻丝「视校回到這個(gè)階段可能需要一段時(shí)間，尤其是在處理大型數(shù)據(jù)集時(shí)葡缰，最佳做法是將對(duì)象另存為本地可輕松加載的文件亏掀。
# Save the split seurat object
saveRDS(split_seurat, "data/split_seurat.rds")
要將
.rds文件加載回您的環(huán)境，您將使用以下代碼：# Load the split seurat object into the environment split_seurat <- readRDS("data/split_seurat.rds")
