# ==== QQ plots ====
library( SuppDists )
library( MASS )
library( Haplin )
## Registered S3 methods overwritten by 'ffbase':
## method from
## [.ff ff
## [.ffdf ff
## [<-.ff ff
## [<-.ffdf ff
library( ggplot2 )
library( ggrepel )
source( "gg_qqplot.R" )
# ==== CL dataset ====
# CL dataset
# ==== GxM ====
# GxM
load( "../GxM_PoOxM_new_betas_mult_CLO.RData" )
gxe.pval.promoter.cl <- t( sapply( gxe.prom.all, function(x){
return( data.frame( pval = x$gxm$pval[ 2 ], df = x$gxm$df[ 2 ] ) )
} ) )
gxe.pval.enhancer.cl <- t( sapply( gxe.enh.all, function(x){
return( data.frame( pval = x$gxm$gxe.test$pval[ 2 ], df = x$gxm$gxe.test$df[ 2 ] ) )
} ) )
gxe.pval.gene.cl <- t( sapply( gxe.gene.all, function(x){
return( data.frame( pval = x$gxm$gxe.test$pval[ 2 ], df = x$gxm$gxe.test$df[ 2 ] ) )
} ) )
# creating the data.frame
cur.snp.names <- as.character( my.SNP.map$Marker[ match( c( rownames( gxe.pval.promoter.cl ),
rownames( gxe.pval.enhancer.cl ), rownames( gxe.pval.gene.cl ) ),
my.SNP.map$gene_name ) ] )
gxe.pval.all.cl <- data.frame( region = c( rep( "promoter", nrow( gxe.pval.promoter.cl ) ),
rep( "enhancer", nrow( gxe.pval.enhancer.cl ) ), rep( "gene", nrow( gxe.pval.gene.cl ) ) ),
gene.name = c( rownames( gxe.pval.promoter.cl ), rownames( gxe.pval.enhancer.cl ),
rownames( gxe.pval.gene.cl ) ),
snp.name = cur.snp.names,
pval = unlist( c( gxe.pval.promoter.cl[ , "pval" ], gxe.pval.enhancer.cl[ , "pval" ],
gxe.pval.gene.cl[ , "pval" ] ) ),
df = unlist( c( gxe.pval.promoter.cl[ , "df" ], gxe.pval.enhancer.cl[ , "df" ],
gxe.pval.gene.cl[ , "df" ] ) ) )
head( gxe.pval.all.cl )
## region gene.name snp.name pval df
## 1 promoter PAX7 rs742071 0.06400049 2
## 2 promoter ABCA4 rs560426 0.42805451 2
## 3 promoter IRF6 rs642961 0.55340321 2
## 4 promoter THADA rs7590268 0.17320790 2
## 5 promoter 8q21.3 rs12543318 0.32332097 2
## 6 promoter 8q24 rs987525 0.42999212 2
dim( gxe.pval.all.cl )
## [1] 25 5
df.all.cl <- unique( gxe.pval.all.cl$df )
df.all.cl
## [1] 2
# plotting
pvals.GxM.CL <- gxe.pval.all.cl$pval
names( pvals.GxM.CL ) <- paste( gxe.pval.all.cl$region, gxe.pval.all.cl$snp.name,
sep = "_" )
pvals.GxM.CL
## promoter_rs742071 promoter_rs560426 promoter_rs642961
## 0.06400049 0.42805451 0.55340321
## promoter_rs7590268 promoter_rs12543318 promoter_rs987525
## 0.17320790 0.32332097 0.42999212
## promoter_rs3758249 promoter_rs8001641 promoter_rs1873147
## 0.87703914 0.53417280 0.39764783
## promoter_rs227731 promoter_rs13041247 enhancer_rs742071
## 0.94003164 0.55666521 0.34534483
## enhancer_rs560426 enhancer_rs7590268 enhancer_rs987525
## 0.90704433 0.10178491 0.41424021
## enhancer_rs7078160 enhancer_rs8001641 gene_rs742071
## 0.65667196 0.63076406 0.91766032
## gene_rs560426 gene_rs642961 gene_rs7590268
## 0.34030217 0.15917225 0.96394185
## gene_rs3758249 gene_rs7078160 gene_rs1873147
## 0.79767473 0.65667196 0.63297345
## gene_rs13041247
## 0.92200517
cur.len <- length( pvals.GxM.CL )
cur.nlabs <- 2
cur.xcoords <- get_x_coords( length( pvals.GxM.CL ), cur.nlabs )
cur.ycoords <- -log10( sort( pvals.GxM.CL ) )[ 1:cur.nlabs ]
cur.labels <- data.frame( x = cur.xcoords, y = cur.ycoords,
label = names( cur.ycoords ) )
gg_qqplot( pvals.GxM.CL ) +
labs( title = "GxMe, CLO dataset" ) +
geom_label_repel( data = cur.labels, aes( x, y, label = label ) ) +
theme( axis.text = element_text( size = 12 ), axis.title = element_text( size = 14 ),
title = element_text( size = 16 ) )
## Warning: Removed 2 rows containing missing values (geom_path).
ggsave( filename = "GxM_CL_mult_betas_gg-qqplot_square.png", width = 14, units = "cm",
height = 14 )
## Warning: Removed 2 rows containing missing values (geom_path).
# ==== POOxM ====
# now, PoOxM for CLO
pooxm.pval.promoter <- t( sapply( gxe.prom.all, function(x){
return( data.frame( pval = x$pooxm$pval[ 2 ], df = x$pooxm$df[ 2 ] ) )
} ) )
pooxm.pval.enhancer <- t( sapply( gxe.enh.all, function(x){
return( data.frame( pval = x$pooxm$gxe.test$pval[ 2 ], df = x$pooxm$gxe.test$df[ 2 ] ) )
} ) )
pooxm.pval.gene <- t( sapply( gxe.gene.all, function(x){
return( data.frame( pval = x$pooxm$gxe.test$pval[ 2 ], df = x$pooxm$gxe.test$df[ 2 ] ) )
} ) )
# creating the data.frame
my.SNP.map$Marker <- tolower( my.SNP.map$Marker )
cur.snp.names <- my.SNP.map$Marker[ match( c( rownames( pooxm.pval.promoter ),
rownames( pooxm.pval.enhancer ), rownames( pooxm.pval.gene ) ),
my.SNP.map$gene_name ) ]
pooxm.pval.all <- data.frame( region = c( rep( "promoter", nrow( pooxm.pval.promoter ) ),
rep( "enhancer", nrow( pooxm.pval.enhancer ) ), rep( "gene", nrow( pooxm.pval.gene ) ) ),
gene.name = c( rownames( pooxm.pval.promoter ), rownames( pooxm.pval.enhancer ), rownames( pooxm.pval.gene ) ),
snp.name = cur.snp.names,
pval = unlist( c( pooxm.pval.promoter[ , "pval" ], pooxm.pval.enhancer[ , "pval" ], pooxm.pval.gene[ , "pval" ] ) ),
df = unlist( c( pooxm.pval.promoter[ , "df" ], pooxm.pval.enhancer[ , "df" ], pooxm.pval.gene[ , "df" ] ) ) )
head( pooxm.pval.all )
## region gene.name snp.name pval df
## 1 promoter PAX7 rs742071 0.92671452 2
## 2 promoter ABCA4 rs560426 0.71507668 2
## 3 promoter IRF6 rs642961 0.10381984 2
## 4 promoter THADA rs7590268 0.60047269 2
## 5 promoter 8q21.3 rs12543318 0.02428055 2
## 6 promoter 8q24 rs987525 0.08685998 2
dim( pooxm.pval.all )
## [1] 25 5
df.all.pooxm <- unique( pooxm.pval.all$df )
df.all.pooxm
## [1] 2
# plotting
pvals.POOxM.CL <- pooxm.pval.all$pval
names( pvals.POOxM.CL ) <- paste( pooxm.pval.all$region, pooxm.pval.all$snp.name, sep = "_" )
pvals.POOxM.CL
## promoter_rs742071 promoter_rs560426 promoter_rs642961
## 0.92671452 0.71507668 0.10381984
## promoter_rs7590268 promoter_rs12543318 promoter_rs987525
## 0.60047269 0.02428055 0.08685998
## promoter_rs3758249 promoter_rs8001641 promoter_rs1873147
## 0.61687631 0.14989054 0.32261930
## promoter_rs227731 promoter_rs13041247 enhancer_rs742071
## 0.01193684 0.31593297 0.54294633
## enhancer_rs560426 enhancer_rs7590268 enhancer_rs987525
## 0.95324526 0.11209232 0.52952510
## enhancer_rs7078160 enhancer_rs8001641 gene_rs742071
## 0.55193606 0.02946621 0.29444331
## gene_rs560426 gene_rs642961 gene_rs7590268
## 0.01350852 0.91256337 0.33880289
## gene_rs3758249 gene_rs7078160 gene_rs1873147
## 0.57162216 0.55193606 0.07496376
## gene_rs13041247
## 0.40997150
cur.len <- length( pvals.POOxM.CL )
cur.nlabs <- 8
cur.xcoords <- get_x_coords( length( pvals.POOxM.CL ), cur.nlabs )
cur.ycoords <- -log10( sort( pvals.POOxM.CL ) )[ 1:cur.nlabs ]
cur.labels <- data.frame( x = cur.xcoords, y = cur.ycoords,
label = names( cur.ycoords ) )
gg_qqplot( pvals.POOxM.CL ) +
labs( title = "PoOxMe, CLO dataset" ) +
geom_label_repel( data = cur.labels, aes( x, y, label = label ), xlim = c( 1,NA ) ) +
theme( axis.text = element_text( size = 12 ), axis.title = element_text( size = 14 ),
title = element_text( size = 16 ) )
## Warning: Removed 2 rows containing missing values (geom_path).
ggsave( filename = "PoOxM_CL_mult_betas_gg-qqplot_square.png", width = 14, units = "cm",
height = 14 )
## Warning: Removed 2 rows containing missing values (geom_path).
# ==== CL/P dataset ====
load( "../GxM_PoOxM_new_betas_mult_CLorP.RData" )
# ==== GxM ====
# GxM
gxe.pval.promoter.clorp <- t( sapply( gxe.prom.all, function(x){
return( data.frame( pval = x$gxm$pval[ 2 ], df = x$gxm$df[ 2 ] ) )
} ) )
gxe.pval.enhancer.clorp <- t( sapply( gxe.enh.all, function(x){
return( data.frame( pval = x$gxm$gxe.test$pval[ 2 ], df = x$gxm$gxe.test$df[ 2 ] ) )
} ) )
gxe.pval.gene.clorp <- t( sapply( gxe.gene.all, function(x){
return( data.frame( pval = x$gxm$gxe.test$pval[ 2 ], df = x$gxm$gxe.test$df[ 2 ] ) )
} ) )
# creating the data.frame
cur.snp.names <- as.character( my.SNP.map$Marker[ match( c( rownames( gxe.pval.promoter.clorp ),
rownames( gxe.pval.enhancer.clorp ), rownames( gxe.pval.gene.clorp ) ),
my.SNP.map$gene_name ) ] )
gxe.pval.all.clorp <- data.frame( region = c( rep( "promoter", nrow( gxe.pval.promoter.clorp ) ),
rep( "enhancer", nrow( gxe.pval.enhancer.clorp ) ), rep( "gene", nrow( gxe.pval.gene.clorp ) ) ),
gene.name = c( rownames( gxe.pval.promoter.clorp ), rownames( gxe.pval.enhancer.clorp ),
rownames( gxe.pval.gene.clorp ) ),
snp.name = cur.snp.names,
pval = unlist( c( gxe.pval.promoter.clorp[ , "pval" ], gxe.pval.enhancer.clorp[ , "pval" ],
gxe.pval.gene.clorp[ , "pval" ] ) ),
df = unlist( c( gxe.pval.promoter.clorp[ , "df" ], gxe.pval.enhancer.clorp[ , "df" ],
gxe.pval.gene.clorp[ , "df" ] ) ) )
head( gxe.pval.all.clorp )
## region gene.name snp.name pval df
## 1 promoter PAX7 rs742071 0.19699287 2
## 2 promoter ABCA4 rs560426 0.78069733 2
## 3 promoter IRF6 rs642961 0.79869026 2
## 4 promoter THADA rs7590268 0.69939504 2
## 5 promoter 8q21.3 rs12543318 0.01231238 2
## 6 promoter 8q24 rs987525 0.11788477 2
dim( gxe.pval.all.clorp )
## [1] 25 5
df.all.clorp <- unique( gxe.pval.all.clorp$df )
df.all.clorp
## [1] 2
# plotting
pvals.GxM.clorp <- gxe.pval.all.clorp$pval
names( pvals.GxM.clorp ) <- paste( gxe.pval.all.clorp$region, gxe.pval.all.clorp$snp.name,
sep = "_" )
pvals.GxM.clorp
## promoter_rs742071 promoter_rs560426 promoter_rs642961
## 0.19699287 0.78069733 0.79869026
## promoter_rs7590268 promoter_rs12543318 promoter_rs987525
## 0.69939504 0.01231238 0.11788477
## promoter_rs3758249 promoter_rs8001641 promoter_rs1873147
## 0.65311965 0.05713693 0.30101309
## promoter_rs227731 promoter_rs13041247 enhancer_rs742071
## 0.78033536 0.87885490 0.59389868
## enhancer_rs560426 enhancer_rs7590268 enhancer_rs987525
## 0.64617122 0.48595842 0.44971466
## enhancer_rs7078160 enhancer_rs8001641 gene_rs742071
## 0.63353369 0.14277888 0.29639614
## gene_rs560426 gene_rs642961 gene_rs7590268
## 0.12167850 0.58935782 0.64241648
## gene_rs3758249 gene_rs7078160 gene_rs1873147
## 0.06864719 0.63353369 0.83821146
## gene_rs13041247
## 0.81390511
cur.len <- length( pvals.GxM.clorp )
cur.nlabs <- 2
cur.xcoords <- get_x_coords( length( pvals.GxM.clorp ), cur.nlabs )
cur.ycoords <- -log10( sort( pvals.GxM.clorp ) )[ 1:cur.nlabs ]
cur.labels <- data.frame( x = cur.xcoords, y = cur.ycoords,
label = names( cur.ycoords ) )
gg_qqplot( pvals.GxM.clorp ) +
labs( title = "GxMe, CL/P dataset" ) +
geom_label_repel( data = cur.labels, aes( x, y, label = label ) ) +
theme( axis.text = element_text( size = 12 ), axis.title = element_text( size = 14 ),
title = element_text( size = 16 ) )
## Warning: Removed 2 rows containing missing values (geom_path).
ggsave( filename = "GxM_CLorP_mult_betas_gg-qqplot_square.png", width = 14, units = "cm",
height = 14 )
## Warning: Removed 2 rows containing missing values (geom_path).
# ==== POOxM ====
# now, PoOxM for CL/P
pooxm.pval.promoter <- t( sapply( gxe.prom.all, function(x){
return( data.frame( pval = x$pooxm$pval[ 2 ], df = x$pooxm$df[ 2 ] ) )
} ) )
pooxm.pval.enhancer <- t( sapply( gxe.enh.all, function(x){
return( data.frame( pval = x$pooxm$gxe.test$pval[ 2 ], df = x$pooxm$gxe.test$df[ 2 ] ) )
} ) )
pooxm.pval.gene <- t( sapply( gxe.gene.all, function(x){
return( data.frame( pval = x$pooxm$gxe.test$pval[ 2 ], df = x$pooxm$gxe.test$df[ 2 ] ) )
} ) )
# creating the data.frame
my.SNP.map$Marker <- tolower( my.SNP.map$Marker )
cur.snp.names <- my.SNP.map$Marker[ match( c( rownames( pooxm.pval.promoter ),
rownames( pooxm.pval.enhancer ), rownames( pooxm.pval.gene ) ),
my.SNP.map$gene_name ) ]
pooxm.pval.all <- data.frame( region = c( rep( "promoter", nrow( pooxm.pval.promoter ) ),
rep( "enhancer", nrow( pooxm.pval.enhancer ) ), rep( "gene", nrow( pooxm.pval.gene ) ) ),
gene.name = c( rownames( pooxm.pval.promoter ), rownames( pooxm.pval.enhancer ), rownames( pooxm.pval.gene ) ),
snp.name = cur.snp.names,
pval = unlist( c( pooxm.pval.promoter[ , "pval" ], pooxm.pval.enhancer[ , "pval" ], pooxm.pval.gene[ , "pval" ] ) ),
df = unlist( c( pooxm.pval.promoter[ , "df" ], pooxm.pval.enhancer[ , "df" ], pooxm.pval.gene[ , "df" ] ) ) )
head( pooxm.pval.all )
## region gene.name snp.name pval df
## 1 promoter PAX7 rs742071 0.4903087 2
## 2 promoter ABCA4 rs560426 0.7649082 2
## 3 promoter IRF6 rs642961 0.3370442 2
## 4 promoter THADA rs7590268 0.4365098 2
## 5 promoter 8q21.3 rs12543318 0.4206556 2
## 6 promoter 8q24 rs987525 0.0584274 2
dim( pooxm.pval.all )
## [1] 25 5
df.all.pooxm <- unique( pooxm.pval.all$df )
df.all.pooxm
## [1] 2
# plotting
pvals.POOxM.clorp <- pooxm.pval.all$pval
names( pvals.POOxM.clorp ) <- paste( pooxm.pval.all$region, pooxm.pval.all$snp.name, sep = "_" )
pvals.POOxM.clorp
## promoter_rs742071 promoter_rs560426 promoter_rs642961
## 0.490308668 0.764908243 0.337044234
## promoter_rs7590268 promoter_rs12543318 promoter_rs987525
## 0.436509753 0.420655613 0.058427399
## promoter_rs3758249 promoter_rs8001641 promoter_rs1873147
## 0.059053176 0.343891913 0.126549691
## promoter_rs227731 promoter_rs13041247 enhancer_rs742071
## 0.003223564 0.105488490 0.434832066
## enhancer_rs560426 enhancer_rs7590268 enhancer_rs987525
## 0.547990224 0.054367326 0.867650691
## enhancer_rs7078160 enhancer_rs8001641 gene_rs742071
## 0.134494271 0.608589146 0.997708266
## gene_rs560426 gene_rs642961 gene_rs7590268
## 0.906781682 0.014342080 0.486847137
## gene_rs3758249 gene_rs7078160 gene_rs1873147
## 0.040261280 0.134494271 0.420966767
## gene_rs13041247
## 0.211418716
cur.len <- length( pvals.POOxM.clorp )
cur.nlabs <- 8
cur.xcoords <- get_x_coords( length( pvals.POOxM.clorp ), cur.nlabs )
cur.ycoords <- -log10( sort( pvals.POOxM.clorp ) )[ 1:cur.nlabs ]
cur.labels <- data.frame( x = cur.xcoords, y = cur.ycoords,
label = names( cur.ycoords ) )
gg_qqplot( pvals.POOxM.clorp ) +
labs( title = "PoOxMe, CL/P dataset" ) +
geom_label_repel( data = cur.labels, aes( x, y, label = label ), #xlim = c( 1.5,NA ),
force = 10, nudge_y = cur.labels$x - 1, nudge_x = 2 - cur.labels$x ) +
theme( axis.text = element_text( size = 12 ), axis.title = element_text( size = 14 ),
title = element_text( size = 16 ) )
ggsave( filename = "PoOxM_CLorP_mult_betas_gg-qqplot_square.png", width = 14, units = "cm",
height = 14 )
# ==== CLP dataset ====
load( "../GxM_PoOxM_new_betas_mult_CLP.RData" )
# ==== GxM ====
# GxM
gxe.pval.promoter.clp <- t( sapply( gxe.prom.all, function(x){
return( data.frame( pval = x$gxm$pval[ 2 ], df = x$gxm$df[ 2 ] ) )
} ) )
gxe.pval.enhancer.clp <- t( sapply( gxe.enh.all, function(x){
return( data.frame( pval = x$gxm$gxe.test$pval[ 2 ], df = x$gxm$gxe.test$df[ 2 ] ) )
} ) )
gxe.pval.gene.clp <- t( sapply( gxe.gene.all, function(x){
return( data.frame( pval = x$gxm$gxe.test$pval[ 2 ], df = x$gxm$gxe.test$df[ 2 ] ) )
} ) )
# creating the data.frame
cur.snp.names <- as.character( my.SNP.map$Marker[ match( c( rownames( gxe.pval.promoter.clp ),
rownames( gxe.pval.enhancer.clp ), rownames( gxe.pval.gene.clp ) ),
my.SNP.map$gene_name ) ] )
gxe.pval.all.clp <- data.frame( region = c( rep( "promoter", nrow( gxe.pval.promoter.clp ) ),
rep( "enhancer", nrow( gxe.pval.enhancer.clp ) ), rep( "gene", nrow( gxe.pval.gene.clp ) ) ),
gene.name = c( rownames( gxe.pval.promoter.clp ), rownames( gxe.pval.enhancer.clp ),
rownames( gxe.pval.gene.clp ) ),
snp.name = cur.snp.names,
pval = unlist( c( gxe.pval.promoter.clp[ , "pval" ], gxe.pval.enhancer.clp[ , "pval" ],
gxe.pval.gene.clp[ , "pval" ] ) ),
df = unlist( c( gxe.pval.promoter.clp[ , "df" ], gxe.pval.enhancer.clp[ , "df" ],
gxe.pval.gene.clp[ , "df" ] ) ) )
head( gxe.pval.all.clp )
## region gene.name snp.name pval df
## 1 promoter PAX7 rs742071 0.2806608 2
## 2 promoter ABCA4 rs560426 0.2512834 2
## 3 promoter IRF6 rs642961 0.8077746 2
## 4 promoter THADA rs7590268 0.6258189 2
## 5 promoter 8q21.3 rs12543318 0.2889103 2
## 6 promoter 8q24 rs987525 0.2100116 2
dim( gxe.pval.all.clp )
## [1] 25 5
df.all.clp <- unique( gxe.pval.all.clp$df )
df.all.clp
## [1] 2
# plotting
pvals.GxM.clp <- gxe.pval.all.clp$pval
names( pvals.GxM.clp ) <- paste( gxe.pval.all.clp$region, gxe.pval.all.clp$snp.name,
sep = "_" )
pvals.GxM.clp
## promoter_rs742071 promoter_rs560426 promoter_rs642961
## 0.280660840 0.251283381 0.807774607
## promoter_rs7590268 promoter_rs12543318 promoter_rs987525
## 0.625818873 0.288910280 0.210011574
## promoter_rs3758249 promoter_rs8001641 promoter_rs1873147
## 0.703687666 0.217081847 0.230290238
## promoter_rs227731 promoter_rs13041247 enhancer_rs742071
## 0.786378069 0.955859512 0.500717570
## enhancer_rs560426 enhancer_rs7590268 enhancer_rs987525
## 0.430706012 0.899082901 0.541993252
## enhancer_rs7078160 enhancer_rs8001641 gene_rs742071
## 0.439422383 0.247811659 0.190955186
## gene_rs560426 gene_rs642961 gene_rs7590268
## 0.686837728 0.885820221 0.975074491
## gene_rs3758249 gene_rs7078160 gene_rs1873147
## 0.001039765 0.439422383 0.385127538
## gene_rs13041247
## 0.842246318
cur.len <- length( pvals.GxM.clp )
cur.nlabs <- 2
cur.xcoords <- get_x_coords( length( pvals.GxM.clp ), cur.nlabs )
cur.ycoords <- -log10( sort( pvals.GxM.clp ) )[ 1:cur.nlabs ]
cur.labels <- data.frame( x = cur.xcoords, y = cur.ycoords,
label = names( cur.ycoords ) )
gg_qqplot( pvals.GxM.clp ) +
labs( title = "GxMe, CLP dataset" ) +
geom_label_repel( data = cur.labels, aes( x, y, label = label ) ) +
theme( axis.text = element_text( size = 12 ), axis.title = element_text( size = 14 ),
title = element_text( size = 16 ) )
ggsave( filename = "GxM_CLP_mult_betas_gg-qqplot_square.png", width = 14, units = "cm",
height = 14 )
# ==== POOxM ====
# now, PoOxM for CLP
pooxm.pval.promoter <- t( sapply( gxe.prom.all, function(x){
return( data.frame( pval = x$pooxm$pval[ 2 ], df = x$pooxm$df[ 2 ] ) )
} ) )
pooxm.pval.enhancer <- t( sapply( gxe.enh.all, function(x){
return( data.frame( pval = x$pooxm$gxe.test$pval[ 2 ], df = x$pooxm$gxe.test$df[ 2 ] ) )
} ) )
pooxm.pval.gene <- t( sapply( gxe.gene.all, function(x){
return( data.frame( pval = x$pooxm$gxe.test$pval[ 2 ], df = x$pooxm$gxe.test$df[ 2 ] ) )
} ) )
# creating the data.frame
my.SNP.map$Marker <- tolower( my.SNP.map$Marker )
cur.snp.names <- my.SNP.map$Marker[ match( c( rownames( pooxm.pval.promoter ),
rownames( pooxm.pval.enhancer ), rownames( pooxm.pval.gene ) ),
my.SNP.map$gene_name ) ]
pooxm.pval.all <- data.frame( region = c( rep( "promoter", nrow( pooxm.pval.promoter ) ),
rep( "enhancer", nrow( pooxm.pval.enhancer ) ), rep( "gene", nrow( pooxm.pval.gene ) ) ),
gene.name = c( rownames( pooxm.pval.promoter ), rownames( pooxm.pval.enhancer ), rownames( pooxm.pval.gene ) ),
snp.name = cur.snp.names,
pval = unlist( c( pooxm.pval.promoter[ , "pval" ], pooxm.pval.enhancer[ , "pval" ], pooxm.pval.gene[ , "pval" ] ) ),
df = unlist( c( pooxm.pval.promoter[ , "df" ], pooxm.pval.enhancer[ , "df" ], pooxm.pval.gene[ , "df" ] ) ) )
head( pooxm.pval.all )
## region gene.name snp.name pval df
## 1 promoter PAX7 rs742071 0.167697016 2
## 2 promoter ABCA4 rs560426 0.986305286 2
## 3 promoter IRF6 rs642961 0.049270673 2
## 4 promoter THADA rs7590268 0.765851843 2
## 5 promoter 8q21.3 rs12543318 0.569978278 2
## 6 promoter 8q24 rs987525 0.003653522 2
dim( pooxm.pval.all )
## [1] 25 5
df.all.pooxm <- unique( pooxm.pval.all$df )
df.all.pooxm
## [1] 2
# plotting
pvals.POOxM.clp <- pooxm.pval.all$pval
names( pvals.POOxM.clp ) <- paste( pooxm.pval.all$region, pooxm.pval.all$snp.name, sep = "_" )
pvals.POOxM.clp
## promoter_rs742071 promoter_rs560426 promoter_rs642961
## 0.167697016 0.986305286 0.049270673
## promoter_rs7590268 promoter_rs12543318 promoter_rs987525
## 0.765851843 0.569978278 0.003653522
## promoter_rs3758249 promoter_rs8001641 promoter_rs1873147
## 0.274574235 0.882450263 0.328849307
## promoter_rs227731 promoter_rs13041247 enhancer_rs742071
## 0.009281047 0.113938408 0.061957648
## enhancer_rs560426 enhancer_rs7590268 enhancer_rs987525
## 0.568982718 0.145883127 0.917629435
## enhancer_rs7078160 enhancer_rs8001641 gene_rs742071
## 0.079368084 0.015638002 0.665827849
## gene_rs560426 gene_rs642961 gene_rs7590268
## 0.372431164 0.014836717 0.682446660
## gene_rs3758249 gene_rs7078160 gene_rs1873147
## 0.184688210 0.079368084 0.214447224
## gene_rs13041247
## 0.144298379
cur.len <- length( pvals.POOxM.clp )
cur.nlabs <- 8
cur.xcoords <- get_x_coords( length( pvals.POOxM.clp ), cur.nlabs )
cur.ycoords <- -log10( sort( pvals.POOxM.clp ) )[ 1:cur.nlabs ]
cur.labels <- data.frame( x = cur.xcoords, y = cur.ycoords,
label = names( cur.ycoords ) )
gg_qqplot( pvals.POOxM.clp ) +
labs( title = "PoOxMe, CLP dataset" ) +
geom_label_repel( data = cur.labels, aes( x, y, label = label ),# xlim = c( 1,NA ) ) +
force = 10, nudge_y = cur.labels$x - 1, nudge_x = 2 - cur.labels$x ) +
theme( axis.text = element_text( size = 12 ), axis.title = element_text( size = 14 ),
title = element_text( size = 16 ) )
ggsave( filename = "PoOxM_CLP_mult_betas_gg-qqplot_square.png", width = 14, units = "cm",
height = 14 )
# ==== CPO dataset ====
load( "../GxM_PoOxM_new_betas_mult_CPO.RData" )
# ==== GxM ====
# GxM
gxe.pval.promoter.cpo <- t( sapply( gxe.prom.all, function(x){
return( data.frame( pval = x$gxm$pval[ 2 ], df = x$gxm$df[ 2 ] ) )
} ) )
gxe.pval.enhancer.cpo <- t( sapply( gxe.enh.all, function(x){
return( data.frame( pval = x$gxm$gxe.test$pval[ 2 ], df = x$gxm$gxe.test$df[ 2 ] ) )
} ) )
gxe.pval.gene.cpo <- t( sapply( gxe.gene.all, function(x){
return( data.frame( pval = x$gxm$gxe.test$pval[ 2 ], df = x$gxm$gxe.test$df[ 2 ] ) )
} ) )
# creating the data.frame
cur.snp.names <- as.character( my.SNP.map$Marker[ match( c( rownames( gxe.pval.promoter.cpo ),
rownames( gxe.pval.enhancer.cpo ), rownames( gxe.pval.gene.cpo ) ),
my.SNP.map$gene_name ) ] )
gxe.pval.all.cpo <- data.frame( region = c( rep( "promoter", nrow( gxe.pval.promoter.cpo ) ),
rep( "enhancer", nrow( gxe.pval.enhancer.cpo ) ), rep( "gene", nrow( gxe.pval.gene.cpo ) ) ),
gene.name = c( rownames( gxe.pval.promoter.cpo ), rownames( gxe.pval.enhancer.cpo ),
rownames( gxe.pval.gene.cpo ) ),
snp.name = cur.snp.names,
pval = unlist( c( gxe.pval.promoter.cpo[ , "pval" ], gxe.pval.enhancer.cpo[ , "pval" ],
gxe.pval.gene.cpo[ , "pval" ] ) ),
df = unlist( c( gxe.pval.promoter.cpo[ , "df" ], gxe.pval.enhancer.cpo[ , "df" ],
gxe.pval.gene.cpo[ , "df" ] ) ) )
head( gxe.pval.all.cpo )
## region gene.name snp.name pval df
## 1 promoter PAX7 rs742071 0.54252058 2
## 2 promoter ABCA4 rs560426 0.51196686 2
## 3 promoter IRF6 rs642961 0.79769113 2
## 4 promoter THADA rs7590268 0.43698267 2
## 5 promoter 8q21.3 rs12543318 0.05562257 2
## 6 promoter 8q24 rs987525 0.03553574 2
dim( gxe.pval.all.cpo )
## [1] 25 5
df.all.cpo <- unique( gxe.pval.all.cpo$df )
df.all.cpo
## [1] 2
# plotting
pvals.GxM.cpo <- gxe.pval.all.cpo$pval
names( pvals.GxM.cpo ) <- paste( gxe.pval.all.cpo$region, gxe.pval.all.cpo$snp.name,
sep = "_" )
pvals.GxM.cpo
## promoter_rs742071 promoter_rs560426 promoter_rs642961
## 0.54252058 0.51196686 0.79769113
## promoter_rs7590268 promoter_rs12543318 promoter_rs987525
## 0.43698267 0.05562257 0.03553574
## promoter_rs3758249 promoter_rs8001641 promoter_rs1873147
## 0.90363749 0.36284283 0.83325550
## promoter_rs227731 promoter_rs13041247 enhancer_rs742071
## 0.12712688 0.67827329 0.55712772
## enhancer_rs560426 enhancer_rs7590268 enhancer_rs987525
## 0.41986313 0.78638193 0.87588013
## enhancer_rs7078160 enhancer_rs8001641 gene_rs742071
## 0.86240975 0.48897261 0.48787983
## gene_rs560426 gene_rs642961 gene_rs7590268
## 0.37242464 0.98266263 0.82984563
## gene_rs3758249 gene_rs7078160 gene_rs1873147
## 0.07646270 0.86240975 0.75762619
## gene_rs13041247
## 0.45940485
cur.len <- length( pvals.GxM.cpo )
cur.nlabs <- 2
cur.xcoords <- get_x_coords( length( pvals.GxM.cpo ), cur.nlabs )
cur.ycoords <- -log10( sort( pvals.GxM.cpo ) )[ 1:cur.nlabs ]
cur.labels <- data.frame( x = cur.xcoords, y = cur.ycoords,
label = names( cur.ycoords ) )
gg_qqplot( pvals.GxM.cpo ) +
labs( title = "GxMe, CPO dataset" ) +
geom_label_repel( data = cur.labels, aes( x, y, label = label ) ) +
theme( axis.text = element_text( size = 12 ), axis.title = element_text( size = 14 ),
title = element_text( size = 16 ) )
## Warning: Removed 2 rows containing missing values (geom_path).
ggsave( filename = "GxM_CPO_mult_betas_gg-qqplot_square.png", width = 14, units = "cm",
height = 14 )
## Warning: Removed 2 rows containing missing values (geom_path).
# ==== POOxM ====
# now, PoOxM for CPO
pooxm.pval.promoter <- t( sapply( gxe.prom.all, function(x){
return( data.frame( pval = x$pooxm$pval[ 2 ], df = x$pooxm$df[ 2 ] ) )
} ) )
pooxm.pval.enhancer <- t( sapply( gxe.enh.all, function(x){
return( data.frame( pval = x$pooxm$gxe.test$pval[ 2 ], df = x$pooxm$gxe.test$df[ 2 ] ) )
} ) )
pooxm.pval.gene <- t( sapply( gxe.gene.all, function(x){
return( data.frame( pval = x$pooxm$gxe.test$pval[ 2 ], df = x$pooxm$gxe.test$df[ 2 ] ) )
} ) )
# creating the data.frame
my.SNP.map$Marker <- tolower( my.SNP.map$Marker )
cur.snp.names <- my.SNP.map$Marker[ match( c( rownames( pooxm.pval.promoter ),
rownames( pooxm.pval.enhancer ), rownames( pooxm.pval.gene ) ),
my.SNP.map$gene_name ) ]
pooxm.pval.all <- data.frame( region = c( rep( "promoter", nrow( pooxm.pval.promoter ) ),
rep( "enhancer", nrow( pooxm.pval.enhancer ) ), rep( "gene", nrow( pooxm.pval.gene ) ) ),
gene.name = c( rownames( pooxm.pval.promoter ), rownames( pooxm.pval.enhancer ), rownames( pooxm.pval.gene ) ),
snp.name = cur.snp.names,
pval = unlist( c( pooxm.pval.promoter[ , "pval" ], pooxm.pval.enhancer[ , "pval" ], pooxm.pval.gene[ , "pval" ] ) ),
df = unlist( c( pooxm.pval.promoter[ , "df" ], pooxm.pval.enhancer[ , "df" ], pooxm.pval.gene[ , "df" ] ) ) )
head( pooxm.pval.all )
## region gene.name snp.name pval df
## 1 promoter PAX7 rs742071 0.1707909 2
## 2 promoter ABCA4 rs560426 0.5706230 2
## 3 promoter IRF6 rs642961 0.8569329 2
## 4 promoter THADA rs7590268 0.3167452 2
## 5 promoter 8q21.3 rs12543318 0.4721891 2
## 6 promoter 8q24 rs987525 0.3090114 2
dim( pooxm.pval.all )
## [1] 25 5
df.all.pooxm <- unique( pooxm.pval.all$df )
df.all.pooxm
## [1] 2
# plotting
pvals.POOxM.cpo <- pooxm.pval.all$pval
names( pvals.POOxM.cpo ) <- paste( pooxm.pval.all$region, pooxm.pval.all$snp.name, sep = "_" )
pvals.POOxM.cpo
## promoter_rs742071 promoter_rs560426 promoter_rs642961
## 0.17079094 0.57062300 0.85693294
## promoter_rs7590268 promoter_rs12543318 promoter_rs987525
## 0.31674519 0.47218915 0.30901137
## promoter_rs3758249 promoter_rs8001641 promoter_rs1873147
## 0.96604798 0.21219071 0.20065785
## promoter_rs227731 promoter_rs13041247 enhancer_rs742071
## 0.35923455 0.02423747 0.22089377
## enhancer_rs560426 enhancer_rs7590268 enhancer_rs987525
## 0.94292065 0.24664514 0.17507179
## enhancer_rs7078160 enhancer_rs8001641 gene_rs742071
## 0.69055131 0.29045258 0.81670411
## gene_rs560426 gene_rs642961 gene_rs7590268
## 0.80351828 0.89775943 0.12314425
## gene_rs3758249 gene_rs7078160 gene_rs1873147
## 0.23027679 0.69055131 0.69751762
## gene_rs13041247
## 0.13585876
cur.len <- length( pvals.POOxM.cpo )
cur.nlabs <- 8
cur.xcoords <- get_x_coords( length( pvals.POOxM.cpo ), cur.nlabs )
cur.ycoords <- -log10( sort( pvals.POOxM.cpo ) )[ 1:cur.nlabs ]
cur.labels <- data.frame( x = cur.xcoords, y = cur.ycoords,
label = names( cur.ycoords ) )
gg_qqplot( pvals.POOxM.cpo ) +
labs( title = "PoOxMe, CPO dataset" ) +
geom_label_repel( data = cur.labels, aes( x, y, label = label ), #xlim = c( 1,NA ) ) +
force = 10, nudge_y = cur.labels$x - 1, nudge_x = 1.5 - cur.labels$x ) +
theme( axis.text = element_text( size = 12 ), axis.title = element_text( size = 14 ),
title = element_text( size = 16 ) )
## Warning: Removed 2 rows containing missing values (geom_path).
ggsave( filename = "PoOxM_CPO_mult_betas_gg-qqplot_square.png", width = 14, units = "cm",
height = 14 )
## Warning: Removed 2 rows containing missing values (geom_path).
# ==== control dataset ====
load( "../GxM_PoOxM_new_betas_mult_cntrl.RData" )
# ==== GxM ====
# GxM
gxe.pval.promoter.cntrl <- t( sapply( gxe.prom.all, function(x){
return( data.frame( pval = x$gxm$pval[ 2 ], df = x$gxm$df[ 2 ] ) )
} ) )
gxe.pval.enhancer.cntrl <- t( sapply( gxe.enh.all, function(x){
return( data.frame( pval = x$gxm$gxe.test$pval[ 2 ], df = x$gxm$gxe.test$df[ 2 ] ) )
} ) )
gxe.pval.gene.cntrl <- t( sapply( gxe.gene.all, function(x){
return( data.frame( pval = x$gxm$gxe.test$pval[ 2 ], df = x$gxm$gxe.test$df[ 2 ] ) )
} ) )
# creating the data.frame
cur.snp.names <- as.character( my.SNP.map$Marker[ match( c( rownames( gxe.pval.promoter.cntrl ),
rownames( gxe.pval.enhancer.cntrl ), rownames( gxe.pval.gene.cntrl ) ),
my.SNP.map$gene_name ) ] )
gxe.pval.all.cntrl <- data.frame( region = c( rep( "promoter", nrow( gxe.pval.promoter.cntrl ) ),
rep( "enhancer", nrow( gxe.pval.enhancer.cntrl ) ), rep( "gene", nrow( gxe.pval.gene.cntrl ) ) ),
gene.name = c( rownames( gxe.pval.promoter.cntrl ), rownames( gxe.pval.enhancer.cntrl ),
rownames( gxe.pval.gene.cntrl ) ),
snp.name = cur.snp.names,
pval = unlist( c( gxe.pval.promoter.cntrl[ , "pval" ], gxe.pval.enhancer.cntrl[ , "pval" ],
gxe.pval.gene.cntrl[ , "pval" ] ) ),
df = unlist( c( gxe.pval.promoter.cntrl[ , "df" ], gxe.pval.enhancer.cntrl[ , "df" ],
gxe.pval.gene.cntrl[ , "df" ] ) ) )
head( gxe.pval.all.cntrl )
## region gene.name snp.name pval df
## 1 promoter PAX7 rs742071 0.71998085 2
## 2 promoter ABCA4 rs560426 0.83397007 2
## 3 promoter IRF6 rs642961 0.12346638 2
## 4 promoter THADA rs7590268 0.53910778 2
## 5 promoter 8q21.3 rs12543318 0.14634632 2
## 6 promoter 8q24 rs987525 0.04779387 2
dim( gxe.pval.all.cntrl )
## [1] 25 5
df.all.cntrl <- unique( gxe.pval.all.cntrl$df )
df.all.cntrl
## [1] 2
# plotting
pvals.GxM.cntrl <- gxe.pval.all.cntrl$pval
names( pvals.GxM.cntrl ) <- paste( gxe.pval.all.cntrl$region, gxe.pval.all.cntrl$snp.name,
sep = "_" )
pvals.GxM.cntrl
## promoter_rs742071 promoter_rs560426 promoter_rs642961
## 7.199808e-01 8.339701e-01 1.234664e-01
## promoter_rs7590268 promoter_rs12543318 promoter_rs987525
## 5.391078e-01 1.463463e-01 4.779387e-02
## promoter_rs3758249 promoter_rs8001641 promoter_rs1873147
## 5.900314e-01 3.211253e-01 5.599778e-05
## promoter_rs227731 promoter_rs13041247 enhancer_rs742071
## 5.739594e-01 2.748482e-01 5.950565e-01
## enhancer_rs560426 enhancer_rs7590268 enhancer_rs987525
## 1.553378e-01 8.476451e-01 2.213102e-01
## enhancer_rs7078160 enhancer_rs8001641 gene_rs742071
## 2.545081e-01 8.419980e-01 3.336203e-01
## gene_rs560426 gene_rs642961 gene_rs7590268
## 2.996933e-02 3.812790e-01 4.578903e-01
## gene_rs3758249 gene_rs7078160 gene_rs1873147
## 8.592355e-05 2.545081e-01 4.092147e-01
## gene_rs13041247
## 8.796321e-01
cur.len <- length( pvals.GxM.cntrl )
cur.nlabs <- 2
cur.xcoords <- get_x_coords( length( pvals.GxM.cntrl ), cur.nlabs )
cur.ycoords <- -log10( sort( pvals.GxM.cntrl ) )[ 1:cur.nlabs ]
cur.labels <- data.frame( x = cur.xcoords, y = cur.ycoords,
label = names( cur.ycoords ) )
gg_qqplot( pvals.GxM.cntrl ) +
labs( title = "GxMe, control dataset" ) +
geom_label_repel( data = cur.labels, aes( x, y, label = label ) ) +
theme( axis.text = element_text( size = 12 ), axis.title = element_text( size = 14 ),
title = element_text( size = 16 ) )
ggsave( filename = "GxM_cntrl_mult_betas_gg-qqplot_square.png", width = 14, units = "cm",
height = 14 )
# ==== POOxM ====
# now, PoOxM for cntrl
pooxm.pval.promoter <- t( sapply( gxe.prom.all, function(x){
return( data.frame( pval = x$pooxm$pval[ 2 ], df = x$pooxm$df[ 2 ] ) )
} ) )
pooxm.pval.enhancer <- t( sapply( gxe.enh.all, function(x){
return( data.frame( pval = x$pooxm$gxe.test$pval[ 2 ], df = x$pooxm$gxe.test$df[ 2 ] ) )
} ) )
pooxm.pval.gene <- t( sapply( gxe.gene.all, function(x){
return( data.frame( pval = x$pooxm$gxe.test$pval[ 2 ], df = x$pooxm$gxe.test$df[ 2 ] ) )
} ) )
# creating the data.frame
my.SNP.map$Marker <- tolower( my.SNP.map$Marker )
cur.snp.names <- my.SNP.map$Marker[ match( c( rownames( pooxm.pval.promoter ),
rownames( pooxm.pval.enhancer ), rownames( pooxm.pval.gene ) ),
my.SNP.map$gene_name ) ]
pooxm.pval.all <- data.frame( region = c( rep( "promoter", nrow( pooxm.pval.promoter ) ),
rep( "enhancer", nrow( pooxm.pval.enhancer ) ), rep( "gene", nrow( pooxm.pval.gene ) ) ),
gene.name = c( rownames( pooxm.pval.promoter ), rownames( pooxm.pval.enhancer ), rownames( pooxm.pval.gene ) ),
snp.name = cur.snp.names,
pval = unlist( c( pooxm.pval.promoter[ , "pval" ], pooxm.pval.enhancer[ , "pval" ], pooxm.pval.gene[ , "pval" ] ) ),
df = unlist( c( pooxm.pval.promoter[ , "df" ], pooxm.pval.enhancer[ , "df" ], pooxm.pval.gene[ , "df" ] ) ) )
head( pooxm.pval.all )
## region gene.name snp.name pval df
## 1 promoter PAX7 rs742071 0.406144153 2
## 2 promoter ABCA4 rs560426 0.103986612 2
## 3 promoter IRF6 rs642961 0.386871106 2
## 4 promoter THADA rs7590268 0.486207252 2
## 5 promoter 8q21.3 rs12543318 0.005556453 2
## 6 promoter 8q24 rs987525 0.080081329 2
dim( pooxm.pval.all )
## [1] 25 5
df.all.pooxm <- unique( pooxm.pval.all$df )
df.all.pooxm
## [1] 2
# plotting
pvals.POOxM.cntrl <- pooxm.pval.all$pval
names( pvals.POOxM.cntrl ) <- paste( pooxm.pval.all$region, pooxm.pval.all$snp.name, sep = "_" )
pvals.POOxM.cntrl
## promoter_rs742071 promoter_rs560426 promoter_rs642961
## 0.406144153 0.103986612 0.386871106
## promoter_rs7590268 promoter_rs12543318 promoter_rs987525
## 0.486207252 0.005556453 0.080081329
## promoter_rs3758249 promoter_rs8001641 promoter_rs1873147
## 0.609387191 0.630839815 0.910796277
## promoter_rs227731 promoter_rs13041247 enhancer_rs742071
## 0.486795993 0.977763642 0.802386461
## enhancer_rs560426 enhancer_rs7590268 enhancer_rs987525
## 0.028611471 0.512101622 0.731971506
## enhancer_rs7078160 enhancer_rs8001641 gene_rs742071
## 0.884796158 0.515237910 0.656457456
## gene_rs560426 gene_rs642961 gene_rs7590268
## 0.027490359 0.457588341 0.393289836
## gene_rs3758249 gene_rs7078160 gene_rs1873147
## 0.572414729 0.884796158 0.195249397
## gene_rs13041247
## 0.530862256
cur.len <- length( pvals.POOxM.cntrl )
cur.nlabs <- 8
cur.xcoords <- get_x_coords( length( pvals.POOxM.cntrl ), cur.nlabs )
cur.ycoords <- -log10( sort( pvals.POOxM.cntrl ) )[ 1:cur.nlabs ]
cur.labels <- data.frame( x = cur.xcoords, y = cur.ycoords,
label = names( cur.ycoords ) )
gg_qqplot( pvals.POOxM.cntrl ) +
labs( title = "PoOxMe, control dataset" ) +
geom_label_repel( data = cur.labels, aes( x, y, label = label ),# xlim = c( 1,NA ) ) +
force = 10, nudge_y = cur.labels$x - 1, nudge_x = 2 - cur.labels$x ) +
theme( axis.text = element_text( size = 12 ), axis.title = element_text( size = 14 ),
title = element_text( size = 16 ) )
ggsave( filename = "PoOxM_cntrl_mult_betas_gg-qqplot_square.png", width = 14, units = "cm",
height = 14 )
# ==== RANDOM SNPS, on CL/P ====
load( "../RandomSNPs/GxM_PoOxM_new_betas_mult_randomSNPs_CLorP.RData" )
# ==== GxM ====
# GxM
gxe.pval.promoter.random <- lapply( gxe.prom.all, function(x){
return( data.frame( pval = x$gxm$pval[ 2 ], df = x$gxm$df[ 2 ] ) )
} )
gxe.pval.promoter.random <- do.call( rbind, gxe.pval.promoter.random )
gxe.pval.enhancer.random <- lapply( gxe.enh.all, function(x){
return( data.frame( pval = x$gxm$gxe.test$pval[ 2 ], df = x$gxm$gxe.test$df[ 2 ] ) )
} )
gxe.pval.enhancer.random <- do.call( rbind, gxe.pval.enhancer.random )
gxe.pval.gene.random <- lapply( gxe.gene.all, function(x){
return( data.frame( pval = x$gxm$gxe.test$pval[ 2 ], df = x$gxm$gxe.test$df[ 2 ] ) )
} )
gxe.pval.gene.random <- do.call( rbind, gxe.pval.gene.random )
# creating the data.frame
cur.snp.names <- as.character( my.SNP.map$Marker[ match( c( rownames( gxe.pval.promoter.random ),
rownames( gxe.pval.enhancer.random ), rownames( gxe.pval.gene.random ) ),
my.SNP.map$Marker ) ] )
gxe.pval.all.random <- data.frame( region = c( rep( "promoter", nrow( gxe.pval.promoter.random ) ),
rep( "enhancer", nrow( gxe.pval.enhancer.random ) ), rep( "gene", nrow( gxe.pval.gene.random ) ) ),
gene.name = c( rownames( gxe.pval.promoter.random ), rownames( gxe.pval.enhancer.random ),
rownames( gxe.pval.gene.random ) ),
snp.name = cur.snp.names,
pval = unlist( c( gxe.pval.promoter.random[ , "pval" ], gxe.pval.enhancer.random[ , "pval" ],
gxe.pval.gene.random[ , "pval" ] ) ),
df = unlist( c( gxe.pval.promoter.random[ , "df" ], gxe.pval.enhancer.random[ , "df" ],
gxe.pval.gene.random[ , "df" ] ) ) )
head( gxe.pval.all.random )
## region gene.name snp.name pval df
## 1 promoter rs4621895 rs4621895 0.1180659 2
## 2 promoter rs1546124 rs1546124 0.1633977 2
## 3 promoter rs3750817 rs3750817 0.5874762 2
## 4 promoter rs2981428 rs2981428 0.4894996 2
## 5 promoter rs10984103 rs10984103 0.9874819 2
## 6 promoter rs9939609 rs9939609 0.2574505 2
dim( gxe.pval.all.random )
## [1] 33 5
df.all.random <- unique( gxe.pval.all.random$df )
df.all.random
## [1] 2
# plotting
pvals.GxM.random <- gxe.pval.all.random$pval
names( pvals.GxM.random ) <- paste( gxe.pval.all.random$region, gxe.pval.all.random$snp.name,
sep = "_" )
pvals.GxM.random
## promoter_rs4621895 promoter_rs1546124 promoter_rs3750817
## 0.11806593 0.16339770 0.58747621
## promoter_rs2981428 promoter_rs10984103 promoter_rs9939609
## 0.48949958 0.98748190 0.25745050
## promoter_rs6010718 promoter_rs3752462 promoter_rs2867125
## 0.97113217 0.78197580 0.67086580
## promoter_rs2912771 promoter_rs6659735 promoter_rs10767664
## 0.98656372 0.70087949 0.24549981
## promoter_rs7138803 promoter_rs713586 promoter_rs1443433
## 0.10721363 0.11450703 0.63696071
## promoter_rs4752566 enhancer_rs3750817 enhancer_rs2981428
## 0.54866549 0.80866325 0.91807293
## enhancer_rs2912771 enhancer_rs6659735 enhancer_rs7138803
## 0.95393044 0.49345554 0.85627120
## enhancer_rs713586 gene_rs4621895 gene_rs1546124
## 0.32166882 0.18940757 0.25943800
## gene_rs3750817 gene_rs2981428 gene_rs9939609
## 0.88149195 0.85487948 0.27799145
## gene_rs6010718 gene_rs3752462 gene_rs2912771
## 0.60463286 0.09274268 0.33032421
## gene_rs6659735 gene_rs10767664 gene_rs4752566
## 0.79403559 0.07807974 0.74079287
cur.len <- length( pvals.GxM.random )
cur.nlabs <- 4
cur.xcoords <- get_x_coords( length( pvals.GxM.random ), cur.nlabs )
cur.ycoords <- -log10( sort( pvals.GxM.random ) )[ 1:cur.nlabs ]
cur.labels <- data.frame( x = cur.xcoords, y = cur.ycoords,
label = names( cur.ycoords ) )
gg_qqplot( pvals.GxM.random ) +
labs( title = "GxMe, random SNPs, CL/P dataset" ) +
geom_label_repel( data = cur.labels, aes( x, y, label = label ) ) +
theme( axis.text = element_text( size = 12 ), axis.title = element_text( size = 14 ),
title = element_text( size = 16 ) )
## Warning: Removed 2 rows containing missing values (geom_path).
ggsave( filename = "GxM_random_CLorP_mult_betas_gg-qqplot_square.png", width = 14, units = "cm",
height = 14 )
## Warning: Removed 2 rows containing missing values (geom_path).
# ==== POOxM ====
# now, PoOxM for CLO
pooxm.pval.promoter <- lapply( gxe.prom.all, function(x){
return( data.frame( pval = x$pooxm$pval[ 2 ], df = x$pooxm$df[ 2 ] ) )
} )
pooxm.pval.promoter <- do.call( rbind, pooxm.pval.promoter )
pooxm.pval.enhancer <- lapply( gxe.enh.all, function(x){
return( data.frame( pval = x$pooxm$gxe.test$pval[ 2 ], df = x$pooxm$gxe.test$df[ 2 ] ) )
} )
pooxm.pval.enhancer <- do.call( rbind, pooxm.pval.enhancer )
pooxm.pval.gene <- lapply( gxe.gene.all, function(x){
return( data.frame( pval = x$pooxm$gxe.test$pval[ 2 ], df = x$pooxm$gxe.test$df[ 2 ] ) )
} )
pooxm.pval.gene <- do.call( rbind, pooxm.pval.gene )
# creating the data.frame
my.SNP.map$Marker <- tolower( my.SNP.map$Marker )
cur.snp.names <- my.SNP.map$Marker[ match( c( rownames( pooxm.pval.promoter ),
rownames( pooxm.pval.enhancer ), rownames( pooxm.pval.gene ) ),
my.SNP.map$Marker ) ]
pooxm.pval.all <- data.frame( region = c( rep( "promoter", nrow( pooxm.pval.promoter ) ),
rep( "enhancer", nrow( pooxm.pval.enhancer ) ), rep( "gene", nrow( pooxm.pval.gene ) ) ),
gene.name = c( rownames( pooxm.pval.promoter ), rownames( pooxm.pval.enhancer ), rownames( pooxm.pval.gene ) ),
snp.name = cur.snp.names,
pval = unlist( c( pooxm.pval.promoter[ , "pval" ], pooxm.pval.enhancer[ , "pval" ], pooxm.pval.gene[ , "pval" ] ) ),
df = unlist( c( pooxm.pval.promoter[ , "df" ], pooxm.pval.enhancer[ , "df" ], pooxm.pval.gene[ , "df" ] ) ) )
head( pooxm.pval.all )
## region gene.name snp.name pval df
## 1 promoter rs4621895 rs4621895 0.72807227 2
## 2 promoter rs1546124 rs1546124 0.53742222 2
## 3 promoter rs3750817 rs3750817 0.56388313 2
## 4 promoter rs2981428 rs2981428 0.84228727 2
## 5 promoter rs10984103 rs10984103 0.01005862 2
## 6 promoter rs9939609 rs9939609 0.21176658 2
dim( pooxm.pval.all )
## [1] 33 5
df.all.pooxm <- unique( pooxm.pval.all$df )
df.all.pooxm
## [1] 2
# plotting
pvals.POOxM.random <- pooxm.pval.all$pval
names( pvals.POOxM.random ) <- paste( pooxm.pval.all$region, pooxm.pval.all$snp.name, sep = "_" )
pvals.POOxM.random
## promoter_rs4621895 promoter_rs1546124 promoter_rs3750817
## 0.72807227 0.53742222 0.56388313
## promoter_rs2981428 promoter_rs10984103 promoter_rs9939609
## 0.84228727 0.01005862 0.21176658
## promoter_rs6010718 promoter_rs3752462 promoter_rs2867125
## 0.34162379 0.05130543 0.61980121
## promoter_rs2912771 promoter_rs6659735 promoter_rs10767664
## 0.06483495 0.61302723 0.70609158
## promoter_rs7138803 promoter_rs713586 promoter_rs1443433
## 0.02368679 0.37280208 0.28315669
## promoter_rs4752566 enhancer_rs3750817 enhancer_rs2981428
## 0.21296224 0.38827109 0.86498538
## enhancer_rs2912771 enhancer_rs6659735 enhancer_rs7138803
## 0.56378499 0.63643147 0.52693061
## enhancer_rs713586 gene_rs4621895 gene_rs1546124
## 0.23641855 0.92409800 0.54578404
## gene_rs3750817 gene_rs2981428 gene_rs9939609
## 0.11088681 0.06325315 0.20866275
## gene_rs6010718 gene_rs3752462 gene_rs2912771
## 0.70173975 0.40685657 0.61717747
## gene_rs6659735 gene_rs10767664 gene_rs4752566
## 0.92562959 0.63967902 0.90081823
cur.len <- length( pvals.POOxM.random )
cur.nlabs <- 4
cur.xcoords <- get_x_coords( length( pvals.POOxM.random ), cur.nlabs )
cur.ycoords <- -log10( sort( pvals.POOxM.random ) )[ 1:cur.nlabs ]
cur.labels <- data.frame( x = cur.xcoords, y = cur.ycoords,
label = names( cur.ycoords ) )
gg_qqplot( pvals.POOxM.random ) +
labs( title = "PoOxMe, random SNPs, CL/P dataset" ) +
geom_label_repel( data = cur.labels, aes( x, y, label = label ), xlim = c( 1,NA ) ) +
theme( axis.text = element_text( size = 12 ), axis.title = element_text( size = 14 ),
title = element_text( size = 16 ) )
## Warning: Removed 2 rows containing missing values (geom_path).
ggsave( filename = "PoOxM_random_CLorP_mult_betas_gg-qqplot_square.png", width = 14, units = "cm",
height = 14 )
## Warning: Removed 2 rows containing missing values (geom_path).
# ==== TOP 20 SNPs FROM PoO SCAN on CLP ====
load( "../PoO_scan/GxM_PoOxM_new_betas_mult_top-SNPs-PoO-scan_CLP.RData" )
# ==== GxM ====
# GxM
gxe.pval.promoter.scan.top20 <- lapply( gxe.prom.all, function(x){
return( data.frame( pval = x$gxm$pval[ 2 ], df = x$gxm$df[ 2 ] ) )
} )
gxe.pval.promoter.scan.top20 <- do.call( rbind, gxe.pval.promoter.scan.top20 )
gxe.pval.enhancer.scan.top20 <- lapply( gxe.enh.all, function(x){
return( data.frame( pval = x$gxm$gxe.test$pval[ 2 ], df = x$gxm$gxe.test$df[ 2 ] ) )
} )
gxe.pval.enhancer.scan.top20 <- do.call( rbind, gxe.pval.enhancer.scan.top20 )
gxe.pval.gene.scan.top20 <- lapply( gxe.gene.all, function(x){
return( data.frame( pval = x$gxm$gxe.test$pval[ 2 ], df = x$gxm$gxe.test$df[ 2 ] ) )
} )
gxe.pval.gene.scan.top20 <- do.call( rbind, gxe.pval.gene.scan.top20 )
# creating the data.frame
cur.snp.names <- as.character( my.SNP.map$Marker[ match( c( rownames( gxe.pval.promoter.scan.top20 ),
rownames( gxe.pval.enhancer.scan.top20 ), rownames( gxe.pval.gene.scan.top20 ) ),
my.SNP.map$Marker ) ] )
gxe.pval.all.scan.top20 <- data.frame(
region = c( rep( "promoter", nrow( gxe.pval.promoter.scan.top20 ) ),
rep( "enhancer", nrow( gxe.pval.enhancer.scan.top20 ) ),
rep( "gene", nrow( gxe.pval.gene.scan.top20 ) ) ),
gene.name = c( rownames( gxe.pval.promoter.scan.top20 ),
rownames( gxe.pval.enhancer.scan.top20 ),
rownames( gxe.pval.gene.scan.top20 ) ),
snp.name = cur.snp.names,
pval = unlist( c( gxe.pval.promoter.scan.top20[ , "pval" ],
gxe.pval.enhancer.scan.top20[ , "pval" ],
gxe.pval.gene.scan.top20[ , "pval" ] ) ),
df = unlist( c( gxe.pval.promoter.scan.top20[ , "df" ],
gxe.pval.enhancer.scan.top20[ , "df" ],
gxe.pval.gene.scan.top20[ , "df" ] ) ) )
head( gxe.pval.all.scan.top20 )
## region gene.name snp.name pval df
## 1 promoter rs29941 rs29941 0.82654721 2
## 2 promoter rs2235371 rs2235371 0.66949676 2
## 3 promoter rs2912760 rs2912760 0.22333362 2
## 4 promoter rs11084753 rs11084753 0.03117635 2
## 5 promoter rs8001641 rs8001641 0.08656249 2
## 6 promoter rs2912771 rs2912771 0.99998991 2
dim( gxe.pval.all.scan.top20 )
## [1] 32 5
df.all.scan.top20 <- unique( gxe.pval.all.scan.top20$df )
df.all.scan.top20
## [1] 2
# plotting
pvals.GxM.scan.top20 <- gxe.pval.all.scan.top20$pval
names( pvals.GxM.scan.top20 ) <- paste( gxe.pval.all.scan.top20$region, gxe.pval.all.scan.top20$snp.name,
sep = "_" )
pvals.GxM.scan.top20
## promoter_rs29941 promoter_rs2235371 promoter_rs2912760
## 0.82654721 0.66949676 0.22333362
## promoter_rs11084753 promoter_rs8001641 promoter_rs2912771
## 0.03117635 0.08656249 0.99998991
## promoter_rs1002246 promoter_rs6701037 promoter_rs766325
## 0.81700052 0.91714460 0.59691883
## promoter_rs4621895 promoter_rs470563 promoter_rs7138803
## 0.32458643 0.15879066 0.97677359
## promoter_rs987525 enhancer_rs4752028 enhancer_rs8001641
## 0.35728352 0.94383682 0.55298564
## enhancer_rs2912771 enhancer_rs7078160 enhancer_rs766325
## 0.99999005 0.60788200 0.97881594
## enhancer_rs470563 enhancer_rs7138803 enhancer_rs987525
## 0.41654277 0.22454267 0.68264461
## enhancer_rs3752075 enhancer_rs7864322 gene_rs2235371
## 0.13738084 0.14358089 0.28213909
## gene_rs4752028 gene_rs2912760 gene_rs2912771
## 0.68510218 0.90893240 0.99999007
## gene_rs1002246 gene_rs7078160 gene_rs4621895
## 0.27991851 0.60788200 0.65649353
## gene_rs470563 gene_rs3752075
## 0.05219897 0.56835826
cur.len <- length( pvals.GxM.scan.top20 )
cur.nlabs <- 4
cur.xcoords <- get_x_coords( length( pvals.GxM.scan.top20 ), cur.nlabs )
cur.ycoords <- -log10( sort( pvals.GxM.scan.top20 ) )[ 1:cur.nlabs ]
cur.labels <- data.frame( x = cur.xcoords, y = cur.ycoords,
label = names( cur.ycoords ) )
gg_qqplot( pvals.GxM.scan.top20 ) +
labs( title = "GxMe, CLP dataset",
subtitle = "top 20 SNPs from PoO scan on CLP dataset" ) +
geom_label_repel( data = cur.labels, aes( x, y, label = label ) ) +
theme( axis.text = element_text( size = 12 ), axis.title = element_text( size = 14 ),
title = element_text( size = 16 ) )
## Warning: Removed 2 rows containing missing values (geom_path).
# ==== POOxM ====
# now, PoOxM for CLP
pooxm.pval.promoter <- lapply( gxe.prom.all, function(x){
return( data.frame( pval = x$pooxm$pval[ 2 ], df = x$pooxm$df[ 2 ] ) )
} )
pooxm.pval.promoter <- do.call( rbind, pooxm.pval.promoter )
pooxm.pval.enhancer <- lapply( gxe.enh.all, function(x){
return( data.frame( pval = x$pooxm$gxe.test$pval[ 2 ], df = x$pooxm$gxe.test$df[ 2 ] ) )
} )
pooxm.pval.enhancer <- do.call( rbind, pooxm.pval.enhancer )
pooxm.pval.gene <- lapply( gxe.gene.all, function(x){
return( data.frame( pval = x$pooxm$gxe.test$pval[ 2 ], df = x$pooxm$gxe.test$df[ 2 ] ) )
} )
pooxm.pval.gene <- do.call( rbind, pooxm.pval.gene )
# creating the data.frame
my.SNP.map$Marker <- tolower( my.SNP.map$Marker )
cur.snp.names <- my.SNP.map$Marker[ match( c( rownames( pooxm.pval.promoter ),
rownames( pooxm.pval.enhancer ), rownames( pooxm.pval.gene ) ),
my.SNP.map$Marker ) ]
pooxm.pval.all <- data.frame( region = c( rep( "promoter", nrow( pooxm.pval.promoter ) ),
rep( "enhancer", nrow( pooxm.pval.enhancer ) ), rep( "gene", nrow( pooxm.pval.gene ) ) ),
gene.name = c( rownames( pooxm.pval.promoter ), rownames( pooxm.pval.enhancer ), rownames( pooxm.pval.gene ) ),
snp.name = cur.snp.names,
pval = unlist( c( pooxm.pval.promoter[ , "pval" ], pooxm.pval.enhancer[ , "pval" ], pooxm.pval.gene[ , "pval" ] ) ),
df = unlist( c( pooxm.pval.promoter[ , "df" ], pooxm.pval.enhancer[ , "df" ], pooxm.pval.gene[ , "df" ] ) ) )
head( pooxm.pval.all )
## region gene.name snp.name pval df
## 1 promoter rs29941 rs29941 0.43921706 2
## 2 promoter rs2235371 rs2235371 0.10088696 2
## 3 promoter rs2912760 rs2912760 0.33704873 2
## 4 promoter rs11084753 rs11084753 0.13590825 2
## 5 promoter rs8001641 rs8001641 0.02341064 2
## 6 promoter rs2912771 rs2912771 0.99999491 2
dim( pooxm.pval.all )
## [1] 32 5
df.all.pooxm <- unique( pooxm.pval.all$df )
df.all.pooxm
## [1] 2
# plotting
pvals.POOxM.scan.top20 <- pooxm.pval.all$pval
names( pvals.POOxM.scan.top20 ) <- paste( pooxm.pval.all$region, pooxm.pval.all$snp.name, sep = "_" )
pvals.POOxM.scan.top20
## promoter_rs29941 promoter_rs2235371 promoter_rs2912760
## 0.4392170565 0.1008869623 0.3370487300
## promoter_rs11084753 promoter_rs8001641 promoter_rs2912771
## 0.1359082533 0.0234106432 0.9999949074
## promoter_rs1002246 promoter_rs6701037 promoter_rs766325
## 0.3174203675 0.9408714893 0.0006464076
## promoter_rs4621895 promoter_rs470563 promoter_rs7138803
## 0.9120701317 0.3110710536 0.7712654563
## promoter_rs987525 enhancer_rs4752028 enhancer_rs8001641
## 0.6038221712 0.5257461092 0.0624808580
## enhancer_rs2912771 enhancer_rs7078160 enhancer_rs766325
## 0.9999949768 0.2489809872 0.8189611982
## enhancer_rs470563 enhancer_rs7138803 enhancer_rs987525
## 0.3622765345 0.4319407359 0.1249379440
## enhancer_rs3752075 enhancer_rs7864322 gene_rs2235371
## 0.6039579732 0.5228906273 0.5629594490
## gene_rs4752028 gene_rs2912760 gene_rs2912771
## 0.0897178994 0.6048225816 0.9999949348
## gene_rs1002246 gene_rs7078160 gene_rs4621895
## 0.4707194915 0.2489809872 0.8953377754
## gene_rs470563 gene_rs3752075
## 0.4853539693 0.0757098844
cur.len <- length( pvals.POOxM.scan.top20 )
cur.nlabs <- 4
cur.xcoords <- get_x_coords( length( pvals.POOxM.scan.top20 ), cur.nlabs )
cur.ycoords <- -log10( sort( pvals.POOxM.scan.top20 ) )[ 1:cur.nlabs ]
cur.labels <- data.frame( x = cur.xcoords, y = cur.ycoords,
label = names( cur.ycoords ) )
gg_qqplot( pvals.POOxM.scan.top20 ) +
labs( title = "PoOxMe, CLP dataset",
subtitle = "top 20 SNPs from PoO scan on CLP dataset" ) +
geom_label_repel( data = cur.labels, aes( x, y, label = label ), xlim = c( 1,NA ) ) +
theme( axis.text = element_text( size = 12 ), axis.title = element_text( size = 14 ),
title = element_text( size = 16 ) )
ggsave( filename = "PoOxM_top20SNPs-PoO-scan_CLP_mult_betas_gg-qqplot_square.png", width = 14, units = "cm",
height = 14 )
# ==== GxM and PoOxM on controls, TOP 20 SNPs FROM PoO SCAN on CLP ====
load( "../PoO_scan/GxM_PoOxM_new_betas_mult_top-SNPs-PoO-scan_cntrl.RData" )
# ==== GxM ====
# GxM
gxe.pval.promoter.cntrl.scan.top20 <- lapply( gxe.prom.all, function(x){
return( data.frame( pval = x$gxm$pval[ 2 ], df = x$gxm$df[ 2 ] ) )
} )
gxe.pval.promoter.cntrl.scan.top20 <- do.call( rbind, gxe.pval.promoter.cntrl.scan.top20 )
gxe.pval.enhancer.cntrl.scan.top20 <- lapply( gxe.enh.all, function(x){
return( data.frame( pval = x$gxm$gxe.test$pval[ 2 ], df = x$gxm$gxe.test$df[ 2 ] ) )
} )
gxe.pval.enhancer.cntrl.scan.top20 <- do.call( rbind, gxe.pval.enhancer.cntrl.scan.top20 )
gxe.pval.gene.cntrl.scan.top20 <- lapply( gxe.gene.all, function(x){
return( data.frame( pval = x$gxm$gxe.test$pval[ 2 ], df = x$gxm$gxe.test$df[ 2 ] ) )
} )
gxe.pval.gene.cntrl.scan.top20 <- do.call( rbind, gxe.pval.gene.cntrl.scan.top20 )
# creating the data.frame
cur.snp.names <- as.character( my.SNP.map$Marker[ match( c( rownames( gxe.pval.promoter.cntrl.scan.top20 ),
rownames( gxe.pval.enhancer.cntrl.scan.top20 ), rownames( gxe.pval.gene.cntrl.scan.top20 ) ),
my.SNP.map$Marker ) ] )
gxe.pval.all.cntrl.scan.top20 <- data.frame(
region = c( rep( "promoter", nrow( gxe.pval.promoter.cntrl.scan.top20 ) ),
rep( "enhancer", nrow( gxe.pval.enhancer.cntrl.scan.top20 ) ),
rep( "gene", nrow( gxe.pval.gene.cntrl.scan.top20 ) ) ),
gene.name = c( rownames( gxe.pval.promoter.cntrl.scan.top20 ),
rownames( gxe.pval.enhancer.cntrl.scan.top20 ),
rownames( gxe.pval.gene.cntrl.scan.top20 ) ),
snp.name = cur.snp.names,
pval = unlist( c( gxe.pval.promoter.cntrl.scan.top20[ , "pval" ],
gxe.pval.enhancer.cntrl.scan.top20[ , "pval" ],
gxe.pval.gene.cntrl.scan.top20[ , "pval" ] ) ),
df = unlist( c( gxe.pval.promoter.cntrl.scan.top20[ , "df" ],
gxe.pval.enhancer.cntrl.scan.top20[ , "df" ],
gxe.pval.gene.cntrl.scan.top20[ , "df" ] ) ) )
head( gxe.pval.all.cntrl.scan.top20 )
## region gene.name snp.name pval df
## 1 promoter rs29941 rs29941 0.1372708 2
## 2 promoter rs2235371 rs2235371 0.6699029 2
## 3 promoter rs2912760 rs2912760 0.3141108 2
## 4 promoter rs11084753 rs11084753 0.3286412 2
## 5 promoter rs8001641 rs8001641 0.2137793 2
## 6 promoter rs2912771 rs2912771 0.8798302 2
dim( gxe.pval.all.cntrl.scan.top20 )
## [1] 32 5
df.all.cntrl.scan.top20 <- unique( gxe.pval.all.cntrl.scan.top20$df )
df.all.cntrl.scan.top20
## [1] 2
# plotting
pvals.GxM.cntrl.scan.top20 <- gxe.pval.all.cntrl.scan.top20$pval
names( pvals.GxM.cntrl.scan.top20 ) <- paste( gxe.pval.all.cntrl.scan.top20$region, gxe.pval.all.cntrl.scan.top20$snp.name,
sep = "_" )
pvals.GxM.cntrl.scan.top20
## promoter_rs29941 promoter_rs2235371 promoter_rs2912760
## 0.137270790 0.669902904 0.314110830
## promoter_rs11084753 promoter_rs8001641 promoter_rs2912771
## 0.328641169 0.213779268 0.879830202
## promoter_rs1002246 promoter_rs6701037 promoter_rs766325
## 0.845813433 0.766016661 0.601984456
## promoter_rs4621895 promoter_rs470563 promoter_rs7138803
## 0.226496387 0.862097102 0.552290095
## promoter_rs987525 enhancer_rs4752028 enhancer_rs8001641
## 0.394732285 0.008949105 0.756633977
## enhancer_rs2912771 enhancer_rs7078160 enhancer_rs766325
## 0.806504578 0.936330519 0.369985147
## enhancer_rs470563 enhancer_rs7138803 enhancer_rs987525
## 0.573976931 0.305741501 0.033988270
## enhancer_rs3752075 enhancer_rs7864322 gene_rs2235371
## 0.513942721 0.023297383 0.128936886
## gene_rs4752028 gene_rs2912760 gene_rs2912771
## 0.209120232 0.948572340 0.776234184
## gene_rs1002246 gene_rs7078160 gene_rs4621895
## 0.201886639 0.936330519 0.146912777
## gene_rs470563 gene_rs3752075
## 0.516586450 0.026807263
cur.len <- length( pvals.GxM.cntrl.scan.top20 )
cur.nlabs <- 4
cur.xcoords <- get_x_coords( length( pvals.GxM.cntrl.scan.top20 ), cur.nlabs )
cur.ycoords <- -log10( sort( pvals.GxM.cntrl.scan.top20 ) )[ 1:cur.nlabs ]
cur.labels <- data.frame( x = cur.xcoords, y = cur.ycoords,
label = names( cur.ycoords ) )
gg_qqplot( pvals.GxM.cntrl.scan.top20 ) +
labs( title = "GxMe, control dataset",
subtitle = "top 20 SNPs from PoO scan on CLP dataset" ) +
geom_label_repel( data = cur.labels, aes( x, y, label = label ) ) +
theme( axis.text = element_text( size = 12 ), axis.title = element_text( size = 14 ),
title = element_text( size = 16 ) )
## Warning: Removed 1 rows containing missing values (geom_path).
# ==== POOxM ====
# now, PoOxM for controls
pooxm.pval.promoter <- lapply( gxe.prom.all, function(x){
return( data.frame( pval = x$pooxm$pval[ 2 ], df = x$pooxm$df[ 2 ] ) )
} )
pooxm.pval.promoter <- do.call( rbind, pooxm.pval.promoter )
pooxm.pval.enhancer <- lapply( gxe.enh.all, function(x){
return( data.frame( pval = x$pooxm$gxe.test$pval[ 2 ], df = x$pooxm$gxe.test$df[ 2 ] ) )
} )
pooxm.pval.enhancer <- do.call( rbind, pooxm.pval.enhancer )
pooxm.pval.gene <- lapply( gxe.gene.all, function(x){
return( data.frame( pval = x$pooxm$gxe.test$pval[ 2 ], df = x$pooxm$gxe.test$df[ 2 ] ) )
} )
pooxm.pval.gene <- do.call( rbind, pooxm.pval.gene )
# creating the data.frame
my.SNP.map$Marker <- tolower( my.SNP.map$Marker )
cur.snp.names <- my.SNP.map$Marker[ match( c( rownames( pooxm.pval.promoter ),
rownames( pooxm.pval.enhancer ), rownames( pooxm.pval.gene ) ),
my.SNP.map$Marker ) ]
pooxm.pval.all <- data.frame( region = c( rep( "promoter", nrow( pooxm.pval.promoter ) ),
rep( "enhancer", nrow( pooxm.pval.enhancer ) ), rep( "gene", nrow( pooxm.pval.gene ) ) ),
gene.name = c( rownames( pooxm.pval.promoter ), rownames( pooxm.pval.enhancer ), rownames( pooxm.pval.gene ) ),
snp.name = cur.snp.names,
pval = unlist( c( pooxm.pval.promoter[ , "pval" ], pooxm.pval.enhancer[ , "pval" ], pooxm.pval.gene[ , "pval" ] ) ),
df = unlist( c( pooxm.pval.promoter[ , "df" ], pooxm.pval.enhancer[ , "df" ], pooxm.pval.gene[ , "df" ] ) ) )
head( pooxm.pval.all )
## region gene.name snp.name pval df
## 1 promoter rs29941 rs29941 0.24115969 2
## 2 promoter rs2235371 rs2235371 0.08800383 2
## 3 promoter rs2912760 rs2912760 0.41365840 2
## 4 promoter rs11084753 rs11084753 0.99976080 2
## 5 promoter rs8001641 rs8001641 0.79689498 2
## 6 promoter rs2912771 rs2912771 0.62263364 2
dim( pooxm.pval.all )
## [1] 32 5
df.all.pooxm <- unique( pooxm.pval.all$df )
df.all.pooxm
## [1] 2
# plotting
pvals.POOxM.cntrl.scan.top20 <- pooxm.pval.all$pval
names( pvals.POOxM.cntrl.scan.top20 ) <- paste( pooxm.pval.all$region, pooxm.pval.all$snp.name, sep = "_" )
pvals.POOxM.cntrl.scan.top20
## promoter_rs29941 promoter_rs2235371 promoter_rs2912760
## 0.241159688 0.088003825 0.413658402
## promoter_rs11084753 promoter_rs8001641 promoter_rs2912771
## 0.999760796 0.796894975 0.622633637
## promoter_rs1002246 promoter_rs6701037 promoter_rs766325
## 0.373770725 0.631668512 0.442064844
## promoter_rs4621895 promoter_rs470563 promoter_rs7138803
## 0.847711773 0.262633849 0.885167826
## promoter_rs987525 enhancer_rs4752028 enhancer_rs8001641
## 0.100198269 0.647343113 0.992102507
## enhancer_rs2912771 enhancer_rs7078160 enhancer_rs766325
## 0.845707298 0.204245362 0.771067726
## enhancer_rs470563 enhancer_rs7138803 enhancer_rs987525
## 0.672004643 0.572507479 0.004942489
## enhancer_rs3752075 enhancer_rs7864322 gene_rs2235371
## 0.079418842 0.340948665 0.292493187
## gene_rs4752028 gene_rs2912760 gene_rs2912771
## 0.961928933 0.252333479 0.796431053
## gene_rs1002246 gene_rs7078160 gene_rs4621895
## 0.719148045 0.204245362 0.824632091
## gene_rs470563 gene_rs3752075
## 0.540884050 0.470202890
cur.len <- length( pvals.POOxM.cntrl.scan.top20 )
cur.nlabs <- 4
cur.xcoords <- get_x_coords( length( pvals.POOxM.cntrl.scan.top20 ), cur.nlabs )
cur.ycoords <- -log10( sort( pvals.POOxM.cntrl.scan.top20 ) )[ 1:cur.nlabs ]
cur.labels <- data.frame( x = cur.xcoords, y = cur.ycoords,
label = names( cur.ycoords ) )
gg_qqplot( pvals.POOxM.cntrl.scan.top20 ) +
labs( title = "PoOxMe, control dataset",
subtitle = "top 20 SNPs from PoO scan on CLP dataset" ) +
geom_label_repel( data = cur.labels, aes( x, y, label = label ), xlim = c( 1.5,NA ) ) +
theme( axis.text = element_text( size = 12 ), axis.title = element_text( size = 14 ),
title = element_text( size = 16 ) )
## Warning: Removed 1 rows containing missing values (geom_path).
ggsave( filename = "PoOxM_top20SNPs-PoO-scan_CLP_cntrl_mult_betas_gg-qqplot_square.png",
width = 14, units = "cm", height = 14 )
## Warning: Removed 1 rows containing missing values (geom_path).