單細胞轉錄組整合分析——seurat包

Seurat是一個分析轉錄組數據的R包，我們之前的推文對其進行過描述：

Seurat 學習筆記

該包於去年新推出了整合功能。文章19年6月份發表於cell雜誌，原文題目為：Comprehensive Integration of Single-Cell Data    被引量超過300次

我們一起來看一下。

該方法的目的是識別不同數據集中存在的共享細胞狀態，即使它們是從不同的個體、實驗條件、技術甚至物種中收集來的。

重點是找到不同數據集中的錨點anchors，這些「錨點」然後用於協調數據集，或將信息從一個數據集傳輸到另一個數據集。

步驟如下：

數據預處理

作者把單細胞數據放在了SeuratData等一系列包中，如果你的網速不行，可以直接到網頁下載數據。

 library(Seurat)
 #devtools::install_github('satijalab/seurat-data')
 library(SeuratData)
 #InstallData("panc8")
 #data("panc8")
 load('panc8.SeuratData/data/panc8.rda')
 
 #To construct a reference, we will identify 『anchors』 between the individual datasets.
 #首先，將組合的數據分成列表，每個數據集是單獨的元素
 pancreas.list <- SplitObject(panc8, split.by = "tech")
 pancreas.list <- pancreas.list[c("celseq", "celseq2", "fluidigmc1", "smartseq2")]
對數據先進行標準化，並識別variable feature。
 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個胰島細胞數據集
整合三個數據集作為參考，並使用FindIntegrationAnchors函數識別錨點。參數默認。
 reference.list <- pancreas.list[c("celseq", "celseq2", "smartseq2")]
 pancreas.anchors <- FindIntegrationAnchors(object.list = reference.list, dims = 1:30)
然後我們將這些錨點傳遞給IntegrateData函數，該函數返回一個Seurat對象。
 pancreas.integrated <- IntegrateData(anchorset = pancreas.anchors, dims = 1:30)
現在我們得到了seurat對象——一個整合後的表達矩陣pancreas.integrated。
然後我們可以使用這個新的表達矩陣進行下遊分析和可視化。
包括進行標準化，運行PCA，並使用UMAP可視化結果。
 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
 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)
 
 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)
左圖按照技術聚類，右圖按照細胞類型聚類。

使用參考數據集進行細胞類型分類
找到錨點之後，我們使用TransferData函數基於參考數據尋找細胞。
 pancreas.query <- pancreas.list[["fluidigmc1"]]
 pancreas.anchors <- FindTransferAnchors(reference = pancreas.integrated, query = pancreas.query,
     dims = 1:30)
 predictions <- TransferData(anchorset = pancreas.anchors, refdata = pancreas.integrated$celltype,
     dims = 1:30)
 pancreas.query <- AddMetaData(pancreas.query, metadata = predictions)
因為我們有來自完整整合分析的原始標籤注釋，所以我們可以評估我們預測的細胞類型注釋與完整參考的匹配程度。在這個例子中，我們發現在細胞類型分類上有很高的一致性，超過97%的細胞被正確標記。
 pancreas.query$prediction.match <- pancreas.query$predicted.id == pancreas.query$celltype
 table(pancreas.query$prediction.match)
 table(pancreas.query$predicted.id)
 VlnPlot(pancreas.query, c("REG1A", "PPY", "SST", "GHRL", "VWF", "SOX10"), group.by = "predicted.id")
可以看到這幾個基因在水平表達量的高低。
未完待續...
文末友情宣傳
強烈建議你推薦給身邊的博士後以及年輕生物學PI，多一點數據認知，讓他們的科研上一個臺階：