在文章单细胞转录组中的pseudocell又是什么中，我们介绍了pseudocell的概念，并且在文章最后贴上了sc-MCA的计算代码。作为一个Seurat的深度用户，我们不禁要想：能不能把这段代码写成可以接受Seurat对象的函数呢？而且要保留Seurat的一般风格。下面，就让我们试试吧。
首先要获得不同assay，即对哪个assay来计算？获得assay之后，assay的哪个slot？已经确定对哪种分群来计算。确定参数后，我们来写代码：
library(purrr)
GatherData <- function(object, ...) {
UseMethod("GatherData")
}
GatherData.Seurat <- function(object,
assay,
slot_use,
...) {
assay <- assay %||% "RNA"
slot_use <- slot_use %||% "data"
obj_data <- GetAssayData(
object = object,
assay = assay,
slot = slot_use
) %>%
as.matrix()
return(obj_data)
}
这段代码来获取assay和slot的表达矩阵，返回一个Seurat的对象。
PseudoCell <- function(object,
assay_use = NULL,
slot_use = NULL,
cluster_use =NULL,
pseudocell.size =20){
message("tips:
Cluster_use : one col in metadata
pseudocell.size : how many cell will be pseudo,should min(cell number of cluster) < pseudocell.size <max(cell number of cluster)")
Inter<- GatherData(object = object,
assay = assay_use,
slot_use = slot_use)
Inter[Inter<0]=0
idd<-object@meta.data
Inter.id<-cbind(rownames(idd),as.vector(idd[,cluster_use]))
rownames(Inter.id)<-rownames(idd)
colnames(Inter.id)<-c("CellID","Celltype")
Inter.id<-as.data.frame(Inter.id)
Inter1<-Inter[,Inter.id$CellID]
Inter<-as.matrix(Inter1)
pseudocell.size = pseudocell.size ## 10 test
new_ids_list = list()
for (i in 1:length(levels(Inter.id$Celltype))) {
cluster_id = levels(Inter.id$Celltype)[i]
cluster_cells <- rownames(Inter.id[Inter.id$Celltype == cluster_id,])
cluster_size <- length(cluster_cells)
pseudo_ids <- floor(seq_along(cluster_cells)/pseudocell.size)
pseudo_ids <- paste0(cluster_id, "_Cell", pseudo_ids)
names(pseudo_ids) <- sample(cluster_cells)
new_ids_list[[i]] <- pseudo_ids
}
new_ids <- unlist(new_ids_list)
new_ids <- as.data.frame(new_ids)
new_ids_length <- table(new_ids)
new_colnames <- rownames(new_ids) ###add
all.data<-Inter[,as.character(new_colnames)] ###add
all.data <- t(all.data)###add
new.data<-aggregate(list(all.data[,1:length(all.data[1,])]),
list(name=new_ids[,1]),FUN=mean)
rownames(new.data)<-new.data$name
new.data<-new.data[,-1]
new_ids_length<-as.matrix(new_ids_length)##
short<-which(new_ids_length<10)##
new_good_ids<-as.matrix(new_ids_length[-short,])##
result<-t(new.data)[,rownames(new_good_ids)]
rownames(result)<-rownames(Inter)
Tool(object) <- list(pseuoData = result,metaData = new_ids)
#object <- LogSeuratCommand(object, return.command = TRUE)
return(object)
}
pseudocell.size 的意思是几个单细胞变成一个pseudocell，注意不要小于最小细胞群的细胞数，也不要大于最大细胞群的细胞数哦。
这个函数完成计算，并用Tool函数来吧计算的结果保存在Seurat的对象中。下面让我们来试试吧
library(Seurat)
head(pbmc_small@meta.data)
orig.ident nCount_RNA nFeature_RNA RNA_snn_res.0.8 letter.idents groups RNA_snn_res.1
ATGCCAGAACGACT SeuratProject 70 47 0 A g2 0
CATGGCCTGTGCAT SeuratProject 85 52 0 A g1 0
GAACCTGATGAACC SeuratProject 87 50 1 B g2 0
TGACTGGATTCTCA SeuratProject 127 56 0 A g2 0
AGTCAGACTGCACA SeuratProject 173 53 0 A g2 0
TCTGATACACGTGT SeuratProject 70 48 0 A g1 0
mypbmc <- PseudoCell(pbmc_small, "RNA","data","RNA_snn_res.1",35)
tips:
Cluster_use : one col in metadata
pseudocell.size : how many cell will be pseudo,should min(cell number of cluster) < pseudocell.size <max(cell number of cluster)
查看我们运行的结果,pseudocell矩阵，这个可以拿来再构建Seurat对象，也可以传给其他Bulk-RNA的分析工具。
(mypbmc@tools$PseudoCell$pseuoData)[1:4,1:3]
0_Cell0 1_Cell0 2_Cell0
MS4A1 0.0000000 0.2626990 2.976558
CD79B 0.5289529 1.1268602 3.274061
CD79A 0.0000000 0.6454348 2.822091
HLA-DRA 1.7836394 5.3341723 6.170035
新旧barcode的对应关系,这个可以AddMetaData到我们的mypbmc的对象中。
head(mypbmc@tools$PseudoCell$metaData)
new_ids
GGCATATGGGGAGT 0_Cell0
ATCATCTGACACCA 0_Cell0
CTTCATGACCGAAT 0_Cell0
AGATATACCCGTAA 0_Cell0
CTAAACCTCTGACA 0_Cell0
CATTACACCAACTG 0_Cell0
注意一下 pseudocell的 命名规则： 0_Cell0。_ 之前是细胞群，Cell之后是该群的第几个pseudocell（从零开始编号）。当然，你可以根据自己的心绪，自行命名。
最后，大家看到函数中注释掉的那一行了嘛：
#object <- LogSeuratCommand(object, return.command = TRUE)
本来打算把这个函数的参数也记录在mypbmc@commands的里面的，可惜简书的编辑器页面太小，写不下了。感兴趣的同学可以尝试一下啊。
这样，我们就为Seurat写了一个函数啦。以后相对自己的scrna数据做什么操作，直接以函数的形式嫁接到Seurat里就可以啦。
Seurat只是一个工具吗？不，它已经变成我们的一部分了。
网友评论