Manhattan plots

email title: "Manhattan plots"

First email: May19th

Location on server: /space/chen-syn01/1/data/cinliu/plots/manhattan_plots/qqman

Initial instructions:

Here are R scripts for Manhattan and Miami plots (eg, male vs female GWAS) for GWAS result visualization.

Dennis uses the data here (/space/chen-syn01/1/data/cmakowski/GWASgclus_UKB/Mar2020_fastGWA/phenotypes/main_analysis/sumstats/*fastGWA).

I use scp to fetch data on the server to my local laptop

scp chc101@ip99.ucsd.edu:/space/chen-syn01/1/data/chen/gwas_gclust/UKB/selection/analysis/summary.txt local_filepath

@Nini, we may ask your help adjust the Manhanttan plots later so please take a look at the code.

Here are an initial R script for plotting and three plots generated by the script. It'd be great to have your help to combine the three plots side-by-side into one single plot.

Initial files attached: 8 files

Initial code of runManhattanPlot.R:

made using ggplot

#! /usr/bin/env Rscript
rm(list=ls())library(data.table)library(ggplot2)
Path <- "/Users/dennis/"setwd(Path)
######################################### # point towards sumstats and set params ##########################################
# ls *3103*area*fastGWA > ChenSumstatsArea.txt# ls *3103*thick*fastGWA > ChenSumstatsThick.txt
# FileNames <- fread("/Users/dennis/ChenSumstatsArea.txt",data.table=F,header = F)# FileNames$V2 <- gsub("UKB39k_Chen_GWAS_DM_310320_ChenA_bh_","",FileNames$V1)# Sets <- gsub("_area_sumstat.fastGWA","",FileNames$V2)
FileNames <- fread("/Users/dennis/ChenSumstatsThick.txt",data.table=F,header = F)FileNames$V2 <- gsub("UKB39k_Chen_GWAS_DM_310320_ChenT_bh_","",FileNames$V1)Sets <- gsub("_thickness_sumstat.fastGWA","",FileNames$V2)
for(i in 1:length(Sets)){#dat <- fread(paste0("/Users/dennis/UKB39k_Chen_GWAS_DM_310320_ChenA_bh_",Sets[i],"_area_sumstat.fastGWA"),data.table=F)[,c("CHR","POS","P")]dat <- fread(paste0("/Users/dennis/UKB39k_Chen_GWAS_DM_310320_ChenT_bh_",Sets[i],"_thickness_sumstat.fastGWA"),data.table=F)[,c("CHR","POS","P")]
cols <- c("grey72","grey42") # set colors for the alternating chromosomespcutoff <- 0.2 # cut away SNPs with P above this value, removing most SNPs, so it runs faster. ymax=63 # Same y-axis range for all plots (-log10 of most sign p-val)
################ Prep inputs ################
colnames(dat) <- c("chr","pos","p")dat$p[dat$p==0] <- NAdat <- dat[complete.cases(dat),]
###### Create df with chr positionschrmaxpos <- data.frame("chr"=1:22,"maxpos"=NA,"genomestartpos"=NA)for(j in 1:22){ chrmaxpos$maxpos[j] <- max(c(dat$pos[dat$chr==j])) if(j==1){chrmaxpos$genomestartpos[j] <- 0 } else { chrmaxpos$genomestartpos[j] <- chrmaxpos$genomestartpos[j-1] + chrmaxpos$maxpos[j-1]} }
## label position chrmaxpos$labpos <- chrmaxpos$genomestartpos + 0.5*chrmaxpos$maxpos
####### Add genomic position to sumstatsdat$gpos <- NAfor(k in 1:22){dat$gpos[dat$chr==k] <- dat$pos[dat$chr==k]+chrmaxpos$genomestartpos[k]}
####### set graph parameters#maxp <- -log10(min(dat$p))plotdat <- dat[dat$p<pcutoff,]col <- rep(cols,22)[plotdat$chr]
######## ggplotManplot <- ggplot(plotdat, aes(x=gpos, y=-log10(p))) + theme_bw() + geom_hline(yintercept = -log10(5e-8),color="navyblue",linetype="dashed",size=1) + geom_hline(yintercept = 0,color="black",linetype="solid",size=0.2) + geom_segment(x = -143500000,xend=-143500000,y=0,yend=ymax,color="black",linetype="solid",size=0.2) + geom_point(col=col) + scale_y_continuous(name="-log10(p-value) \n",limits = c(-2,ymax)) + geom_text(data=chrmaxpos,aes(x=labpos,y=-1),label=as.character(1:22),size=4,inherit.aes = F) + #ggtitle(Sets[i]) + theme(plot.title=element_text(size=28,hjust = 0.5,margin = margin(b = -20)), axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.title.y = element_text(size = 24), axis.text.y = element_text(size = 22), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), plot.background = element_blank(), panel.border = element_blank(), panel.background = element_rect(colour = NA, fill = 'transparent')) #ggsave(paste0(Path,"Dropbox/Chen/Plots/Manhattan_",Sets[i],"_Area.png"), plot = Manplot, width = 15, height = 7,bg="transparent")ggsave(paste0(Path,"Dropbox/Chen/Plots/Manhattan_",Sets[i],"_Thick.png"), plot = Manplot, width = 15, height = 7,bg="transparent")
}

Second email: Jun 17-Jun 30

email title: Manhattan plots for GWAS project

Instructions:

We want to update five of the 24 regions to have a "break" in the y-axis.

I have attached:

1) R file: manhattan plot script

2) text files: that include names of area and thickness regions that is called in the area and thickness loop, respectively

3) png file: an example of one output plot (with a y-max of 60).

Specifically, we can re-run to set ymax to 30 for all plots (I've already updated that in the script). You can test the script out as is to make sure that you're able to generate the plots.

Then for five area regions, we would need to create a break in the y-axis manually because they have a log-pval greater than 30 for the following five regions:

parsopercularis_area

inferiorparietal_area

occipital_area

precuneus_area

motorpremotor_area

Summary statistics are in: /space/chen-syn01/1/data/cmakowski/GWASgclus_UKB/Mar2020_fastGWA/phenotypes/main_analysis/sumstats/*fastGWA

What I did:

Plots had to be remade using qqman instead of ggplot because we need axis breaks

Scrip rewritten

qqman.R

bold parts are the major edits

#! /usr/bin/env Rscript
rm(list=ls())library(data.table)library(ggplot2)library(plotrix)library(qqman)
setwd("/Users/niniliu/Desktop/2020lab/ukb/")
######################################### # point towards sumstats and set params ###########################################run separately for area, then thickness
#####1. AREAFileNames <- fread("area_mmil_features.txt",data.table=F,header = F)
for(i in 1:12){ dat <- fread(paste0("UKB39k_Chen_GWAS_DM_310320_",FileNames[i,],"_sumstat.fastGWA"),data.table=F)[,c("CHR","POS","P")]
cols <- c("grey72","grey42") # set colors for the alternating chromosomespcutoff <- 0.2 # cut away SNPs with P above this value, removing most SNPs, so it runs faster.
#------CM COMMENT: could also test just having ymax=30 for all plots with the exception of the select few regions that have very sig hitsymax=30 # Same y-axis range for all plots (-log10 of most sign p-val), then will need to re-run plots for five area regions with log pvals above 30 (parsopercularis, inferior parietal, occipital, precuneus, motorpremotor)
################ Prep inputs ################
colnames(dat) <- c("chr","pos","p")dat$p[dat$p==0] <- NAdat <- dat[complete.cases(dat),]
###### Create df with chr positionschrmaxpos <- data.frame("chr"=1:22,"maxpos"=NA,"genomestartpos"=NA)for(j in 1:22){ chrmaxpos$maxpos[j] <- max(c(dat$pos[dat$chr==j])) if(j==1){chrmaxpos$genomestartpos[j] <- 0 } else { chrmaxpos$genomestartpos[j] <- chrmaxpos$genomestartpos[j-1] + chrmaxpos$maxpos[j-1]} }
## label position chrmaxpos$labpos <- chrmaxpos$genomestartpos + 0.5*chrmaxpos$maxpos
####### Add genomic position to sumstatsdat$gpos <- NAfor(k in 1:22){dat$gpos[dat$chr==k] <- dat$pos[dat$chr==k]+chrmaxpos$genomestartpos[k]}
####### set graph parameters#maxp <- -log10(min(dat$p))plotdat <- dat[dat$p<pcutoff,]col <- rep(cols,22)[plotdat$chr]
######## ggplot
png(paste0(FileNames[i,], "_example_plot.png"),width=3000,height=1430) #raiot is 2.1par(mar=c(6,10,2,0))#par(mar=c(b,l,t,r))manhattan(plotdat, chr = "chr", bp = "pos", p = "p", snp = "gpos", ylim= c(0,30), cex.axis=0.1,suggestiveline =F,genomewideline = F, cex=4, xlab="", ylab="", tck=0)abline(h=-log10(5e-8), col="blue3", lwd=10, lty=2)#abline(h=30, col="blue3", lwd=10)axis(2,at=c(0,10,20,30),labels=c(0,10, 20, 30),las = 2, cex.axis=5, lwd=4) axis(1, at= chrmaxpos$labpos, labels=1:22, pos = 0, cex.axis=1.9, font=2)title(ylab="-log10(p-value)", cex.lab=4 , mgp = c(7, 3, 0))#mtext("xlab", side=1, cex = 4, line=6)
dev.off()

###change where you save output
}

#####2. THICKNESSFileNames <- fread("thickness_mmil_features.txt",data.table=F,header = F)
for(i in 1:12){ dat <- fread(paste0("UKB39k_Chen_GWAS_DM_310320_",FileNames[i,],"_sumstat.fastGWA"),data.table=F)[,c("CHR","POS","P")] cols <- c("grey72","grey42") # set colors for the alternating chromosomes pcutoff <- 0.2 # cut away SNPs with P above this value, removing most SNPs, so it runs faster. #------CM COMMENT: could also test just having ymax=30 for all plots with the exception of the select few regions that have very sig hits ymax=30 # Same y-axis range for all plots (-log10 of most sign p-val) ############### # Prep inputs # ############### colnames(dat) <- c("chr","pos","p") dat$p[dat$p==0] <- NA dat <- dat[complete.cases(dat),] ###### Create df with chr positions chrmaxpos <- data.frame("chr"=1:22,"maxpos"=NA,"genomestartpos"=NA) for(j in 1:22){ chrmaxpos$maxpos[j] <- max(c(dat$pos[dat$chr==j])) if(j==1){chrmaxpos$genomestartpos[j] <- 0 } else { chrmaxpos$genomestartpos[j] <- chrmaxpos$genomestartpos[j-1] + chrmaxpos$maxpos[j-1] } } ## label position chrmaxpos$labpos <- chrmaxpos$genomestartpos + 0.5*chrmaxpos$maxpos ####### Add genomic position to sumstats dat$gpos <- NA for(k in 1:22){dat$gpos[dat$chr==k] <- dat$pos[dat$chr==k]+chrmaxpos$genomestartpos[k]} ####### set graph parameters #maxp <- -log10(min(dat$p)) plotdat <- dat[dat$p<pcutoff,] col <- rep(cols,22)[plotdat$chr] ######## ggplot png(paste0(FileNames[i,], "_example_plot.png"),width=3000,height=1430) #raiot is 2.1 par(mar=c(6,10,2,0)) #par(mar=c(b,l,t,r)) manhattan(plotdat, chr = "chr", bp = "pos", p = "p", snp = "gpos", ylim= c(0,30), cex.axis=0.1,suggestiveline =F,genomewideline = F, cex=4, xlab="", ylab="", tck=0) abline(h=-log10(5e-8), col="blue3", lwd=10, lty=2) #abline(h=30, col="blue3", lwd=10) axis(2,at=c(0,10,20,30),labels=c(0,10, 20, 30),las = 2, cex.axis=5, lwd=4) axis(1, at= chrmaxpos$labpos, labels=1:22, pos = 0, cex.axis=1.9, font=2) title(ylab="-log10(p-value)", cex.lab=4 , mgp = c(7, 3, 0)) #mtext("xlab", side=1, cex = 4, line=6) dev.off()}

Plots with axis breaks are seperatly made. One example presented below, the rest is attached in the google drive folder.

#! /usr/bin/env Rscriptrm(list=ls())library(data.table)library(qqman) library(plotrix)
setwd('/Users/niniliu/Desktop/2020lab/ukb/') #update
######################################### # point towards sumstats and set params ###########################################run separately for area, then thickness
#####1. AREAdat <- fread(paste0("UKB39k_Chen_GWAS_DM_310320_ChenA_bh_inferiorparietal_area_sumstat.fastGWA"),data.table=F)[,c("CHR","POS","P")] cols <- c("grey72","grey42") # set colors for the alternating chromosomes pcutoff <- 0.2 # cut away SNPs with P above this value, removing most SNPs, so it runs faster. #------CM COMMENT: could also test just having ymax=30 for all plots with the exception of the select few regions that have very sig hits ymax=30 # Same y-axis range for all plots (-log10 of most sign p-val), then will need to re-run plots for five area regions with log pvals above 30 (parsopercularis, inferior parietal, occipital, precuneus, motorpremotor) ############### # Prep inputs # ############### colnames(dat) <- c("chr","pos","p") dat$p[dat$p==0] <- NA dat <- dat[complete.cases(dat),] ###### Create df with chr positions chrmaxpos <- data.frame("chr"=1:22,"maxpos"=NA,"genomestartpos"=NA) for(j in 1:22){ chrmaxpos$maxpos[j] <- max(c(dat$pos[dat$chr==j])) if(j==1){chrmaxpos$genomestartpos[j] <- 0 } else { chrmaxpos$genomestartpos[j] <- chrmaxpos$genomestartpos[j-1] + chrmaxpos$maxpos[j-1] } } ## label position chrmaxpos$labpos <- chrmaxpos$genomestartpos + 0.5*chrmaxpos$maxpos ####### Add genomic position to sumstats dat$gpos <- NA for(k in 1:22){dat$gpos[dat$chr==k] <- dat$pos[dat$chr==k]+chrmaxpos$genomestartpos[k]} ####### set graph parameters #maxp <- -log10(min(dat$p)) plotdat <- dat[dat$p<pcutoff,]#add neg logpval to a line for double checcking (can delete this line if you want) plotdat$neglogPval <- -log10(plotdat$p) col <- rep(cols,22)[plotdat$chr]
##only run once, offsetingidx <- which((-log10(plotdat$p)) > 30) view_dat_above30 <- plotdat[idx,] #viewing the data mover 30 for double chcking view_dat_above30#viewing the data above 30 before offesting for double chcking
###CUSTOM PLOT inferioparietal########only run once, offseting parsopercularis_areaplotdat[-log10(plotdat$p)>31,"p"] = plotdat[-log10(plotdat$p)>31,'p']/1e-3 #offest over the 40 rangeplotdat[-log10(plotdat$p)>40,"p"] = plotdat[-log10(plotdat$p)>40,'p']/1e-19 #offset over 50 range
png(file="ChenA_bh_inferiorparietal_area_example_plot.png",width=3000,height=1430) #raiot is 2.1par(mar=c(6,10,2,0)) #reset if change the size #par(mar=c(b,l,t,r))manhattan(plotdat, chr = "chr", bp = "pos", p = "p", snp = "gpos", ylim= c(0,37), cex.axis=0.1,suggestiveline =F,genomewideline = F, cex=4, xlab="", ylab="", tck=0,lwd=6)axis(2,at=c(0,10,20,30,34,37),labels=c(0,10, 20, 30,40,60),las = 2, cex.axis=5, lwd=6)axis(1, at= chrmaxpos$labpos, labels=1:22, pos = 0, cex.axis=1.9, font=2)title(ylab="-log10(p-value)", cex.lab=4 , mgp = c(7, 3, 0))abline(h=-log10(5e-8), col="blue3", lwd=10, lty=2)axis.break(2,31,style="slash", brw=0.01)axis.break(2,35,style="slash", brw=0.01)
dev.off()

Third email: Jul 15-Jul 17

email title: Manhattan plot adjustments

Additional adjustments made in terms of color and font size

The plots look something like this:

Full links to plots:

runManhattanPlot_cm1_qqman-1.R

The general plot:

#! /usr/bin/env Rscript
rm(list=ls())library(data.table)library(ggplot2)library(plotrix)library(qqman)
setwd("/Users/niniliu/Desktop/2020lab/ukb/")
######################################### # point towards sumstats and set params ###########################################run separately for area, then thickness
#####1. AREAFileNames <- fread("area_mmil_features.txt",data.table=F,header = F)
for(i in 1:12){ dat <- fread(paste0("UKB39k_Chen_GWAS_DM_310320_",FileNames[i,],"_sumstat.fastGWA"),data.table=F)[,c("CHR","POS","P")]
cols <- c("grey72","grey42") # set colors for the alternating chromosomespcutoff <- 0.2 # cut away SNPs with P above this value, removing most SNPs, so it runs faster.
#------CM COMMENT: could also test just having ymax=30 for all plots with the exception of the select few regions that have very sig hitsymax=30 # Same y-axis range for all plots (-log10 of most sign p-val), then will need to re-run plots for five area regions with log pvals above 30 (parsopercularis, inferior parietal, occipital, precuneus, motorpremotor)
################ Prep inputs ################
colnames(dat) <- c("chr","pos","p")dat$p[dat$p==0] <- NAdat <- dat[complete.cases(dat),]
###### Create df with chr positionschrmaxpos <- data.frame("chr"=1:22,"maxpos"=NA,"genomestartpos"=NA)for(j in 1:22){ chrmaxpos$maxpos[j] <- max(c(dat$pos[dat$chr==j])) if(j==1){chrmaxpos$genomestartpos[j] <- 0 } else { chrmaxpos$genomestartpos[j] <- chrmaxpos$genomestartpos[j-1] + chrmaxpos$maxpos[j-1]} }
## label position chrmaxpos$labpos <- chrmaxpos$genomestartpos + 0.5*chrmaxpos$maxpos
####### Add genomic position to sumstatsdat$gpos <- NAfor(k in 1:22){dat$gpos[dat$chr==k] <- dat$pos[dat$chr==k]+chrmaxpos$genomestartpos[k]}
####### set graph parameters#maxp <- -log10(min(dat$p))plotdat <- dat[dat$p<pcutoff,]col <- rep(cols,22)[plotdat$chr]
######## ggplot
png(paste0(FileNames[i,], "_example_plot.png"),width=3000,height=1430) #raiot is 2.1par(mar=c(6,10,2,0))#par(mar=c(b,l,t,r))manhattan(plotdat, chr = "chr", bp = "pos", p = "p", snp = "gpos", ylim= c(0,30), cex.axis=0.1,suggestiveline =F,genomewideline = F, cex=4, xlab="", ylab="", tck=0)abline(h=-log10(5e-8), col="blue3", lwd=10, lty=2)#abline(h=30, col="blue3", lwd=10)axis(2,at=c(0,10,20,30),labels=c(0,10, 20, 30),las = 2, cex.axis=5, lwd=4) axis(1, at= chrmaxpos$labpos, labels=1:22, pos = 0, cex.axis=1.9, font=2)title(ylab="-log10(p-value)", cex.lab=4 , mgp = c(7, 3, 0))#mtext("xlab", side=1, cex = 4, line=6)
dev.off()

###change where you save output
}

#####2. THICKNESSFileNames <- fread("thickness_mmil_features.txt",data.table=F,header = F)
for(i in 1:12){ dat <- fread(paste0("UKB39k_Chen_GWAS_DM_310320_",FileNames[i,],"_sumstat.fastGWA"),data.table=F)[,c("CHR","POS","P")] cols <- c("grey72","grey42") # set colors for the alternating chromosomes pcutoff <- 0.2 # cut away SNPs with P above this value, removing most SNPs, so it runs faster. #------CM COMMENT: could also test just having ymax=30 for all plots with the exception of the select few regions that have very sig hits ymax=30 # Same y-axis range for all plots (-log10 of most sign p-val) ############### # Prep inputs # ############### colnames(dat) <- c("chr","pos","p") dat$p[dat$p==0] <- NA dat <- dat[complete.cases(dat),] ###### Create df with chr positions chrmaxpos <- data.frame("chr"=1:22,"maxpos"=NA,"genomestartpos"=NA) for(j in 1:22){ chrmaxpos$maxpos[j] <- max(c(dat$pos[dat$chr==j])) if(j==1){chrmaxpos$genomestartpos[j] <- 0 } else { chrmaxpos$genomestartpos[j] <- chrmaxpos$genomestartpos[j-1] + chrmaxpos$maxpos[j-1] } } ## label position chrmaxpos$labpos <- chrmaxpos$genomestartpos + 0.5*chrmaxpos$maxpos ####### Add genomic position to sumstats dat$gpos <- NA for(k in 1:22){dat$gpos[dat$chr==k] <- dat$pos[dat$chr==k]+chrmaxpos$genomestartpos[k]} ####### set graph parameters #maxp <- -log10(min(dat$p)) plotdat <- dat[dat$p<pcutoff,] col <- rep(cols,22)[plotdat$chr] ######## ggplot png(paste0(FileNames[i,], "_example_plot.png"),width=3000,height=1430) #raiot is 2.1 par(mar=c(6,10,2,0)) #par(mar=c(b,l,t,r)) manhattan(plotdat, chr = "chr", bp = "pos", p = "p", snp = "gpos", ylim= c(0,30), cex.axis=0.1,suggestiveline =F,genomewideline = F, cex=4, xlab="", ylab="", tck=0) abline(h=-log10(5e-8), col="blue3", lwd=10, lty=2) #abline(h=30, col="blue3", lwd=10) axis(2,at=c(0,10,20,30),labels=c(0,10, 20, 30),las = 2, cex.axis=5, lwd=4) axis(1, at= chrmaxpos$labpos, labels=1:22, pos = 0, cex.axis=1.9, font=2) title(ylab="-log10(p-value)", cex.lab=4 , mgp = c(7, 3, 0)) #mtext("xlab", side=1, cex = 4, line=6) dev.off()}

motorpremotor.R

Sample axis break code

#! /usr/bin/env Rscriptrm(list=ls())library(data.table)library(qqman) library(plotrix)
setwd('/Users/niniliu/Desktop/2020lab/ukb/')
######################################### # point towards sumstats and set params ###########################################run separately for area, then thickness
#####1. AREAdat <- fread(paste0("UKB39k_Chen_GWAS_DM_310320_ChenA_bh_motorpremotor_area_sumstat.fastGWA"),data.table=F)[,c("CHR","POS","P")] cols <- c("grey72","grey42") # set colors for the alternating chromosomes pcutoff <- 0.2 # cut away SNPs with P above this value, removing most SNPs, so it runs faster. #------CM COMMENT: could also test just having ymax=30 for all plots with the exception of the select few regions that have very sig hits ymax=30 # Same y-axis range for all plots (-log10 of most sign p-val), then will need to re-run plots for five area regions with log pvals above 30 (parsopercularis, inferior parietal, occipital, precuneus, motorpremotor) ############### # Prep inputs # ############### colnames(dat) <- c("chr","pos","p") dat$p[dat$p==0] <- NA dat <- dat[complete.cases(dat),] ###### Create df with chr positions chrmaxpos <- data.frame("chr"=1:22,"maxpos"=NA,"genomestartpos"=NA) for(j in 1:22){ chrmaxpos$maxpos[j] <- max(c(dat$pos[dat$chr==j])) if(j==1){chrmaxpos$genomestartpos[j] <- 0 } else { chrmaxpos$genomestartpos[j] <- chrmaxpos$genomestartpos[j-1] + chrmaxpos$maxpos[j-1] } } ## label position chrmaxpos$labpos <- chrmaxpos$genomestartpos + 0.5*chrmaxpos$maxpos ####### Add genomic position to sumstats dat$gpos <- NA for(k in 1:22){dat$gpos[dat$chr==k] <- dat$pos[dat$chr==k]+chrmaxpos$genomestartpos[k]} ####### set graph parameters #maxp <- -log10(min(dat$p)) plotdat <- dat[dat$p<pcutoff,]
plotdat$neglogPval <- -log10(plotdat$p) #add neg logpval to a line (for double checking, can remove this line if desired) col <- rep(cols,22)[plotdat$chr]
##only run once, offsetingidx <- which((-log10(plotdat$p)) > 30) view_dat_above30 <- plotdat[idx,] #viewing the data mover 30 for double chcking view_dat_above30#viewing the data above 30 before offesting for double chcking
###CUSTOM PLOT motorpremotor_area########only run once, offseting parsopercularis_areaplotdat[-log10(plotdat$p)>31,"p"] = plotdat[-log10(plotdat$p)>31,'p']/1e-3 #offest over the 40 range
png(file="ChenA_bh_motorpremotor_area_example_plot.png",width=3000,height=1430) #raiot is 2.1par(mar=c(6,10,2,0)) #reset if change the size #par(mar=c(b,l,t,r))manhattan(plotdat, chr = "chr", bp = "pos", p = "p", snp = "gpos", ylim= c(0,35), cex.axis=0.1,suggestiveline =F,genomewideline = F, cex=4, xlab="", ylab="", tck=0)abline(h=-log10(5e-8), col="blue3", lwd=10, lty=2)axis(2,at=c(0,10,20,30,35),labels=c(0,10, 20, 30,40),las = 2, cex.axis=5, lwd=4) axis(1, at= chrmaxpos$labpos, labels=1:22, pos = 0, cex.axis=1.9, font=2)#mtext("xlab", side=1, cex = 4, line=6)title(ylab="-log10(p-value)", cex.lab=4 , mgp = c(7, 3, 0))axis.break(2,30.5,style="slash", brw=0.01)dev.off()

Page updated

Report abuse