Estimating genetic (co)variance and heritability

line_means %>% select(contains(" ")) %>% cor() %>% kable(digits =3) %>% kable_styling()

run_command(glue("{gcta} --bfile dgrp2_QC_all_lines",
                 " --make-grm",
                 " --out gcta_GRM"), path = "/data/derived/")

do_univariate_greml <- function(trait1){

  # First make 'focal_data', a 2-column data frame with the line and the focal phenotype value
  focal_data <- line_means %>%
    select(line, !! trait1)
  names(focal_data)[2] <- c("focal_pheno")
  
  # Write a file with the line and phenotype data called phenotype.txt, for gcta64
  pheno_data <- focal_data %>%
    # filter(!(line %in% lines_to_prune)) %>% 
    mutate(line_copy = line) %>%
    select(line, line_copy, focal_pheno) %>%
    as.matrix()
  
  pheno_data %>%
    write.table(row.names = FALSE, col.names = FALSE,
                file = "data/derived/phenotype.txt", quote = FALSE)
  
  console_output <- suppressWarnings(capture.output(
    run_command(glue("{gcta}  --reml --reml-maxit 10000 --grm gcta_GRM  --pheno phenotype.txt  --out delete_me"), 
                path = "/data/derived/")))
  
  errors <- console_output[str_detect(console_output, "[Ee]rror")]
  if(length(errors) == 0) errors <- NA
  if(length(errors) > 1) errors <- paste(errors, sep = "; ")

  rbind(read.delim("data/derived/delete_me.hsq", sep = "\t") , 
        data.frame(Source = "Errors", Variance = errors, SE = NA)) %>% 
    mutate(Trait = trait1) %>% select(Trait, everything())
}

bind_rows(do_univariate_greml("Female fitness early"),
          do_univariate_greml("Female fitness late"),
          do_univariate_greml("Male fitness early"),
          do_univariate_greml("Male fitness late")) %>% 
  filter(Source %in% c("V(G)", "V(e)", "Vp", "V(G)/Vp")) %>% 
  mutate(Variance = as.numeric(Variance)) %>% 
  rename(Parameter = Source) %>% 
  kable(digits = 3) %>% 
  kable_styling(full_width = FALSE)

do_bivar_greml <- function(trait1, trait2){

  # First make 'focal_data', a 2-column data frame with the line and the focal phenotype value
  focal_data <- line_means %>%
    select(line, !! trait1, !! trait2)
  names(focal_data)[2:3] <- c("focal_pheno1", "focal_pheno2")
  
  # Write a file with the line and phenotype data called phenotype.txt, for gcta64
  pheno_data <- focal_data %>%
    # filter(!(line %in% lines_to_prune)) %>% 
    mutate(line_copy = line) %>%
    select(line, line_copy, focal_pheno1, focal_pheno2) %>%
    as.matrix()
  
  pheno_data %>%
    write.table(row.names = FALSE, col.names = FALSE,
                file = "data/derived/phenotype_pair.txt", quote = FALSE)
  
  # Note that the warnings come when gcta could not fit the model properly. I save these errors especially, so no need to print them.
  console_output <- suppressWarnings(capture.output(
    run_command(glue("{gcta}  --reml-bivar --reml-bivar-lrt-rg 0 --reml-maxit 10000 --grm gcta_GRM  --pheno phenotype_pair.txt  --out delete_me"), 
                path = "/data/derived/")))
  
  errors <- console_output[str_detect(console_output, "[Ee]rror")]
  if(length(errors) == 0) errors <- NA
  if(length(errors) > 1) errors <- paste(errors, sep = "; ")
  traits <- sort(c(trait1, trait2))
  
  rbind(read.delim("data/derived/delete_me.hsq", sep = "\t") , 
        data.frame(Source = "Errors", Variance = errors, SE = NA)) %>% 
    mutate(`Trait 1` = traits[1], `Trait 2` = traits[2]) %>% select(`Trait 1`, `Trait 2`, everything())
}

results <- bind_rows(do_bivar_greml("Female fitness early", "Female fitness late"),
          do_bivar_greml("Male fitness early", "Male fitness late"),
          do_bivar_greml("Female fitness early", "Male fitness early"),
          do_bivar_greml("Female fitness late", "Male fitness late"),
          do_bivar_greml("Female fitness early", "Male fitness late"),
          do_bivar_greml("Female fitness late", "Male fitness early"))

unlink(c("data/derived/delete_me.log", "data/derived/phenotype_pair.txt", "data/derived/delete_me.hsq"))
unlink(list.files("data/derived/", pattern = "gcta_GRM", full.names = TRUE))

results %>% 
  filter(Source %in% c("V(G)_tr1", "V(e)_tr1", "Vp_tr1", "V(G)/Vp_tr1",
                       "V(G)_tr2", "V(e)_tr2", "Vp_tr2", "V(G)/Vp_tr2",
                       "C(G)_tr12", "rG")) %>% 
  mutate(Variance = as.numeric(Variance)) %>% 
  rename(Parameter = Source) %>% 
  kable(digit = 3) %>% 
  kable_styling(full_width = FALSE)