在過去的一周里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:
而單獨使用基因名是可以的：
WhichCells(seurat_object, expression = KIF22 > 0)
為什么呢超凳？為什么基因名作為循環(huán)變量就不行了呢？我們發(fā)現(xiàn)循環(huán)的時候基因名是一個字符串耀态，而單獨寫的時候并不是字符串（沒加引號）轮傍，GitHub的回答是：
Much like R's
subsetfunction,subset.Seuratis 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 andFetchData:
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分析多個樣本。今天選擇這個問題的原因是它讓我們明白在使用工具的時候仙逻，對工具原理的了解是多么重要驰吓。在運用概念時揍魂，對概念之前的關(guān)系理解是多么重要。
提問者拿出nFeatures_RNA 的箱型圖棚瘟，問這是不是批次效應(yīng)造成的现斋，這能不能說明樣本間有批次效應(yīng)？
這是一個經(jīng)典的問題偎蘸，類似地我們可以拿出兩個樣本間任何有差異的東西來做出這樣的提問: 這是不是批次效應(yīng)造成的庄蹋，這能不能說明樣本間有批次影響？那么迷雪，批次效應(yīng)真的是萬惡之源了限书。其實，在觀察到差異的時候章咧，我們首先要問的是：為什么倦西。然后觀察造成這種差異的原因，也許批次效應(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ì)胞盾饮，一個一個看是不切實際的采桃。那么，我的細(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
Seuratv3 objects natively. It's available on CRAN and can be installed with a simpleinstall.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ì)性也提煉出來單獨作為一群了
- 上一層級分群不足，還包含了不該有的異質(zhì)性沙廉。
這里就帶來靈魂拷問了拘荡，就拿B細(xì)胞來說吧，它本身也是有異質(zhì)性的啊撬陵，那么他的異質(zhì)性是如何的呢珊皿？我們知道，某一類細(xì)胞內(nèi)的異質(zhì)性一般是要小于細(xì)胞群之間的異質(zhì)性的巨税。所以蟋定，拿到這個圖我們就可以根據(jù)自己帶著生物學(xué)意義的期望來做一個判斷了。
其實草添，我們也知道分群終究是非監(jiān)督的驶兜，只是數(shù)據(jù)驅(qū)動的，并不摻雜著數(shù)據(jù)（表達(dá)譜）以外的生物學(xué)意義远寸。如果拋開這些生物學(xué)意義抄淑，其實是有一些辦法來評價分群結(jié)果的：
這些方法也是在做群內(nèi)和群之間的比較蝇狼，得出類似群純度的度量單位來評價分群結(jié)果。在不久前張澤民老師團(tuán)隊的一篇文章中提到過一種方法：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é)果纽哥，看你敢用不敢用钠乏。
