Volcano plots

bootstrap.R

# bootstrap.R
cat("Loading packages\n")
library(tidyverse)
library(RefFreeEWAS)
library(qvalue)

cat("Loading data\n")
load("pdat.merge.rdata")
load("betas.rdata")

cat("Filtering data\n")
phenotypes <- pdat.merge[, -10]  %>% # drop extra 'Sex' column
    rownames_to_column("patient") %>% # add patient column
    filter(Tissue == "BUFFY") %>% # only blood samples
    filter(!is.na(PSAS1) & !is.na(Education) & !is.na(income))
save(phenotypes, files="phenotypes.rdata")

measurements <- betas[, phenotypes$patient] # this is not necessary
cat("Using", nrow(phenotypes), "patients and", nrow(measurements), "markers\n")
save(measurements, files="measurements.rdata")

cat("Determening degrees of freedom\n")
rmt <- EstDimRMT(measurements)
cat("Estimated dof for markers", rmt$dim, "\n")

cat("Building model\n");
m <- model.matrix(patient ~ (PSAS1 * Sex) + MALE_RecentEvent_Score + income + Education + Ethnicity,
		data = phenotypes)

rfm <- RefFreeEwasModel(measurements, m, K = rmt$dim)

n_bootstrap <- 100
cat("Bootstrapping n=", n_bootstrap, "\n")
b <- BootRefFreeEwasModel(rfm, nboot = n_bootstrap)

cat("Summarizing bootstraps\n")
res <- summary(b)
save(res, file = paste0("bootstrap_", n_bootstrap,".rdata"))
cat("Done\n")

volcano.R

volcano_plot <- function(bootstrap_summary, variable) {
    mean_beta <- res[, variable, "B", "mean"]
    sd_beta <- res[, variable, "B", "sd"]
    p_value <- pt(mean_beta / sd_beta, dof1)
    q_value <- qvalue(p_value)

    jpeg(paste0("volcano_", variable, ".jpg"), width = 1000, height = 1000)
    
    print(data_frame(
        gene = rownames(measurements), 
        variable = mean_beta, 
        p_value = p_value,
        q_value = q_value$qvalues,
        logq = -log10(q_value),
        significant = logq > 2) %>%
        ggplot(aes(x = variable, y = logq)) + 
        geom_point(aes(color = significant)) +
        geom_text(aes(label = ifelse(significant, gene, '')), vjust = 0, hjust = 0, angle = 45) +
        scale_color_brewer(palette = "Set1"))
    dev.off()
}

# load `res` from bootstraping:
load("bootstrap_100.rdata")
# load other data
load("phenotypes.rdata")
load("measurements.rdata")

volcano_plot(res, "PSAS1:Sex")