【原】学徒作业完成！多个单细胞亚群各自差异分析后如何汇总可视化

健明 2023-09-06 发布于广东

展开全文

这个是目前为止最快被完成的学徒作业，三天前才发布的练习题：
多个单细胞亚群各自差异分析后如何汇总可视化

数据获取

首先拿到ifnb的数据

#加载需要的R包
library(SeuratData)
library(ifnb.SeuratData)
library(ReactomePA)
library(org.Hs.eg.db)
library(ggplot2)
library(Seurat)
library(pheatmap)
library(reshape2)
library(tidyr)
library(dplyr)
library(tibble)
library(RColorBrewer)
library(ggsci)
library(ggrepel)
library(ggplotify)
library(cowplot)

#获取ifnb数据

data("ifnb")
table(ifnb$orig.ident)
sce = ifnb
sce$celltype=sce$seurat_annotations
Idents(sce) <- sce$celltype 
sce$group=sce$stim
gplots::balloonplot(table(sce$group,sce$celltype )) 
Idents(sce) = sce$group
table(Idents(sce))
degs = lapply(unique(sce$celltype), function(x){
  FindMarkers(sce[,sce$celltype==x],ident.1 = 'STIM',
              ident.2 = 'CTRL')
})
names(degs) <- unique(sce$celltype)

第一种可视化：热图绘制

先找出各细胞类型上下调的gene,然后拿到gene-cell type的表达矩阵，将其分为上调的和下调的

#获取上下调基因
up <- data.frame()
down <- data.frame()
for (i in 1:length(degs)){
  dat <- degs[[i]]
  dat$cell_type <- rep(names(degs)[[i]],times=nrow(dat))
  dat$gene <- rownames(dat)
  up_regulate_df <- dat[dat$avg_log2FC>0 & dat$p_val_adj<0.05,c('avg_log2FC','cell_type','gene')]
  
  down_regulate_df <- dat[dat$avg_log2FC<0 & dat$p_val_adj<0.05,c('avg_log2FC','cell_type','gene')]
  up <- rbind(up,up_regulate_df)
  down <- rbind(down,down_regulate_df)
}

#整理上下调表达矩阵信息
up_df <- pivot_wider(up,names_from = cell_type,values_from = avg_log2FC)
up_df <- column_to_rownames(up_df,'gene')
down_df <- pivot_wider(down,names_from = cell_type,values_from = avg_log2FC)
down_df <- column_to_rownames(down_df,'gene')

再分别找到每种cell type独有的gene，得到specific gene-cell type表达矩阵

#上调
up_specific_gene <- c()
for (i in colnames(up_df)){
  target_cell_gene_up <- rownames(up_df[!is.na(up_df[,i]),])
  other_cell_df_up <- up_df[target_cell_gene_up,setdiff(colnames(up_df),i)]
  other_cell_df_up$count_na <- rowSums(!is.na(other_cell_df_up))
  other_cell_df_up <- other_cell_df_up[other_cell_df_up$count_na==0,]
  other_cell_gene_up <- rownames(other_cell_df_up)
  up_specific_gene <- c(up_specific_gene,other_cell_gene_up)
}

up_specific_gene <- up_specific_gene[!is.na(up_specific_gene)]
up_specific <- up_df[up_specific_gene,]
# 填充空值
up_specific[is.na(up_specific)] <- 0

#绘制图形
p1 <- pheatmap(as.matrix(up_specific),scale = 'row',show_rownames = F,cluster_rows = F,
         cluster_cols = F)

#下调
down_specific_gene <- c()
for (i in colnames(down_df)){
  target_cell_gene_down <- rownames(down_df[!is.na(down_df[,i]),])
  other_cell_df_down <- down_df[target_cell_gene_down,setdiff(colnames(down_df),i)]
  other_cell_df_down$count_na <- rowSums(!is.na(other_cell_df_down))
  other_cell_df_down <- other_cell_df_down[other_cell_df_down$count_na==0,]
  other_cell_gene_down <- rownames(other_cell_df_down)
  down_specific_gene <- c(down_specific_gene,other_cell_gene_down)
}
down_specific_gene <- down_specific_gene[!is.na(down_specific_gene)]
down_specific <- down_df[down_specific_gene,]
# 填充空值
down_specific[is.na(down_specific)] <- 0
pheatmap(as.matrix(down_specific),scale = 'row',show_rownames = F,cluster_rows = F,
         cluster_cols = F)

然后是shared gene的可视化热图

up_df$count_na <- rowSums(is.na(up_df))
up_shared_df <- up_df[up_df$count_na<(ncol(up_df)-2),]

down_df$count_na <- rowSums(is.na(down_df))
down_shared_df <- down_df[down_df$count_na<(ncol(down_df)-2),]

p2 <- pheatmap(as.matrix(up_shared_df[,-ncol(up_shared_df)]),
         show_rownames = F,cluster_rows = F,na_col = 'white',cluster_cols = F,
         color = colorRampPalette(rev(brewer.pal(n = 7, name = "Reds")))(100))

pheatmap(as.matrix(down_shared_df[,-ncol(down_shared_df)]),
         show_rownames = F,cluster_rows = F,na_col = 'white',cluster_cols = F,
         color = colorRampPalette(rev(brewer.pal(n = 7, name = "Blues")))(100))

最后我们可以把两个图拼起来展示

p1 <- pheatmap(as.matrix(up_specific),scale = 'row',show_rownames = F,cluster_rows = F,
         cluster_cols = F,legend=F,show_colnames = F)
p2 <- pheatmap(as.matrix(up_shared_df[,-ncol(up_shared_df)]),legend = F,
         show_rownames = F,cluster_rows = F,na_col = 'white',cluster_cols = F,
         color = colorRampPalette(rev(brewer.pal(n = 7, name = "Reds")))(100))
plot_grid(as.ggplot(p1),as.ggplot(p2),nrow = 2)

第二种可视化方式

先把所有cell type的差异gene合并成一张表，还需要获取上下调topn的gene的列表，方便后面进行文本标记

df2 <- data.frame()
top_n_df <- data.frame()
for (i in 1:length(degs)){
  dat <- degs[[i]]
  dat$cell_type <- names(degs)[[i]]
  dat$gene <- rownames(dat)
  dat$change <- ifelse(dat$avg_log2FC>0 & dat$p_val_adj<0.05,'up',
                       ifelse(dat$avg_log2FC<0 & dat$p_val_adj<0.05,'down','nochange'))
  dat_sig <- dat[dat$change!='nochange',]
  up_top2 <- top_n(dat_sig[dat_sig$change=='up',],3,avg_log2FC)[,'gene']
  down_top2 <- top_n(dat_sig[dat_sig$change=='down',],-3,avg_log2FC)[,'gene']
  label_gene <- c(up_top2,down_top2)
  top_n_df <- rbind(top_n_df,dat[label_gene,]) 
  df2 <- rbind(df2,dat)
}

rownames(df2) <- 1:nrow(df2)

df3 <- df2
df2 <- df2[df2$change!='nochange',]

p1 <- ggplot()+geom_jitter(data = df2,aes(x=cell_type,y=avg_log2FC,color=change),width=0.2,size=1.5)
p1

添加相应的文本标记

p2 <- p1+geom_jitter(data = top_n_df,aes(x=cell_type,y=avg_log2FC,color=change),width=0.2,size=1.5)+
  geom_text_repel(data=top_n_df,aes(x=cell_type,y=avg_log2FC,label=gene),size=3)+
  theme(axis.text.x = element_text(angle = 90,hjust = 1))

画出基本的柱子

up_bar <- top_n(group_by(df2,cell_type),1,avg_log2FC)
down_bar <- top_n(group_by(df2,cell_type),-1,avg_log2FC)
bar_df <- data.frame(row.names=up_bar$cell_type,x=up_bar$cell_type,
                     up=up_bar$avg_log2FC+0.3,
                     down=down_bar$avg_log2FC-0.3)
bar_df <- bar_df[sort(unique(df2$cell_type)),]
# 有些细胞没有下调gene
bar_df$down <- ifelse(bar_df$down>=0,NA,bar_df$down)

ggplot()+geom_col(data = bar_df,aes(x=x,y=up),alpha = 0.2)+
  geom_col(data = bar_df,aes(x=x,y=down),alpha = 0.2)

p3 <- p2+geom_col(data = bar_df,aes(x=x,y=up),alpha = 0.2)+
  geom_col(data = bar_df,aes(x=x,y=down),alpha = 0.2)
p3

再插入中间的色块

box_df <- data.frame(row.names=up_bar$cell_type,x=up_bar$cell_type,y=0)

p4 <- p3+geom_tile(data = box_df,aes(x=x,y=y),height=0.8,fill=pal_d3('category20')(13))
p4

第三种可视化方式

其实就是把上面的图翻转一下

ggplot()+geom_jitter(data = df3[df3$change!='nochange',],aes(x=avg_log2FC,y=cell_type,color=cell_type))+
  geom_jitter(data = df3[df3$change=='nochange',],aes(x=avg_log2FC,y=cell_type,color=cell_type),alpha=0.1)+
  geom_vline(xintercept = 0,linewidth=0.5)+theme_classic()+
  theme(axis.title.y = element_blank(),axis.ticks.y = element_blank(),axis.line.y = element_blank())