學(xué)習(xí)目標(biāo)：

• 評估是否存在聚類過程產(chǎn)生的技術(shù)誤差
• 使用 PCA 和 UMAP 圖確定聚類質(zhì)量，并了解何時重新聚類
• 評估已知的細胞類型標(biāo)記與假設(shè)簇的細胞類型同一性

單細胞 RNA-seq 聚類分析

現(xiàn)在我們已經(jīng)進行了整合，我們想知道我們的細胞群中存在哪些不同細胞類型。

image

目標(biāo)：

• *生成特定于細胞類型的簇，并使用已知的標(biāo)記來確定簇的身份。
• 確定分群是否代表真實的細胞類型或由于生物或技術(shù)差異而形成的群集 ，例如處于細胞周期S期的細胞群集、特定批次的群集或高線粒體含量的細胞。

挑戰(zhàn)：

• 識別每個分群的細胞類型
• 保持耐心，因為這可能是聚類和標(biāo)記識別之間的高度迭代過程（有時甚至回到 QC 過濾或標(biāo)準(zhǔn)化）

建議：

• 對存在的細胞類型和這些細胞類型的幾個標(biāo)記基因有一個很好的預(yù)期。了解你所期望的是低復(fù)雜性的細胞類型還是高線粒體含量的細胞類型，以及這些細胞是否正在分化。
• 確定要刪除的任何垃圾群。垃圾群可能是包含高線粒體含量和低UMIs/基因的那些。
• 如果未將所有細胞類型檢測為單獨的群 ，請嘗試更改UMAP分辨率、用于分群的PC數(shù)量或使用的可變基因數(shù)量。

探討質(zhì)制指標(biāo)

為了確定我們的分群是否可能是由于細胞周期階段或線粒體表達等人為因素造成的，可視化探索這些指標(biāo)以查看是否有任何簇表現(xiàn)出富集或與其他簇不同，這是很有用的。但是，如果觀察到特定簇的富集或差異，它可以用細胞類型來解釋，那就可以不必?fù)?dān)憂。
為了探索和可視化各種質(zhì)量指標(biāo)，我們將使用來自Seurat的通用的DimPlot() 和FeaturePlot() 函數(shù)。

按樣本分離聚類

我們可以從每個樣本中每個簇的細胞分布開始：

# Extract identity and sample information from seurat object to determine the number of cells per cluster per sample
n_cells <- FetchData(seurat_integrated, 
                     vars = c("ident", "orig.ident")) %>%
        dplyr::count(ident, orig.ident) %>%
        tidyr::spread(ident, n)

# View table
View(n_cells)

image

我們可以使用 UMAP 對每個樣本的每個簇的細胞進行可視化：

# UMAP of cells in each cluster by sample
DimPlot(seurat_integrated, 
        label = TRUE, 
        split.by = "sample")  + NoLegend()

image

一般來說，我們希望在所有條件下，大多數(shù)細胞類型簇都存在；然而，根據(jù)實驗，我們可能會看到一些特定條件的細胞類型存在。這些簇在不同條件下看起來非常相似，這很好，因為我們預(yù)計在對照和刺激條件下都存在類似的細胞類型。

按細胞周期分離簇

接下來，我們將探索細胞是否會因不同的細胞周期階段而出現(xiàn)聚集。當(dāng)我們對無意義的變異源進行SCTransform歸一化和回歸時，并沒有因為細胞周期階段而使變異消退。如果我們的細胞簇在線粒體表達上表現(xiàn)出很大的差異，這預(yù)示著我們要重新運行SCTransform，并將 S.Score 和 G2M.Score 添加到我們的變量中以進行回歸，然后重新運行其余步驟。
接下來，我們將探討細胞是否會因不同的細胞周期出現(xiàn)階段聚集。

# Explore whether clusters segregate by cell cycle phase
DimPlot(seurat_integrated,
        label = TRUE, 
        split.by = "Phase")  + NoLegend()

image

我們看不到細胞周期評分的太多聚類，因此我們可以繼續(xù)進行 QC。

按照各種意義的變化來源分離聚類

接下來，我們將探索其他指標(biāo)，例如每個細胞的 UMI 和基因數(shù)量、S 期和 G2M 期標(biāo)記，以及 UMAP 的線粒體基因表達。查看各自的 S 和 G2M 分?jǐn)?shù)可以為我們提供額外的信息來檢查細胞周期因素的影響，就像我們之前所做的那樣。

# Determine metrics to plot present in seurat_integrated@meta.data
metrics <-  c("nUMI", "nGene", "S.Score", "G2M.Score", "mitoRatio")

FeaturePlot(seurat_integrated, 
            reduction = "umap", 
            features = metrics,
            pt.size = 0.4, 
            order = TRUE,
            min.cutoff = 'q10',
            label = TRUE)

image

注：該order參數(shù)將在低表達細胞上方繪制高表達細胞，而min.cutoff參數(shù)將確定著色的閾值。 min.cutoff的 q10意味著基因表達最低的10%的細胞，不會表現(xiàn)出任何紫色陰影（完全灰色）

指標(biāo)在整個群集中相對均勻，但nUMI和nGene在群集3、9、14和15中顯示出更高的值，也許還包括群集17。我們將密切關(guān)注這些群集，看看細胞類型是否可以解釋這種增加。

如果我們在此時看到與這些指標(biāo)對應(yīng)的差異，那么我們通常會注意到它們，然后在確定細胞類型標(biāo)識后決定是否采取進一步的行動。

探索驅(qū)動不同集群的PC

我們還可以探索群集通過不同的PC分開的情況。我們希望定義的PC能夠很好地區(qū)分細胞類型。要可視化這個信息，我們需要提取細胞的UMAP坐標(biāo)信息及其對應(yīng)的每個PC分?jǐn)?shù)，以便通過UMAP查看每個PC。
首先，我們確定要從Seurat對象提取的信息，然后，可以使用 FetchData() 函數(shù)提取它。

# Defining the information in the seurat object of interest
columns <- c(paste0("PC_", 1:16),
            "ident",
            "UMAP_1", "UMAP_2")

# Extracting this data from the seurat object
pc_data <- FetchData(seurat_integrated, 
                     vars = columns)

注意：我們?nèi)绾沃涝?code>FetchData()函數(shù)中包含UMAP_1以獲取 UMAP 坐標(biāo)？所述seurat的cheatsheet描述了該函數(shù)為能夠從表達式矩陣，細胞的嵌入，或元數(shù)據(jù)中提取任何數(shù)據(jù)。

例如，如果您瀏覽seurat_integrated@reductions列表對象，第一個組件用于 PCA，并包含一個用于cell.embeddings. 我們可以使用列名（PC_1、PC_2、PC_3等）為每個 PC 提取與每個細胞對應(yīng)的坐標(biāo)或 PC 分?jǐn)?shù)。

我們可以為 UMAP 做同樣的事情：

# Extract the UMAP coordinates for the first 10 cells
seurat_integrated@reductions$umap@cell.embeddings[1:10, 1:2]

該FetchData()函數(shù)可以讓我們更容易地提取數(shù)據(jù)。

在下面的 UMAP 圖中，細胞按每個主成分的 PC 分?jǐn)?shù)著色。

讓我們?yōu)g覽一下排名前 16 的 PC：

# Adding cluster label to center of cluster on UMAP
umap_label <- FetchData(seurat_integrated, 
                        vars = c("ident", "UMAP_1", "UMAP_2"))  %>%
  group_by(ident) %>%
  summarise(x=mean(UMAP_1), y=mean(UMAP_2))

# Plotting a UMAP plot for each of the PCs
map(paste0("PC_", 1:16), function(pc){
        ggplot(pc_data, 
               aes(UMAP_1, UMAP_2)) +
                geom_point(aes_string(color=pc), 
                           alpha = 0.7) +
                scale_color_gradient(guide = FALSE, 
                                     low = "grey90", 
                                     high = "blue")  +
                geom_text(data=umap_label, 
                          aes(label=ident, x, y)) +
                ggtitle(pc)
}) %>% 
        plot_grid(plotlist = .)

image

我們可以看到如何由不同的 PC 表示的集群。例如，PC_2基因在簇 6、11 和 17 中表現(xiàn)出更高的表達（在 15 中也可能更高一些）。我們可以回顧一下驅(qū)動這臺 PC 的基因，以了解細胞類型可能是什么：

# Examine PCA results 
print(seurat_integrated[["pca"]], dims = 1:5, nfeatures = 5)

image

以 CD79A 和 CD74 基因以及 HLA 基因作為 PC_2的陽性標(biāo)記，我們可以假設(shè)簇6、11 和17 對應(yīng)于 B 細胞。這只是暗示了集群身份可能是什么，集群的身份是通過 PC 的組合來確定的。

為了真正確定集群的身份以及resolution是否合適，探索一些已知的預(yù)期細胞類型的基因標(biāo)記是有幫助的。

探索已知的細胞類型標(biāo)記

隨著細胞的聚集，我們可以通過尋找已知的標(biāo)記來探索細胞類型的身份。顯示了帶有簇標(biāo)記的 UMAP 圖，然后是預(yù)期的不同細胞類型。

DimPlot(object = seurat_integrated, 
        reduction = "umap", 
        label = TRUE) + NoLegend()

image

seurat的FeaturePlot()函數(shù)可以使用存儲在 Seurat 對象中的基因 ID 輕松可視化少數(shù)基因。我們可以在 UMAP 可視化之上輕松探索已知基因標(biāo)記的表達。讓我們通過并確定集群的身份。為了訪問所有基因的標(biāo)準(zhǔn)化表達水平，我們可以使用存儲在RNA檢測槽中的標(biāo)準(zhǔn)化計數(shù)數(shù)據(jù)。

注意： SCTransform 歸一化僅對 3000 個最易變的基因執(zhí)行，因此我們感興趣的許多基因可能不存在于該數(shù)據(jù)中。

# Select the RNA counts slot to be the default assay
DefaultAssay(seurat_integrated) <- "RNA"

# Normalize RNA data for visualization purposes
seurat_integrated <- NormalizeData(seurat_integrated, verbose = FALSE)

注意： Assay 是在 Seurat 對象中定義的一個插槽，其中有多個插槽。在給定的檢測中，counts槽存儲非標(biāo)準(zhǔn)化的原始計數(shù)，data槽存儲標(biāo)準(zhǔn)化的表達數(shù)據(jù)。因此，當(dāng)我們運行NormalizeData()上述代碼中的函數(shù)時，歸一化數(shù)據(jù)將存儲在dataRNA 分析的槽中，而counts槽將保持不變。

根據(jù)我們感興趣的標(biāo)記，它們可能是特定細胞類型的陽性或陰性標(biāo)記。我們選擇的少數(shù)標(biāo)記的組合表達應(yīng)該讓我們了解一個集群是否對應(yīng)于特定的細胞類型。

對于此處使用的標(biāo)記，我們正在尋找陽性標(biāo)記和標(biāo)記在整個集群中的表達一致性。例如，如果一種細胞類型有兩個標(biāo)記，而一個簇中只有其中一個被表達——那么我們不能可靠地將該簇分配給細胞類型。

CD14+ 單核細胞標(biāo)志物

FeaturePlot(seurat_integrated, 
            reduction = "umap", 
            features = c("CD14", "LYZ"), 
            order = TRUE,
            min.cutoff = 'q10', 
            label = TRUE)

image

CD14+ 單核細胞似乎對應(yīng)于簇 1、3 和 14。我們不會包括簇 9 和 15，因為它們不高度表達這兩種標(biāo)記。

FCGR3A+ 單核細胞標(biāo)志物

FeaturePlot(seurat_integrated, 
            reduction = "umap", 
            features = c("FCGR3A", "MS4A7"), 
            order = TRUE,
            min.cutoff = 'q10', 
            label = TRUE)

image

FCGR3A+ 單核細胞標(biāo)記明顯突出了簇 9，盡管我們確實在簇 1、3 和 14 中看到了一些不錯的表達。我們希望在執(zhí)行標(biāo)記識別時看到 FCGR3A+ 細胞的其他標(biāo)記出現(xiàn)。

巨噬細胞

FeaturePlot(seurat_integrated, 
            reduction = "umap", 
            features = c("MARCO", "ITGAM", "ADGRE1"), 
            order = TRUE,
            min.cutoff = 'q10', 
            label = TRUE)

image

我們沒有看到我們的標(biāo)記有太多重疊，因此似乎沒有與巨噬細胞對應(yīng)的簇；也許細胞培養(yǎng)條件對巨噬細胞進行了負(fù)面選擇（更高的粘附性）。

常規(guī)樹突狀細胞標(biāo)志物

FeaturePlot(seurat_integrated, 
            reduction = "umap", 
            features = c("FCER1A", "CST3"), 
            order = TRUE,
            min.cutoff = 'q10', 
            label = TRUE)

image

對應(yīng)于常規(guī)樹突細胞的標(biāo)記物識別簇 15（兩個標(biāo)記物始終顯示表達）。

漿細胞樣樹突狀細胞標(biāo)志物

FeaturePlot(seurat_integrated, 
            reduction = "umap", 
            features = c("IL3RA", "GZMB", "SERPINF1", "ITM2C"), 
            order = TRUE,
            min.cutoff = 'q10', 
            label = TRUE)

image

漿細胞樣樹突狀細胞代表簇 19。雖然這些標(biāo)記的表達存在很多差異，但我們看到簇 19 始終強烈表達。

鍛煉

假設(shè)表中每個不同集群對應(yīng)的集群：

注意：如果任何集群似乎包含兩種不同的細胞類型，則提高我們的集群分辨率以正確地對集群進行子集化是有幫助的?；蛘?，如果我們?nèi)匀粺o法使用增加的分辨率分離出簇，那么可能我們使用的主成分太少，以至于我們沒有分離出這些感興趣的細胞類型。為了告知我們對 PC 的選擇，我們可以查看與 UMAP 圖重疊的 PC 基因表達，并確定我們的細胞群是否被包含的 PC 分開。

現(xiàn)在我們對大多數(shù)集群對應(yīng)的細胞類型有了一個不錯的想法，但仍然存在一些問題：

1. 集群 7 和 20 的細胞類型標(biāo)識是什么？
2. 對應(yīng)于相同細胞類型的簇是否具有生物學(xué)意義的差異？是否存在這些細胞類型的亞群？
3. 通過識別這些簇的其他標(biāo)記基因，我們能否對這些細胞類型身份獲得更高的可信度？

標(biāo)記識別分析可以幫助我們解決所有這些問題??！

下一步將是進行標(biāo)記識別分析，這將輸出聚類之間表達顯著差異的基因。使用這些基因，我們可以確定或提高對集群/亞群群識別的信心。