scRNA-seq整合簡介
對兩個或兩個以上單細(xì)胞數(shù)據(jù)集的整合分析提出了獨特的挑戰(zhàn)。特別是，在標(biāo)準(zhǔn)工作流下，識別存在于多個數(shù)據(jù)集中的基因可能存在問題蛇更。Seurat v4 包括一組方法，以匹配（或"對齊"）跨數(shù)據(jù)集共同的基因。這些方法首先識別處于匹配生物狀態(tài)的交叉數(shù)據(jù)集對（"錨點"）派任，既可用于糾正數(shù)據(jù)集之間的技術(shù)差異（即批量效應(yīng)校正）砸逊，又可用于對整個實驗條件進行比較scRNA-seq分析。
下面掌逛，我們演示了 ScRNA-seq 集成的方法师逸，在ctrl或干擾素刺激狀態(tài)下對人體免疫細(xì)胞 （PBMC） 進行比較分析。
整合目標(biāo)
以下教程旨在為您概述使用 Seurat 集成程序可能的復(fù)雜細(xì)胞類型的比較分析類型豆混。在這里篓像，我們討論幾個關(guān)鍵目標(biāo)：
- 創(chuàng)建"整合"數(shù)據(jù)分析，用于下游分析
- 識別兩個數(shù)據(jù)集中存在的細(xì)胞類型
- 獲取在對照和刺激細(xì)胞中保守的細(xì)胞類型標(biāo)記
- 比較數(shù)據(jù)集皿伺，找到細(xì)胞型對刺激的特定反應(yīng)
設(shè)置seurat對象
為了方便起見员辩，我們通過我們的SeuratData包分發(fā)此數(shù)據(jù)集。
library(Seurat)
library(SeuratData)
library(patchwork)
# install dataset
InstallData("ifnb")
# load dataset
LoadData("ifnb")
# split the dataset into a list of two seurat objects (stim and CTRL)
ifnb.list <- SplitObject(ifnb, split.by = "stim")
# normalize and identify variable features for each dataset independently
ifnb.list <- lapply(X = ifnb.list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})
# select features that are repeatedly variable across datasets for integration
features <- SelectIntegrationFeatures(object.list = ifnb.list)
執(zhí)行整合
然后鸵鸥，我們使用該功能識別錨點奠滑，以 Seurat 對象列表作為輸入，并使用這些錨點將兩個數(shù)據(jù)集集成在一起妒穴。
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, anchor.features = features)
# this command creates an 'integrated' data assay
immune.combined <- IntegrateData(anchorset = immune.anchors)
執(zhí)行整合分析
現(xiàn)在宋税，我們可以對所有細(xì)胞進行單次整合分析！
# specify that we will perform downstream analysis on the corrected data note that the original
# unmodified data still resides in the 'RNA' assay
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)
immune.combined <- RunUMAP(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindNeighbors(immune.combined, reduction = "pca", dims = 1:30)
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, repel = TRUE)
p1 + p2
為了并排可視化這兩個條件讼油，我們可以使用參數(shù)按群著色每個條件杰赛。
DimPlot(immune.combined, reduction = "umap", split.by = "stim")
識別保守的細(xì)胞類型標(biāo)記
為了識別在條件下保守的細(xì)胞類型標(biāo)記基因，我們提供該功能汁讼。該功能為每個數(shù)據(jù)集/組執(zhí)行基因表達(dá)測試淆攻，并使用 MetaDE R 包中的元分析方法組合 p 值阔墩。例如嘿架，我們可以計算出無論刺激條件如何，第6組（NK細(xì)胞）中保守標(biāo)記的基因啸箫。
# For performing differential expression after integration, we switch back to the original data
DefaultAssay(immune.combined) <- "RNA"
nk.markers <- FindConservedMarkers(immune.combined, ident.1 = 6, grouping.var = "stim", verbose = FALSE)
head(nk.markers)
## CTRL_p_val CTRL_avg_log2FC CTRL_pct.1 CTRL_pct.2 CTRL_p_val_adj
## GNLY 0 6.006422 0.944 0.045 0
## FGFBP2 0 3.223246 0.503 0.020 0
## CLIC3 0 3.466418 0.599 0.024 0
## PRF1 0 2.654683 0.424 0.017 0
## CTSW 0 2.991829 0.533 0.029 0
## KLRD1 0 2.781453 0.497 0.019 0
## STIM_p_val STIM_avg_log2FC STIM_pct.1 STIM_pct.2 STIM_p_val_adj
## GNLY 0.000000e+00 5.853573 0.956 0.060 0.000000e+00
## FGFBP2 7.275492e-161 2.200379 0.260 0.016 1.022425e-156
## CLIC3 0.000000e+00 3.549919 0.627 0.031 0.000000e+00
## PRF1 0.000000e+00 4.102686 0.862 0.057 0.000000e+00
## CTSW 0.000000e+00 3.139620 0.596 0.035 0.000000e+00
## KLRD1 0.000000e+00 2.880055 0.558 0.027 0.000000e+00
## max_pval minimump_p_val
## GNLY 0.000000e+00 0
## FGFBP2 7.275492e-161 0
## CLIC3 0.000000e+00 0
## PRF1 0.000000e+00 0
## CTSW 0.000000e+00 0
## KLRD1 0.000000e+00 0
我們可以探索每個群的這些標(biāo)記基因耸彪，并使用它們來注釋我們的簇為特定的細(xì)胞類型。
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` = "NK", `7` = "T activated", `8` = "DC", `9` = "B Activated",
`10` = "Mk", `11` = "pDC", `12` = "Eryth", `13` = "Mono/Mk Doublets", `14` = "HSPC")
DimPlot(immune.combined, label = TRUE)
參數(shù)的功能可用于查看跨條件下的保存細(xì)胞類型標(biāo)記忘苛，顯示表示任何給定基因在群中的表達(dá)水平和細(xì)胞百分比蝉娜。在這里，我們?yōu)?4個cluster中的每一個繪制2-3個強標(biāo)記基因扎唾。
Idents(immune.combined) <- factor(Idents(immune.combined), levels = c("HSPC", "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", "PRSS57")
DotPlot(immune.combined, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8, split.by = "stim") +
RotatedAxis()
識別不同條件的差異基因
現(xiàn)在召川，我們已經(jīng)對齊了刺激和對照細(xì)胞，我們可以開始做比較分析胸遇，看看刺激引起的差異荧呐。觀察這些變化的一種方法是繪制受刺激細(xì)胞和對照細(xì)胞的平均表達(dá)，并尋找在散點圖上離群的基因。在這里倍阐，我們采取刺激和對照組的幼稚T細(xì)胞和CD14+單核細(xì)胞群的平均表達(dá)概疆，并產(chǎn)生散點圖，突出顯示干擾素刺激劇烈反應(yīng)的基因峰搪。
library(ggplot2)
library(cowplot)
theme_set(theme_cowplot())
t.cells <- subset(immune.combined, idents = "CD4 Naive T")
Idents(t.cells) <- "stim"
avg.t.cells <- as.data.frame(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 <- as.data.frame(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)
p1 + p2
正如你所看到的岔冀，許多相同的基因在這兩種細(xì)胞類型中都得到了調(diào)節(jié)，并且可能代表一個保守的干擾素反應(yīng)途徑概耻。
我們有信心在不同條件下識別出常見的細(xì)胞類型使套，我們可以查看同類型細(xì)胞的基因在不同條件下有什么變化。首先咐蚯，我們在meta.data 插槽中創(chuàng)建一個列童漩，以同時保存細(xì)胞類型和刺激信息，并將當(dāng)前標(biāo)識切換到該列春锋。然后矫膨，我們用來尋找受刺激和對照B細(xì)胞之間的不同基因。請注意期奔，此處顯示的許多top基因與我們之前繪制的核心干擾素響應(yīng)基因相同侧馅。此外，我們所看到的CXCL10等基因是單核細(xì)胞和B細(xì)胞干擾素反應(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)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## ISG15 5.398167e-155 4.5889194 0.998 0.240 7.586044e-151
## IFIT3 2.209577e-150 4.5032297 0.964 0.052 3.105118e-146
## IFI6 7.060888e-150 4.2375542 0.969 0.080 9.922666e-146
## ISG20 7.147214e-146 2.9387415 1.000 0.672 1.004398e-141
## IFIT1 7.650201e-138 4.1295888 0.914 0.032 1.075083e-133
## MX1 1.124186e-120 3.2883709 0.905 0.115 1.579819e-116
## LY6E 2.504364e-118 3.1297866 0.900 0.152 3.519383e-114
## TNFSF10 9.454398e-110 3.7783774 0.791 0.025 1.328627e-105
## IFIT2 1.672384e-105 3.6569980 0.783 0.035 2.350201e-101
## B2M 5.564362e-96 0.6100242 1.000 1.000 7.819599e-92
## PLSCR1 1.128239e-93 2.8205802 0.796 0.117 1.585514e-89
## IRF7 6.602529e-92 2.5832239 0.838 0.190 9.278534e-88
## CXCL10 4.402118e-82 5.2406913 0.639 0.010 6.186297e-78
## UBE2L6 2.995453e-81 2.1487435 0.852 0.300 4.209510e-77
## PSMB9 1.755809e-76 1.6379482 0.940 0.573 2.467438e-72
可視化基因表達(dá)變化的另一個有用的方法。顯示給定基因列表的特征圖肺孤，按分組變量（此處為刺激對照）進行拆分罗晕。CD3D和GNLY等基因是細(xì)胞類型標(biāo)記（用于T細(xì)胞和NK/CD8 T細(xì)胞），它們幾乎不受干擾素刺激的影響赠堵，并在對照和刺激組中顯示類似的基因表達(dá)模式小渊。另一方面，IFI6和ISG15是核心干擾素反應(yīng)基因茫叭，在所有細(xì)胞類型中都受到相應(yīng)的調(diào)節(jié)酬屉。最后，CD14和CXCL10是顯示細(xì)胞類型特定干擾素反應(yīng)的基因揍愁。CD14表達(dá)在CD14單核細(xì)胞刺激后減少呐萨，這可能導(dǎo)致監(jiān)督分析框架中的誤分類，強調(diào)整合分析的價值莽囤。CXCL10顯示干擾素刺激后單核細(xì)胞和B細(xì)胞的明顯上升調(diào)節(jié)谬擦，但其他細(xì)胞類型不顯示。
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)
wrap_plots(plots = plots, ncol = 1)
使用sctransform整合數(shù)據(jù)
介紹一個改進的整合方法朽缎，該方法被命名為"sctransform"惨远。
下面蔚舀，我們演示如何修改 Seurat 整合工作流，以實現(xiàn)用sctransform工作流標(biāo)準(zhǔn)化數(shù)據(jù)集锨络。命令大致相似赌躺，有幾個關(guān)鍵差異：
- 通過SCTransform()單獨（而不是在整合之前）實現(xiàn)數(shù)據(jù)集標(biāo)準(zhǔn)化
- 我們使用 3，000 個或更多基因來進行sctransform 下游分析羡儿。
- 在識別錨點之前運行該功能[PrepSCTIntegration()]
- 運行[FindIntegrationAnchors()]和[IntegrateData()]時礼患，設(shè)置參數(shù)到normalization.method為
SCT - 當(dāng)運行基于 sctransform 的工作流（包括整合）時，不運行該功能
[ScaleData()]
LoadData("ifnb")
ifnb.list <- SplitObject(ifnb, split.by = "stim")
ifnb.list <- lapply(X = ifnb.list, FUN = SCTransform)
features <- SelectIntegrationFeatures(object.list = ifnb.list, nfeatures = 3000)
ifnb.list <- PrepSCTIntegration(object.list = ifnb.list, anchor.features = features)
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, normalization.method = "SCT",
anchor.features = features)
immune.combined.sct <- IntegrateData(anchorset = immune.anchors, normalization.method = "SCT")
immune.combined.sct <- RunPCA(immune.combined.sct, verbose = FALSE)
immune.combined.sct <- RunUMAP(immune.combined.sct, reduction = "pca", dims = 1:30)
p1 <- DimPlot(immune.combined.sct, reduction = "umap", group.by = "stim")
p2 <- DimPlot(immune.combined.sct, reduction = "umap", group.by = "seurat_annotations", label = TRUE,
repel = TRUE)
p1 + p2
現(xiàn)在掠归，數(shù)據(jù)集已經(jīng)整合缅叠，您可以按照此教程中的先前步驟識別細(xì)胞類型和特定于細(xì)胞類型的響應(yīng)。
