[R] how to change manhattan plot code to get a different color per chromosome

Jim Lemon drj|m|emon @end|ng |rom gm@||@com
Tue Jun 16 11:00:06 CEST 2020


Hi Ana,
Your attached image seems to have bailed out before landing in the
list. Here's how to do it using a simple manhattan plot with the data
from:

https://reneshbedre.github.io//assets/posts/mhat/gwas_res_sim.csv

gwas<-read.csv("gwas_res_sim.csv",stringsAsFactors=FALSE)
# get the data into chromosome order
gwas<-gwas[order(gwas$chr,gwas$SNP),]
gwas$BP<-rep(NA,dim(gwas)[1]
# fake some base positions
for(chromosome in 1:20) {
 snporder<-order(as.numeric(gsub("rs","",gwas[gwas$chr == chromosome,"SNP"])))
 gwas$BP[gwas$chr == chromosome]<-
  as.numeric(paste(chromosome,snporder,sep="."))
}
library(plotrix)
# set the chromosome colors here - be more creative than me
chrcol<-color.scale(1:10,extremes=c("red","blue"))

# simple manhattan plot function
manhattan<-function(x,CHR="CHR",BP="BP",SNP="SNP",p="p",
 main="Manhattan Plot",xlab="Chromosome",ylab="-log10(p)",
 pch=".",cex=1,siglevel=0.00001,sigcol="green",sigcex=2,
 chrlab=NULL,chrcol=NULL,annotate=FALSE) {

 par(xaxs="i",yaxs="i")
 nchr<-length(unique(x[,CHR]))
 if(is.null(chrlab)) chrlab<-1:nchr
 ypos<--log10(x[,p])
 plot(x[,BP],ypos,ylim=c(0,max(ypos,na.rm=TRUE)*1.05),
  main=main,xlab=xlab,ylab=ylab,
  pch=pch,col=chrcol[x[,CHR]],cex=4,xaxt="n")
 abline(h=-log10(siglevel),lty=2)
 staxlab(1,at=(1:nchr) + 0.5,labels=chrlab)
 sigindx<-which(ypos >= -log10(siglevel))
 sigSNP<-unique(x[sigindx,SNP])
 cat(length(sigindx),"sig\n")
 points(x[sigindx,BP],ypos[sigindx],
  pch=pch,col=sigcol,cex=sigcex)
 if(annotate) text(x[sigindx,BP],ypos[sigindx]*1.03,x[sigindx,SNP])
 return(sigSNP=sigSNP)
}

manhattan(gwas,CHR="chr",p="pvalue",cex=4,sigcex=6,
 chrcol=chrcol,annotate=TRUE)

Jim

On Tue, Jun 16, 2020 at 8:06 AM Ana Marija <sokovic.anamarija using gmail.com> wrote:
>
> Hello,
>
> Is there is a way to set colors in this plot to look like this one in
> attach (different color for each CHR-there is 22 of them)?
>
>
> library(qqman)
> results_log <- read.table("meta_p_pos_chr.F", head=TRUE,stringsAsFactors=FALSE)
> png("META.png")
> manhattan(results_log,chr="CHR",bp="POS",p="META_pval",snp="MARKER",ylim
> = c(0, 10))
> dev.off()
>
> Thanks
> Ana
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.



More information about the R-help mailing list