什么,需要5Tb内存啊!

众所周知,单细胞数据分析对计算机资源的要求有一点高,尤其是10X单细胞转录组数据的多样本合并那个步骤!

最近我就接到一个粉丝咨询,说他想处理一个公共数据集,只有8个原位肿瘤+3个转移肿瘤的10X单细胞转录组样品,但是数据处理的过程发现系统提示说需要5Tb内存,虽然说他自己有一个512G内存的服务器,但是也承受不起5Tb内存,问我有没有渠道!

额,给他配置一个5Tb内存服务器倒是简单,我自己就有2.5T内存的服务器,不就是加倍嘛!不过,我注意到他就是11个10X转录组样品,理论上不可能是需要5Tb内存的,所以让他把代码发过来我检查看看.

首先是读取多个10x单细胞转录组文件夹,需要保证每个文件夹下面都是有3个文件哦:

# For output from CellRanger < 3.0
#  Should show barcodes.tsv, genes.tsv, and matrix.mtx
# For output from CellRanger >= 3.0 with multiple data types
 # Should show barcodes.tsv.gz, features.tsv.gz, and matrix.mtx.gz
sceList = lapply(samples,function(pro){ 
  #pro=samples[1]
  folder=file.path('outputs',pro ) 
  print(pro)
  print(folder)
  print(list.files(folder))
  print( Sys.time() )
  sce=CreateSeuratObject(counts = Read10X(folder), 
                     project = pro )
  print( Sys.time() )
  return(sce)
})

然后是直接走CCA整合:

sceList
for (i in 1:length(sceList)) {
  print( Sys.time() )
  sceList[[i]] <- NormalizeData(sceList[[i]], verbose = FALSE)
  sceList[[i]] <- FindVariableFeatures(sceList[[i]], selection.method = "vst"
                                             nfeatures = 2000, verbose = FALSE)
}
sceList
Sys.time()
sce.anchors <- FindIntegrationAnchors(object.list = sceList, dims = 1:30)
# load('sce.anchors.RDATA')
Sys.time()
sce.integrated <- IntegrateData(anchorset = sce.anchors, dims = 1:30)
Sys.time()

看起来中规中矩,就是我一直教学使用的代码。没办法,我只好让他把数据发过来了。

我自己也读取看看,让我留意到了它居然每个样品有70万个细胞!!!所以我猜测应该是他的10X的3个文件里面并没有过滤,把全部的barcode输出了,我就给他加上了一个简单的检查代码,以及两个标准过滤:

lapply(sceList, function(x) print(x))
# retaining cells that had unique
# molecular identifiers (UMIs) greater than 400, 
# expressed 100 and 8000 genes
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
sceList = lapply(sceList, function(x) {
  subset(x, subset = nFeature_RNA > 200 & nFeature_RNA < 8000 & nCount_RNA   > 400)
})

lapply(sceList, function(x) print(x))


# had mitochondrial content less than 10 percent.
 
sceList =  sceList = lapply(sceList, function(x) {
  rownames(x)[grepl('^mt-',rownames(x),ignore.case = T)]
  x[["percent.mt"]] <- PercentageFeatureSet(x, pattern = "^MT-")
  subset(x, subset =  percent.mt < 10 )
})
lapply(sceList, function(x) print(x))

接下来使用我的32G内存的Mac,都可以走这个CCA流程啦!

太有意思了,为什么我想讲解这个故事呢,因为在很多交流群都看到有粉丝问内存不够,实际上很多情况下,内存不够是因为你代码学的很差。

如果你连512G内存都没有呢?

临时使用的话,可以考虑我们的共享云哦!见:共享云服务器又又又又又来了(只有50个名额) 

生物医学科研方法

李兰娟院士就被 Nature 评为“2020年度十大科学人物”发表声明

2021-1-15 5:14:35

生物医学科研方法

真开心一分钟赚100块钱

2021-1-15 5:24:51