MEQTLvsGEMMA.R

library(dplyr)
library(ggplot2)
library(tidyr)
library(data.table)
METS<-as.data.frame(fread("zcat ~/mets_analysis/meqtl/combined_pop/METS_FDR0.05_PC0_PF0.txt.gz", colClasses = "character"))
colnames(METS)<-c("snps", "gene", "statistic", "pvalue", "FDR", "beta")

gemma<-as.data.frame(fread("zcat ~/mets_analysis/for_gemma/gemma_whole_genome_FDR0.05.txt.gz", header = F, sep = "\t"))
colnames(gemma)<-c("chr", "rs", "ps", "n_miss", "allele1", "allele0", "af", "beta", "se", "l_remle", "l_mle", "p_wald", "p_lrt", "p_score", "gene_inf", "FDR")
info<-separate(gemma, gene_inf, into = c("gene1","gene2",NA,NA,NA,NA), sep = "\\.")
gem2<-cbind(gemma,paste(info$gene1,".",info$gene2, sep=""))
colnames(gem2)<-c("chr", "rs", "ps", "n_miss", "allele1", "allele0", "af", "beta", "se", "l_remle", "l_mle", "p_wald", "p_lrt", "p_score", "gene_inf", "FDR", "gene")
comparison<-inner_join(METS,gem2, by = c("snps"="rs"),c("gene"="gene"))
str(head(comparison))
p1 <- ggplot(data = comparison, aes(x = as.numeric(FDR.x),y=as.numeric(FDR.y)))
png("~/mets_analysis/meqtl/meqtl_gemma_FDR.png")
p1 + xlab("MEQTL FDR") + ylab("GEMMA FDR") + geom_point() +geom_smooth(method="lm") + ggtitle("meqt FDR vs gemma p_wald")
dev.off()
png("~/mets_analysis/meqtl/meqtl_gemma_FDR_log10.png")
p1 + xlab("MEQTL FDR log10") + ylab("GEMMA FDR log10") + geom_point() +geom_smooth(method="lm") + scale_x_log10() + scale_y_log10() + ggtitle("meqt FDR vs gemma FDR, log10")
dev.off()

#p2 <- ggplot(data = comparison, aes(x = "FDR",y="p_wald"))
#png("~/mets_analysis/meqtl/meqtl_gemma_fdr.png")
#p + xlab("MEQTL fdr") + ylab("GEMMA fdr") + geom_point() +geom_smooth(method="lm") + scale_x_log10() + scale_y_log10() + ggtitle("meqt fdr vs gemma fdr")
#dev.off()
#png("~/mets_analysis/meqtl/meqtl_gemma_fdr_log10.png")
#p + xlab("MEQTL pvalue") + ylab("GEMMA p_wald") + geom_point() +geom_smooth(method="lm") + scale_x_log10() + scale_y_log10() + ggtitle("meqt pvalue vs gemma p_wald, log10")
#dev.off()