R/run_exwas_full.R

In groovearmada/exwas:

run_exwas_full <- function(input_data, variable_list, sdq, subcohort = NULL, file_name) {

  data_sel <- input_data %>%
    tidyr::gather(key = "Term", value = "concentration", -c(HelixID, h_sdq_external, h_sdq_internal, cohort, smoking:sex)) %>%
    tidyr::gather(key = "SDQ", value = "score", h_sdq_external:h_sdq_internal) %>%
    group_by(SDQ, Term) %>%
    group_modify(~ broom::tidy(MASS::glm.nb(score ~ concentration + cohort + smoking + mother_bmi + mother_age + mother_edu + parity + breastf_cat + mother_work + child_age + sex, data = .), conf.int = TRUE, exponentiate = TRUE)) %>% # model outputs as data frame
    dplyr::filter(term == "concentration")

  # transform rownames into 1st row
  tidy_fit <- data_sel %>%
    dplyr::mutate(CI = stringr::str_c(round(conf.low, 2), round(conf.high, 2), sep = "; "),
                  estCI = stringr::str_c(round(estimate, 2), " (", CI, ")", sep = ""),
                  p_value = round(p.value, 3)) %>%
    dplyr::select(-c(term, std.error, statistic, p.value)) %>%
    dplyr::rename(Estimate = estimate, conf_low = conf.low, conf_high = conf.high)

  # Add exposure names
  tidy_exwas_result <- tidy_fit %>%

    # add variable descriptions
    dplyr::left_join(variable_list[, c(1:3)], ., by = "Term") %>%

    # Order the table by Exposure
    dplyr::arrange(SDQ, Family, Exposure)

  # remove unnecessary variables
  pooled_fits_plot <- dplyr::select(tidy_exwas_result, SDQ, Family, Exposure, Estimate, conf_low, conf_high)

  exwas_sel <- dplyr::select(tidy_exwas_result, SDQ, Family, Exposure, estCI, p_value)

  # export table result to a .csv file
  write.csv(pooled_fits_plot, file = paste0(file_name, "_plot.csv"))

  # export table result to a .csv file
  write.csv(exwas_sel, file = paste0(file_name, ".csv"))

  return(data_sel)
}