Plots and models of variant effect sizes

make_evidence_ratio_plot <- function(ERs, ymax, ylab){
  # Argument needs to be a dataframe of loci or transcripts, each identified by a column called "identifier"
  # There should be a col called "evidence_ratio", calculated from the LFSR from ED mashr. The "trait" col says which type of ER it is.
  
  ERs$trait[ERs$trait == "inter_sex_early"] <- "Inter-sex (early life)"
  ERs$trait[ERs$trait == "inter_sex_late"] <- "Inter-sex (late life)"
  ERs$trait[ERs$trait == "inter_age_females"] <- "Inter-age (females)"
  ERs$trait[ERs$trait == "inter_age_males"] <- "Inter-age (males)"
  
  antagonism_evidence_ratios <- ERs %>%
    mutate(trait = factor(trait, c("Inter-sex (early life)",
                                   "Inter-sex (late life)",
                                   "Inter-age (females)",
                                   "Inter-age (males)")))
  
  antagonism_evidence_ratios %>%
    ggplot(aes(log2(evidence_ratio))) + 
    geom_histogram(data=subset(antagonism_evidence_ratios, evidence_ratio < 1), 
                   bins = 500, fill = "#FF635C") + 
    geom_histogram(data=subset(antagonism_evidence_ratios, evidence_ratio > 1), 
                   bins = 500, fill = "#5B8AFD") +
    coord_cartesian(xlim = c(-10, 10), ylim = c(0, ymax)) + 
    scale_x_continuous(breaks = c(-10,  -6,  -2,   2,   6,  10), 
                       labels = c(paste("1/",2 ^ c(10,  6,  2), sep = ""), 2 ^ c(2,6,10))) + 
    facet_wrap(~ trait) + 
    xlab("Evidence ratio (log2 scale)") + ylab(ylab) + 
    theme_bw() + 
    theme(panel.border = element_rect(size = 0.8),
          text = element_text(family = "Raleway", size = 12),
          axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
          strip.background = element_blank())
}

GWAS_ratios_plot <- antagonism_evidence_ratios_gwas %>% 
  select(SNP_clump,  starts_with("inter")) %>%
  distinct() %>% # Include SNPs in 100% LD a single time, such that these clumps are the unit of replication
  gather(trait, evidence_ratio, -SNP_clump) %>% 
  rename(identifier = SNP_clump) %>% 
  make_evidence_ratio_plot(ymax = 5000, ylab = "Number of loci")

# For counting the number of transcripts with 50x more evidence for antagonism 
# antagonism_evidence_ratios_twas %>% 
#   select(FBID,  starts_with("inter")) %>%
#   gather(trait, evidence_ratio, -FBID) %>% 
#   rename(identifier = FBID) %>% filter(evidence_ratio < 1/50) %>% summarise(n())

TWAS_ratios_plot <- antagonism_evidence_ratios_twas %>% 
  select(FBID,  starts_with("inter")) %>%
  gather(trait, evidence_ratio, -FBID) %>% 
  rename(identifier = FBID) %>% 
  make_evidence_ratio_plot(ymax = 500, ylab = "Number of transcripts")

pp <- plot_grid(GWAS_ratios_plot, TWAS_ratios_plot, 
                   labels = c('A', 'B'), label_size = 12)

ggsave("figures/fig5_antagonism_ratios.pdf", pp, width = 8.5, height = 4.9)

pp

# gwas models

# Annotate it with MAF, site, class, etc. 
dat <- left_join(antagonism_evidence_ratios_gwas,
          collect(tbl(db, "variants")), by = "SNP") %>% 
  distinct() %>% mutate(site.class = str_to_title(site.class))

# Focus only on the commonest site classes:
dat <- dat %>% 
  filter(site.class %in% c("Intron", "Intergenic", "Downstream", "Upstream",
                           "Synonymous_coding", "Non_synonymous_coding", "Utr_3_prime",
                           "Exon", "Utr_5_prime")) 

# Remove chromosome 4 (too few sites)
dat <- dat %>% filter(chr != "4") 

# If there are multiple site classes for a SNP, or multiple SNPs 
# in the same 100% LD clump, pick one SNP and/or 1 site class at random
set.seed(1)
dat <- dat[sample(nrow(dat), nrow(dat)), ] %>%
  split(.$SNP_clump) %>%
  map_df(~ .x[1, ])

dat <- dat %>% arrange(site.class)
dat$site.class <- relevel(factor(dat$site.class), ref = "Synonymous_coding")
dat$chr <- relevel(factor(dat$chr), ref = "X")

n_loci <- prettyNum(nrow(dat), big.mark = ",", scientific = FALSE)


intersex_early_model_gwas <- glm(antagonistic ~ chr + MAF + site.class, 
             data = dat %>% 
               mutate(antagonistic = as.numeric(inter_sex_early <= 
                        quantile(dat$inter_sex_early, probs=0.01))), 
    family = "binomial")
car::Anova(intersex_early_model_gwas, type = 3)

intersex_early_model_gwas <- update(intersex_early_model_gwas,  ~ . -site.class)
car::Anova(intersex_early_model_gwas, type = 3)

summary(intersex_early_model_gwas)

intersex_late_model_gwas <- glm(antagonistic ~ chr + MAF + site.class, 
             data = dat %>% 
               mutate(antagonistic = inter_sex_late <= 
                        quantile(dat$inter_sex_late, probs=0.01)), 
    family = "binomial")
car::Anova(intersex_late_model_gwas, type = 3)

intersex_late_model_gwas <- update(intersex_late_model_gwas,  ~ . -site.class)
car::Anova(intersex_late_model_gwas, type = 3)

summary(intersex_late_model_gwas)

# twas models
n_transcripts <- prettyNum(nrow(antagonism_evidence_ratios_twas), big.mark = ",", scientific = FALSE)

twas_model_data <- antagonism_evidence_ratios_twas %>% 
  mutate(antagonistic = inter_sex_early <= 
           quantile(antagonism_evidence_ratios_twas$inter_sex_early, probs=0.01),
         absolute_sex_bias_in_expression = abs(male_bias_in_expression),
         mean_expression_level = AveExpr)


intersex_early_twas <- glm(antagonistic ~ chromosome + absolute_sex_bias_in_expression + mean_expression_level, 
             data = twas_model_data, 
    family = "binomial")
car::Anova(intersex_early_twas, type = 3)

summary(intersex_early_twas)

intersex_late_model_twas <- glm(antagonistic ~ chromosome + absolute_sex_bias_in_expression + mean_expression_level, 
             data = antagonism_evidence_ratios_twas %>% 
  mutate(antagonistic = inter_sex_late <= 
           quantile(antagonism_evidence_ratios_twas$inter_sex_late, probs=0.01),
         absolute_sex_bias_in_expression = abs(male_bias_in_expression),
         mean_expression_level = AveExpr), 
    family = "binomial")
car::Anova(intersex_late_model_twas, type = 3)

intersex_late_model_twas <- update(intersex_late_model_twas,  ~ . -absolute_sex_bias_in_expression -chromosome)
car::Anova(intersex_late_model_twas, type = 3)

summary(intersex_late_model_twas)

# plotting models
new <- data.frame(MAF = seq(from = 0.05, to = 0.5, by = 0.01), chr = "X")
inv_logit <- function(x) exp(x)/(1+exp(x))
preds <- data.frame(new, prop = predict(intersex_early_model_gwas, newdata = new, se.fit=T, type = "link")) %>% 
  mutate(y = prop.fit, lwr = prop.fit - 1.96*prop.se.fit, upr = prop.fit + 1.96*prop.se.fit,
         y = inv_logit(y), lwr = inv_logit(lwr), upr = inv_logit(upr))
p1 <- preds %>% 
  mutate(panel = "A.") %>% 
  ggplot(aes(MAF,  y)) + geom_line() + 
  geom_hline(yintercept = .01, linetype = 3) + 
  geom_line(aes(y = upr), linetype = 2) + 
  geom_line(aes(y = lwr), linetype = 2) + 
  scale_x_continuous(limits = c(0, 0.5)) + 
  scale_y_continuous(limits = c(0, max(preds$y))) + 
  facet_wrap( ~ panel) + 
  theme_bw() + 
  theme(text = element_text(family = "Raleway", size = 12)) +
  theme(strip.background = element_blank(), strip.text = element_text(hjust = 0)) +
  ylab("Probability of being in the top 1%\nsexually antagonistic loci (\u00B1 95% CIs)") +
  xlab("Minor allele frequency")

new2 <- data.frame(chr = levels(dat$chr), MAF = 0.3)
inv_logit <- function(x) exp(x)/(1+exp(x))
preds2 <- data.frame(new2, prop = predict(intersex_early_model_gwas, newdata = new2, se.fit=T, type = "link")) %>% 
  mutate(y = prop.fit, lwr = prop.fit - 1.96*prop.se.fit, upr = prop.fit + 1.96*prop.se.fit,
         y = inv_logit(y), lwr = inv_logit(lwr), upr = inv_logit(upr))
p2 <- preds2 %>% 
  mutate(panel = "B.") %>%
  ggplot(aes(chr,  y)) + 
  geom_point(size = 3) + 
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0) + 
  facet_wrap( ~ panel) + 
  theme_bw() +
  theme(text = element_text(family = "Raleway", size = 12)) +
  theme(strip.background = element_blank(), strip.text = element_text(hjust = 0)) +
  ylab(" \n ") +
  xlab("Chromosome arm")

new <- data.frame(
  mean_expression_level = seq(from = min(twas_model_data$mean_expression_level), 
                                to = max(twas_model_data$mean_expression_level), by = 0.01), 
                  chromosome = "X", 
  absolute_sex_bias_in_expression = 0,
                  mean_expression_level = mean(twas_model_data$absolute_sex_bias_in_expression))
inv_logit <- function(x) exp(x)/(1+exp(x))
preds <- data.frame(new, prop = predict(intersex_early_twas, newdata = new, se.fit=T, type = "link")) %>% 
  mutate(y = prop.fit, lwr = prop.fit - 1.96*prop.se.fit, upr = prop.fit + 1.96*prop.se.fit,
         y = inv_logit(y), lwr = inv_logit(lwr), upr = inv_logit(upr))
p3 <- preds %>% 
  mutate(panel = "C.") %>%
  ggplot(aes(mean_expression_level,  y)) + geom_line() + 
  geom_hline(yintercept = .01, linetype = 3) + 
  geom_line(aes(y = upr), linetype = 2) + 
  geom_line(aes(y = lwr), linetype = 2) + 
  ylab("Probability of being in the top 1%\nsexually antagonistic transcripts (\u00B1 95% CIs)") +
  xlab("Average expression level") + 
  facet_wrap( ~ panel) + 
  theme_bw() + 
  theme(text = element_text(family = "Raleway", size = 12)) +
  theme(strip.background = element_blank(), strip.text = element_text(hjust = 0)) 


blank <- ggplot() + theme_void()

grid.arrange(arrangeGrob(p1, p2, ncol = 2), 
             arrangeGrob(blank, p3, blank, nrow = 1, widths = c(0.2,0.6,0.2)))

pdf("figures/fig6_models.pdf", height = 7, width = 7)
grid.arrange(arrangeGrob(p1, p2, ncol = 2), 
             arrangeGrob(blank, p3, blank, nrow = 1, widths = c(0.2,0.6,0.2)))
dev.off()