Skip to content

Instantly share code, notes, and snippets.

@ivan-krukov
Last active June 5, 2017 19:23
Show Gist options
  • Save ivan-krukov/0fb0622a864e3592f022fdd4b314346f to your computer and use it in GitHub Desktop.
Save ivan-krukov/0fb0622a864e3592f022fdd4b314346f to your computer and use it in GitHub Desktop.
Volcano plots
# 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_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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment