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

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

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

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

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

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

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

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

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

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

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

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

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

3. 細(xì)胞周期相關(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對(duì)象導(dǎo)入Monocle

其實(shí)這個(gè)問(wèn)題我也遇到了，并且已經(jīng)有人給出了解決方案。(生信寶典注*：官方最開(kāi)始沒(méi)支持Seurat3，我寫了個(gè)簡(jiǎn)易的，在**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語(yǔ)言 - 熱圖繪制 (heatmap)

• R語(yǔ)言 - 熱圖簡(jiǎn)化

• R語(yǔ)言 - 熱圖美化

• 利用ComplexHeatmap繪制熱圖(一)

• 獲取pheatmap聚類后和標(biāo)準(zhǔn)化后的結(jié)果

6. Seurat對(duì)象轉(zhuǎn)化為.h5ad對(duì)象

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)換確實(shí)是一個(gè)問(wèn)題，好在易生信提供了好的解決方案。來(lái)這試試？**易生信線下課推遲，線上課程免費(fèi)看*

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

也可以提取特定Cluster，進(jìn)行進(jìn)一步細(xì)分。

# 提取特定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è)置

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

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

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

• 重磅綜述：三萬(wàn)字長(zhǎng)文讀懂單細(xì)胞RNA測(cè)序分析的最佳實(shí)踐教程 （原理、代碼和評(píng)述）

陷阱和建議：

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

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

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

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

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

9. RunTSNE不是在聚類

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

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

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

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

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

• 還在用PCA降維？快學(xué)學(xué)大牛最愛(ài)的t-SNE算法吧, 附Python/R代碼

10. Seurat與線粒體基因

這是一個(gè)簡(jiǎn)單的操作問(wèn)題，讀懂代碼就好解決了。送您一句話：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个问题答疑精选！的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

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