NGS系列文章包括NGS基礎(chǔ)、轉(zhuǎn)錄組分析?（Nature重磅綜述|關(guān)于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細胞測序分析?(重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程 （原理、代碼和評述）)、GEO數(shù)據(jù)挖掘（典型醫(yī)學設(shè)計實驗GEO數(shù)據(jù)分析 (step-by-step) - Limma差異分析、火山圖、功能富集）等內(nèi)容。

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

于是我看見了https://satijalab.org/seurat/，Seurat是一個R軟件包，設(shè)計用于QC、分析和探索單細胞RNA-seq數(shù)據(jù)。Seurat旨在幫助用戶能夠識別和解釋單細胞轉(zhuǎn)錄組學中的的異質(zhì)性來源，并通過整合各種類型的單細胞數(shù)據(jù)，能夠在單個細胞層面上進行系統(tǒng)分析。

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

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

1. 有關(guān)merge函數(shù)的問題

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

2. 有關(guān)PC的選擇

Seurat應(yīng)用JackStraw隨機抽樣構(gòu)建一個特征基因與主成分相關(guān)性值的背景分布，選擇富集特征基因相關(guān)性顯著的主成分用于后續(xù)分析。對大的數(shù)據(jù)集，這一步計算會比較慢，有時也可能不會找到合適的臨界點。

建議通過ElbowPlot來選，找到拐點或使得所選PC包含足夠大的variation了 (80%以上)，便合適。然后再可以在這個數(shù)目上下都選幾個值試試，最好測試的時候往下游測試些，越下游越好，看看對結(jié)果是否有影響。

另外一個方式可以是根據(jù)與各個主成分相關(guān)的基因的GSEA功能富集分析選擇合適的主成分?(一文掌握GSEA，超詳細教程)。

一般選擇7-12都合適。實際分析時，可以嘗試選擇不同的值如10,?15或所有主成分，結(jié)果通常差別不大。但如果選擇的主成分比較少如5等，結(jié)果可能會有一些變化。

另外要注意：用了這么多年的PCA可視化竟然是錯的！！！

3. 細胞周期相關(guān)基因

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

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

seurat_data <- newimport(SeuratObject)###重新定義了import函數(shù)，稱為newimport newimport <- function(otherCDS, import_all = FALSE) {if(class(otherCDS)[1] == 'Seurat') {requireNamespace("Seurat")data <- otherCDS@assays$RNA@countsif(class(data) == "data.frame") {data <- as(as.matrix(data), "sparseMatrix")}pd <- tryCatch( {pd <- new("AnnotatedDataFrame", data = otherCDS@meta.data)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 Seuratif(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 <- 0if(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 <- NULLotherCDS@misc$Monocle@auxClusteringData$scran <- NULLmonocle_cds <- otherCDS@misc$Monoclemist_list <- otherCDS} else {# mist_list <- list(ident = ident)mist_list <- otherCDS}} else {mist_list <- list()}if(1==1) {var.genes <- setOrderingFilter(monocle_cds, otherCDS@assays$RNA@var.features)}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@logExprsOffsetfd <- otherCDS@featureDatapd <- otherCDS@phenoDataexperimentData = otherCDS@experimentDataif("is.expr" %in% slotNames(otherCDS))lowerDetectionLimit <- otherCDS@is.exprelselowerDetectionLimit <- 1if(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_listmonocle_cds@auxOrderingData$scran <- mist_list} else {stop('the object type you want to export to is not supported yet')}return(monocle_cds) }

5. 不同條件下畫熱圖

• R語言 - 熱圖繪制 (heatmap)

• R語言 - 熱圖簡化

• R語言 - 熱圖美化

• 利用ComplexHeatmap繪制熱圖(一)

• 獲取pheatmap聚類后和標準化后的結(jié)果

6. Seurat對象轉(zhuǎn)化為.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.

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

7. 根據(jù)基因取細胞子集

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

# 提取特定cluster，繼續(xù)后續(xù)分析。 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. 篩選條件的設(shè)置

這是最開始讀入的初始設(shè)置，后面質(zhì)控還有更多篩選方式。

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

• 對一篇單細胞RNA綜述的評述：細胞和基因質(zhì)控參數(shù)的選擇

• 重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程 （原理、代碼和評述）

陷阱和建議：

• 通過可視化檢測到的基因數(shù)量、總分子數(shù)和線粒體基因的表達比例的分布中的異常峰來執(zhí)行QC。

聯(lián)合考慮這3個變量，而不是單獨考慮一個變量。

• 設(shè)置寬松的QC閾值；

如果下游聚類無法解釋時再重新設(shè)定嚴格的QC閾值。

• 如果樣品之間的QC變量分布不同（存在多個強峰），則需要考慮樣品質(zhì)量差異，應(yīng)按照 Plasschaert et al. (2018)的方法為每個樣品分別確定QC閾值。

9. RunTSNE不是在聚類

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

• 聚類是直接基于距離矩陣的經(jīng)典無監(jiān)督機器學習問題。

  通過最小化簇內(nèi)距離或在降維后的表達空間中鑒定密集區(qū)域，將細胞分配給簇。

  方法有層級聚類、K-menas、基于圖的聚類等 （基因共表達聚類分析及可視化，基因表達聚類分析之初探SOM）。

• tSNE和UMAP目前主要用于可視化。用于可視化的降維必然涉及信息丟失并改變細胞之間的距離。因此tSNE/UMAP圖應(yīng)僅只用于解釋或傳達基于更精確的、更多維度的定量分析結(jié)果。這樣可以保證分析充分利用了壓縮到二維空間時丟失的信息。假如二維圖上呈現(xiàn)的細胞分布與使用更多數(shù)目的PC進行聚類獲得的結(jié)果之間存在差異，應(yīng)傾向于相信后者（聚類）的結(jié)果。(如何使用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

總結(jié)

以上是生活随笔為你收集整理的单细胞分析Seurat使用相关的10个问题答疑精选！的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。