Seurat分組隨機選取細(xì)胞數(shù)實戰(zhàn)（隨機采樣后找亞群DEG） 2022-06-01

關(guān)鍵詞

  • 隨機取樣細(xì)胞
  • Downsample cells
  • 分組隨機選取細(xì)胞

適用背景

之前的博客提到，R語言處理大數(shù)據(jù)效率較低，耗時長，一種解決方案是可以轉(zhuǎn)用Python語言流程，但如果對Python語言比較陌生，任務(wù)又急，那可以采用另一種方案——分組隨機取樣。
盡管Seurat這個軟件包功能極其強大，但是當(dāng)細(xì)胞數(shù)達到幾十萬甚至上百萬時，把常規(guī)流程跑一遍少則幾天，多則幾周，實在是極其消耗時間。而且有時吧，只是單純想測試一下某些參數(shù)或者流程是否可用，如果用全數(shù)據(jù)集來測試實在有點浪費時間，所有可用考慮分組隨機選取細(xì)胞數(shù)進行分析。

主函數(shù)

這里封裝了一個函數(shù)sample_seob，以下是參數(shù)解釋：

  • obj Seurat對象
  • group.by 分組名，默認(rèn)使用聚類結(jié)果seurat_clusters
  • sp.size 取樣大小，也就是對分組里的每一個類別選取的細(xì)胞數(shù)，例如設(shè)置為100，將對cluster 1取100個細(xì)胞，cluster 2也取100個細(xì)胞，以此類推，如果某個cluster的細(xì)胞數(shù)不足100個將選取這個cluster的所有細(xì)胞。
  • diet 是否使用 DietSeurat函數(shù)對Seurat對象進行瘦身，默認(rèn)為true，因為如果Seurat對象包含scale.data等信息會很耗內(nèi)存，瘦身后能減少內(nèi)存并加快分析速度。
  • sp.total 選取的總細(xì)胞數(shù)，為可選項，有時候不知道每個亞群取多少細(xì)胞數(shù)合適，只想大概取到一定的細(xì)胞數(shù)，例如1000，就可以用sp.total參數(shù)，注意，只有sp.size沒有賦值的時候，sp.total參數(shù)才會生效，即這兩個參數(shù)是二選一即可。
sample_seob <- function(obj,group.by="seurat_clusters",sp.size=NULL,diet="true",sp.total=1000) {
all <- obj
if (diet=="true") {
all <- DietSeurat(all)
}
if (is.null(sp.size)) {
nlen <- length(unique(all@meta.data[,group.by]))
sp.size <- ceiling(sp.total/nlen)
}
seob_list <- list()
i <- 1
for (sc in unique(all@meta.data[,group.by])){
cellist <- colnames(all)[which(all@meta.data[,group.by] == sc)]
ob <- subset(all, cells=cellist)
if (length(colnames(ob)) > sp.size) {
ob <- subset(ob,cells=sample(colnames(ob), sp.size))
}
seob_list[[i]] <- ob
i <- i+1
}
all <- Reduce(merge,seob_list)
return(all)
}

腳本示例

數(shù)據(jù)集采用Seurat內(nèi)置數(shù)據(jù)集pbmc_small，有80個細(xì)胞，按RNA_snn_res.1分組有3種類型

> pbmc_small
An object of class Seurat
230 features across 80 samples within 1 assay
Active assay: RNA (230 features, 20 variable features)
 2 dimensional reductions calculated: pca, tsne
> unique(pbmc_small$RNA_snn_res.1)
[1] 0 2 1
Levels: 0 1 2

使用sp.size=10，按RNA_snn_res.1分組，每種類型取10個細(xì)胞。

> all <- sample_seob(pbmc_small,sp.size=10,group.by='RNA_snn_res.1')
> all
An object of class Seurat
230 features across 30 samples within 1 assay
Active assay: RNA (230 features, 0 variable features)

使用sp.total=50，按RNA_snn_res.1分組，取大約50個細(xì)胞。

> all <- sample_seob(pbmc_small,sp.total=50,group.by='RNA_snn_res.1')
> all
An object of class Seurat
230 features across 51 samples within 1 assay
Active assay: RNA (230 features, 0 variable features)

分析實戰(zhàn)——分組隨機采樣后找每個亞群的DEG

在實際的分析中，我發(fā)現(xiàn)FindAllMarkers經(jīng)常跑著跑著就斷了，出現(xiàn)以下報錯:

Warning: When testing 15 versus all:
        The total size of the 5 globals that need to be exported for the future expression (‘FUN()’) is 2.06 GiB. This exceeds the maximum allowed size of 500.00 MiB (option 'future.globals.maxSize'). The three largest globals are ‘data.use’ (2.06 GiB of class ‘S4’), ‘group.info’ (2.56 MiB of class ‘list’) and ‘FUN’ (9.76 KiB of class ‘function’).

這種情況一般是因為數(shù)據(jù)集太大了，由運行內(nèi)存不足導(dǎo)致，這種情況目前找到的解決辦法是隨機取樣后找DEG，結(jié)果也比較可靠，運行起來也快很多。以下是示例封裝的函數(shù)**subset_deg **，部分參數(shù)與上面一致，其它參數(shù)與FindAllMarkers一致。

subset_deg <- function(obj,group.by="seurat_clusters",sp.size=NULL,output="./",
                       min.pct=0.25,logfc.threshold=0.25,only.pos=F,assays ="RNA",order=F) {
all <- DietSeurat(obj)
if (!is.null(sp.size)) {
seob_list <- list()
i <- 1
for (sc in unique(all@meta.data[,group.by])){
cellist <- colnames(all)[which(all@meta.data[,group.by] == sc)]
ob <- subset(all, cells=cellist)
if (length(colnames(ob)) > sp.size) {
ob <- subset(ob,cells=sample(colnames(ob), sp.size))
}
seob_list[[i]] <- ob
i <- i+1
}
all <- Reduce(merge,seob_list)
}
all_markers <- FindAllMarkers(all, only.pos = only.pos, min.pct = min.pct, logfc.threshold = logfc.threshold, verbose = F,assays = assays,order=order)
write.table(all_markers, paste0(output, "deg_sample",sp.size,".xls"),sep="\t",quote=F)
}

也可以簡寫成下面的形式：

subset_deg <- function(obj,group.by="seurat_clusters",sp.size=NULL,output="./",
                       min.pct=0.25,logfc.threshold=0.25,only.pos=F,assays ="RNA",order=F,diet="true",sp.total=1000) {
all <- sample_seob(obj, sp.size=sp.size,group.by=group.by,diet=diet,sp.total=sp.total)
all_markers <- FindAllMarkers(all, only.pos = only.pos, min.pct = min.pct, logfc.threshold = logfc.threshold, verbose = F,assays = assays,order=order)
write.table(all_markers, paste0(output, "deg_sample",sp.size,".xls"),sep="\t",quote=F)
}

升級版

因為上面的版本循環(huán)里每次都要subset，太麻煩了，最后還要merge一下，比較耗時間，因此修改代碼為下面的結(jié)果，只需要subset一次，也不用merge，能較快得到結(jié)果。

Sample_seob <- function(obj,group.by="seurat_clusters",sp.size=NULL,diet="true",sp.total=1000) {
all <- obj
if (diet=="true") {
all <- DietSeurat(all,dimreducs = c('pca','umap'))
}

if (is.null(sp.size)) {
nlen <- length(unique(all@meta.data[,group.by]))
sp.size <- ceiling(sp.total/nlen)

}
ncellist <- c()
for (sc in unique(all@meta.data[,group.by])){
cellist <- colnames(all)[which(all@meta.data[,group.by] == sc)]
if (length(cellist) > sp.size) {
cellist=sample(cellist, sp.size)
}
ncellist <- c(ncellist,cellist)
}
all <- subset(all,cells=ncellist)
return(all)
}

極簡版

突然看到官網(wǎng)有函數(shù)是能簡便實現(xiàn)這種分組隨機取細(xì)胞數(shù)的，轉(zhuǎn)換Idents后就能根據(jù)Idents信息分類別隨機選取細(xì)胞。

Idents(pbmc) <- "orig.ident"
# Downsample the number of cells per identity class
ob1 <- subset(x = pbmc, downsample = 100)

小結(jié)與補充

不多說了，搬磚去了。祝大家六一兒童節(jié)快樂！

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容