加载所需的R包
library(ggplot2)
library(reshape2)
设置工作路径
setwd("/Users/liang/Desktop")
data <- read.table("test_ks.txt",header = T,sep="\t")
data <- melt(data,variable.name="Species",value.name = "Ks")
#去除缺失的行
data = na.omit(data)
head(data)
## Species Ks
## 1 SpeciesA_SpeciesB 0.0915
## 2 SpeciesA_SpeciesB 0.2535
## 3 SpeciesA_SpeciesB 0.0386
## 4 SpeciesA_SpeciesB 0.1385
## 5 SpeciesA_SpeciesB 0.1125
## 6 SpeciesA_SpeciesB 0.1960
提取不同类型的数据
Sab <- data[data$Species=="SpeciesA_SpeciesB",]
Saa <- data[data$Species=="SpeciesA_SpeciesA",]
Sbb <- data[data$Species=="SpeciesB_SpeciesB",]
定义一个寻找密度图峰值的函数
densFindPeak <- function(x){
td <- density(x)
maxDens <- which.max(td$y)
list(x=td$x[maxDens],y=td$y[maxDens])
}
限定取值范围
Sab_limit <- Sab$Ks[Sab$Ks>=0 & Sab$Ks<=1]
Saa_limit <- Saa$Ks[Saa$Ks>=0 & Saa$Ks<=1]
Sbb_limit <- Sbb$Ks[Sbb$Ks>=0 & Sbb$Ks<=1]
获得峰值
abPeak = densFindPeak(Sab_limit)
aaPeak = densFindPeak(Saa_limit)
bbPeak = densFindPeak(Sbb_limit)
使用ggplot2绘图
p1 <- ggplot(data,aes(Ks,fill=Species,color=Species,alpha=0.8)) +
geom_density() + xlim(0,1) + theme_classic() + geom_rug() +
geom_vline(xintercept = abPeak[[1]],colour="red",linetype="dashed") +
geom_text(aes(x=abPeak[[1]], y= abPeak[[2]]+0.3,label=paste("Ks =",round(abPeak[[1]],4))),color="black") +
geom_vline(xintercept = aaPeak[[1]],colour="red",linetype="dashed") +
geom_text(aes(x=aaPeak[[1]]-0.02, y= aaPeak[[2]]+0.3,label=paste("Ks =",round(aaPeak[[1]],4))),color="black") +
geom_vline(xintercept = bbPeak[[1]],colour="red",linetype="dashed") +
geom_text(aes(x=bbPeak[[1]]+0.02, y= bbPeak[[2]]+0.7,label=paste("Ks =",round(bbPeak[[1]],4))),color="black") +
guides(alpha=FALSE) +
theme(axis.title = element_text(size=16),axis.text=element_text(size=16),legend.position = "top")
p1
image
Find the density plot peak and second peak value
library(ggplot2)
head(faithful)
## eruptions waiting
## 1 3.600 79
## 2 1.800 54
## 3 3.333 74
## 4 2.283 62
## 5 4.533 85
## 6 2.883 55
p1 <- ggplot(faithful, aes(waiting)) + geom_density()
image
density(faithful$waiting)
##
## Call:
## density.default(x = faithful$waiting)
##
## Data: faithful$waiting (272 obs.); Bandwidth 'bw' = 3.988
##
## x y
## Min. : 31.04 Min. :5.880e-06
## 1st Qu.: 50.27 1st Qu.:1.949e-03
## Median : 69.50 Median :1.224e-02
## Mean : 69.50 Mean :1.299e-02
## 3rd Qu.: 88.73 3rd Qu.:1.909e-02
## Max. :107.96 Max. :3.660e-02
head(density(faithful$waiting)$y)
## [1] 8.951193e-06 1.017743e-05 1.152920e-05 1.301184e-05 1.474573e-05
## [6] 1.666005e-05
head(density(faithful$waiting)$x)
## [1] 31.03732 31.18786 31.33840 31.48894 31.63948 31.79002
MaxY1_index <- which.max(density(faithful$waiting)$y)
MaxY1_index
## [1] 326
MaxX1 <- density(faithful$waiting)$x[MaxY1_index]
MaxX1
## [1] 79.96245
p2 <- ggplot(faithful, aes(waiting)) + geom_density() +
geom_vline(xintercept = density(faithful$waiting)$x[MaxY1_index],linetype="dashed")
p2
image
MaxY2_index <- which.max(density(faithful$waiting)$y[density(faithful$waiting)$x < 60])
MaxY2_index
## [1] 151
MaxX2 <- density(faithful$waiting)$x[MaxY2_index]
MaxX2
## [1] 53.61815
p3 <- ggplot(faithful, aes(waiting)) + geom_density() +
geom_vline(xintercept = density(faithful$waiting)$x[MaxY2_index],linetype="dashed")
p3
image
参考来源：http://ianmadd.github.io/pages/PeakDensityDistribution.html
sessionInfo()
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: OS X El Capitan 10.11.3
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] zh_CN.UTF-8/zh_CN.UTF-8/zh_CN.UTF-8/C/zh_CN.UTF-8/zh_CN.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggpubr_0.1.7.999 magrittr_1.5 reshape2_1.4.3 ggplot2_3.0.0
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.18 rstudioapi_0.7 bindr_0.1.1 knitr_1.20
## [5] tidyselect_0.2.4 munsell_0.5.0 colorspace_1.3-2 R6_2.2.2
## [9] rlang_0.2.2 stringr_1.3.1 plyr_1.8.4 dplyr_0.7.6
## [13] tools_3.5.1 grid_3.5.1 gtable_0.2.0 withr_2.1.2
## [17] htmltools_0.3.6 assertthat_0.2.0 yaml_2.2.0 lazyeval_0.2.1
## [21] rprojroot_1.3-2 digest_0.6.16 tibble_1.4.2 crayon_1.3.4
## [25] bindrcpp_0.2.2 purrr_0.2.5 glue_1.3.0 evaluate_0.11
## [29] rmarkdown_1.10 labeling_0.3 stringi_1.2.4 compiler_3.5.1
## [33] pillar_1.3.0 scales_1.0.0 backports_1.1.2 pkgconfig_2.0.2</pre>
