Embedded Files
Location on server: /space/chen-syn01/1/data/cinliu/plots/manhattan_plots/genomicSEM_manhattan
/space/chen-syn01/1/data/cinliu/plots/manhattan_plots/genomicSEM_manhattan/genomicSEM_manhattan
/space/chen-syn01/1/data/cinliu/plots/manhattan_plots/genomicSEM_manhattan/genomicSEM_manhattan.zip
Initial instructions:
Initial instructions:
@Nini, could you make four manhattan plots using four files in this directory
/space/chen-syn01/1/data/chen/gwas_gclust/genomicSEM/results
The column name with "Pval_Estimate" is used for y-axis (logP) in the manhattan plots.
Goal plot:
genomicSEM_qqma_cmedit.R
#genomicSEM manhattan plot rm(list=ls())setwd("/Users/nini/Desktop/2021lab/CH/genomicSEM_manhattan/")#files from server /space/chen-syn01/1/data/chen/gwas_gclust/genomicSEM/results#The column name with "Pval_Estimate" is used for y-axis (logP) in the manhattan plots.
library(qqman)#FileNames <- list.files("/Users/nini/Desktop/2021lab/CH/genomicSEM_manhattan/results/")areaFile <- c("F1_area_sumstats.txt","F2_area_sumstats.txt")
#####area#####for (i in 1:2) {sumstats <- read.csv(paste0("~/Desktop/2021lab/CH/genomicSEM_manhattan/results/",areaFile[i]), sep="")#F1_thick_sumstats <- read.csv("~/Desktop/2021lab/CH/genomicSEM_manhattan/results/F1_thick_sumstats.txt", sep="")#F2_area_sumstats <- read.csv("~/Desktop/2021lab/CH/genomicSEM_manhattan/results/F2_area_sumstats.txt", sep="")#F2_thick_sumstats <- read.csv("~/Desktop/2021lab/CH/genomicSEM_manhattan/results/F2_thick_sumstats.txt", sep="")sumstats$CHR <- as.numeric(sumstats$CHR)sumstats$BP <- as.numeric(sumstats$BP)sumstats$Pval_Estimate <- as.numeric(sumstats$Pval_Estimate)sumstats$Z_Estimate <- as.numeric(sumstats$Z_Estimate)
#CM EDIT: keep only rows with error = 0, this should filter out about 1000 SNPsdf<-sumstats[sumstats$error==0,]summary(-log10(df$Z_Estimate)) #max is 55, 13,33,14#summary(df)#####1. png(paste0(areaFile[i], "_60max.png"),width=1600,height=800) #raiot is 2.1par(mar=c(4,8,2,0))manhattan(df, ylim=c(0,60), chr = "CHR", bp = "BP", snp = "SNP",p = "Z_Estimate", suggestiveline =F,genomewideline = F, cex=3,cex.axis=3, ylab= "", xlab = "",tck=0)abline(h=-log10(5e-8), col="blue3", lwd=10, lty=2)title(ylab="-log10(p-value)", cex.lab=3 , mgp = c(5, 3, 0))dev.off()}
####thickness #####Thick_file <- c("F1_thick_sumstats.txt","F2_thick_sumstats.txt")for (i in 1:2) { sumstats <- read.csv(paste0("~/Desktop/2021lab/CH/genomicSEM_manhattan/results/",Thick_file[i]), sep="") #F1_thick_sumstats <- read.csv("~/Desktop/2021lab/CH/genomicSEM_manhattan/results/F1_thick_sumstats.txt", sep="") #F2_area_sumstats <- read.csv("~/Desktop/2021lab/CH/genomicSEM_manhattan/results/F2_area_sumstats.txt", sep="") #F2_thick_sumstats <- read.csv("~/Desktop/2021lab/CH/genomicSEM_manhattan/results/F2_thick_sumstats.txt", sep="") sumstats$CHR <- as.numeric(sumstats$CHR) sumstats$BP <- as.numeric(sumstats$BP) sumstats$Pval_Estimate <- as.numeric(sumstats$Pval_Estimate) sumstats$Z_Estimate <- as.numeric(sumstats$Z_Estimate) #CM EDIT: keep only rows with error = 0, this should filter out about 1000 SNPs df<-sumstats[sumstats$error==0,] summary(-log10(df$Z_Estimate)) #max is 55, 13,33,14 #summary(df) #####1. png(paste0(Thick_file[i], "_20max.png"),width=1600,height=800) #raiot is 2.1 par(mar=c(4,8,2,0)) manhattan(df, ylim=c(0,20), chr = "CHR", bp = "BP", snp = "SNP",p = "Z_Estimate", suggestiveline =F,genomewideline = F, cex=3,cex.axis=3, ylab= "", xlab = "",tck=0) abline(h=-log10(5e-8), col="blue3", lwd=10, lty=2) title(ylab="-log10(p-value)", cex.lab=3 , mgp = c(5, 3, 0)) dev.off() }
library(qqman)#FileNames <- list.files("/Users/nini/Desktop/2021lab/CH/genomicSEM_manhattan/results/")areaFile <- c("F1_area_sumstats.txt","F2_area_sumstats.txt")
#####area#####for (i in 1:2) {sumstats <- read.csv(paste0("~/Desktop/2021lab/CH/genomicSEM_manhattan/results/",areaFile[i]), sep="")#F1_thick_sumstats <- read.csv("~/Desktop/2021lab/CH/genomicSEM_manhattan/results/F1_thick_sumstats.txt", sep="")#F2_area_sumstats <- read.csv("~/Desktop/2021lab/CH/genomicSEM_manhattan/results/F2_area_sumstats.txt", sep="")#F2_thick_sumstats <- read.csv("~/Desktop/2021lab/CH/genomicSEM_manhattan/results/F2_thick_sumstats.txt", sep="")sumstats$CHR <- as.numeric(sumstats$CHR)sumstats$BP <- as.numeric(sumstats$BP)sumstats$Pval_Estimate <- as.numeric(sumstats$Pval_Estimate)sumstats$Z_Estimate <- as.numeric(sumstats$Z_Estimate)
#CM EDIT: keep only rows with error = 0, this should filter out about 1000 SNPsdf<-sumstats[sumstats$error==0,]summary(-log10(df$Z_Estimate)) #max is 55, 13,33,14#summary(df)#####1. png(paste0(areaFile[i], "_60max.png"),width=1600,height=800) #raiot is 2.1par(mar=c(4,8,2,0))manhattan(df, ylim=c(0,60), chr = "CHR", bp = "BP", snp = "SNP",p = "Z_Estimate", suggestiveline =F,genomewideline = F, cex=3,cex.axis=3, ylab= "", xlab = "",tck=0)abline(h=-log10(5e-8), col="blue3", lwd=10, lty=2)title(ylab="-log10(p-value)", cex.lab=3 , mgp = c(5, 3, 0))dev.off()}
####thickness #####Thick_file <- c("F1_thick_sumstats.txt","F2_thick_sumstats.txt")for (i in 1:2) { sumstats <- read.csv(paste0("~/Desktop/2021lab/CH/genomicSEM_manhattan/results/",Thick_file[i]), sep="") #F1_thick_sumstats <- read.csv("~/Desktop/2021lab/CH/genomicSEM_manhattan/results/F1_thick_sumstats.txt", sep="") #F2_area_sumstats <- read.csv("~/Desktop/2021lab/CH/genomicSEM_manhattan/results/F2_area_sumstats.txt", sep="") #F2_thick_sumstats <- read.csv("~/Desktop/2021lab/CH/genomicSEM_manhattan/results/F2_thick_sumstats.txt", sep="") sumstats$CHR <- as.numeric(sumstats$CHR) sumstats$BP <- as.numeric(sumstats$BP) sumstats$Pval_Estimate <- as.numeric(sumstats$Pval_Estimate) sumstats$Z_Estimate <- as.numeric(sumstats$Z_Estimate) #CM EDIT: keep only rows with error = 0, this should filter out about 1000 SNPs df<-sumstats[sumstats$error==0,] summary(-log10(df$Z_Estimate)) #max is 55, 13,33,14 #summary(df) #####1. png(paste0(Thick_file[i], "_20max.png"),width=1600,height=800) #raiot is 2.1 par(mar=c(4,8,2,0)) manhattan(df, ylim=c(0,20), chr = "CHR", bp = "BP", snp = "SNP",p = "Z_Estimate", suggestiveline =F,genomewideline = F, cex=3,cex.axis=3, ylab= "", xlab = "",tck=0) abline(h=-log10(5e-8), col="blue3", lwd=10, lty=2) title(ylab="-log10(p-value)", cex.lab=3 , mgp = c(5, 3, 0)) dev.off() }
The 4 plots are placed together using powerpoint. Power point file can be found in the server
/space/chen-syn01/1/data/cinliu/plots/manhattan_plots/genomicSEM_manhattan/genomicSEM_manhattan
/space/chen-syn01/1/data/cinliu/plots/manhattan_plots/genomicSEM_manhattan/genomicSEM_manhattan.zip
Page updated
Report abuse