单细胞差异基因的分析很多人一直使用的是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分析结果中的pct1和pct2.

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)

对比下结果：set1是Findmarkers，set2是muscat，可以看出差异基因共同的不少，但是不相同的也挺多。

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科研分享与服务”微信公众号

如有侵权，请联系本站删除！