前言
一般而言蹈丸，我們做完pathway富集分析后成黄，會做下氣泡圖或bar圖來進行展示，這樣實際上只考慮了富集因子和Pvalue逻杖。如果我們不關(guān)注這兩個因素奋岁，而是在乎樣本本身的pathway豐度呢？
對于KEGG熱圖繪制荸百，大部分是做到KO層級闻伶，因為基因/蛋白和KO的絕大部分都是一對一的對應(yīng)關(guān)系，十分方便地得到想要的結(jié)果够话。如果一定要做Pathway的豐度熱圖呢蓝翰？一般的方法是將該通路中的基因/蛋白的豐度進行累加來表示該pathway的豐度。
好了女嘲，現(xiàn)在我們來計算并繪制熱圖吧畜份。
數(shù)據(jù)處理
得到pathway富集分析結(jié)果文件一般是這樣的：
Proteins字段中的基因/蛋白是用分號隔開的。
> colnames(path)
[1] "X.Pathway" "Sample1..1113." "Sample2..15327." "Pvalue" "Pathway.ID" "Level1"
[7] "Level2" "Proteins" "KOs"
除此之外欣尼，我們還需要一個基因表達矩陣：
這個數(shù)據(jù)有四組樣本爆雹，每組3個重復(fù)，共12個樣本愕鼓。
我們的目標(biāo)就是整理成這樣的table钙态，用來繪制熱圖：
從兩個表可知，數(shù)據(jù)處理關(guān)鍵就是pathway中的蛋白豐度求和菇晃。把pathway中對應(yīng)的各蛋白展開册倒，再匹配到表達矩陣上，最后歸并求和就好了磺送，思路清晰了就動手吧剩失。
library(tidyverse)
path2 <- path %>% dplyr::select(X.Pathway,Level1,Level2,Proteins)
#下面這一步最關(guān)鍵，dplyr中為我們提供了一個有用的函數(shù)unnest
path3 <- path2 %>% mutate(ProteinID = strsplit(Proteins, ";")) %>% unnest()
colnames(path3)[1] <- "Pathway"
#如果不熟悉册着，這一步也可用Map函數(shù)配合do.call來完成：
out <- do.call(rbind, Map(cbind, path2$X.Pathway,path2$Level1,path2$Level2,strsplit(path2$Proteins, ";")))
out <- as.data.frame(out)
colnames(out) <- colnames(path2)
處理后得到的結(jié)果是這樣的：
Proteins列中的蛋白都一一和Pathway對應(yīng)起來，后面就好辦了脾歧，直接貼代碼：
#sum scale
ibaq2 <- sweep(ibaq,2,apply(ibaq, 2, sum),FUN = "/")
#caculate each group mean value
group <- 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參數(shù)
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)
現(xiàn)在已經(jīng)得到想要的數(shù)據(jù)了。
繪圖
這個就不用多解釋了鞭执。
library(pheatmap)
Exp_log2=plotdat #實際上我中間還進行了其他處理司顿，這里便于繪圖直接賦值
colnames(Exp_log2)
exp_plot <- select(Exp_log2,S11SC,S12SC,S12CC,S01CC)
rownames(exp_plot) <- Exp_log2$Pathway
annotation_row <- select(Exp_log2,Level2,Level1)
rownames(annotation_row) <- Exp_log2$Pathway
pheatmap(exp_plot,cluster_rows = F,cluster_cols = F,scale = "row",
annotation_row = annotation_row,
border_color = NA,
#angle_col=45,
color = colorRampPalette(c("blue","white","red"))(50))
圖片大概成這樣：
根據(jù)自己需要挑選一些pathway展示吧芒粹，太多不好看。
Ref: https://stackoverflow.com/questions/28719088/r-semicolon-delimited-a-column-into-rows
