投稿问答最小化  关闭

万维书刊APP下载

单细胞差异基因分析:pseudobulk(muscat包)

2023/1/29 9:08:27  阅读:290 发布者:

单细胞差异基因的分析很多人一直使用的是Seurat自带函数Findmarkers。但是也有人提出这个分析结果不准确,之前有篇Nature communications专门研究过不同的分析方法,文章链接:https://www.nature.com/articles/s41467-021-25960-2

众所周知,我们做单细胞差异基因的目的是为了分析两组间细胞的差异,而不是单个细胞本身的差异。这里我们使用muscat包的pseudobulk分析,muscat可以进行多组分析,我们看一看差异分析效果,并与Findmarkers比较一下结果!

muscat的官方说明和流程见如下链接:

http://www.bioconductor.org/packages/release/bioc/vignettes/muscat/inst/doc/analysis.html

However,无论是哪种方法,都可以任人使用,很多高分文章中也在使用Findmarkers,不能说它不好吧。接下来,我们跟随muscat的标准流程做一下分析。详细注释代码和示例数据已上传群文件,请自行下载使用!!!

安装加载包:

library(Seurat)

library(muscat)

library(SingleCellExperiment)

library(dplyr)

human_data <- readRDS("D:/KS项目/公众号文章/human_data.rds")

human_sce <- subset(human_data, celltype=='Mast', invert=T)

human_sce$celltype = droplevels(human_sce$celltype,

                                exclude = setdiff(levels(human_sce$celltype),

                                                  unique(human_sce$celltype)))

转化对象:

DefaultAssay(human_sce) <- "RNA"

human_sce <- as.SingleCellExperiment(human_sce)

human_sce <- prepSCE(human_sce,

                     kid = "celltype",

                     gid = "group",#group——id

                     sid = "orig.ident",#sample_id,

                     drop=T)

human_sce$group_id <- factor(human_sce$group_id, levels=c("BM", "GM"))

nk <- length(kids <- levels(human_sce$cluster_id))

ns <- length(sids <- levels(human_sce$sample_id))

names(kids) <- kids

names(sids) <- sids

t(table(human_sce$cluster_id, human_sce$sample_id))

将单细胞数据聚合为pseudobulk data,看看分群!

pb <- aggregateData(human_sce,

                    assay = "counts",

                    fun = "sum",

                    by = c("cluster_id", "sample_id"))

(pb_mds <- pbMDS(pb))

差异分析:

pb$group_id <- factor(pb$group_id, levels=c("BM", "GM"))

res <- pbDS(pb, verbose = FALSE)#pseudobulk DS analysis

tmp <- human_sce

counts(tmp) <- as.matrix(counts(tmp))

result_table <- resDS(tmp, res, bind = "row", frq = FALSE, cpm = FALSE)

rm(tmp)

将表达某一基因的BM/GM细胞群中的比例添加到结果中,类似于Findmarker分析结果中的pct1pct2.

human_sce <- subset(human_data, celltype=='Mast', invert=T)

count_mat <- as.matrix(human_sce[["RNA"]]@data) > 0

cluster_list <- unique(result_table$cluster_id)

result_table$BM.frq <- 0

BM_cells <- colnames(human_sce)[human_sce$group == "BM"]

for(i in 1:length(cluster_list)){

  cluster_cells <- colnames(human_sce)[human_sce$celltype == cluster_list[i]]

  test_cells <- intersect(BM_cells, cluster_cells)

  row_ind <- which(result_table$cluster_id == cluster_list[i])

  frq <- rowSums(count_mat[result_table$gene[row_ind],test_cells]) / length(test_cells)

  result_table$BM.frq[row_ind] <- frq

}

GM_cells <- colnames(human_sce)[human_sce$group == "GM"]

result_table$GM.frq <- 0

for(i in 1:length(cluster_list)){

  cluster_cells <- colnames(human_sce)[human_sce$celltype == cluster_list[i]]

  test_cells <- intersect(GM_cells, cluster_cells)

  row_ind <- which(result_table$cluster_id == cluster_list[i])

  frq <- rowSums(count_mat[result_table$gene[row_ind],test_cells]) / length(test_cells)

  result_table$GM.frq[row_ind] <- frq

}

write.csv(result_table, file = "result_table.csv", row.names = F)

对比下结果:set1Findmarkersset2muscat,可以看出差异基因共同的不少,但是不相同的也挺多。

Mac <- subset(human_data, celltype=="Macrophage")

Mac_DEGs <- FindMarkers(human_data,

                        min.pct = 0,

                        logfc.threshold = 0.25,

                        group.by = 'group',

                        ident.1 ="GM",

                        ident.2="BM")

Mac_DEGs_Findmarker_sig <- Mac_DEGs[Mac_DEGs$p_val < 0.05 &

                                      abs(Mac_DEGs$avg_log2FC) > 0.25, ]

####muscat

Mac_muscat_DEGs <- subset(result_table, cluster_id=='Macrophage')

Mac_muscat_DEGs_sig <- Mac_muscat_DEGs[Mac_muscat_DEGs$p_val < 0.05 &

                                         abs(Mac_muscat_DEGs$logFC) > 0.25&

                                         Mac_muscat_DEGs$BM.frq>0.1&

                                         Mac_muscat_DEGs$GM.frq>0.1, ]

rownames(Mac_muscat_DEGs_sig) <- Mac_muscat_DEGs_sig$gene

library(devtools)

install_github("js229/Vennerable")

library(Vennerable)

Set1 <- as.list(rownames(Mac_DEGs_Findmarker_sig))

Set2 <- as.list(Mac_muscat_DEGs_sig$gene)

example <-list(Set1=Set1,Set2=Set2)

Veenplot <- Venn(example)

Veenplot<-Veenplot[, c("Set1", "Set2")]

plot(Veenplot, doWeights = TRUE)

做一下共同1099个基因的相关性。

same_gene <- Veenplot@IntersectionSets$`11`

data <- cbind(Mac_DEGs_Findmarker_sig[same_gene,][,2],

              Mac_muscat_DEGs_sig[same_gene,][,3])

colnames(data) <- c("Seurat", "Muscat")

data <- as.data.frame(data)

library(ggpubr)

ggscatter(data,x="Seurat",y="Muscat",

          add = "reg.line",

          conf.int = T,

          color = '#0f8096')+

  stat_cor(label.x = 0.2, label.y = 0)

可以看出,两种方法共同的差异基因上下调是一致的,只有一个基因相反。还是那句话,方法无所谓好坏,都可以使用,主要是讲顺自己的故事。好了这就是所有内容了,关于单细胞差异基因的分析方法还有很多,这里的muscat只是一个,其他的感兴趣的可以去探索。

转自:KS科研分享与服务”微信公众号

如有侵权,请联系本站删除!


  • 万维QQ投稿交流群    招募志愿者

    版权所有 Copyright@2009-2015豫ICP证合字09037080号

     纯自助论文投稿平台    E-mail:eshukan@163.com