区区20万个单细胞居然超内存了

我有一个备用的小mac电脑，只有64G内存，32线程，4T硬盘。最近搞活动，详见：春节期间单细胞转录组数据分析全免费，收到了上百个需求，就把它拿出来了避免吃灰。

其中一个数据集是2020发在NC的肺癌单细胞文章：《Single-cell RNA sequencing demonstrates the molecular and cellular reprogramming of metastatic lung adenocarcinoma》，是44个肺癌病人的58个10x技术的单细胞转录组样品，因为是两三年前的10x单细胞技术那个时候都比较辣鸡所以细胞数量就平均每个样品是3千，这样的话合计是20万个单细胞。

作者提供的是 GSE131907_Lung_Cancer_raw_UMI_matrix.txt.gz 文件，我们来演示一下读取它，并且进行单细胞降维聚类分群哈：

代码语言：javascript复制a=fread('inputs/GSE131907_Lung_Cancer_raw_UMI_matrix.txt.gz',

data.table = F)

a[1:4,1:4]

rownames(a)=a[,1]

a=a[,-1]

# 208,506 cells derived from 58 lung adenocarcinomas from 44 patients,

# 11 tumour, 11 distant normal lung, 10 normal lymph node, and 10 metastatic brain tissue samples

# were obtained from LUAD patients conserving surgery without prior treatment.

# 7 metastatic lymph node and 4 lung tumour tissue samples from advanced stage LUAD patient

# endobronchial ultrasound (EBUS) and bronchoscopy (BRONCHO).

kp1 = grepl('LN',colnames(a))

table(kp1)

kp2 = grepl('T',colnames(a))

table(kp2 )

a[1:4,1:4]

dim(a)

因为数据量很小，所以读取起来完全没有压力，有了两万个基因在20万个单细胞里面的表达量矩阵，就可以直接构建单细胞对象了，但是CreateSeuratObject(counts = ct ) 会报错，居然是内存问题。

虽然我这个小mac电脑，只有64G内存，32线程，4T硬盘仅仅是备用电脑平时吃灰的， 但是不应该是区区20万个单细胞都搞不定啊，我思考了一下，大概是因为在R里面的两万个基因在20万个单细胞里面的表达量矩阵存储格式不经济，我简单看了看，居然占用了20G内存，怪不得呢。

本来是想把它转为稀疏矩阵，但是也失败了，代码如下所示：

代码语言：javascript复制

if(F){

format(object.size(a) ,units='Gb')

A <- as( a , "sparseMatrix")

B <- Matrix(as.matrix(a), sparse = TRUE)

library(Matrix)

library(mltools)

A = mltools::sparsify(a)

format(object.size(A) ,units='Gb')

}

然后灵机一动，虽然这个两万个基因在20万个单细胞里面的表达量矩阵失败，但是它本来就是44个肺癌病人的58个10x技术的单细胞转录组样品，我只需要拆分构建Seurat对象，最后合并这些对象，也是同样的操作啊；

代码语言：javascript复制

sampl = paste0(phe[,2],phe[,3])

sceList = lapply(unique(sampl) ,function(pro){

# pro=unique(sampl)[1]

print(pro)

sce =CreateSeuratObject(counts = a[,sampl %in% pro],

project = gsub( '_matrix','',pro) ,

min.cells = 5,

min.features = 500 )

return(sce)

})

names(sceList)

sce.all=merge(x=sceList[[1]],

y=sceList[ -1 ] )

as.data.frame(sce.all@assays$RNA@counts[1:10, 1:2])

head(sce.all@meta.data, 10)

table(sce.all@meta.data$orig.ident)

这次果然就不报错了，我们有了对象，后面的质量控制和降维聚类分群就顺理成章的事情了。

然后质量控制并且降维聚类分群即可质量控制其实每个数据集不一样的，取决于单细胞转录组来源的技术，特殊情况下需要个性化的修改 scRNA_scripts/qc.R 里面的代码：

代码语言：javascript复制###### step2:QC质控 ######

dir.create("./1-QC")

setwd("./1-QC")

# 如果过滤的太狠,就需要去修改这个过滤代码

source('../scRNA_scripts/qc.R')

sce.all.filt = basic_qc(sce.all)

print(dim(sce.all))

print(dim(sce.all.filt))

setwd('../')

但是多样品整合，以及后续的降维聚类分群这个代码在任意单细胞转录组数据集都是大概率不需要修改的：

代码语言：javascript复制sp='human'

###### step3: harmony整合多个单细胞样品 ######

dir.create("2-harmony")

getwd()

setwd("2-harmony")

source('../scRNA_scripts/harmony.R')

# 默认 ScaleData 没有添加"nCount_RNA", "nFeature_RNA"

# 默认的

sce.all.int = run_harmony(sce.all.filt)

setwd('../')

###### step4: 降维聚类分群和看标记基因库 ######

# 原则上分辨率是需要自己肉眼判断,取决于个人经验

# 为了省力,我们直接看 0.1和0.8即可

table(Idents(sce.all.int))

table(sce.all.int$seurat_clusters)

table(sce.all.int$RNA_snn_res.0.1)

table(sce.all.int$RNA_snn_res.0.8)

getwd()

dir.create('check-by-0.1')

setwd('check-by-0.1')

sel.clust = "RNA_snn_res.0.1"

sce.all.int <- SetIdent(sce.all.int, value = sel.clust)

table(sce.all.int@active.ident)

source('../scRNA_scripts/check-all-markers.R')

setwd('../')

getwd()

dir.create('check-by-0.8')

setwd('check-by-0.8')

sel.clust = "RNA_snn_res.0.8"

sce.all.int <- SetIdent(sce.all.int, value = sel.clust)

table(sce.all.int@active.ident)

source('../scRNA_scripts/check-all-markers.R')

setwd('../')

getwd()

然后就是合理的生物学命名了我给大家准备的基因列表主要是针对肿瘤单细胞的第一层次降维聚类分群 ， 是：

immune (CD45+,PTPRC),epithelial/cancer (EpCAM+,EPCAM),stromal (CD10+,MME,fibo or CD31+,PECAM1,endo)参考我前面介绍过 CNS图表复现08—肿瘤单细胞数据第一次分群通用规则，这3大单细胞亚群构成了肿瘤免疫微环境的复杂。绝大部分文章都是抓住免疫细胞亚群进行细分，包括淋巴系（T,B,NK细胞）和髓系（单核，树突，巨噬，粒细胞）的两大类作为第二次细分亚群。但是也有不少文章是抓住stromal 里面的fibo 和endo进行细分，并且编造生物学故事的。

我这里给出来了的单细胞亚群的生物学名字是：

代码语言：javascript复制 celltype=data.frame(ClusterID=0:13 ,

celltype= 0:13)

#定义细胞亚群

celltype[celltype$ClusterID %in% c( 0, 11,13 ),2]='Tcells'

celltype[celltype$ClusterID %in% c( 6 ),2]='Bcells'

celltype[celltype$ClusterID %in% c( 5 ),2]='mast'

celltype[celltype$ClusterID %in% c( 3 ),2]='cycle'

celltype[celltype$ClusterID %in% c( 4 ),2]='fibro'

celltype[celltype$ClusterID %in% c( 7 ),2]='endo'

celltype[celltype$ClusterID %in% c( 9 ),2]='pDC'

celltype[celltype$ClusterID %in% c( 1 ),2]='myeloids'

celltype[celltype$ClusterID %in% c( 2,12,8 ),2]='epi'

最后的降维聚类分群如下所示：

降维聚类分群

当然了，降维聚类分群仅仅是故事的开始，后面还可以根据实验设计进行多种多样的分析。