目的：計算密度曲線對應(yīng)的峰值并作圖。
library(ggplot2)
library(reshape2)
mytheme2 <- theme_bw() + theme(legend.title=element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
axis.text.x=element_text(size=6,angle=0), legend.position="top")
argv<-commandArgs(TRUE)
data <- read.table(argv[1], header = TRUE, sep="\t")
densFindPeak <- function(x){
td <- density(x)
maxDens <- which.max(td$y)
list(x=td$x[maxDens],y=td$y[maxDens])
}
peak = densFindPeak(data$depth)
peak
p1 <- ggplot(data,aes(depth, alpha=0.8)) +
geom_density() +
geom_rug() +
geom_vline(xintercept = peak[[1]],colour="red",linetype="dashed") +
geom_text(aes(x=peak[[1]], y= peak[[2]]+0.3,label=paste("peak =",round(peak[[1]],4))), color="black") +
ylab("Density") +
mytheme2
ggsave(file=argv[2], plot=p1, width = 5, height = 4)
