R语言中绘制曼哈顿图--线段图

 

 

001 完整代码

dat <- read.table("result.fst", header = T)
head(dat)
chrlen <- vector()

for (i in unique(dat$CHR)) {
  temp <- dat[dat$CHR == i, ]
  temp <- temp[order(temp[,3], decreasing = T),]
  chrlen <- c(chrlen, temp[,3][1])
}
gap <- 0.005 * sum(chrlen)

chrpos <- cumsum(chrlen[-length(chrlen)] + gap)
chrpos <- c(0, chrpos)

names(chrpos) <- unique(dat$CHR)
dat$POS <- dat$POS + chrpos[dat$CHR]

colsource <- c("red", "black", "blue", "green", "cyan", "purple")
colchr <- rep(colsource, 10)[1:length(chrlen)]
chr <- unique(dat$CHR)
idx <- match(dat$CHR, chr)
dat$COL <- colchr[idx]

chrloc <- vector()

for (i in unique(dat$CHR)) {
  chrloc <- c(chrloc, chrpos[i] + round(1/2 * chrlen[i]))
}

par(mai = c(1, 1, 0.6, 0.6), mgp = c(3, 1, 0))
plot(dat$POS, dat$FST, ylim = c(0, 0.8),type = "n",ann = F, axes = F, yaxs = "i", xaxs = "i",bty = "l")
for (i in 1:nrow(dat)) {
  segments(dat$POS[i],0, dat$POS[i], dat$FST[i], col = dat$COL[i])  
}
axis(2,at = c(0,0.2,0.4,0.6, 0.8),labels = c(0,0.2,0.4,0.6, 0.8),
     font=2,tck=-0.015,lwd=3,las = 1,lwd.ticks =3,cex.axis= 1,line = 0)

axis(1,at = chrloc,labels = unique(dat$CHR),
     font=2,tck=-0.015,lwd=3,las = 1,lwd.ticks =3,cex.axis= 1,line = 0)
axis(1,at = c(0, dat$POS[nrow(dat)]),labels = NA,
     font=2,tck=0,lwd=3,las = 1,lwd.ticks =3,cex.axis= 1,line = 0)

threshold <- quantile(dat$FST, 0.99)
abline(h = threshold, lwd = 3, col = "grey")

mtext(side = 2,expression(Fst),adj = 0.5,line = 2.8,font = 2,cex=1.8)
mtext(side = 1,"Chromosome",adj = 0.5,line = 2.8,font = 2,cex=1.7)

 

002 运行过程：

> dat <- read.table("result.fst", header = T)
> head(dat)          ## 数据结构
  CHR      SNP    POS NMISS      FST
1   1 snp00001  55910   100 0.256630
2   1 snp00002  85204   100 0.000000
3   1 snp00003 122948   100 0.101524
4   1 snp00004 203750   100 0.181254
5   1 snp00005 312707   100 0.052965
6   1 snp00006 400518   100 0.000000
> chrlen <- vector()
> for (i in unique(dat$CHR)) {
+   temp <- dat[dat$CHR == i, ]
+   temp <- temp[order(temp[,3], decreasing = T),]
+   chrlen <- c(chrlen, temp[,3][1])
+ }
> gap <- 0.005 * sum(chrlen)
> chrpos <- cumsum(chrlen[-length(chrlen)] + gap)
> chrpos <- c(0, chrpos)
> names(chrpos) <- unique(dat$CHR)
> dat$POS <- dat$POS + chrpos[dat$CHR]
> colsource <- c("red", "black", "blue", "green", "cyan", "purple")
> colchr <- rep(colsource, 10)[1:length(chrlen)]
> chr <- unique(dat$CHR)
> idx <- match(dat$CHR, chr)
> dat$COL <- colchr[idx]
> chrloc <- vector()
> for (i in unique(dat$CHR)) {
+   chrloc <- c(chrloc, chrpos[i] + round(1/2 * chrlen[i]))
+ }
> par(mai = c(1, 1, 0.6, 0.6), mgp = c(3, 1, 0))
> plot(dat$POS, dat$FST, ylim = c(0, 0.8),type = "n",ann = F, axes = F, yaxs = "i", xaxs = "i",bty = "l")
> for (i in 1:nrow(dat)) {
+   segments(dat$POS[i],0, dat$POS[i], dat$FST[i], col = dat$COL[i])  
+ }
> axis(2,at = c(0,0.2,0.4,0.6, 0.8),labels = c(0,0.2,0.4,0.6, 0.8),
+      font=2,tck=-0.015,lwd=3,las = 1,lwd.ticks =3,cex.axis= 1,line = 0)
> axis(1,at = chrloc,labels = unique(dat$CHR),
+      font=2,tck=-0.015,lwd=3,las = 1,lwd.ticks =3,cex.axis= 1,line = 0)
> axis(1,at = c(0, dat$POS[nrow(dat)]),labels = NA,
+      font=2,tck=0,lwd=3,las = 1,lwd.ticks =3,cex.axis= 1,line = 0)
> threshold <- quantile(dat$FST, 0.99)
> abline(h = threshold, lwd = 3, col = "grey")
> mtext(side = 2,expression(Fst),adj = 0.5,line = 2.8,font = 2,cex=1.8)
> mtext(side = 1,"Chromosome",adj = 0.5,line = 2.8,font = 2,cex=1.7)

 

003 绘图结果