單細胞分析Seurat使用相關的10個問題答疑精選！

2020 年 2 月 26 日
筆記

NGS系列

作為一個剛剛開始進行單細胞轉錄組分析的菜鳥，R語言底子沒有，有時候除了會copy外，如果你讓我寫個for循環，我只能cross my fingers。。。。

於是我看見了https://satijalab.org/seurat/，Seurat是一個R軟體包，設計用於QC、分析和探索單細胞RNA-seq數據。Seurat旨在幫助用戶能夠識別和解釋單細胞轉錄組學中的的異質性來源，並通過整合各種類型的單細胞數據，能夠在單個細胞層面上進行系統分析。

裡面非常詳細的介紹了這個單細胞轉錄組測序的workflow,包括添加了很多的其他功能，如細胞周期 (Seurat亮點之細胞周期評分和回歸)等。

但裡面有蠻多的程式碼的原理其實我並不太清楚 (讀完這個，還不懂，來找我，重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程（原理、程式碼和評述）)，這次我就介紹一下裡面讓我曾經困惑的幾個問題以及比較nice的解答！（令我萬萬沒想到，找答案比自己寫答案確實困難的多。。。）

1. 有關merge函數的問題

merge只是放在一起，fastMNN才是真正的整合分析。

2. 有關PC的選擇

Seurat應用JackStraw隨機抽樣構建一個特徵基因與主成分相關性值的背景分布，選擇富集特徵基因相關性顯著的主成分用於後續分析。對大的數據集，這一步計算會比較慢，有時也可能不會找到合適的臨界點。
建議通過ElbowPlot來選，找到拐點或使得所選PC包含足夠大的variation了 (80%以上)，便合適。然後再可以在這個數目上下都選幾個值試試，最好測試的時候往下游測試些，越下游越好，看看對結果是否有影響。
另外一個方式可以是根據與各個主成分相關的基因的GSEA功能富集分析選擇合適的主成分 (一文掌握GSEA，超詳細教程)。
一般選擇7-12都合適。實際分析時，可以嘗試選擇不同的值如10, 15或所有主成分，結果通常差別不大。但如果選擇的主成分比較少如5等，結果可能會有一些變化。
另外要注意：用了這麼多年的PCA可視化竟然是錯的！！！

3. 細胞周期相關基因

Pairs of genes (A, B) are identified from a training data set where in each pair, the fraction of cells in phase G1 with expression of A > B (based on expression values in training.data) and the fraction with B > A in each other phase exceeds fraction. These pairs are defined as the marker pairs for G1. This is repeated for each phase to obtain a separate marker pair set.

4. Seurat對象導入Monocle

其實這個問題我也遇到了，並且已經有人給出了解決方案。(生信寶典註：官方最開始沒支援Seurat3，我寫了個簡易的，在https://github.com/cole-trapnell-lab/monocle-release/issues/311，現在應該更新全面支援Seurat3了。)

seurat_data <- newimport(SeuratObject)    ###重新定義了import函數，稱為newimport  newimport <- function(otherCDS, import_all = FALSE) {    if(class(otherCDS)[1] == 'Seurat') {      requireNamespace("Seurat")      data <- otherCDS@assays$RNA@counts        if(class(data) == "data.frame") {        data <- as(as.matrix(data), "sparseMatrix")      }        pd <- tryCatch( {        pd <- new("AnnotatedDataFrame", data = [email protected])        pd      },      #warning = function(w) { },      error = function(e) {        pData <- data.frame(cell_id = colnames(data), row.names = colnames(data))        pd <- new("AnnotatedDataFrame", data = pData)          message("This Seurat object doesn't provide any meta data");        pd      })        # remove filtered cells from Seurat      if(length(setdiff(colnames(data), rownames(pd))) > 0) {        data <- data[, rownames(pd)]      }        fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))      fd <- new("AnnotatedDataFrame", data = fData)      lowerDetectionLimit <- 0        if(all(data == floor(data))) {        expressionFamily <- negbinomial.size()      } else if(any(data < 0)){        expressionFamily <- uninormal()      } else {        expressionFamily <- tobit()      }        valid_data <- data[, row.names(pd)]        monocle_cds <- newCellDataSet(data,                                    phenoData = pd,                                    featureData = fd,                                    lowerDetectionLimit=lowerDetectionLimit,                                    expressionFamily=expressionFamily)        if(import_all) {        if("Monocle" %in% names(otherCDS@misc)) {          otherCDS@misc$Monocle@auxClusteringData$seurat <- NULL          otherCDS@misc$Monocle@auxClusteringData$scran <- NULL            monocle_cds <- otherCDS@misc$Monocle          mist_list <- otherCDS          } else {          # mist_list <- list(ident = ident)          mist_list <- otherCDS        }      } else {        mist_list <- list()      }        if(1==1) {        var.genes <- setOrderingFilter(monocle_cds, otherCDS@[email protected])        }      monocle_cds@auxClusteringData$seurat <- mist_list      } else if (class(otherCDS)[1] == 'SCESet') {      requireNamespace("scater")        message('Converting the exprs data in log scale back to original scale ...')      data <- 2^otherCDS@assayData$exprs - otherCDS@logExprsOffset        fd <- otherCDS@featureData      pd <- otherCDS@phenoData      experimentData = otherCDS@experimentData      if("is.expr" %in% slotNames(otherCDS))        lowerDetectionLimit <- [email protected]      else        lowerDetectionLimit <- 1        if(all(data == floor(data))) {        expressionFamily <- negbinomial.size()      } else if(any(data < 0)){        expressionFamily <- uninormal()      } else {        expressionFamily <- tobit()      }        if(import_all) {        # mist_list <- list(iotherCDS@sc3,        #                   otherCDS@reducedDimension)        mist_list <- otherCDS        } else {        mist_list <- list()      }        monocle_cds <- newCellDataSet(data,                                    phenoData = pd,                                    featureData = fd,                                    lowerDetectionLimit=lowerDetectionLimit,                                    expressionFamily=expressionFamily)      # monocle_cds@auxClusteringData$sc3 <- otherCDS@sc3      # monocle_cds@auxOrderingData$scran <- mist_list        monocle_cds@auxOrderingData$scran <- mist_list      } else {      stop('the object type you want to export to is not supported yet')    }      return(monocle_cds)  }

5. 不同條件下畫熱圖

6. Seurat對象轉化為.h5ad對象

Answer：Looks like the way to do it is to write to loom format via loomR(https://github.com/mojaveazure/loomR), then read that into anndata (https://anndata.readthedocs.io/en/latest/api.html) to be written as an .h5ad file.

Single cell folks need to a pick a file format/structure and stick with it.

生信寶典註：不同文件類型和格式之間的轉換確實是一個問題，好在易生信提供了好的解決方案。來這試試？易生信線下課推遲，線上課程免費看

7. 根據基因取細胞子集

也可以提取特定Cluster，進行進一步細分。

# 提取特定cluster，繼續後續分析。  ident_df <- data.frame(cell=names(Idents(pbmc)), cluster=Idents(pbmc))  pbmc_subcluster <- subset(pbmc, cells=as.vector(ident_df[ident_df$cluster=="Naive CD4 T",1]))

8. 篩選條件的設置

這是最開始讀入的初始設置，後面質控還有更多篩選方式。

質控是用於保證下游分析時數據品質足夠好。由於不能先驗地確定什麼是足夠好的數據品質，所以只能基於下游分析結果（例如，細胞簇注釋）來對其進行判斷。因此，在分析數據時可能需要反覆調整參數進行品質控制 (生信寶典註：反覆分析，多次分析，是做生信的基本要求。這也是為啥需要上易生信培訓班而不是單純依賴公司的分析。)。從寬鬆的QC閾值開始並研究這些閾值的影響，然後再執行更嚴格的QC，通常是有益的。這種方法對於細胞種類混雜的數據集特別有用，在這些數據集中，特定細胞類型或特定細胞狀態可能會被誤解為低品質的異常細胞。在低品質數據集中，可能需要嚴格的品質控制閾值。具體見：

陷阱和建議：

通過可視化檢測到的基因數量、總分子數和線粒體基因的表達比例的分布中的異常峰來執行QC。聯合考慮這3個變數，而不是單獨考慮一個變數。
儘可能設置寬鬆的QC閾值；如果下游聚類無法解釋時再重新設定嚴格的QC閾值。
如果樣品之間的QC變數分布不同（存在多個強峰），則需要考慮樣品品質差異，應按照Plasschaert et al. (2018)的方法為每個樣品分別確定QC閾值。

9. RunTSNE不是在聚類

區分好聚類 (FindClusters)和降維 (PCA，tSNE，UMAP)。

聚類是直接基於距離矩陣的經典無監督機器學習問題。通過最小化簇內距離或在降維後的表達空間中鑒定密集區域，將細胞分配給簇。方法有層級聚類、K-menas、基於圖的聚類等（基因共表達聚類分析及可視化，基因表達聚類分析之初探SOM）。
tSNE和UMAP目前主要用於可視化。用於可視化的降維必然涉及資訊丟失並改變細胞之間的距離。因此tSNE/UMAP圖應僅只用於解釋或傳達基於更精確的、更多維度的定量分析結果。這樣可以保證分析充分利用了壓縮到二維空間時丟失的資訊。假如二維圖上呈現的細胞分布與使用更多數目的PC進行聚類獲得的結果之間存在差異，應傾向於相信後者（聚類）的結果。(如何使用Bioconductor進行單細胞分析？)
還在用PCA降維？快學學大牛最愛的t-SNE演算法吧, 附Python/R程式碼

10. Seurat與線粒體基因

這是一個簡單的操作問題，讀懂程式碼就好解決了。送您一句話：Some years back when I was less experienced not being sure if I did an analysis right used to keep me up at night … I am still not sure if I do it right now but I sleep better 😉

圖源：Giphy

參考來源

撰文/Tiger 編輯/生信寶典