導(dǎo)讀
撰寫本文的主要目的是：整合 處理與對照后的 PBMC（Human peripheral blood mononuclear cell，人外周血單個核細胞） 數(shù)據(jù)集以了解細胞類型特異性反應(yīng)和整合的作用着饥。
本教程介紹了來自 Kang 等人犀农，2017 年 的兩組 PBMC 的比對情況。在這個實驗中宰掉，PBMCs 被分成處理組（刺激組）和對照組呵哨，處理組使用 干擾素β 處理。對干擾素的應(yīng)激導(dǎo)致細胞類型特異性基因表達發(fā)生變化轨奄，這使得對所有數(shù)據(jù)的聯(lián)合分析變得困難孟害。在這里，我們展示了我們的整合策略挪拟，如 Stuart 和 Butler 等人挨务，2018 年 所述，執(zhí)行整合分析以促進常見細胞類型的識別并進行比較分析玉组。雖然此示例只演示了兩個數(shù)據(jù)集（條件）的整合谎柄，但這個方法可以擴展到多個數(shù)據(jù)集。
1. 目的
以下教程旨在為您概述使用 Seurat 整合后對復(fù)雜細胞類型進行的各種比較分析惯雳。
在這里朝巫，我們有以下三個目標：
- 識別兩個數(shù)據(jù)集中都存在的細胞類型
- 獲得對照組和處理組中都保守的細胞類型標記（
markers） - 通過比較數(shù)據(jù)集來尋找對刺激處理產(chǎn)生特異性反應(yīng)的細胞類型
2. 創(chuàng)建對象
基因表達矩陣可以在文末的鏈接找到下載地址，或 點我石景。我們首先讀入兩個計數(shù)矩陣并創(chuàng)建 Seurat 對象劈猿。
# 加載包
library(Seurat)
library(cowplot)
# 讀取數(shù)據(jù)
ctrl.data <- read.table(file = "../data/immune_control_expression_matrix.txt.gz", sep = "\t")
stim.data <- read.table(file = "../data/immune_stimulated_expression_matrix.txt.gz", sep = "\t")
# 創(chuàng)建 對照組 對象
ctrl <- CreateSeuratObject(counts = ctrl.data, project = "IMMUNE_CTRL", min.cells = 5)
ctrl$stim <- "CTRL"
ctrl <- subset(ctrl, subset = nFeature_RNA > 500)
ctrl <- NormalizeData(ctrl, verbose = FALSE)
ctrl <- FindVariableFeatures(ctrl, selection.method = "vst", nfeatures = 2000)
# 創(chuàng)建 處理組 對象
stim <- CreateSeuratObject(counts = stim.data, project = "IMMUNE_STIM", min.cells = 5)
stim$stim <- "STIM"
stim <- subset(stim, subset = nFeature_RNA > 500)
stim <- NormalizeData(stim, verbose = FALSE)
stim <- FindVariableFeatures(stim, selection.method = "vst", nfeatures = 2000)
3. 整合
然后，我們使用 FindIntegrationAnchors 函數(shù)識別 anchors（錨點）潮孽，該函數(shù)將 Seurat 對象列表作為輸入揪荣，并使用這些錨點，利用 IntegrateData 函數(shù)將兩個數(shù)據(jù)集整合在一起往史。
# 識別 anchors
immune.anchors <- FindIntegrationAnchors(object.list = list(ctrl, stim), dims = 1:20)
# 整合
immune.combined <- IntegrateData(anchorset = immune.anchors, dims = 1:20)
4. 整體分析
現(xiàn)在我們可以對所有細胞進行一個綜合分析变逃！
DefaultAssay(immune.combined) <- "integrated"
# 運行可視化和聚類的標準工作流程
immune.combined <- ScaleData(immune.combined, verbose = FALSE)
immune.combined <- RunPCA(immune.combined, npcs = 30, verbose = FALSE)
# t-SNE and 聚類
immune.combined <- RunUMAP(immune.combined, reduction = "pca", dims = 1:20)
immune.combined <- FindNeighbors(immune.combined, reduction = "pca", dims = 1:20)
immune.combined <- FindClusters(immune.combined, resolution = 0.5)
# 可視化
p1 <- DimPlot(immune.combined, reduction = "umap", group.by = "stim")
p2 <- DimPlot(immune.combined, reduction = "umap", label = TRUE)
plot_grid(p1, p2)
為了并排可視化這兩個條件，我們可以使用 split.by 參數(shù)來顯示按簇著色的每個條件怠堪。
DimPlot(immune.combined, reduction = "umap", split.by = "stim")
5. 保守marker鑒定
為了識別跨條件保守的細胞類型標記基因揽乱，Seurat 提供了 FindConservedMarkers 函數(shù)。此函數(shù)對每個數(shù)據(jù)集或組執(zhí)行差異基因表達 test粟矿，并使用 MetaDE R 包中的元分析（Meta-analysis）方法組合 p 值凰棉。例如，我們可以計算簇 7（NK 細胞）中無論條件如何陌粹，都是保守的標記基因撒犀。
DefaultAssay(immune.combined) <- "RNA"
nk.markers <- FindConservedMarkers(immune.combined, ident.1 = 7, grouping.var = "stim", verbose = FALSE)
head(nk.markers)
我們可以探索每個簇的這些標記基因，并使用它們掏秩，將我們的簇注釋為特定的細胞類型或舞。
FeaturePlot(immune.combined, features = c("CD3D", "SELL", "CREM", "CD8A", "GNLY", "CD79A", "FCGR3A",
"CCL2", "PPBP"), min.cutoff = "q9")
immune.combined <- RenameIdents(immune.combined, `0` = "CD14 Mono", `1` = "CD4 Naive T", `2` = "CD4 Memory T",
`3` = "CD16 Mono", `4` = "B", `5` = "CD8 T", `6` = "T activated", `7` = "NK", `8` = "DC", `9` = "B Activated",
`10` = "Mk", `11` = "pDC", `12` = "Eryth", `13` = "Mono/Mk Doublets")
DimPlot(immune.combined, label = TRUE)
帶有 split.by 參數(shù)的 DotPlot 函數(shù)可用于查看跨條件的保守細胞類型標記，顯示表達水平和表達任何給定基因的簇中細胞的百分比蒙幻。在這里映凳，我們?yōu)?13 個聚類中的每一個繪制了 2-3 個強標記基因。
Idents(immune.combined) <- factor(Idents(immune.combined), levels = c("Mono/Mk Doublets", "pDC",
"Eryth", "Mk", "DC", "CD14 Mono", "CD16 Mono", "B Activated", "B", "CD8 T", "NK", "T activated",
"CD4 Naive T", "CD4 Memory T"))
markers.to.plot <- c("CD3D", "CREM", "HSPH1", "SELL", "GIMAP5", "CACYBP", "GNLY", "NKG7", "CCL5",
"CD8A", "MS4A1", "CD79A", "MIR155HG", "NME1", "FCGR3A", "VMO1", "CCL2", "S100A9", "HLA-DQA1",
"GPR183", "PPBP", "GNG11", "HBA2", "HBB", "TSPAN13", "IL3RA", "IGJ")
DotPlot(immune.combined, features = rev(markers.to.plot), cols = c("blue", "red"), dot.scale = 8,
split.by = "stim") + RotatedAxis()
6. 跨條件識別差異表達基因
現(xiàn)在我們已經(jīng)對齊了受刺激組（處理組）細胞和對照組細胞邮破，我們可以開始進行比較分析并查看刺激引起的差異诈豌。觀察這些變化的一種方法是繪制受刺激細胞和對照細胞的平均表達，并在散點圖上尋找異常值的基因抒和。在這里矫渔，我們?nèi)∈艽碳ず蛯φ?naive T 細胞和 CD14 單核細胞群的平均表達，并生成散點圖摧莽，突出顯示對干擾素刺激有顯著反應(yīng)的基因庙洼。
t.cells <- subset(immune.combined, idents = "CD4 Naive T")
Idents(t.cells) <- "stim"
avg.t.cells <- log1p(AverageExpression(t.cells, verbose = FALSE)$RNA)
avg.t.cells$gene <- rownames(avg.t.cells)
cd14.mono <- subset(immune.combined, idents = "CD14 Mono")
Idents(cd14.mono) <- "stim"
avg.cd14.mono <- log1p(AverageExpression(cd14.mono, verbose = FALSE)$RNA)
avg.cd14.mono$gene <- rownames(avg.cd14.mono)
genes.to.label = c("ISG15", "LY6E", "IFI6", "ISG20", "MX1", "IFIT2", "IFIT1", "CXCL10", "CCL8")
p1 <- ggplot(avg.t.cells, aes(CTRL, STIM)) + geom_point() + ggtitle("CD4 Naive T Cells")
p1 <- LabelPoints(plot = p1, points = genes.to.label, repel = TRUE)
p2 <- ggplot(avg.cd14.mono, aes(CTRL, STIM)) + geom_point() + ggtitle("CD14 Monocytes")
p2 <- LabelPoints(plot = p2, points = genes.to.label, repel = TRUE)
plot_grid(p1, p2)
正如您所看到的，許多相同的基因在這兩種細胞類型中都上調(diào)镊辕，并且可能代表了一種保守的干擾素反應(yīng)途徑油够。
因為我們有信心在不同條件下識別出常見的細胞類型，所以我們可以查看相同類型細胞在不同條件下哪些基因會發(fā)生變化丑蛤。首先叠聋，我們在 meta.data 中創(chuàng)建一個列來保存細胞類型和處理信息，并將當(dāng)前標識切換到該列受裹。然后使用 FindMarkers 來查找受刺激 B 細胞和對照 B 細胞之間不同的基因碌补。請注意，此處顯示的許多 Top 基因與我們之前繪制的核心干擾素反應(yīng)基因相同棉饶。此外厦章，我們看到的 CXCL10 等特定于單核細胞和 B 細胞干擾素反應(yīng)的基因在此列表中也顯示出非常重要的意義。
immune.combined$celltype.stim <- paste(Idents(immune.combined), immune.combined$stim, sep = "_")
immune.combined$celltype <- Idents(immune.combined)
Idents(immune.combined) <- "celltype.stim"
b.interferon.response <- FindMarkers(immune.combined, ident.1 = "B_STIM", ident.2 = "B_CTRL", verbose = FALSE)
head(b.interferon.response, n = 15)
另一種可視化基因表達變化的方法是使用 FeaturePlot 或 VlnPlot 函數(shù)的 split.by 選項照藻。這將顯示給定基因列表的特征圖袜啃，按分組變量（此處為刺激條件）拆分。CD3D 和 GNLY 等基因是典型的細胞類型標記（用于 T 細胞和 NK/CD8 T 細胞）幸缕，它們幾乎不受干擾素刺激的影響群发，并且在對照組和受刺激組中顯示出相似的基因表達模式晰韵。另一方面，IFI6 和 ISG15 是核心干擾素反應(yīng)基因熟妓，并在所有細胞類型中上調(diào)雪猪。最后，CD14 和 CXCL10 是顯示細胞類型特異性干擾素反應(yīng)的基因起愈。刺激 CD14 單核細胞后 CD14 表達降低只恨，這可能導(dǎo)致監(jiān)督分析框架中發(fā)生錯誤分類，強調(diào)了整合分析的價值抬虽。CXCL10 在干擾素刺激后在單核細胞和 B 細胞中顯示出明顯的上調(diào)官觅，但在其他細胞類型中則沒有。
FeaturePlot(immune.combined, features = c("CD3D", "GNLY", "IFI6"), split.by = "stim", max.cutoff = 3,
cols = c("grey", "red"))
plots <- VlnPlot(immune.combined, features = c("LYZ", "ISG15", "CXCL10"), split.by = "stim", group.by = "celltype",
pt.size = 0, combine = FALSE)
CombinePlots(plots = plots, ncol = 1)
本文由mdnice多平臺發(fā)布
