单细胞差异基因的分析很多人一直使用的是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科研分享与服务
