本教程介紹了Kang等人（2017）的兩組PBMC的對齊方式杀饵。在該實驗中，將PBMC分為刺激組和對照組谬擦，并用干擾素β治療刺激組切距。對干擾素的反應引起細胞類型特異性基因表達的變化，這使得對所有數(shù)據(jù)進行聯(lián)合分析變得困難惨远，并且細胞按刺激條件和細胞類型聚類谜悟。在這里，我們證明了我們的整合策略北秽，如Stuart和Butler等人（2018年）所述葡幸，用于執(zhí)行整合分析以促進常見細胞類型的鑒定并進行比較分析。盡管此示例演示了兩個數(shù)據(jù)集（條件）的集成贺氓，但是這些方法已擴展到多個數(shù)據(jù)集蔚叨。這個工作流程 提供了整合四個胰島數(shù)據(jù)集的示例。
整合目標
以下教程旨在概述使用Seurat集成過程可能進行的復雜細胞類型的比較分析辙培。在這里蔑水，我們解決了三個主要目標：
- 識別兩個數(shù)據(jù)集中都存在的單元格類型
- 獲得在對照和刺激細胞中均保守的細胞類型標記
- 比較數(shù)據(jù)集以查找對刺激的細胞類型特異性反應
設置Seurat對象
基因表達矩陣可以在這里找到。我們首先讀取兩個計數(shù)矩陣并設置Seurat對象扬蕊。
library(Seurat)
library(cowplot)
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")
# Set up control object
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)
# Set up stimulated object
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)
執(zhí)行整合
然后搀别，我們使用FindIntegrationAnchors函數(shù)識別錨點，該函數(shù)將Seurat對象的列表作為輸入尾抑，并使用這些錨點將兩個數(shù)據(jù)集與集成在一起IntegrateData歇父。
immune.anchors <- FindIntegrationAnchors(object.list = list(ctrl, stim), dims = 1:20)
immune.combined <- IntegrateData(anchorset = immune.anchors, dims = 1:20)
進行綜合分析
現(xiàn)在，我們可以在所有單元上運行單個集成分析蛮穿！
DefaultAssay(immune.combined) <- "integrated"
# Run the standard workflow for visualization and clustering
immune.combined <- ScaleData(immune.combined, verbose = FALSE)
immune.combined <- RunPCA(immune.combined, npcs = 30, verbose = FALSE)
# t-SNE and Clustering
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)
# Visualization
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")
識別保守的細胞類型標記
為了確定跨條件保守的規(guī)范細胞類型標記基因践磅，我們提供了該FindConservedMarkers功能。此功能對每個數(shù)據(jù)集/組執(zhí)行差異基因表達測試灸异，并使用MetaDE R軟件包中的薈萃分析方法組合p值府适。例如羔飞，無論簇7中的刺激條件如何，我們都可以計算出保守標記的基因（NK細胞）檐春。
DefaultAssay(immune.combined) <- "RNA"
nk.markers <- FindConservedMarkers(immune.combined, ident.1 = 7, grouping.var = "stim", verbose = FALSE)
head(nk.markers)
## CTRL_p_val CTRL_avg_logFC CTRL_pct.1 CTRL_pct.2 CTRL_p_val_adj
## GNLY 0 4.186117 0.943 0.046 0
## NKG7 0 3.164712 0.953 0.085 0
## GZMB 0 2.915692 0.839 0.044 0
## CLIC3 0 2.407695 0.601 0.024 0
## FGFBP2 0 2.241968 0.500 0.021 0
## CTSW 0 2.088278 0.537 0.030 0
## STIM_p_val STIM_avg_logFC STIM_pct.1 STIM_pct.2 STIM_p_val_adj
## GNLY 0.000000e+00 4.033650 0.955 0.061 0.000000e+00
## NKG7 0.000000e+00 2.914724 0.952 0.082 0.000000e+00
## GZMB 0.000000e+00 3.142391 0.898 0.061 0.000000e+00
## CLIC3 0.000000e+00 2.470769 0.629 0.031 0.000000e+00
## FGFBP2 9.524349e-156 1.483922 0.259 0.016 1.338457e-151
## CTSW 0.000000e+00 2.196390 0.604 0.035 0.000000e+00
## max_pval minimump_p_val
## GNLY 0.000000e+00 0
## NKG7 0.000000e+00 0
## GZMB 0.000000e+00 0
## CLIC3 0.000000e+00 0
## FGFBP2 9.524349e-156 0
## CTSW 0.000000e+00 0
我們可以為每個簇探索這些標記基因逻淌，并使用它們將我們的簇注釋為特定的細胞類型。
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)
DotPlot帶有split.by參數(shù)的函數(shù)可用于跨條件查看保守的細胞類型標記疟暖，顯示表達水平和表達任何給定基因的簇中細胞的百分比卡儒。在這里，我們?yōu)?3個簇中的每一個繪制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()
確定跨條件的差異表達基因
現(xiàn)在骨望，我們已經(jīng)將刺激細胞與對照細胞對齊，我們可以開始進行比較分析欣舵，并觀察刺激引起的差異擎鸠。廣泛觀察這些變化的一種方法是繪制受刺激細胞和對照細胞的平均表達，并在散點圖上尋找視覺異常值的基因缘圈。在這里劣光，我們采用受刺激的和對照的原始T細胞和CD14單核細胞群體的平均表達，并生成散點圖糟把，突出顯示對干擾素刺激表現(xiàn)出戲劇性反應的基因绢涡。
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)垂寥，并且可能代表保守的干擾素應答途徑。
因為我們有信心確定出跨條件的常見細胞類型另锋，所以我們可以詢問相同條件下不同條件下哪些基因會發(fā)生變化滞项。首先，我們在meta.data插槽中創(chuàng)建一列夭坪，以保存單元格類型和刺激信息文判，并將當前標識切換到該列。然后室梅，我們使用它FindMarkers來找到受刺激的B細胞和對照B細胞之間不同的基因戏仓。請注意，此處顯示的許多頂級基因與我們之前繪制的核心干擾素反應基因相同亡鼠。此外赏殃，我們看到的諸如CXCL10的基因?qū)魏思毎虰細胞干擾素的反應也具有特異性，在該列表中也顯示出很高的意義间涵。
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)
## p_val avg_logFC pct.1 pct.2 p_val_adj
## ISG15 8.611499e-155 3.1934171 0.998 0.236 1.210174e-150
## IFIT3 1.319470e-150 3.1195144 0.965 0.053 1.854251e-146
## IFI6 4.716672e-148 2.9264004 0.964 0.078 6.628339e-144
## ISG20 1.061563e-145 2.0390802 1.000 0.664 1.491814e-141
## IFIT1 1.830963e-136 2.8706318 0.909 0.030 2.573053e-132
## MX1 1.775606e-120 2.2540787 0.909 0.118 2.495259e-116
## LY6E 2.824749e-116 2.1460522 0.896 0.153 3.969620e-112
## TNFSF10 4.227184e-109 2.6372382 0.785 0.020 5.940461e-105
## IFIT2 4.627440e-106 2.5102230 0.789 0.038 6.502941e-102
## B2M 1.344345e-94 0.4193618 1.000 1.000 1.889208e-90
## PLSCR1 5.170871e-94 1.9769476 0.794 0.113 7.266624e-90
## IRF7 1.451494e-92 1.7994058 0.838 0.190 2.039785e-88
## CXCL10 6.201621e-84 3.6906104 0.650 0.010 8.715138e-80
## UBE2L6 1.324818e-81 1.4879509 0.854 0.301 1.861767e-77
## PSMB9 1.098134e-76 1.1378896 0.940 0.571 1.543208e-72
可視化基因表達中這些變化的另一種有用方法是split.by選擇FeaturePlot或VlnPlot功能仁热。這將顯示給定基因列表的FeaturePlots，并按分組變量（此處為刺激條件）進行劃分勾哩。諸如CD3D和GNLY的基因是典型的細胞類型標記（對于T細胞和NK / CD8 T細胞）抗蠢，實際上不受干擾素刺激的影響举哟，并且在對照組和受刺激組中顯示出相似的基因表達模式。另一方面迅矛，IFI6和ISG15是核心干擾素反應基因妨猩，因此在所有細胞類型中均被上調(diào)。最后秽褒，CD14和CXCL10是顯示細胞類型特異性干擾素應答的基因壶硅。CD14單核細胞受刺激后，CD14表達下降销斟，這可能導致在有監(jiān)督的分析框架中進行錯誤分類庐椒，從而強調(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)
