GWAS曼哈顿图总结

一页多GWAS图

文件如下：

内部格式如下：

贴出代码

# setwd("路径")
library(dplyr)
library(stringr)
library(tidyverse)
library(qqman)
pc <- 10
#添加画布
#pdf(paste0("mafplink_miami.pdf"),width=40, height=10)
tiff(filename = "CML333-eigen_res.tiff",width = 5000,height = 1500,res = 140)
layout(matrix(1:10,2,5,byrow = T))
par(oma = c(0,0,3,0))
for(ipc in 1:pc){
  # ipc=1
  d = read.table(paste0("CML333-eigen_res_",ipc,".txt"),as.is = T,header = T,sep = "\t",quote = "",na.strings = "Inf")
  d$chrom <- as.integer(str_replace(d$chrom,"chr",""))
  # d$chrom[74237:nrow(d)] <- as.integer(d$Chr[74237:nrow(d)])
  d <- d[-6]
  d <- na.omit(d)
  colnames(d)[2] = "BP"
  a <- c(37,37,37,37,37,37,37,37,313,37)
  #设置前提条件
  P1="PGC"
  P2="Fst"
  Log1=TRUE
  Log2=FALSE
  colors = c("#EE7600", "#458B00")
  suggestiveline=-log10(1e-5)
  title="" 
  colnames(d)[1] = "CHR"
  pidx1=which(colnames(d) == P1)
  pidx2 = which(colnames(d) == P2)
  
  if (!("CHR" %in% names(d) & "BP" %in% names(d) & P1 %in% names(d) & P2 %in% names(d))) stop("Make sure your data frame contains columns CHR, BP, and P")
  d=subset(na.omit(d[order(d$CHR, d$BP), ]), (!is.na(d[,pidx1]) & !is.na(d[,pidx2]))) # remove na's, sort, and keep only 0<P<=1
  d <- na.omit(d)
  d[d$CHR==11,]$BP <- d[d$CHR==11,]$BP*3000
  if (Log1){
    d$logp = -log10(d[,pidx1])
  } else {
    d$logp = d[,pidx1]
  }
  
  if (Log2){
    d$logp2 = log10(d[,pidx2])
  } else {
    d$logp2 = -1*d[,pidx2]
  }
  
  ytick = c(ceiling(max(abs(d$logp2))), 0, a[ipc])
  
  if (max(d$logp) > max(abs(d$logp2))){
    d$logp2 = d$logp2 * max(d$logp)/max(abs(d$logp2))
  } else {
    d$logp = d$logp * max(abs(d$logp2))/max(d$logp)
  }
  
  d$pos=NA
  ticks=NULL
  lastbase=0
  colors <- rep(colors,max(d$CHR))[1:max(d$CHR)]
  ymax = ceiling(max(d$logp))
  ymin = floor(min(d$logp2))
  
  numchroms=length(unique(d$CHR))
  
  #change the position
  Uchr=unique(d$CHR)
  # Uchr=sort(Uchr)
  for (i in 1:length(Uchr)) {
    if (i==1) {
      d[d$CHR==Uchr[i], ]$pos=d[d$CHR==Uchr[i], ]$BP
    } else {
      lastbase=lastbase+tail(subset(d,CHR==Uchr[i-1])$BP, 1)
      d[d$CHR==Uchr[i], ]$pos=d[d$CHR==Uchr[i], ]$BP+lastbase
    }
    ticks=c(ticks, d[d$CHR==Uchr[i], ]$pos[floor(length(d[d$CHR==Uchr[i], ]$pos)/2)+1])
  }
  
  with(d, plot(main=title, pos, logp, ylim=c(ymin, ymax), ylab=expression(paste(paste(italic("F")[italic("ST")])," ",-log[10](italic(p)))), xlab="Chromosome", axes=F, type="n" ))
  axis(1, at=ticks, labels = Uchr)
  axis(2, at=c(-1*ceiling(max(abs(d$logp2))), 0, ceiling(max(abs(d$logp)))), lab=ytick,ylab='n')
  icol=1
  for (i in 1:length(Uchr)) {
    with(d[d$CHR==Uchr[i], ],points(pos, logp, col=colors[icol],cex=1,pch=16, bty="l"))
    with(d[d$CHR==Uchr[i], ],points(pos, logp2, col=colors[icol],cex=1,pch=16, bty="l"))
    icol=icol+1
  }
  
  #添加虚线   
  abline(h=-log10(0.05/nrow(d)), col="gray",lty = 2)
  abline(h=log10(0.05/nrow(d)), col="gray",lty = 2)
  abline(h=0, col="white", lwd=2, lty=2)
  
  
  title(paste("",ipc,sep = ""),adj=0,line = 0)
}
mtext("CML333-GWAS",cex = 3, side = 3, line = 0,outer = TRUE)
dev.off()

结果大概是这个样子：

继续补充。