18年1月29日宏基因組轉(zhuǎn)載了中科院生態(tài)中心鄧曄組的文章《土壤細(xì)菌定量方法結(jié)合相對豐度分析揭示種群的真實(shí)變化》。其中的圖3基于堆疊柱狀圖，添加組間各成分連線，可以更容易的觀察和比較組間的變化。如下圖：
image.png
于是我寫了一個R腳本來實(shí)現(xiàn)這個圖的畫法，腳本如下：
library(tidyverse)
library(optparse)
library(ggplot2)
option_list=list(
make_option(c("-f","--file"),type = "character",default = FALSE,
help = "The input file"),
make_option(c("-c","--colname"),type = "character",default = "None",
help="Input a file that contains the column names use to redefine the column order of the data farme,
The order is not changed by default"),
make_option(c("-o","--out"),type = "character",default = FALSE,
help = "the out put file name")
)
opt = parse_args(OptionParser(option_list = option_list, usage = "This Script is use for Stacked histogram"))
out_name=paste(opt$out,"pdf",sep = ".")
df=read.table(opt$file,sep = "\t",header = T)
if(opt$colname=="None"){
cat("The order is",colnames(df))
colname=colnames(df)
}else{
colname=read.table(opt$colname,header = F)
colname=colname$V1
df=df[,colname]
cat("The order is",colname[-1])
}
df.long <- df %>% gather(group, abundance, -Phylum)###(data,header_name,value_name,dataFarm)
group=colname[-1]
df.long$group=factor(df.long$group,levels = group,ordered = T)
link_dat <- df %>%
arrange(by=desc(Phylum)) %>%
mutate_if(is.numeric, cumsum)
bar.width <- 0.7
link_dat <- link_dat[, c(1,2,rep(3:(ncol(link_dat)-1),each=2), ncol(link_dat))]
link_dat <- data.frame(y=t(matrix(t(link_dat[,-1]), nrow=2)))
link_dat$x.1 <- 1:(ncol(df)-2)+bar.width/2
link_dat$x.2 <- 1:(ncol(df)-2)+(1-bar.width/2)
p=ggplot(df.long, aes(x=group, y=abundance, fill=Phylum)) +
geom_bar(stat = "identity", width=bar.width, col='black') +
geom_segment(data=link_dat,
aes(x=x.1, xend=x.2, y=y.1, yend=y.2), inherit.aes = F)
ggsave(p,file=out_name,width =8.3 ,height =5.8 )
腳本有三個參數(shù)
-f是指要輸入的文件，格式如下：
image.png
第一列是物種。后面四列是物種的豐度值，第一行是列名
-c是畫圖時X軸的分組排列順序，以文件的形式輸入，默認(rèn)是輸入數(shù)據(jù)的順序，文件格式如下：
image.png
-o是輸出文件的前綴
使用示例：
Rscript Stacked_histogram.R -f test.txt -c colname.txt -o 123
輸出結(jié)果如下：
image.png
