Created
October 14, 2015 08:45
-
-
Save chr1swallace/49dd3e056ba84a9f8f93 to your computer and use it in GitHub Desktop.
Create a Manhattan plot of T1D GWAS statistics downloaded from T1DBase
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
## download https://www.t1dbase.org/downloads/protected_data/GWAS_Data/hg19_gwas_t1d_barrett_4_18_0.gz | |
## you will need to create a login! | |
## Then... | |
x <- read.table("hg19_gwas_t1d_barrett_4_18_0.gz",sep="\t",skip=1,as.is=TRUE) | |
colnames(x) <- c("snp","Chr","Pos","PValue","OR","kk") | |
## libraries | |
library(ggplot2) | |
library(cowplot) # nicer ggplot theme | |
x$Chr <- as.numeric(x$Chr) | |
x <- x[order(x$Chr,x$Pos),] | |
x$newchr <- c(0,diff(x$chr)) | |
x$dPos <- c(x$Pos[1],diff(x$Pos)) | |
x$dPos[ x$dPos < 0 ] <- 1000 | |
x$cpos <- cumsum(x$dPos) | |
x$logp <- pmin(-log10(x$PValue),20) | |
## one way to get a two tone image is to substitute *Chr* in the line | |
## marked *** below with *odd* | |
x$odd <- as.numeric(x$Chr) %% 2 | |
x$odd <- as.character(x$odd) | |
x$chrlabel <- ifelse(c(1,diff(x$Chr))==0,"",x$Chr) | |
x$chrx <- ifelse(x$chrlabel=="",NA,x$cpos) | |
g <- ggplot(x,aes(x=cpos,y=logp,col=as.factor(Chr))) + # *** | |
geom_point(size=1) + | |
scale_x_continuous("Chromosome",breaks=x$cpos[x$chrlabel!=""], | |
labels=x$chrlabel[x$chrlabel!=""]) + | |
theme(legend.position="none",axis.text.x=element_text(angle=90), | |
text=element_text(size=16)) + | |
xlab("") + ylab("Evidence for association") | |
g | |
ggsave(g, file="T1DManhattan-JDRF.jpg",height=4.67,width=12.8) | |
message("DONE") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment