加载所需的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>

▼更多精彩推荐，请关注我们▼

image