Plots and models of variant effect sizes

do_ppcheck <- function(rr, title){
    pp_check(mut_model, resp = rr, ndraws = 50) + 
    scale_y_continuous(limits = c(0, 0.55)) +
    labs(subtitle=title) + 
    theme_bw() + 
    theme(legend.position = "none") 
}

grid.arrange(
  do_ppcheck("Femalefitnessearly", "Female fitness early"),
  do_ppcheck("Femalefitnesslate", "Female fitness late"),
  do_ppcheck("Malefitnessearly", "Male fitness early"),
  do_ppcheck("Malefitnesslate", "Male fitness late"))

chisq_table <- as.table(rbind(c(129+103, 208987 - (129+103)), c(93+48, 208987 - (93+48))))
dimnames(chisq_table) <- list(life_stage = c("Early_life", "Late_life"),
                              antagonism_status = c("Sex_antag","Non_sex_antag"))
chisq_table

chisq.test(chisq_table)

chisq_table <- as.table(rbind(c(691+689, 14286 - (691+689)), c(677 + 664, 14286 - (677 + 664))))
dimnames(chisq_table) <- list(life_stage = c("Early_life", "Late_life"),
                              antagonism_status = c("Sex_antag","Non_sex_antag"))
chisq_table

chisq.test(chisq_table)

get_antagonism_ratios <- function(dat){
  dat %>%
    
    # Convert the LFSR to the probability that the effect size is positive
    mutate(pp_female_early = ifelse(pheno1 > 0, LFSR1, 1 - LFSR1),
           pp_female_late  = ifelse(pheno2 > 0, LFSR2, 1 - LFSR2),
           pp_male_early   = ifelse(pheno3 > 0, LFSR3, 1 - LFSR3),
           pp_male_late    = ifelse(pheno4 > 0, LFSR4, 1 - LFSR4)) %>%
    
    # Calculate the probabilities that beta_i and beta_j have the same/opposite signs
    mutate(p_sex_concord_early = pp_female_early * pp_male_early + 
             (1 - pp_female_early) * (1 - pp_male_early),
           p_sex_antag_early = pp_female_early * (1 - pp_male_early) + 
             (1 - pp_female_early) * pp_male_early,
           p_sex_concord_late  = pp_female_late * pp_male_late + 
             (1 - pp_female_late) * (1 - pp_male_late),
           p_sex_antag_late = pp_female_late * (1 - pp_male_late) + 
             (1 - pp_female_late) * pp_male_late,
           p_age_concord_females = pp_female_early * pp_female_late + 
             (1 - pp_female_early) * (1 - pp_female_late),
           p_age_antag_females = pp_female_early * (1 - pp_female_late) + 
             (1 - pp_female_early) * pp_female_late,
           p_age_concord_males = pp_male_early * pp_male_late + (1 - pp_male_early) * (1 - pp_male_late),
           p_age_antag_males = pp_male_early * (1 - pp_male_late) + (1 - pp_male_early) * pp_male_late) %>%
    
    # Find the ratios of some of these two probabilities (i.e. the "evidence ratios")
    mutate(inter_sex_early = p_sex_concord_early / p_sex_antag_early,
           inter_sex_late = p_sex_concord_late / p_sex_antag_late,
           inter_age_females = p_age_concord_females / p_age_antag_females,
           inter_age_males = p_age_concord_males / p_age_antag_males) 
}

antagonism_evidence_ratios_gwas <- all_snps %>%
  rename(pheno1 = beta_female_early_mashr_ED, 
         pheno2 = beta_female_late_mashr_ED, 
         pheno3 = beta_male_early_mashr_ED, 
         pheno4 = beta_male_late_mashr_ED,
         LFSR1 = LFSR_female_early_mashr_ED, 
         LFSR2 = LFSR_female_late_mashr_ED, 
         LFSR3 = LFSR_male_early_mashr_ED, 
         LFSR4 = LFSR_male_late_mashr_ED) %>% 
  select(SNP, SNP_clump, starts_with("pheno"), LFSR1, LFSR2, LFSR3, LFSR4) %>% 
  filter(!is.na(pheno1)) %>% 
  collect(n=Inf) %>% 
  distinct() %>% 
  get_antagonism_ratios() %>%
  select(SNP, SNP_clump, starts_with("inter")) 

TWAS_ED <- readRDS("data/derived/TWAS/TWAS_ED.rds")

twas <- data.frame(
  FBID = read_csv("data/derived/TWAS/TWAS_results.csv")$FBID,
      as.data.frame(get_pm(TWAS_ED)),
      as.data.frame(get_lfsr(TWAS_ED)) %>% 
        rename_all(~str_replace_all(., "beta", "LFSR"))) %>% 
      as_tibble() %>%
  left_join(tbl(db, "genes") %>% select(FBID, chromosome) %>% collect(), by = "FBID") %>%
  left_join(read_csv("data/derived/gene_expression_by_sex.csv"), by = "FBID") %>%
  filter(chromosome %in% c("2L", "2R", "3L", "3R", "X")) %>%
  rename(pheno1 = beta_FE, 
         pheno2 = beta_FL, 
         pheno3 = beta_ME, 
         pheno4 = beta_ML,
         LFSR1 = LFSR_FE, 
         LFSR2 = LFSR_FL, 
         LFSR3 = LFSR_ME, 
         LFSR4 = LFSR_ML)

antagonism_evidence_ratios_twas <- twas %>% 
  get_antagonism_ratios() %>%
  select(FBID, chromosome, male_bias_in_expression, AveExpr, starts_with("inter")) %>% 
  mutate(chromosome = relevel(factor(chromosome), ref = "X"))

# Annotate data with MAF, site, class, etc. 
dat <- left_join(antagonism_evidence_ratios_gwas,
          collect(tbl(db, "variants")), by = "SNP", multiple = "all") %>% 
  distinct() %>% rename(chromosome = chr)

# Remove chromosome 4 (too few sites)
dat <- dat %>% filter(chromosome != "4") 
dat$chromosome <- relevel(factor(dat$chromosome), ref = "X")

# Classify mutations as being major effect or not
dat <- dat %>% 
  mutate(major_effect = "No",
         major_effect = replace(major_effect, site.class %in% major_effect_types, "Yes"))

# For mutations with multiple site classes, record whether any site class is major or not,
# ensuring that there is one row in 'dat' for each SNP
dat <- dat %>% 
  group_by(SNP_clump) %>% 
  summarise(inter_sex_early = inter_sex_early[1],
            inter_sex_late = inter_sex_late[1],
            MAF = MAF[1],
            chromosome = chromosome[1],
            major_effect = ifelse(any(major_effect == "Yes"), "Yes", "No"))

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

intersex_early_model_gwas <- glm(antagonistic ~ chromosome + MAF + major_effect, 
             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,  ~ . -major_effect)
car::Anova(intersex_early_model_gwas, type = 3)

summary(intersex_early_model_gwas)

intersex_late_model_gwas <- glm(antagonistic ~ chromosome + MAF + major_effect, 
             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,  ~ . -major_effect)
car::Anova(intersex_late_model_gwas, type = 3)

summary(intersex_late_model_gwas)

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)

intersex_early_twas <- update(
  intersex_early_twas,  
  ~ . -absolute_sex_bias_in_expression -chromosome)
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), chromosome = "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(chromosome = levels(dat$chromosome), 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(chromosome,  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))

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/fig7_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)))
invisible(dev.off())