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 绘图结果