Age of menarche index of suspicion

ios_aom.R

library(TwoSampleMR)
library(MRInstruments)
library(ggplot2)
library(dplyr)
library(reshape2)
library(classInt)


calculate_ios <- function(exposure_dat, outcome_dat)
{
	require(dplyr)
	require(reshape2)
	require(classInt)
	ios <- dplyr::group_by(outcome_dat, SNP) %>%
		dplyr::summarise(
			ios_mean = sum(vgu, na.rm=TRUE),
			ios_sd = sd(vgu, na.rm=TRUE),
			ios_iqr = quantile(vgu, 0.75, na.rm=TRUE) - quantile(vgu, 0.25, na.rm=TRUE),
			ios_median = median(vgu, na.rm=TRUE),
			ios_95 = quantile(vgu, 0.95, na.rm=TRUE),
			ios_max = max(vgu, na.rm=TRUE)
		) %>% melt
	ios <- merge(ios, subset(exposure_dat, select=c(SNP, gene.exposure, vgx)), by="SNP")
	ios$ios_ratio <- ios$value / ios$vgx

	W <- sum(outcome_dat$se.outcome^2)

	ios <- group_by(ios, variable) %>%
		do({
			x <- as.data.frame(., stringsAsFactors=FALSE)
			x$value <- x$value / sum(x$value) * W
			x$ios_ratio <- x$ios_ratio / sum(x$ios_ratio) * W
			return(x)
		})

	ios <- dplyr::group_by(ios, variable) %>%
		dplyr::mutate(
			jenks = findCols(classIntervals(value, n=2, style="jenks")),
			jenks_ratio = findCols(classIntervals(ios_ratio, n=2, style="jenks"))
		)

	ios$lab <- ios$gene.exposure
	ios <- group_by(ios, variable) %>%
	do({
		.$lab[.$value < quantile(.$value, 0.9) & .$vgx < quantile(.$vgx, 0.9)] <- NA
		return(.)
	})
	return(ios)
}

####

# Find 122 AOM SNPs
data(gwas_catalog)
a <- subset(gwas_catalog, grepl("menarche", Phenotype, ignore.case=TRUE) & PubmedID == "25231870")
aom_exposure <- format_data(a)

# Get 452 metabolite traits from Shin et al
ao <- available_outcomes()
shin <- subset(ao, grepl("Shin", author))

# Extracting 122 SNPs from each of 452 metabolutes
# This takes about an hour or so
aom_shin <- extract_outcome_data(aom_exposure$SNP, shin$id)

# Calculate Vg for each SNP-trait association
aom_exposure$vgx <- aom_exposure$beta.exposure^2 * 2 * aom_exposure$eaf.exposure * (1 - aom_exposure$eaf.exposure)
aom_shin$vgu <- aom_shin$beta.outcome^2 * 2 * aom_shin$eaf.outcome * (1 - aom_shin$eaf.outcome)


save(aom_exposure, aom_shin, file="aom_shin.rdata")


####

load("aom_shin.rdata")


aom_ios <- calculate_ios(aom_exposure, aom_shin)

a <- group_by(aom_shin, SNP) %>% summarise(vgu=median(vgu))
b <- merge(a, aom_exposure, by="SNP")
ggplot(b, aes(x=gene.exposure, y=vgu)) + geom_point() + labs(x="Urate gene", y="Index of suspicion (metabolite source)") + theme(axis.text.x=element_text(angle=90))

ggplot(aom_ios, aes(x=gene.exposure, y=value)) +
geom_point(aes(colour=variable)) + theme(axis.text.x=element_text(angle=90, hjust=0.5, vjust=0.5)) +
labs(x="Age of menarche gene", y="Index of suspicion (metabolite source", colour="IoS method")