Contents

library(ggbio)
library(dplyr)
library(IHW)
library(fdrtool)
library(cowplot)
library(tidyr)
library(scales)
library(latex2exp)

Let us start by loading in the data:

file_loc <- system.file("extdata","real_data", "hqtl_chrom1_chrom2", package = "IHWpaper")

First the two tables with the p-values corresponding to the two chromosomes. Note that only p-values <= 1e-4 are stored in these.

chr1_df <- readRDS(file.path(file_loc, "chr1_subset.Rds"))
chr2_df <- readRDS(file.path(file_loc, "chr2_subset.Rds"))
pval_threshold <- 10^(-4)

Also recall each hypothesis corresponds to a peak (which we call gene below) and a SNP. Hence let us load files about each of the SNPs and peaks:

snp_chr1 <- readRDS(file.path(file_loc, "snppos_chr1.Rds"))
snp_chr2 <-  readRDS(file.path(file_loc, "snppos_chr2.Rds"))

all_peaks <- readRDS(file.path(file_loc, "peak_locations.Rds"))
peaks_chr1 <- dplyr::filter(all_peaks, chr=="chr1")
peaks_chr2 <- dplyr::filter(all_peaks, chr=="chr2")

We can use these both to infer how many hypotheses were conducted in total (or at a given distance), but also to calculate our covariates which are a function of SNP and peak (their distance).

Now let us attach the new column with the covariate (distance) to the data frames.

chr1_df <- left_join(chr1_df, select(snp_chr1, snp, pos), by=(c("SNP"="snp"))) %>%
       left_join(peaks_chr1, by=(c("gene"="id"))) %>%
       mutate( dist = pmin( abs(pos-start), abs(pos-end)))

chr2_df <- left_join(chr2_df, select(snp_chr2, snp, pos), by=(c("SNP"="snp"))) %>%
       left_join(peaks_chr2, by=(c("gene"="id"))) %>%
       mutate( dist = pmin( abs(pos-start), abs(pos-end)))

Now let us convert the distance to a categorical covariate by binning:

my_breaks <- c(-1, 
                 seq(from=10000,to=290000, by=10000) , 
                 seq(from=300000, to=0.9*10^6, by=100000),
                 seq(from=10^6, to=25.1*10^7, by=10^7))
myf1 <- cut(chr1_df$dist, my_breaks)
myf2 <- cut(chr2_df$dist, my_breaks)

To apply our method despite the fact that only small p-values are available, we will count how many hypotheses there are in each of the bins. The above code is not very efficient, so we have precomputed these and do not run the below chunk.

cnt = 0
ms <- ms*0
pb = txtProgressBar(min = 0, max = nrow(peaks_chr1), initial = 0) 

for (i in 1:nrow(peaks_chr1)){
  setTxtProgressBar(pb,i)
  start_pos <- peaks_chr1$start[i]
  end_pos <- peaks_chr1$end[i]
  dist_vec <- pmin( abs(snp_chr1$pos - start_pos), abs(snp_chr1$pos - end_pos) )
  ms <- ms + table(cut(dist_vec, my_breaks))
}

saveRDS( ms, file = "m_groups_chr1.Rds" )


cnt = 0
ms_chr2 <- table(myf2)*0
pb = txtProgressBar(min = 0, max = nrow(peaks_chr2), initial = 0) 

for (i in 1:nrow(peaks_chr2)){
  setTxtProgressBar(pb,i)
  start_pos <- peaks_chr1$start[i]
  end_pos <- peaks_chr1$end[i]
  dist_vec <- pmin( abs(snp_chr2$pos - start_pos), abs(snp_chr2$pos - end_pos) )
  ms_chr2 <- ms_chr2 + table(cut(dist_vec, my_breaks))
}

saveRDS( ms_chr2, file = "m_groups_chr2.Rds" )

Let us load the result from the above execution:

ms_chr1 <- readRDS(file.path(file_loc, "m_groups_chr1.Rds"))
ms_chr2 <- readRDS(file.path(file_loc, "m_groups_chr2.Rds"))

Let us put the data for the two chromosomes together:

chr1_chr2_df <- rbind(chr1_df, chr2_df)
chr1_chr2_groups <- as.factor(c(myf1,myf2))
folds_vec <- as.factor(c(rep(1, nrow(chr1_df)), rep(2, nrow(chr2_df))))
m_groups <- cbind(ms_chr1, ms_chr2)
m <- sum(m_groups) #total number of hypotheses
m
## [1] 15725016812

1 Histogram plots

Get our colors:

beyonce_colors <- c("#b72da0", "#7c5bd2", "#0097ed","#00c6c3",
                   "#9cd78a", "#f7f7a7", "#ebab5f", "#e24344",
                   "#04738d")#,"#d8cdc9")
beyonce_colors[6] <- c("#dbcb09") # thicker yellow
pretty_colors <- beyonce_colors[c(2,1,3:5)]
set.seed(1)
n_subsample <- 10000
chr1_chr2_df_sub <- sample_n(filter(chr1_chr2_df, pvalue <= 1e-7),n_subsample)
qs <- c(0.25,0.5, 0.75)
cutoffs <- c(0, quantile(chr1_chr2_df_sub$dist,qs), Inf)
cov_scatter_gg <- ggplot(chr1_chr2_df_sub, aes(x=rank(dist),y=-log10(pvalue))) +
     geom_point(alpha=0.2, col=pretty_colors[1]) + 
     geom_vline(xintercept=qs*n_subsample, linetype="dashed") + 
     ylab(expression(paste(-log[10],"(p-value)"))) +
     xlab(expression(paste("Rank of distance")))
cov_scatter_gg