R语言绘制曼哈顿图——散点图
作者：互联网
1、
完整代码：
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 <- sum(chrlen) * 0.005
chrpos <- cumsum(chrlen[-length(chrlen)] + gap)
chrpos <- c(0, chrpos)
chrpos
names(chrpos) <- unique(dat[,1])
dat$POS <- dat$POS + chrpos[dat$CHR]
head(dat)
col_source <- c("red", "black", "blue", "green", "purple", "cyan", "yellow", "magenta")
length(col_source)
chr <- unique(dat$CHR)
chr_col <- rep(col_source, 10)[1:length(chr)]
col_idx <- match(dat$CHR, chr)
dat$COL <- chr_col[col_idx]
chr_loc <- vector()
for (i in 1:length(chrpos)) {
chr_loc <- c(chr_loc, chrpos[i] + round(1/2 * chrlen[i]))
}
par(mai = c(1, 1, 0.8, 0.8) , mgp = c(3, 1, 1))
plot(dat$POS, dat$FST, ylim = c(0,1),col = dat$COL, pch = 16, ann = F, bty = "l", axes = F, yaxs="i",xaxs="i")
axis(2,at = c(0,0.2,0.4,0.6, 0.8, 1.0),labels = c(0,0.2,0.4,0.6, 0.8, 1.0),
font=2,tck=-0.015,lwd=3,las = 1,lwd.ticks =3,cex.axis= 1,line = 0)
axis(1,at = chr_loc,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 = "brown")
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)
运行过程：
> 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 <- sum(chrlen) * 0.005
> chrpos <- cumsum(chrlen[-length(chrlen)] + gap)
> chrpos <- c(0, chrpos)
> chrpos
[1] 0 283289573 540089908 771981780 899208687 1014025094
[7] 1138579096 1246568806 1345089729 1447636880 1541994479 1611963821
> names(chrpos) <- unique(dat[,1])
> dat$POS <- dat$POS + chrpos[dat$CHR]
> 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
> col_source <- c("red", "black", "blue", "green", "purple", "cyan", "yellow", "magenta")
> length(col_source)
[1] 8
> chr <- unique(dat$CHR)
> chr_col <- rep(col_source, 10)[1:length(chr)]
> col_idx <- match(dat$CHR, chr)
> dat$COL <- chr_col[col_idx]
> chr_loc <- vector()
> for (i in 1:length(chrpos)) {
+ chr_loc <- c(chr_loc, chrpos[i] + round(1/2 * chrlen[i]))
+ }
> par(mai = c(1, 1, 0.8, 0.8) , mgp = c(3, 1, 1))
> plot(dat$POS, dat$FST, ylim = c(0,1),col = dat$COL, pch = 16, ann = F, bty = "l", axes = F, yaxs="i",xaxs="i")
> axis(2,at = c(0,0.2,0.4,0.6, 0.8, 1.0),labels = c(0,0.2,0.4,0.6, 0.8, 1.0),
+ font=2,tck=-0.015,lwd=3,las = 1,lwd.ticks =3,cex.axis= 1,line = 0)
> axis(1,at = chr_loc,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 = "brown")
> 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)
绘图效果：
标签：曼哈顿,散点图,dat,100,font,绘制,col,lwd,axis 来源： https://www.cnblogs.com/liujiaxin2018/p/16279321.html