好久没更了,最近梅雨季节,天天下雨,真的好难过啊。☔️
到暑假了,自己也是忙得不行,再加上拖延症晚期,真的是事事都没有进展。😭
今天继续更一下hdWGCNA
吧。😘
rm(list = ls())
library(Seurat)
library(tidyverse)
library(cowplot)
library(patchwork)
library(WGCNA)
library(hdWGCNA)
记得设置一下哦!~😂
theme_set(theme_cowplot())
set.seed(123)
enableWGCNAThreads(nThreads = 8)
load('hdWGCNA_object.Rdata')
这里我们比较来自INH
的细胞,不同性别分组间的差异。😘
group1 <- seurat_obj@meta.data %>% subset(cell_type == 'INH' & msex == 0) %>% rownames
group2 <- seurat_obj@meta.data %>% subset(cell_type == 'INH' & msex != 0) %>% rownames
head(group1)
head(group2)
利用FindDMEs
函数进行差异分析。🧐
DMEs <- FindDMEs(
seurat_obj,
barcodes1 = group1,
barcodes2 = group2,
test.use='wilcox',
wgcna_name='tutorial'
)
head(DMEs)
可视化一下结果吧!~🙃
棒棒糖图展示。🍭
PlotDMEsLollipop(
seurat_obj,
DMEs,
wgcna_name='tutorial',
pvalue = "p_val_adj"
)
火山图展示一下!~😜
PlotDMEsVolcano(
seurat_obj,
DMEs,
wgcna_name = 'tutorial'
)
我们可以写个loop
来运行DME
分析,在多个Cluster
中进行比较。😘
在这里,我们将对INH
的亚群进行一下FindDMEs
差异分析。😏
clusters <- c("INH1 VIP+", "INH2 PVALB+", "INH3 SST+", 'INH4 LAMP5+', "INH5 PVALB+")
DMEs <- data.frame()
for(cur_cluster in clusters){
# identify barcodes for group1 and group2 in eadh cluster
group1 <- seurat_obj@meta.data %>% subset(annotation == cur_cluster & msex == 0) %>% rownames
group2 <- seurat_obj@meta.data %>% subset(annotation == cur_cluster & msex != 0) %>% rownames
# run the DME test
cur_DMEs <- FindDMEs(
seurat_obj,
barcodes1 = group1,
barcodes2 = group2,
test.use='wilcox',
pseudocount.use=0.01, # we can also change the pseudocount with this param
wgcna_name = 'tutorial'
)
# add the cluster info to the table
cur_DMEs$cluster <- cur_cluster
# append the table
DMEs <- rbind(DMEs, cur_DMEs)
}
head(DMEs)
热图展示下结果吧!~🙊
# get the modules table:
modules <- GetModules(seurat_obj)
mods <- levels(modules$module); mods <- mods[mods != 'grey']
# make a copy of the DME table for plotting
plot_df <- DMEs
# set the factor level for the modules so they plot in the right order:
plot_df$module <- factor(as.character(plot_df$module), levels=mods)
# set a min/max threshold for plotting
maxval <- 0.5; minval <- -0.5
plot_df$avg_log2FC <- ifelse(plot_df$avg_log2FC > maxval, maxval, plot_df$avg_log2FC)
plot_df$avg_log2FC <- ifelse(plot_df$avg_log2FC < minval, minval, plot_df$avg_log2FC)
# add significance levels
plot_df$Significance <- gtools::stars.pval(plot_df$p_val_adj)
# change the text color to make it easier to see
plot_df$textcolor <- ifelse(plot_df$avg_log2FC > 0.2, 'black', 'white')
# make the heatmap with geom_tile
p <- plot_df %>%
ggplot(aes(y=cluster, x=module, fill=avg_log2FC)) +
geom_tile()
# add the significance levels
p <- p +
geom_text(label=plot_df$Significance, color=plot_df$textcolor)
# customize the color and theme of the plot
p <- p +
scale_fill_viridis_b()+
RotatedAxis() +
theme(
panel.border = element_rect(fill=NA, color='black', size=1),
axis.line.x = element_blank(),
axis.line.y = element_blank(),
plot.margin=margin(0,0,0,0)
) + xlab('') + ylab('') +
coord_equal()
p
在这里,我们将group.by
对每种细胞类型进行One-versus-all
分析。🦀
group.by = 'cell_type'
DMEs_all <- FindAllDMEs(
seurat_obj,
group.by = 'cell_type',
wgcna_name = 'tutorial'
)
head(DMEs_all)
p <- PlotDMEsVolcano(
seurat_obj,
DMEs_all,
wgcna_name = 'tutorial',
plot_labels=F,
show_cutoff=F
)
# facet wrap by each cell type
p + facet_wrap(~group, ncol=3)
我们还可以使用module expression scores
而不是module eigengenes
进行差异分析。🧐
但要先使用函数ModuleExprScore
计算。🧮
我们以UCell
为例。🙊
seurat_obj <- ModuleExprScore(
seurat_obj,
n_genes = 25,
method='UCell'
)
group1 <- seurat_obj@meta.data %>% subset(cell_type == 'INH' & msex == 0) %>% rownames
group2 <- seurat_obj@meta.data %>% subset(cell_type == 'INH' & msex != 0) %>% rownames
DMEs_scores <- FindDMEs(
seurat_obj,
features = 'ModuleScores',
barcodes1 = group1,
barcodes2 = group2,
test.use='wilcox',
wgcna_name='tutorial'
)
PlotDMEsLollipop(
seurat_obj,
DMEs_scores,
wgcna_name = 'tutorial',
pvalue = "p_val_adj"
)