單細胞轉錄組整合分析——seurat包
- 2020 年 3 月 30 日
- 筆記
Seurat是一個分析轉錄組數據的R包，我們之前的推文對其進行過描述：
該包於去年新推出了整合功能。文章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")
可以看到這幾個基因在水平表達量的高低。
未完待續…
