搜索
查看: 923|回复: 0

绘制Ks密度分布图

[复制链接]

4

主题

8

帖子

237

积分

中级会员

Rank: 3Rank: 3

积分
237
发表于 2018-8-26 13:25:25 | 显示全部楼层 |阅读模式
[AppleScript] 纯文本查看 复制代码
library(ggplot2)
library(reshape2)
setwd("/Users/Davey/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)
# 提取不同类型的数据
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.2,label=paste("Ks =",abPeak[[1]])),color="black") +
  geom_vline(xintercept = aaPeak[[1]],colour="red",linetype="dashed") +
  geom_text(aes(x=aaPeak[[1]]-0.02, y= aaPeak[[2]]+0.2,label=paste("Ks =", aaPeak[[1]])),color="black") +
  geom_vline(xintercept = bbPeak[[1]],colour="red",linetype="dashed") +
  geom_text(aes(x=bbPeak[[1]]+0.02, y= bbPeak[[2]]+0.4,label=paste("Ks =",bbPeak[[1]])),color="black") +
  guides(alpha=FALSE) +
  theme(axis.title = element_text(size=16),axis.text=element_text(size=16),legend.position = "top")
p1

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:1KG数据库的妙用:利用连锁和重组补全缺失基因型
下一篇:我是来试试怎么发帖的我一会就老实删掉
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-22 02:14 , Processed in 0.029381 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.