Seurat Weekly NO.0 || 開刊詞

在過去的一周里Seurat社區(qū)在github總提問數(shù)由上周的3066上升到3090，當(dāng)然有同一問題反復(fù)討論的情況，也有之前的問題還有人再問的情況，總的來說上周其在github中的issues往來郵件一共182封.本次主要和大家分享一下幾個，本期的封面故事是：

Resolution parameter in Seurat's FindClusters function for larger cell numbers

即，在做細(xì)胞分群的時候如何確定分群數(shù)據(jù)量？在這之前我們先來看幾個有趣的問答吧。

Error using WhichCells/subset in a for loop

https://github.com/satijalab/seurat/issues/1053
https://github.com/satijalab/seurat/issues/1839

當(dāng)我們需要循環(huán)地對某個基因列表用WhichCells/subset做條件篩選的時候：

genes <- c("KIF22", "RPGRIP1L")
v <- c(1:length(genes))
for (i in v){
  cells <- WhichCells(seurat_object, expression = genes[i] > 0)
  print(length(cells))
}

我們發(fā)現(xiàn)：

Error in FetchData(object = object, vars = expr.char[vars.use], cells = cells, : None of the requested variables were found:

而單獨(dú)使用基因名是可以的：

WhichCells(seurat_object, expression = KIF22 > 0) 

為什么呢？為什么基因名作為循環(huán)變量就不行了呢？我們發(fā)現(xiàn)循環(huán)的時候基因名是一個字符串，而單獨(dú)寫的時候并不是字符串（沒加引號），GitHub的回答是：

Much like R's subset function, subset.Seurat is designed for interactive use only. While we currently don't offer a programmatic way to subset Seurat objects based on feature expression, this can be accomplished relatively easily using which and FetchData:

var1 <- "GeneName"
expr <- FetchData(object = object, vars = var1)
object[, which(x = expr > low & expr < high)]

s所以在我們設(shè)計循環(huán)體的時候，不能直接用WhichCells/subset,可以用上面的方式：

library(Seurat)
gene_cells = list()
items <-c( 'LGALS3', 'S100A11')
for (gene1 in items)
{
    if (gene1 %in% rownames(pbmc_small))
        
    {
        expr <- FetchData(object = pbmc_small, vars = gene1)
        gene_cells[[gene1]] <- pbmc_small[, which(x = expr >= 0.5)]
        
    }
    
}

gene_cells

$LGALS3
An object of class Seurat 
230 features across 24 samples within 1 assay 
Active assay: RNA (230 features)
 2 dimensional reductions calculated: pca, tsne

$S100A11
An object of class Seurat 
230 features across 45 samples within 1 assay 
Active assay: RNA (230 features)
 2 dimensional reductions calculated: pca, tsne

關(guān)于Seurat的subset還有一個詭異的地方：

pbmc_small[VariableFeatures(object = pbmc_small), ]
pbmc_small[, 1:10]

subset(x = pbmc_small, subset = MS4A1 > 4)
subset(x = pbmc_small, subset = `DLGAP1-AS1` > 2)
subset(x = pbmc_small, idents = '0', invert = TRUE)
subset(x = pbmc_small, subset = MS4A1 > 3, slot = 'counts')
subset(x = pbmc_small, features = VariableFeatures(object = pbmc_small))

當(dāng)基因名中有-的時候，需要用反單引號引起來才能行。

Fixing batch effect with seurat integrate

biostars : https://www.biostars.org/p/390189/

上一期我們討論了如何在Seurat分析多個樣本。今天選擇這個問題的原因是它讓我們明白在使用工具的時候，對工具原理的了解是多么重要。在運(yùn)用概念時，對概念之前的關(guān)系理解是多么重要。

提問者拿出nFeatures_RNA 的箱型圖，問這是不是批次效應(yīng)造成的，這能不能說明樣本間有批次效應(yīng)？

這是一個經(jīng)典的問題，類似地我們可以拿出兩個樣本間任何有差異的東西來做出這樣的提問: 這是不是批次效應(yīng)造成的，這能不能說明樣本間有批次影響？那么，批次效應(yīng)真的是萬惡之源了。其實(shí)，在觀察到差異的時候，我們首先要問的是：為什么。然后觀察造成這種差異的原因，也許批次效應(yīng)是最后應(yīng)該考慮的。就拿這個問題來說吧，箱型圖說明不同處理間細(xì)胞捕獲的基因數(shù)量不同，那么是哪些基因不同呢可不可以做一個overlap，看看共有的基因是哪些，哪些又是缺失的？只有闡明了差異的來源才能想辦法消除。

Resolution parameter in Seurat's FindClusters function for larger cell numbers

bioinformatics : https://bioinformatics.stackexchange.com/questions/4297/resolution-parameter-in-seurats-findclusters-function-for-larger-cell-numbers
github ： https://github.com/satijalab/seurat/issues/1565

我的細(xì)胞到底分多少個群是合適的？這是一個廣泛而經(jīng)典問題。就單細(xì)胞技術(shù)而言，我們常說每個細(xì)胞都是不同的，也就是說你總可以分到最細(xì)以單細(xì)胞為單位，但是這樣就失去高通量的意義了。在低通量下，我們可以著眼于單個細(xì)胞，現(xiàn)在成千上萬的細(xì)胞，一個一個看是不切實(shí)際的。那么，我的細(xì)胞到底分多少個群是合適的？

這個問題表現(xiàn)在Seurat中就是：Finding optimal cluster resolution in Seurat 3? 我們知道，不同的resolution參數(shù)會帶來不同的分群結(jié)果。先看一下github上面的回答：

While Seurat doesn't have tools for comparing cluster resolutions, there is a tool called clustree designed for this task and works on Seurat v3 objects natively. It's available on CRAN and can be installed with a simple install.packages('clustree')

clustree我們之前講過，可以全局地查看不同分群結(jié)果：

#先執(zhí)行不同resolution 下的分群
library(Seurat)
pbmc_small <- FindClusters(
  object = pbmc_small,
  resolution = c(seq(.4,1.6,.2))
)
clustree(pbmc_small@meta.data, prefix = "RNA_snn_res.")

在clustree的圖中我們看到不同resolution的取值情況下分群的關(guān)系。既然我們最終是以群為單位來分析的，我們肯定是希望每個群是比較純的。如圖可以看到在倒數(shù)第二層級有個亞群來自不同的分群，這有可能是：

• 分群過度，把原來分群的中應(yīng)有的異質(zhì)性也提煉出來單獨(dú)作為一群了
• 上一層級分群不足，還包含了不該有的異質(zhì)性。

這里就帶來靈魂拷問了，就拿B細(xì)胞來說吧，它本身也是有異質(zhì)性的啊，那么他的異質(zhì)性是如何的呢？我們知道，某一類細(xì)胞內(nèi)的異質(zhì)性一般是要小于細(xì)胞群之間的異質(zhì)性的。所以，拿到這個圖我們就可以根據(jù)自己帶著生物學(xué)意義的期望來做一個判斷了。

其實(shí)，我們也知道分群終究是非監(jiān)督的，只是數(shù)據(jù)驅(qū)動的，并不摻雜著數(shù)據(jù)（表達(dá)譜）以外的生物學(xué)意義。如果拋開這些生物學(xué)意義，其實(shí)是有一些辦法來評價分群結(jié)果的：

要方法？一本書夠不夠？

這些方法也是在做群內(nèi)和群之間的比較，得出類似群純度的度量單位來評價分群結(jié)果。在不久前張澤民老師團(tuán)隊(duì)的一篇文章中提到過一種方法：ROGUE: an entropy-based universal metric for assessing the purity of single cell population。

該方法已被封裝為一個R包： https://github.com/PaulingLiu/ROGUE

我們看到已經(jīng)有不少的方法來做分群的評估了，還有：IKAP—Identifying K mAjor cell Population groups in single-cell RNA-sequencing analysis ：

以上這些方法大同小異，核心的問題是，或者研究者真正關(guān)心的是：

哪種分群結(jié)果的生物解釋性高？

正所謂：分析總會有結(jié)果，看你敢用不敢用。