GWAS for fitness in Drosophila

Define functions for prediction

load_and_tidy_fitness_data <- function(){
  female_fitness <- read.csv("data/input/female_fitness.csv", stringsAsFactors = FALSE) %>% 
    arrange(line) %>%
    mutate(vial.ID = paste("F", block, replicate, line, sep = "_"),
           line = paste("line_", line, sep = "")) %>% 
    select(-num.F.late1, -num.F.late2, -num.laying.F.early, -replicate) 
  female_fitness$female.fitness.early[is.na(female_fitness$female.fitness.early)] <- 0
  female_fitness$female.fitness.late[is.na(female_fitness$female.fitness.late)] <- 0
  
  male_fitness <- read.csv("data/input/male_fitness.csv", stringsAsFactors = FALSE ) %>% 
    arrange(line) %>%
    mutate(vial.ID = paste("M", block, replicate, line, sep = "_"),
           line = paste("line_", line, sep = "")) %>%
    select(-num.DGRP.males, -replicate)
  
  # Save these tidy files, for data archiving
  write.csv(female_fitness, file = "data/input/female_fitness_CLEANED.csv", row.names = FALSE)
  write.csv(male_fitness, file = "data/input/male_fitness_CLEANED.csv", row.names = FALSE)
  
  list(female_fitness = female_fitness, male_fitness = male_fitness)
}


# This function calculates the means for each line
get.line.means <- function(){
  left_join(
    load_and_tidy_fitness_data()[[1]] %>% select(line, female.fitness.early, female.fitness.late) %>% 
      group_by(line) %>%
      summarise(female.fitness.early = mean(female.fitness.early), 
                female.fitness.late = mean(female.fitness.late)),
    load_and_tidy_fitness_data()[[2]] %>% 
      group_by(line) %>%
      summarise(early.male.rival = sum(early.male.rival),
                early.male.focal = sum(early.male.focal),
                late.male.rival = sum(late.male.rival),
                late.male.focal = sum(late.male.focal),
                male.fitness.early = early.male.focal / (early.male.focal + early.male.rival),
                male.fitness.late = late.male.focal / (late.male.focal + late.male.rival)) %>%
      select(line, male.fitness.early, male.fitness.late),
    by = "line") %>% as.data.frame()
}

# A function to calculate each line's expected fitness, using brms models to adjust for block + vial random effects
# For the male late life assay, we also adjust for the number of competitor males that died 
# (this assumes that their deaths happened randomly with respect to which line we were measuring)
get_predicted_line_means <- function(){
  
  # If it's already calculated, just load up the file. Otherwise, make the file
  out.file <- "data/derived/predicted_line_means.csv"
  if(file.exists(out.file)) return(read.csv(out.file, stringsAsFactors = FALSE))
  
  ####### Load the individual-level data
  # Load the female data and rearrange to long format
  female_fitness <- load_and_tidy_fitness_data()[[1]] %>% 
    gather(age, fitness, female.fitness.early, female.fitness.late)
  
  # Load the male data and make 5 imputed datasets using the package 'mice'
  # (there are 2 vials with no early-male fitness, and 1 vial with no count of the number of GFP males)
  # See this link for explanation: https://cran.r-project.org/web/packages/brms/vignettes/brms_missings.html
  male_fitness <- load_and_tidy_fitness_data()[[2]]
  male_fitness_imputed <- mice(male_fitness, m = 5, print = FALSE, seed = 1)
  male_fitness_imputed <- lapply(1:5, function(i) complete(male_fitness_imputed, i))
  male_fitness_imputed <- lapply(male_fitness_imputed, function(focal) {
    focal %>% mutate(total1 = early.male.rival + early.male.focal,
                     total2 = late.male.rival + late.male.focal)})
  
  ####### Run model on female data 
  female_model <- brm(fitness ~ age + (age|line) + (age|block) + (1|vial.ID), 
                      family = poisson, data = female_fitness,
                      iter = 4000, chains = 4, cores = 4,
                      seed = 1, control = list(adapt_delta = 0.999, max_treedepth = 15)) %>% 
    add_ic(ic = "R2")
  writeLines(capture.output(summary(female_model)), con = "data/derived/model_summary_females.txt")
  saveRDS(female_model, "data/derived/female_fitness_prediction_model.rds") # maybe delete later - probably not needed?
  
  
  ### Run model on male data 
  # Fit a non-linear smooth for num.GFP.males, but only for the late life data
  # Inspection of the raw data and model results illustrates the non-linear smooth improves model fit appreciably
  # For example, the 'wiggliness' of the smooth term is significantly non-zero
  # See: https://www.fromthebottomoftheheap.net/2018/04/21/fitting-gams-with-brms/
  male1 <- bf(early.male.focal | trials(total1) ~ (1|line) + (1|block) + (1|vial.ID))
  male2 <- bf(late.male.focal  | trials(total2) ~ s(num.GFP.males) + (1|line) + (1|block) + (1|vial.ID))
  
  # Run 5 models, on each of the 5 imputed datasets, and combine the results
  male_model <- brm_multiple(male1 + male2,
                                 family = binomial, data = male_fitness_imputed,
                                 iter = 4000, chains = 4, cores = 4,
                                 seed = 1, control = list(adapt_delta = 0.999, max_treedepth = 15)) %>%
    add_ic(ic = "R2")
  writeLines(capture.output(summary(male_model)), con = "data/derived/model_summary_males.txt")
  saveRDS(male_model, "data/derived/male_fitness_prediction_model.rds") # maybe delete later - probably not needed?
  male_model <- readRDS("data/derived/male_fitness_prediction_model.rds") 
  
  ####### Define new data for prediction
  new.data.female <- female_fitness %>% 
    select(age, line) %>% distinct() %>% arrange(line, age)
  
  new.data.male <- male_fitness %>% 
    select(line) %>% distinct() %>% arrange(line) %>% 
    mutate(total1 = 100, # I checked, and changing this number does not affect the line means, as I expected
           total2 = 100, 
           num.GFP.males = male_fitness %>% .$num.GFP.males %>% mean(na.rm=TRUE))
  
  # Get predicted lines means for early and late life female fitness on linear predictor scale (i.e. normal)
  female_line_preds <- data.frame(new.data.female, 
                                  fitted(female_model, 
                                         newdata = new.data.female, 
                                         re_formula = "~ (age|line)", 
                                         scale = "linear")) %>%
    select(age, line, Estimate) %>% 
    spread(age, Estimate) %>%
    mutate(female.fitness.early = as.numeric(scale(female.fitness.early)),
           female.fitness.late = as.numeric(scale(female.fitness.late)))
  
  # Get predicted lines means for early and late life male fitness on linear predictor scale (i.e. normal)
  male_line_preds <- data.frame(new.data.male, 
                              fitted(male_model_NEW, 
                                     newdata = new.data.male, 
                                     re_formula = "~ (1|line)", 
                                     scale = "linear")) %>%
  select(line, Estimate.earlymalefocal, Estimate.latemalefocal) %>% 
  mutate(male.fitness.early = as.numeric(scale(Estimate.earlymalefocal)),
         male.fitness.late = as.numeric(scale(Estimate.latemalefocal))) %>%
  select(line, male.fitness.early, male.fitness.late)

  # Calculate the proportion of variance explained by line for females
  
  ## Get the posterior predicted values for each data point
  predicted <- data.frame(female_fitness[complete.cases(female_fitness), ],
                          t(fitted(female_model, summary = FALSE)))
  
  ## Define a function to get the ICC for line, for each posterior sample of predictions
  get_ICC_line <- function(predictions, age_class = NULL){
    if(!is.null(age_class)) predictions <- predictions %>% filter(age == age_class) 
    first.col <- which(str_detect(names(predictions), "X"))[1] 
    sapply(0:7999, function(i){
       ICCbare(predictions$line, predictions[, i+first.col])
    })
  }
  
  # Calculate ICC for each sample, separately for old and young females
  Line_ICC_female_early <- get_ICC_line(predicted, "female.fitness.early")
  Line_ICC_female_late <- get_ICC_line(predicted, "female.fitness.late")
  
  # Do the same for males
  male_preds_early <-  data.frame(line = male_fitness$line,
                                  t(fitted(male_model, summary = FALSE)[,,1]))
  male_preds_late <-  data.frame(line = male_fitness$line,
                                 t(fitted(male_model, summary = FALSE)[,,1]))

  Line_ICC_male_early <- get_ICC_line(male_preds_early)
  Line_ICC_male_late <- get_ICC_line(male_preds_late)
  
  
  ICC <- posterior_summary(as.mcmc(data.frame(Line_ICC_female_early, 
                                              Line_ICC_female_late,
                                              Line_ICC_male_early,
                                              Line_ICC_male_late)))
  write.csv(ICC, file = "data/derived/ICC.csv", row.names = TRUE)
  

  # Record which block(s) each line was measured in. Generally one block, except for ref line (all blocks)
  line_blocks <- rbind(female_fitness %>% select(block, line), 
                       male_fitness   %>% select(block, line)) %>% 
    mutate(block = replace(block, line == "352", "reference_line")) %>% 
    distinct() %>% 
    group_by(line) %>% summarise(block = paste0(block, collapse = ", ")) 
  
  predicted_line_means <- female_line_preds %>%
    left_join(male_line_preds, by = "line") %>% 
    left_join(line_blocks, by = "line") %>%
    arrange(line) %>%
    mutate(line = paste("line_", line, sep = "")) 
  
  write.csv(predicted_line_means, out.file, row.names = FALSE)
  predicted_line_means %>% mutate(line = gsub("line_line_", "line_", line))
}

Predict line means and save the results

The object line_means is the raw line means: unscaled, arithmetic means and in the original units.

The object predicted_line_means is the estimated line means from the Bayesian models. These means are on the scale of the linear predictor (i.e. Poisson for females, binomial for males) rather than on the scale of the original response. They have been scaled to have a mean of 0 and a variance of 1.

line_means <- get.line.means()
predicted_line_means <- get_predicted_line_means() 

Summary of model used to predict FEMALE line means:

cat(readLines('data/derived/model_summary_females.txt'), sep = '\n')

 Family: poisson 
  Links: mu = log 
Formula: fitness ~ age + (age | line) + (age | block) + (1 | vial.ID) 
   Data: female_fitness (Number of observations: 1500) 
Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1; 
         total post-warmup samples = 8000
    ICs: LOO = NA; WAIC = NA; R2 = 0.94
 
Group-Level Effects: 
~block (Number of levels: 9) 
                                      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)                             0.32      0.12     0.16     0.62       5698 1.00
sd(agefemale.fitness.late)                1.29      0.43     0.74     2.35       4667 1.00
cor(Intercept,agefemale.fitness.late)    -0.41      0.32    -0.88     0.33       8000 1.00

~line (Number of levels: 115) 
                                      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)                             0.54      0.04     0.46     0.63       4294 1.00
sd(agefemale.fitness.late)                1.58      0.13     1.34     1.87       4293 1.00
cor(Intercept,agefemale.fitness.late)     0.41      0.10     0.21     0.58       3559 1.00

~vial.ID (Number of levels: 750) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     0.48      0.02     0.45     0.51       3005 1.00

Population-Level Effects: 
                       Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept                  4.40      0.13     4.15     4.65       8000 1.00
agefemale.fitness.late    -2.96      0.49    -3.96    -2.00       8000 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

Summary of model used to predict MALE line means:

cat(readLines('data/derived/model_summary_males.txt'), sep = '\n')

 Family: MV(binomial, binomial) 
  Links: mu = logit
         mu = logit 
Formula: early.male.focal | trials(total1) ~ (1 | line) + (1 | block) + (1 | vial.ID) 
         late.male.focal | trials(total2) ~ s(num.GFP.males) + (1 | line) + (1 | block) + (1 | vial.ID) 
   Data: male_fitness_imputed (Number of observations: 750) 
Samples: 20 chains, each with iter = 4000; warmup = 2000; thin = 5; 
         total post-warmup samples = 8000
    ICs: LOO = NA; WAIC = NA; R2 = 0.98
 
Smooth Terms: 
                                    Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sds(latemalefocal_snum.GFP.males_1)     2.80      2.11     0.27     8.30       7683 1.00

Group-Level Effects: 
~block (Number of levels: 9) 
                             Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(earlymalefocal_Intercept)     0.73      0.24     0.41     1.34       7026 1.00
sd(latemalefocal_Intercept)      1.00      0.38     0.48     1.94       7293 1.00

~line (Number of levels: 115) 
                             Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(earlymalefocal_Intercept)     0.57      0.06     0.47     0.68       7060 1.00
sd(latemalefocal_Intercept)      1.37      0.15     1.09     1.67       7202 1.00

~vial.ID (Number of levels: 750) 
                             Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(earlymalefocal_Intercept)     0.81      0.03     0.76     0.87       5630 1.01
sd(latemalefocal_Intercept)      1.91      0.08     1.76     2.08       6832 1.00

Population-Level Effects: 
                               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
earlymalefocal_Intercept           0.30      0.27    -0.25     0.81       7744 1.00
latemalefocal_Intercept           -1.17      0.39    -1.94    -0.37       7624 1.00
latemalefocal_snum.GFP.males_1     0.08      0.60    -1.08     1.43       7941 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

Proportion of variance explained by ‘DGRP line’ in these two models:

read.csv('data/derived/ICC.csv') %>% 
  rename(Parameter = X, Lower.95.CI = X2.5.ile, Upper.95.CI = X97.5.ile) %>% pander()