前情回顧：

• sc-RAN-seq 數(shù)據(jù)分析||Seurat新版教程:Guided Clustering Tutorial
• sc-RAN-seq 數(shù)據(jù)分析||Seurat新版教程: Integrating datasets to learn cell-type specific responses
• sc-RAN-seq 數(shù)據(jù)分析||Seurat新版教程: Using sctransform in Seurat

如Stuart， Butler等Comprehensive Integration of Single-Cell Data所述疮蹦。Seurat v3引入了集成多個單細(xì)胞數(shù)據(jù)集的新方法饰序。這些方法的目的是識別存在于不同數(shù)據(jù)集中的共享細(xì)胞狀態(tài)(shared cell states)俺驶，即使它們是從不同的個體、實(shí)驗(yàn)條件睁搭、技術(shù)甚至物種中收集來的站玄。

我們的方法旨在首先識別數(shù)據(jù)集對之間的“錨(anchors)”。這些代表了個體細(xì)胞之間成對的對應(yīng)關(guān)系(每個數(shù)據(jù)集中有一個)快骗，我們假設(shè)它們源自相同的生物狀態(tài)娜庇。然后塔次，這些“錨”用于協(xié)調(diào)數(shù)據(jù)集，或?qū)⑿畔囊粋€數(shù)據(jù)集傳輸?shù)搅硪粋€數(shù)據(jù)集名秀。下面励负，我們將演示集成分析的多種應(yīng)用，并介紹2019年文章的中描述的之外的新功能（cell匕得，Comprehensive Integration of Single-Cell Data）继榆。為了幫助指引用戶巾表，我們簡要介紹以下功能:

• 標(biāo)準(zhǔn)工作流程

描述標(biāo)準(zhǔn)Seurat v3集成工作流，并將其應(yīng)用于集成(跨不同技術(shù))收集的多個人類胰島數(shù)據(jù)集略吨。我們還演示了如何使用Seurat v3作為分類器集币，將集群標(biāo)簽傳輸?shù)叫率占臄?shù)據(jù)集中。我們向新用戶推薦這個翠忠。

• SCTransform

描述v3集成工作流的一個修改鞠苟，以便應(yīng)用于使用我們的新規(guī)范化方法SCTransform進(jìn)行規(guī)范化的數(shù)據(jù)集。我們將此方法應(yīng)用于與前面描述的相同的胰島數(shù)據(jù)集秽之，并集成來自8種不同技術(shù)的人類PBMC數(shù)據(jù)集 eight different technologies当娱，作為人類細(xì)胞圖譜的系統(tǒng)技術(shù)基準(zhǔn)。

我們向熟悉SCTransform歸一化方法的高級用戶推薦這個方法考榨。您可以在我們最近的預(yù)印本（preprint）中閱讀有關(guān)SCTransform的更多信息跨细，并了解如何將其應(yīng)用于單獨(dú)的方法中的單個數(shù)據(jù)集（ vignette）。

• Reference-based （基于參照）

描述v3集成工作流的修改河质，其中將數(shù)據(jù)集的子集(或單個數(shù)據(jù)集)作為“參考”（reference）列出冀惭。這種方法可以顯著提高速度，特別是在需要集成大量數(shù)據(jù)集的情況下掀鹅。我們將此方法應(yīng)用于上面描述的8個PBMC數(shù)據(jù)集云头，并觀察到相同的結(jié)果，盡管處理時間大大減少淫半。
對于正在集成許多數(shù)據(jù)集并希望提高速度的用戶溃槐，我們推薦使用這個方法。

• Reciprocal PCA
  描述v3集成工作流的一個修改科吭，其中使用互反PCA（reciprocal PCA）代替正則相關(guān)分析來減少錨點(diǎn)查找中使用的維度昏滴。這種方法可以提高處理大型數(shù)據(jù)集的速度和效率。
  我們建議用戶在處理大量數(shù)據(jù)集或細(xì)胞時对人，如使用許多實(shí)驗(yàn)條件谣殊、復(fù)制或患者的實(shí)驗(yàn)設(shè)計，尋找速度/內(nèi)存方面的改進(jìn)牺弄。然而姻几，這個工作流可能很難對齊高度不同的樣本(例如跨物種，或跨模式势告，集成)蛇捌。對于“turbo”模式，可以考慮結(jié)合“基于參照”的方式整合咱台，如下所示络拌。

標(biāo)準(zhǔn)工作流程

在這個例子工作流中，我們演示了我們最近在論文中引入的兩種新方法回溺，Comprehensive Integration of Single Cell Data:

• 將多個不同的scna -seq數(shù)據(jù)集組裝到一個參考數(shù)據(jù)集中
• 將細(xì)胞類型標(biāo)簽從參考數(shù)據(jù)集轉(zhuǎn)移到新的查詢數(shù)據(jù)集

在本例中春贸，我們選擇了通過四種技術(shù)(CelSeq (GSE81076) CelSeq2 (GSE85241)混萝、Fluidigm C1 (GSE86469)和SMART-Seq2 (E-MTAB-5061)生成的人類胰島細(xì)胞數(shù)據(jù)集。為了方便起見萍恕，我們通過SeuratData包分發(fā)這個數(shù)據(jù)集逸嘀。

新方法的代碼在Seurat v3中實(shí)現(xiàn)。您可以使用install.packages從CRAN下載和安裝允粤。

除了新的方法之外厘熟，Seurat v3還包含了許多旨在改進(jìn)Seurat對象和用戶交互的改進(jìn)。為了幫助用戶熟悉這些更改维哈，我們?yōu)槌Ｒ娙蝿?wù)編寫了一個命令備忘單（ command cheat sheet）绳姨。

載入包和數(shù)據(jù)：

library(Seurat)
packageVersion("Seurat")

>[1] ‘3.1.0’

library(SeuratData)
#InstallData("pbmc3k")
#InstallData("panc8")
#install.packages("C:/Users/Administrator/Desktop/panc8.SeuratData_3.0.2.tar.gz", repos = NULL, type = "source")
#install.packages("C:/Users/Administrator/Desktop/pbmcsca.SeuratData_3.0.0.tar.gz", repos = NULL, type = "source")

> library(panc8.SeuratData)
> data("panc8")
> panc8
An object of class Seurat 
34363 features across 14890 samples within 1 assay 
Active assay: RNA (34363 features)
> summary(panc8)
Length  Class   Mode 
     1 Seurat     S4 
> citation('panc8.SeuratData')

在出版物中使用程序包時引用‘panc8.SeuratData’:

  Satija Lab (2019). panc8.SeuratData: Eight Pancreas
  Datasets Across Five Technologies. R package version
  3.0.2.

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {panc8.SeuratData: Eight Pancreas Datasets Across Five Technologies},
    author = {Satija Lab},
    year = {2019},
    note = {R package version 3.0.2},
  }

Warning message:
In citation("panc8.SeuratData") :
  程序包‘panc8.SeuratData’里的DESCRIPTION文件中沒有日期這一域

要構(gòu)造參考數(shù)據(jù)集，我們將在各個數(shù)據(jù)集之間標(biāo)識“錨”阔挠。首先飘庄，我們將組合的對象拆分為一個列表，每個數(shù)據(jù)集作為一個元素购撼。

> pancreas.list <- SplitObject(panc8, split.by = "tech")
> summary(pancreas.list)
           Length Class  Mode
celseq     1      Seurat S4  
celseq2    1      Seurat S4  
smartseq2  1      Seurat S4  
fluidigmc1 1      Seurat S4  
indrop     1      Seurat S4  

在找到錨之前跪削，我們執(zhí)行標(biāo)準(zhǔn)的預(yù)處理(log-normalization)，并為每個錨分別識別變量特性迂求。注意碾盐，Seurat v3實(shí)現(xiàn)了一種改進(jìn)的基于方差穩(wěn)定轉(zhuǎn)換(“vst”)的變量特征選擇方法。

pancreas.list <- pancreas.list[c("celseq", "celseq2", "fluidigmc1", "smartseq2")]
for (i in 1:length(pancreas.list)) {
  pancreas.list[[i]] <- NormalizeData(pancreas.list[[i]], verbose = FALSE)
  pancreas.list[[i]] <- FindVariableFeatures(pancreas.list[[i]], selection.method = "vst",  nfeatures = 2000, verbose = FALSE)
}

整合3個胰島細(xì)胞數(shù)據(jù)集

接下來揩局，我們使用findintegrationanchor函數(shù)來標(biāo)識錨毫玖，該函數(shù)接受Seurat對象的列表（list）作為輸入。在這里凌盯，我們將三個對象集成到一個參考數(shù)據(jù)集中(稍后我們將使用第四個對象)付枫。

> reference.list <- pancreas.list[c("celseq", "celseq2", "smartseq2")]
> pancreas.anchors <- FindIntegrationAnchors(object.list = reference.list, dims = 1:30)
Computing 2000 integration features
Scaling features for provided objects
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed = 02s
Finding all pairwise anchors
  |                                                  | 0 % ~calculating  Running CCA
Merging objects
Finding neighborhoods
Finding anchors
    Found 3514 anchors
Filtering anchors
    Retained 2761 anchors
Extracting within-dataset neighbors
  |+++++++++++++++++                                 | 33% ~25s          Running CCA
Merging objects
Finding neighborhoods
Finding anchors
    Found 3500 anchors
Filtering anchors
    Retained 2728 anchors
Extracting within-dataset neighbors
  |++++++++++++++++++++++++++++++++++                | 67% ~12s          Running CCA
Merging objects
Finding neighborhoods
Finding anchors
    Found 6174 anchors
Filtering anchors
    Retained 4561 anchors
Extracting within-dataset neighbors
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed = 49s

我們在這里使用所有默認(rèn)參數(shù)來標(biāo)識錨，包括數(shù)據(jù)集的“維數(shù)”(30;您可以嘗試在較大范圍內(nèi)更改此參數(shù)驰怎，例如在10到50之間)阐滩。

然后我們將這些錨傳遞給IntegrateData函數(shù)，該函數(shù)返回一個Seurat對象县忌。
返回的對象將包含一個新的Assay掂榔，其中包含一個針對所有細(xì)胞的完整(或“批量校正”)表達(dá)矩陣，使它們能夠被聯(lián)合分析症杏。

> pancreas.integrated <- IntegrateData(anchorset = pancreas.anchors, dims = 1:30)
Merging dataset 1 into 2
Extracting anchors for merged samples
Finding integration vectors
Finding integration vector weights
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Integrating data
Merging dataset 3 into 2 1
Extracting anchors for merged samples
Finding integration vectors
Finding integration vector weights
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Integrating data

> head(pancreas.integrated@meta.data)
        orig.ident nCount_RNA nFeature_RNA   tech replicate
D101_5        D101   4615.810         1986 celseq    celseq
D101_7        D101  29001.563         4209 celseq    celseq
D101_10       D101   6707.857         2408 celseq    celseq
D101_13       D101   8797.224         2964 celseq    celseq
D101_14       D101   5032.558         2264 celseq    celseq
D101_17       D101  13474.866         3982 celseq    celseq
        assigned_cluster celltype dataset
D101_5              <NA>    gamma  celseq
D101_7              <NA>   acinar  celseq
D101_10             <NA>    alpha  celseq
D101_13             <NA>    delta  celseq
D101_14             <NA>     beta  celseq
D101_17             <NA>   ductal  celseq

運(yùn)行IntegrateData后装获，Seurat對象將包含一個使用集成表達(dá)矩陣的Assay 。注意鸳慈，原始值(未更正的值)仍然存儲在“RNA”測試中的對象中饱溢，因此您可以來回切換喧伞。
然后我們可以使用這個新的矩陣進(jìn)行下游分析和可視化走芋。在這里绩郎，我們縮放集成的數(shù)據(jù)，運(yùn)行PCA翁逞，并使用UMAP可視化結(jié)果肋杖。集成的數(shù)據(jù)集按細(xì)胞類型聚類，而不是按技術(shù)分群挖函。

library(ggplot2)
library(cowplot)
# switch to integrated assay. The variable features of this assay are automatically
# set during IntegrateData
DefaultAssay(pancreas.integrated) <- "integrated"

# Run the standard workflow for visualization and clustering
?ScaleData
pancreas.integrated <- ScaleData(pancreas.integrated, verbose = FALSE)

pancreas.integrated <- RunPCA(pancreas.integrated, npcs = 30, verbose = FALSE)
pancreas.integrated <- RunUMAP(pancreas.integrated, reduction = "pca", dims = 1:30)

Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
20:55:25 Read 5683 rows and found 30 numeric columns
20:55:25 Using Annoy for neighbor search, n_neighbors = 30
20:55:26 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:55:28 Writing NN index file to temp file C:\Users\ADMINI~1\AppData\Local\Temp\Rtmpo5X81p\file2870430b3367
20:55:28 Searching Annoy index using 1 thread, search_k = 3000
20:55:30 Annoy recall = 100%
20:55:31 Commencing smooth kNN distance calibration using 1 thread
20:55:31 Found 2 connected components, falling back to 'spca' initialization with init_sdev = 1
20:55:31 Initializing from PCA
20:55:31 PCA: 2 components explained 44.22% variance
20:55:32 Commencing optimization for 500 epochs, with 252460 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:55:59 Optimization finished

p1 <- DimPlot(pancreas.integrated, reduction = "umap", group.by = "tech")
p2 <- DimPlot(pancreas.integrated, reduction = "umap", group.by = "celltype", label = TRUE, 
              repel = TRUE) + NoLegend()
plot_grid(p1, p2)

使用參考數(shù)據(jù)集對細(xì)胞類型進(jìn)行分類

Seurat v3還支持將參考數(shù)據(jù)(或元數(shù)據(jù))投影到查詢對象上状植。雖然許多方法是一致的(這兩個過程都是從識別錨開始的)，但數(shù)據(jù)映射（data transfer）和整合之間有兩個重要的區(qū)別:

• 在數(shù)據(jù)映射中怨喘，Seurat不糾正或修改查詢表達(dá)式數(shù)據(jù)津畸。
• 在數(shù)據(jù)映射中，Seurat有一個選項(默認(rèn)設(shè)置)來將參考的PCA結(jié)構(gòu)投射到查詢上必怜，而不是學(xué)習(xí)與CCA的聯(lián)合結(jié)構(gòu)肉拓。我們通常建議在scna -seq數(shù)據(jù)集之間投影數(shù)據(jù)時使用此選項。

找到錨之后梳庆，我們使用TransferData函數(shù)根據(jù)參考數(shù)據(jù)(參考數(shù)據(jù)的細(xì)胞類型標(biāo)簽向量)對查詢數(shù)據(jù)集的細(xì)胞進(jìn)行分類暖途。TransferData返回一個帶有預(yù)測id和預(yù)測分?jǐn)?shù)的矩陣，我們可以將其添加到查詢元數(shù)據(jù)中（query metadata）

> pancreas.query <- pancreas.list[["fluidigmc1"]]
> pancreas.anchors <- FindTransferAnchors(reference = pancreas.integrated, query = pancreas.query, 
+                                         dims = 1:30)
Performing PCA on the provided reference using 2000 features as input.
Projecting PCA
Finding neighborhoods
Finding anchors
    Found 953 anchors
Filtering anchors
    Retained 885 anchors
Extracting within-dataset neighbors
> predictions <- TransferData(anchorset = pancreas.anchors, refdata = pancreas.integrated$celltype, 
+                             dims = 1:30)
Finding integration vectors
Finding integration vector weights
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Predicting cell labels
> pancreas.query <- AddMetaData(pancreas.query, metadata = predictions)

因?yàn)槲覀冇衼碜酝暾煞治龅脑紭?biāo)簽注釋膏执，所以我們可以評估預(yù)測的細(xì)胞類型注釋與完整引用的匹配程度驻售。在這個例子中，我們發(fā)現(xiàn)在細(xì)胞類型分類上有很高的一致性更米，超過97%的細(xì)胞被正確標(biāo)記欺栗。

pancreas.query$prediction.match <- pancreas.query$predicted.id == pancreas.query$celltype
table(pancreas.query$prediction.match)

FALSE  TRUE 
   18   620 

為了進(jìn)一步驗(yàn)證這一點(diǎn)，我們可以檢查特定胰島細(xì)胞群的一些典型細(xì)胞類型標(biāo)記（cell type markers）征峦。請注意纸巷，即使這些細(xì)胞類型中的一些只由一兩個細(xì)胞(例如epsilon細(xì)胞)，我們?nèi)匀荒軌蛘_地對它們進(jìn)行分類眶痰。

> table(pancreas.query$predicted.id)

            acinar activated_stellate              alpha               beta              delta 
                21                 17                248                258                 22 
            ductal        endothelial            epsilon              gamma         macrophage 
                33                 13                  1                 17                  1 
              mast            schwann 
                 2                  5 
> VlnPlot(pancreas.query, c("REG1A", "PPY", "SST", "GHRL", "VWF", "SOX10"), group.by = "predicted.id")

SCTransform

在前面的選項卡中瘤旨，我們演示了如何在使用標(biāo)準(zhǔn)規(guī)范化（log-normalization）對每個數(shù)據(jù)集進(jìn)行預(yù)處理之后集成數(shù)據(jù)集。在這里竖伯，我們修改工作流存哲，以利用改進(jìn)的預(yù)處理和規(guī)范化工作流:SCTransform。您可以在我們最近的預(yù)印本中閱讀有關(guān)SCTransform的更多信息七婴，并了解如何將其應(yīng)用于單個數(shù)據(jù)集祟偷。

• 創(chuàng)建要集成的Seurat對象列表
• 對每個數(shù)據(jù)集分別執(zhí)行SCTransform規(guī)范化
• 在對象列表上運(yùn)行PrepSCTIntegration函數(shù)
• 集成數(shù)據(jù)集，并進(jìn)行聯(lián)合分析

首先打厘，設(shè)置Seurat對象列表修肠，分別對每個對象運(yùn)行SCTransform:

m(list=ls())
library(Seurat)
library(SeuratData)
library(ggplot2)
options(future.globals.maxSize = 4000 * 1024^2)
data("panc8")
pancreas.list <- SplitObject(panc8, split.by = "tech")
pancreas.list <- pancreas.list[c("celseq", "celseq2", "fluidigmc1", "smartseq2")]
for (i in 1:length(pancreas.list)) {
  pancreas.list[[i]] <- SCTransform(pancreas.list[[i]], verbose = FALSE)
}

接下來，選擇下游集成的特性户盯，并運(yùn)行PrepSCTIntegration嵌施，這將確保計算了所有必需的Pearson殘差饲化。

?SelectIntegrationFeatures
pancreas.features <- SelectIntegrationFeatures(object.list = pancreas.list, nfeatures = 3000)
pancreas.list <- PrepSCTIntegration(object.list = pancreas.list, anchor.features = pancreas.features,  verbose = FALSE)

接下來，識別錨并集成數(shù)據(jù)集吗伤。命令與標(biāo)準(zhǔn)工作流相同吃靠，但請確保設(shè)置了normalization.method = 'SCT':

pancreas.anchors <- FindIntegrationAnchors(object.list = pancreas.list, normalization.method = "SCT", 
                                           anchor.features = pancreas.features, verbose = FALSE)
pancreas.integrated <- IntegrateData(anchorset = pancreas.anchors, normalization.method = "SCT", 
                                     verbose = FALSE)

現(xiàn)在繼續(xù)對集成數(shù)據(jù)集進(jìn)行下游分析(即可視化、聚類)足淆。命令與標(biāo)準(zhǔn)工作流相同巢块，但在集成之后不運(yùn)行ScaleData函數(shù)。您可以看到巧号，在集成之后族奢，細(xì)胞按照它們的生物細(xì)胞類型(已經(jīng)預(yù)先注釋過)分組，而不是按照它們的技術(shù)分組丹鸿。

pancreas.integrated <- RunPCA(pancreas.integrated, verbose = FALSE)
pancreas.integrated <- RunUMAP(pancreas.integrated, dims = 1:30)


plots <- DimPlot(pancreas.integrated, group.by = c("tech", "celltype"), combine = FALSE)
plots <- lapply(X = plots, FUN = function(x) x + theme(legend.position = "top") + guides(color = guide_legend(nrow = 3, 
                                                                                                              byrow = TRUE, override.aes = list(size = 3))))
CombinePlots(plots)

為了進(jìn)一步演示這個工作流歹鱼，我們接下來將它應(yīng)用于來自8種不同技術(shù)的一系列人類PBMC數(shù)據(jù)集，這些數(shù)據(jù)集由human Cell Atlas作為系統(tǒng)的技術(shù)基準(zhǔn)生成卜高。數(shù)據(jù)是單細(xì)胞門戶下載的弥姻。為了方便起見，我們通過SeuratData包分發(fā)了這個數(shù)據(jù)集掺涛。

data("pbmcsca")
pbmc.list <- SplitObject(pbmcsca, split.by = "Method")
for (i in names(pbmc.list)) {
    pbmc.list[[i]] <- SCTransform(pbmc.list[[i]], verbose = FALSE)
}
pbmc.features <- SelectIntegrationFeatures(object.list = pbmc.list, nfeatures = 3000)
pbmc.list <- PrepSCTIntegration(object.list = pbmc.list, anchor.features = pbmc.features)
pbmc.anchors <- FindIntegrationAnchors(object.list = pbmc.list, normalization.method = "SCT", 
    anchor.features = pbmc.features)
pbmc.integrated <- IntegrateData(anchorset = pbmc.anchors, normalization.method = "SCT")

pbmc.integrated <- RunPCA(object = pbmc.integrated, verbose = FALSE)
pbmc.integrated <- RunUMAP(object = pbmc.integrated, dims = 1:30)

plots <- DimPlot(pbmc.integrated, group.by = c("Method", "CellType"), combine = FALSE)
plots <- lapply(X = plots, FUN = function(x) x + theme(legend.position = "top") + guides(color = guide_legend(nrow = 4, 
    byrow = TRUE, override.aes = list(size = 2.5))))
CombinePlots(plots)

請再次注意庭敦，細(xì)胞是按照共享的生物亞型分組的，而不是按照技術(shù)分組的薪缆。在這種情況下秧廉，作者分別在每個數(shù)據(jù)集中標(biāo)注“CellType”，你可以看到這與我們的綜合分析是一致的:

DimPlot(pbmc.integrated, group.by = "CellType", split.by = "Method", ncol = 3)

然而拣帽，我們的綜合分析揭示了這些廣泛的單元類型之間更精細(xì)的細(xì)分疼电，因?yàn)榻M合的數(shù)據(jù)集具有更大的統(tǒng)計能力。我們沒有在集成的數(shù)據(jù)集上進(jìn)行完整的聚類减拭，但是請注意蔽豺，我們可以觀察CD4+ (Memory/Naive)、CD8(多個細(xì)胞毒性群體)和B細(xì)胞(多個發(fā)育階段)的亞群拧粪。下面修陡，我們在集成可視化上可視化原始度量的表達(dá)式(注意，您可以添加split.by = 'Method‘ 為每種技術(shù)單獨(dú)繪制這些圖)"

DefaultAssay(pbmc.integrated) <- "RNA"
# Normalize RNA data for visualization purposes
pbmc.integrated <- NormalizeData(pbmc.integrated, verbose = FALSE)
FeaturePlot(pbmc.integrated, c("CCR7", "S100A4", "GZMB", "GZMK", "GZMH", "TCL1A"))

Reference-based

在這里可霎，我們對Seurat集成工作流進(jìn)行了額外的修改魄鸦，我們將其稱為“基于參照”集成。在前面的工作流中癣朗，我們在所有對數(shù)據(jù)集之間標(biāo)識錨拾因。雖然這使得數(shù)據(jù)集在下游集成中具有同等的權(quán)重，但它也可能成為計算密集型的。例如绢记，當(dāng)集成10個不同的數(shù)據(jù)集時扁达，我們執(zhí)行45種不同的成對比較。

作為一種替代方法庭惜，我們在這里介紹了指定一個或多個數(shù)據(jù)集作為集成分析的“參照”的可能性罩驻，其余的指定為“查詢”數(shù)據(jù)集穗酥。在此工作流中护赊，我們不確定查詢數(shù)據(jù)集對之間的錨，從而減少了比較的數(shù)量砾跃。例如骏啰，當(dāng)將10個數(shù)據(jù)集與一個指定為參照的數(shù)據(jù)集集成時，我們只執(zhí)行9個比較抽高∨懈基于參照的集成可以應(yīng)用于log-normalized 或sctransform規(guī)范化數(shù)據(jù)集。

在實(shí)踐中翘骂，我們經(jīng)常在兩種方法之間觀察到非常相似的結(jié)果壁熄，但是計算時間大大減少。我們使用與前一個選項卡相同的PBMC分析來演示這一點(diǎn)碳竟，輸出高度一致草丧。因此，我們向集成大量數(shù)據(jù)集并希望提高計算效率的用戶推薦這種工作流莹桅。

這種方法要求用戶預(yù)先選擇數(shù)據(jù)集作為參照昌执。在本例中，我們使用了10X v3數(shù)據(jù)集诈泼，因?yàn)閿?shù)據(jù)集同時包含高細(xì)胞數(shù)和高靈敏度懂拾。用戶可以在選擇引用時應(yīng)用類似的標(biāo)準(zhǔn)。此外铐达，您可以在下面的工作流中更改此設(shè)置(例如岖赋，嘗試選擇10X v2或Drop-seq數(shù)據(jù)集作為參考)，在最終結(jié)果中只有非常小的差異瓮孙。

工作流程如下所示贾节。命令與前一個例子相同，只是我們在findintegrationanchor函數(shù)中指定了一個參照參數(shù)衷畦。數(shù)據(jù)最初是從 the Broad Single Cell Portal下載的栗涂，但是為了方便起見，我們通過SeuratData包分發(fā)該數(shù)據(jù)集祈争。

library(Seurat)
library(SeuratData)
library(ggplot2)
InstallData("pbmcsca")
data("pbmcsca")
pbmc.list <- SplitObject(pbmcsca, split.by = "Method")
for (i in names(pbmc.list)) {
    pbmc.list[[i]] <- SCTransform(pbmc.list[[i]], verbose = FALSE)
}
pbmc.features <- SelectIntegrationFeatures(object.list = pbmc.list, nfeatures = 3000)
pbmc.list <- PrepSCTIntegration(object.list = pbmc.list, anchor.features = pbmc.features)

# This command returns dataset 5.  We can also specify multiple refs. (i.e. c(5,6))
reference_dataset <- which(names(pbmc.list) == "10x Chromium (v3)")

pbmc.anchors <- FindIntegrationAnchors(object.list = pbmc.list, normalization.method = "SCT", 
    anchor.features = pbmc.features, reference = reference_dataset)
pbmc.integrated <- IntegrateData(anchorset = pbmc.anchors, normalization.method = "SCT")

pbmc.integrated <- RunPCA(object = pbmc.integrated, verbose = FALSE)
pbmc.integrated <- RunUMAP(object = pbmc.integrated, dims = 1:30)

plots <- DimPlot(pbmc.integrated, group.by = c("Method", "CellType"), combine = FALSE)
plots <- lapply(X = plots, FUN = function(x) x + theme(legend.position = "top") + guides(color = guide_legend(nrow = 4, 
    byrow = TRUE, override.aes = list(size = 2.5))))
CombinePlots(plots)

DefaultAssay(pbmc.integrated) <- "RNA"
# Normalize RNA data for visualization purposes
pbmc.integrated <- NormalizeData(pbmc.integrated, verbose = FALSE)
FeaturePlot(pbmc.integrated, c("CCR7", "S100A4", "GZMB", "GZMK", "GZMH", "TCL1A"))

Reciprocal PCA

對于非常大的數(shù)據(jù)集斤程，CCA（Canonical Correlation Analysis ）有時可能成為計算上的卡點(diǎn)步驟。在這個工作流中，我們引入了另一種降維形式——reciprocal PCA忿墅，以確定尋找錨的有效空間扁藕。在使用reciprocal PCA確定任意兩個數(shù)據(jù)集之間的錨點(diǎn)時，我們將每個數(shù)據(jù)集投影到另一個數(shù)據(jù)集的主元分析空間中疚脐，并根據(jù)相同的互鄰域要求約束錨點(diǎn)亿柑。所有下游集成步驟保持不變，我們能夠“糾正”(或協(xié)調(diào))數(shù)據(jù)集棍弄。

如下所示望薄，工作流程包括以下步驟:

• 創(chuàng)建要集成的Seurat對象列表
• 對每個數(shù)據(jù)集分別執(zhí)行規(guī)范化、特征選擇和縮放
• 對列表中的每個對象運(yùn)行PCA
• 集成數(shù)據(jù)集呼畸，并進(jìn)行聯(lián)合分析

總的來說痕支，我們在這兩個工作流中觀察到驚人的相似結(jié)果，但是當(dāng)使用reciprocal PCA時蛮原，計算時間和內(nèi)存顯著減少卧须。然而，如果數(shù)據(jù)集高度分散(例如儒陨，跨模態(tài)映射或跨物種映射)花嘶，則只能使用一小部分特性來促進(jìn)集成，并且您可以使用CCA觀察到更好的結(jié)果蹦漠。

對于包含許多數(shù)據(jù)集的大型研究椭员，我們還建議將reciprocal PCA與基于引用的集成或SCTransform標(biāo)準(zhǔn)化相結(jié)合(請參閱前面的詳細(xì)信息)。
對于這些例子津辩，我們將使用來自人類細(xì)胞圖譜的“免疫細(xì)胞圖譜”數(shù)據(jù)拆撼，可以在這里找到。首先喘沿，我們將集成這個數(shù)據(jù)集的隨機(jī)下采樣版本闸度，然后是整合的數(shù)據(jù)集。隨機(jī)下采樣的版本包含40K個細(xì)胞(每個供體5K個)蚜印，可以通過我們新的SeuratData包獲得莺禁。

library(Seurat)
library(SeuratData)
# Get data from seurat-data
InstallData(pkgs = "bm40k.SeuratData")

在獲取數(shù)據(jù)后，首先進(jìn)行標(biāo)準(zhǔn)歸一化和變量特征選擇窄赋。

library(ggplot2)
# set up future for parallelization
library(future)
library(future.apply)
plan("multiprocess", workers = 4)
options(future.globals.maxSize = 20000 * 1024^2)
# load in the data
data(hcabm40k)
bm40k.list <- SplitObject(hcabm40k, split.by = "orig.ident")
bm40k.list <- lapply(X = bm40k.list, FUN = function(x) {
    x <- NormalizeData(x, verbose = FALSE)
    x <- FindVariableFeatures(x, verbose = FALSE)
})

接下來哟冬，選擇用于下游集成的特性，并在列表中的每個對象上運(yùn)行PCA忆绰，這是運(yùn)行對等的PCA工作流所必需的浩峡。

features <- SelectIntegrationFeatures(object.list = bm40k.list)
bm40k.list <- lapply(X = bm40k.list, FUN = function(x) {
    x <- ScaleData(x, features = features, verbose = FALSE)
    x <- RunPCA(x, features = features, verbose = FALSE)
})

由于這個數(shù)據(jù)集同時包含男性和女性，我們將選擇一個男性和一個女性(BM1和BM2)在基于引用的工作流中使用错敢。我們通過檢測XIST基因的表達(dá)來確定供體性別翰灾。接下來，識別錨并集成數(shù)據(jù)集。命令與標(biāo)準(zhǔn)工作流相同纸淮，但請確保設(shè)置reduce = 'rpca':

anchors <- FindIntegrationAnchors(object.list = bm40k.list, reference = c(1, 2), reduction = "rpca", 
    dims = 1:30)
bm40k.integrated <- IntegrateData(anchorset = anchors, dims = 1:30)

現(xiàn)在繼續(xù)對集成數(shù)據(jù)集進(jìn)行下游分析(即可視化平斩、聚類)。命令與標(biāo)準(zhǔn)工作流相同咽块。

bm40k.integrated <- ScaleData(bm40k.integrated, verbose = FALSE)
bm40k.integrated <- RunPCA(bm40k.integrated, verbose = FALSE)
bm40k.integrated <- RunUMAP(bm40k.integrated, dims = 1:30)
DimPlot(bm40k.integrated, group.by = "orig.ident")

為了進(jìn)一步演示這個工作流绘面，我們接下來將它應(yīng)用于完整的數(shù)據(jù)集。雖然我們不通過SeuratData提供此功能侈沪，但您可以在HCA數(shù)據(jù)門戶網(wǎng)站上找到counts矩陣揭璃。下面的命令與40k集成相同，只是我們將維度增加到了50峭竣，以反映增加的細(xì)胞數(shù)量和多樣性塘辅。

bm280k.data <- Read10X_h5("../data/ica_bone_marrow_h5.h5")
bm280k <- CreateSeuratObject(counts = bm280k.data, min.cells = 100, min.features = 500)

bm280k.list <- SplitObject(bm280k, split.by = "orig.ident")
bm280k.list <- future_lapply(X = bm280k.list, FUN = function(x) {
    x <- NormalizeData(x, verbose = FALSE)
    x <- FindVariableFeatures(x, verbose = FALSE)
})
features <- SelectIntegrationFeatures(object.list = bm280k.list)
bm280k.list <- future_lapply(X = bm280k.list, FUN = function(x) {
    x <- ScaleData(x, features = features, verbose = FALSE)
    x <- RunPCA(x, features = features, verbose = FALSE)
})

anchors <- FindIntegrationAnchors(object.list = bm280k.list, reference = c(1, 2), reduction = "rpca", 
    dims = 1:50)
bm280k.integrated <- IntegrateData(anchorset = anchors, dims = 1:50)
bm280k.integrated <- ScaleData(bm280k.integrated, verbose = FALSE)
bm280k.integrated <- RunPCA(bm280k.integrated, verbose = FALSE)
bm280k.integrated <- RunUMAP(bm280k.integrated, dims = 1:50)
DimPlot(bm280k.integrated, group.by = "orig.ident")