NGS系列文章包括NGS基礎(chǔ)、轉(zhuǎn)錄組分析 （Nature重磅綜述|關(guān)于RNA-seq你想知道的全在這）、ChIP-seq分析 （ChIP-seq基本分析流程）、單細(xì)胞測(cè)序分析 (重磅綜述：三萬字長(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è)剛剛開始進(jìn)行單細(xì)胞轉(zhuǎn)錄組分析的菜鳥，R語言底子沒有，有時(shí)候除了會(huì)copy外，如果你讓我寫個(gè)for循環(huán)，我只能cross my fingers。。。。
于是我看見了https://satijalab.org/seurat/，Seurat是一個(gè)R軟件包，設(shè)計(jì)用于QC、分析和探索單細(xì)胞RNA-seq數(shù)據(jù)。Seurat旨在幫助用戶能夠識(shí)別和解釋單細(xì)胞轉(zhuǎn)錄組學(xué)中的的異質(zhì)性來源，并通過整合各種類型的單細(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è)，還不懂，來找我，重磅綜述：三萬字長(zhǎng)文讀懂單細(xì)胞RNA測(cè)序分析的最佳實(shí)踐教程 （原理、代碼和評(píng)述）)，這次我就介紹一下里面讓我曾經(jīng)困惑的幾個(gè)問題以及比較nice的解答?。钗胰f萬沒想到，找答案比自己寫答案確實(shí)困難的多。。。）
1. 有關(guān)merge函數(shù)的問題
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)。建議通過
ElbowPlot來選，找到拐點(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ì)有一些變化。
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è)問題我也遇到了，并且已經(jīng)有人給出了解決方案。(生信寶典注：官方最開始沒支持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@counts
if(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 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@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@logExprsOffset
fd <- otherCDS@featureData
pd <- otherCDS@phenoData
experimentData = otherCDS@experimentData
if("is.expr" %in% slotNames(otherCDS))
lowerDetectionLimit <- otherCDS@is.expr
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對(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è)問題，好在易生信提供了好的解決方案。來這試試？易生信線下課推遲，線上課程免費(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è)置
這是最開始讀入的初始設(shè)置，后面質(zhì)控還有更多篩選方式。
質(zhì)控是用于保證下游分析時(shí)數(shù)據(jù)質(zhì)量足夠好。由于不能先驗(yàn)地確定什么是足夠好的數(shù)據(jù)質(zhì)量，所以只能基于下游分析結(jié)果（例如，細(xì)胞簇注釋）來對(duì)其進(jìn)行判斷。因此，在分析數(shù)據(jù)時(shí)可能需要反復(fù)調(diào)整參數(shù)進(jìn)行質(zhì)量控制 (生信寶典注：反復(fù)分析，多次分析，是做生信的基本要求。這也是為啥需要上易生信培訓(xùn)班而不是單純依賴公司的分析。)。從寬松的QC閾值開始并研究這些閾值的影響，然后再執(zhí)行更嚴(yán)格的QC，通常是有益的。這種方法對(duì)于細(xì)胞種類混雜的數(shù)據(jù)集特別有用，在這些數(shù)據(jù)集中，特定細(xì)胞類型或特定細(xì)胞狀態(tài)可能會(huì)被誤解為低質(zhì)量的異常細(xì)胞。在低質(zhì)量數(shù)據(jù)集中，可能需要嚴(yán)格的質(zhì)量控制閾值。具體見：
陷阱和建議：
- 通過可視化檢測(cè)到的基因數(shù)量、總分子數(shù)和線粒體基因的表達(dá)比例的分布中的異常峰來執(zhí)行QC。
聯(lián)合考慮這3個(gè)變量，而不是單獨(dú)考慮一個(gè)變量。
- 設(shè)置寬松的QC閾值；
如果下游聚類無法解釋時(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)典無監(jiān)督機(jī)器學(xué)習(xí)問題。
通過最小化簇內(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ì)胞分析？)
10. Seurat與線粒體基因
這是一個(gè)簡(jiǎ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
