Introduction to scRNA-seq integration

Introduction to scRNA-seq integration

두 개 이상의 싱글셀 데이터세트를 같이 분석하는 것은 unique 한 challenges 가 있습니다. 특히 standard workflows 에서는 여러 데이터세트 걸쳐 존재하는 세포 population 을 식별하는것이 문제가 될 수 있습니다. Seurat v4 에는 여러 데이터세트에서 공유된 세포 population 을 match 시키는 (또는 align하는) 여러 방법이 들어있습니다. 이러한 방법은 먼저 biological state 가 match 되는 세포들의 cross-dataset pairs 를 찾는데 (anchors), 이는 데이터 세트 간의 기술적 차이 (i.e. batch effect correction)를 고치고 비교 scRNA-seq 분석을 하는데에 사용하실 수 있습니다.

아래에서는 Stuart*, Butler* et al, 2019에 설명된 대로 scRNA-seq integration 을 위한 방법을 보여주고, resting 하고 있는 면역세포 (PBMC) 와 interferon-stimulated state 에서의 세포들을 비교하는 분석을 합니다.

Integration goals

이 튜토리얼은 Seurat 의 integration procedure 를 사용해서 할 수 있는 complex cell types 에서의 comparative analysis 의 overview 를 제공하기 위해 만들어졌습니다. 중요한 목표는 아래와 같습니다.

  • Downstream analysis 에 쓰일 수 있는 integrated 한 data assay 만들기
  • 두 데이터세트에 있는 세포 type 식별하기
  • control 과 stimulated 세포에서 다 나오는 세포 type markers 찾기
  • 데이터 세트를 비교해서 stimulation 에 대한 cell-type specific response 를 찾기

Setup the Seurat objects

쉽게 튜토리얼을 하실 수 있도록 데이터 세트는 SeuratData package 로 공유드립니다.

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)

Perform integration

다음으론 FindIntegrationAnchors() 를 사용해서 anchors 를 식별합니다. 이 함수는 Seurat 객체 리스트를 input 으로 사용하고, 이러한 anchor 를 사용하여 두 데이터 세트를 integrateData() 로 통합합니다.

immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, anchor.features = features)
# this command creates an 'integrated' data assay
immune.combined <- IntegrateData(anchorset = immune.anchors)

Perform an integrated analysis

이제는 모든 세포에서 하나의 integrated analysis 를 실행할 수 있습니다!

# 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
https://satijalab.org/seurat/articles/integration_introduction_files/figure-html/viz-1.png

두 conditions 를 옆에서 보려면 split.by argument 를 사용하실 수 있습니다.

DimPlot(immune.combined, reduction = "umap", split.by = "stim"
https://satijalab.org/seurat/articles/integration_introduction_files/figure-html/split.dim-1.png

Identify conserved cell type markers

여러 condition 에서 보존되는 표준 cell-type marker 유전자를 알기위해서는 FindConservedMarkers() 함수가 쓰입니다. 이 함수는 각 데이터세트/그룹에 differential gene expression 테스트를 하고 MetaDeR package 의 meta-analysis method 을 사용해서 p-values 를 조합합니다. 예를 들어, cluster 6 (NK 세포) 에선 stimulation condition 에 상관 없이 보존된 markers 유전자들을 계산할 수 있습니다.

# 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

각 cluster 에 대해 이러한 marker 유전자를 찾아서 세포 type 으로 cluster 에 annotate 할 수 있습니다

FeaturePlot(immune.combined, features = c("CD3D", "SELL", "CREM", "CD8A", "GNLY", "CD79A", "FCGR3A",
    "CCL2", "PPBP"), min.cutoff = "q9")
https://satijalab.org/seurat/articles/integration_introduction_files/figure-html/annotate-1.png
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)
https://satijalab.org/seurat/articles/integration_introduction_files/figure-html/annotate-2.png

split.by parameter 를 사용한 DotPlot() 함수는 expression level 과 cluster 에서의 세포 percentage 를 보여주고 여러 conditions 에서 나오는 cell type markers 를 보는 데 유용합니다. 여기서는 14개의 cluster 각각에 대해 2~3개의 strong marker genes 를 plot 합니다.

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()
https://satijalab.org/seurat/articles/integration_introduction_files/figure-html/splitdotplot-1.png

Identify differential expressed genes across conditions

이제 stimulated 세포와 control 세포를 정렬했으니까 comparative analysis 를 하고 stimulation 에 의해 유발되는 차이를 살펴볼 수 있습니다. 이러한 변화를 폭넓게 보는 한 가지 방법은 stimulated 와 control 세포의 average expression 을 plot 하고 이 scatter plot 에서 보이는 visual outliers 를 찾는 것 입니다. 여기서는 stimulated 와 control 에서 naive T cell 과 CD14 monocyte population 의 average expression 으로 scatter plots 를 만들어서, interferon stimulation 시 많이 변하는 유전자를 볼 수 있습니다.

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
https://satijalab.org/seurat/articles/integration_introduction_files/figure-html/scatterplots-1.png

보시다시피, 두 세포 type 에서 동일한 유전자가 upregulated 되었고 이는 보존된 interferon response pathway 를 나타냅니다. 두 condition 에서 공통적인 cell type 을 확실히 확인했기때문에, 다른 condition 의 같은 세포에서 어떻게 유전자가 다른지 물어볼 수 있습니다. 먼저 meta.data slot 에 열을 만들어 cell type 과 stimulation information 을 저장하고 현재 ident 를 그 열로 전환합니다. 그리고 FindMarkers() 를 이용해서 stimulated 와 control 의 B cell 에서 다르게 나오는 유전자들을 찾습니다. 여기에서 알 수 있는 것은, top gene 은 위에 나온 interferon response 유전자들과 같다는 것 입니다. 추가적으로, CSCL10 와 같은 유전자는 monocyte 에서만 나오고 B cell interferon response 도 이 리스트에서 매우 중요합니다.

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

유전자 발현의 이러한 변화를 visualize 하는 다른 유용한 방법은 FeaturePlot() 또는 VlnPlot() 함수를 사용하는 것 입니다. 그러면 grouping variable (여기서는 stimulation condition) 으로 나눠진 특정 유전자 리스트에 대한 FeaturePlots 이 나옵니다. CD3D 와 GNLY 같은 canonical cell type markers (T cell 그리고 NK/CD8 세포) 는 interferon stimulation 에 영향을 받지 않고 control 이나 stimulated 그룹에서 비슷한 유전자 발현이 나옵니다. 반면 IFI6 와 ISG15 는 interferon response 유전자들이고, 따라서 모든 세포 type 에서 upregulated 되었습니다. 마지막으로 CD14와 CXCL10은 cell type specific 한 interferon response 를 보여주는 유전자입니다. CD14 발현은 CD14 monocytes 의 stimulation 후 줄어들었고, 이는 integrated analysis 의 값을 underscoring 해서 supervised analysis framework 에서 misclassification 이 될 수 있습니다. CXCL10 은 interferon stimulation 후 monocytes 와 B cell 에서 확실히 upregulation 을 보여줍니다.

FeaturePlot(immune.combined, features = c("CD3D", "GNLY", "IFI6"), split.by = "stim", max.cutoff = 3,
    cols = c("grey", "red"))
https://satijalab.org/seurat/articles/integration_introduction_files/figure-html/feature.heatmaps-1.png
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)
https://satijalab.org/seurat/articles/integration_introduction_files/figure-html/splitvln-1.png

Performing integration on datasets normalized with SCTransform

Hafemeister and Satija, 2019 에서는 regulated negative binomial regression 을 기반으로 한 improved 된 scRNA-seq 의 normalization method 을 소개합니다. 이 방법은 ‘sctransform’ 이라고 불리고 pseudocount 와 log-transformation 을 포함해서 standard normalization workflow 의 일부 pitfalls 를 방지합니다. Sctransform 에 대한 더 자세한 내용은 이 manuscript 또는 SCTransform vignette 에서 확인하실 수 있습니다.

아래에서는 sctransform workflow 로 normalize 된 데이터세트에 대해 Seurat integration workflow 를 수정하는 방법을 보여줍니다. 이 commands들은 몇 가지 차이점 외에는 크게 비슷합니다.

  • Integration 하기 전에 NormalizeData() 하는것 보다, SCTransform() 으로 데이터세트를 개별로 normalize 합니다
  • SCTransform vignette 에서 더 설명되는것 처럼, 일반적으론 3,000개의 features 를 사영하여 sctransform 의 downstream analysis 를 수행합니다
  • anchor 를 찾기 전에 PrepSCTIntegration() 함수를 실행합니다
  • FindIntegrationAchords()IntegrateData() 를 실행할땐 normalization.method 의 parameter 값을 SCT 로 설정합니다
  • Integration 을 포함한 sctransform-based workflows 를 실행할때는 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
https://satijalab.org/seurat/articles/integration_introduction_files/figure-html/immunesca.cca.sct.split.dims-1.png

이제 데이터세트가 integrate 되었으니까 이 vignette 의 전 단계에 따라 cell type 과 cell type-specific response 를 식별하실 수 있습니다.

Leave a Reply

Your email address will not be published. Required fields are marked *