🤩 hdWGCNA | 单细胞数据怎么做WGCNA呢！？（四）（差异分析）

生信漫卷

发布于 2024-07-12 16:07:55

900

发布于 2024-07-12 16:07:55

文章被收录于专栏：R语言及实用科研软件

写在前面

好久没更了，最近梅雨季节，天天下雨，真的好难过啊。☔️

到暑假了，自己也是忙得不行，再加上拖延症晚期，真的是事事都没有进展。😭

今天继续更一下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进行多个clusters的差异分析

我们可以写个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

One-versus-all DME差异分析

在这里，我们将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"
)