前言一般而言,我们做完pathway富集分析,就做下气泡图或bar图来进行展示,但它们实际上只考虑了富集因子和Pvalue。如果我们不关注这两个因素,而是在乎样本本身的pathway丰度呢? 对于KEGG热图绘制,大部分是做到KO层级,因为基因/蛋白和KO的绝大部分都是一对一的对应关系。如果一定要做Pathway的丰度热图呢?一般的方法是将该通路中的基因/蛋白的丰度进行累加来表示该pathway的丰度。 好了,现在我们来计算并绘制热图吧。 数据处理得到pathway富集分析结果文件一般是这样的: > colnames(path)[1] 'X.Pathway' 'Sample1..1113.' 'Sample2..15327.' 'Pvalue' 'Pathway.ID' 'Level1' [7] 'Level2' 'Proteins' 'KOs' 除此之外,我们还需要一个基因表达矩阵: 我们的目标就是整理成这样的table,用来绘制热图:
得到的结果是这样的: Proteins列中的蛋白都一一和Pathway对应起来了。后面就好办了,直接贴代码: #sum scaleibaq2 <- sweep(ibaq,2,apply(ibaq, 2, sum),FUN = '/')#caculate each group mean valuegroup <- factor(rep(c('S01CC','S11SC','S12CC','S12SC'),each=3),levels = c('S11SC','S12SC','S12CC','S01CC'))out <- apply(ibaq2,1,function(x){ dat <- data.frame(group=group,value=x) dat_mean <- dat %>% group_by(group) %>% summarise(mean=mean(value)) %>% select(mean)}) #注意此处计算均值未用na.rm参数out[[1]]out2 <- as.data.frame(t(do.call(cbind,out)))colnames(out2) <- levels(group)rownames(out2) <- rownames(ibaq2)exp <- data.frame(ProteinID=rownames(out2),out2)data1 <- left_join(path3,exp,by='ProteinID') %>% dplyr::select(1:3,6:9) %>% gather(Sample,Abundance,-c(Pathway,Level1,Level2)) %>% group_by(Pathway,Sample) %>% summarise(Sum=sum(Abundance)) %>% spread(Sample,Sum)tmp <- path3[1:3]annotation <- tmp[!duplicated(tmp),]length(intersect(data1$Pathway,annotation$Pathway))#先按pathway排序,再按level2,level1排序plotdat <- left_join(annotation,data1,by='Pathway') %>% arrange(Pathway) %>% arrange(Level2) %>% arrange(Level1) 现在已经得到想要的数据了。 绘图这个就不用多解释了。
图片大概成这样:
|
|
来自: 爱读书的云朵 > 《数据获取、统计及处理》