Tables summarising the TWAS results

workflowr

• Summary
• Checks
• Past versions

Last updated: 2023-05-17

Checks: 7 0

Knit directory: fitnessGWAS/

This reproducible R Markdown analysis was created with workflowr (version 1.7.0.4). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

---

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20180914)

The command set.seed(20180914) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: b407d11

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version b407d11. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:

Ignored files:
    Ignored:    .DS_Store
    Ignored:    .Rapp.history
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    .httr-oauth
    Ignored:    .pversion
    Ignored:    analysis/.DS_Store
    Ignored:    code/.DS_Store
    Ignored:    code/Drosophila_GWAS.Rmd
    Ignored:    data/.DS_Store
    Ignored:    data/derived/
    Ignored:    data/input/.DS_Store
    Ignored:    data/input/.pversion
    Ignored:    data/input/dgrp.fb557.annot.txt
    Ignored:    data/input/dgrp2.bed
    Ignored:    data/input/dgrp2.bim
    Ignored:    data/input/dgrp2.fam
    Ignored:    data/input/huang_transcriptome/
    Ignored:    figures/.DS_Store

Untracked files:
    Untracked:  old_analyses/

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

---

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/TWAS_tables.Rmd) and HTML (docs/TWAS_tables.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	b407d11	lukeholman	2023-05-17	wflow_publish("analysis/TWAS_tables.Rmd")
html	d14fe65	lukeholman	2023-04-17	Build site.
html	7b73b2b	lukeholman	2022-07-29	Build site.
Rmd	98d6d09	lukeholman	2022-07-29	wflow_publish("analysis/TWAS_tables.Rmd")
Rmd	68414e9	lukeholman	2022-02-22	Commit in Feb 2022
html	68414e9	lukeholman	2022-02-22	Commit in Feb 2022
html	ca70247	lukeholman	2021-11-10	Build site.
html	7449a90	lukeholman	2021-10-01	Build site.
Rmd	01226ab	lukeholman	2021-10-01	wflow_publish("analysis/*")
html	8d14298	lukeholman	2021-09-26	Build site.
Rmd	af15dd6	lukeholman	2021-09-26	Commit Sept 2021
html	372da3e	lukeholman	2021-03-05	Build site.
Rmd	ac8ee22	lukeholman	2021-03-05	wflow_publish("analysis/TWAS_tables.Rmd")
Rmd	3855d33	lukeholman	2021-03-04	big fist commit 2021
html	3855d33	lukeholman	2021-03-04	big fist commit 2021

---

library(tidyverse)
library(DT)
library(kableExtra)

kable_table <- function(df) { # cool tables
  kable(df, "html") %>%
  kable_styling(full_width = FALSE) %>%
  scroll_box(height = "300px")
}

my_data_table <- function(df){ # Make html tables:
  datatable(
    df, rownames=FALSE,
    autoHideNavigation = TRUE,
    extensions = c("Scroller",  "Buttons"),
    options = list(
      dom = 'Bfrtip',
      deferRender=TRUE,
      scrollX=TRUE, scrollY=400,
      scrollCollapse=TRUE,
      buttons = 
        list('pageLength', 'colvis', 'csv', list(
          extend = 'pdf',
          pageSize = 'A4',
          orientation = 'landscape',
          filename = 'TWAS_sig_genes')),
      pageLength = 50
    )
  ) %>%
    formatStyle(
        columns = c("Female early effect", "Male early effect", "Female late effect", "Male late effect"), 
        color = styleInterval(cuts = 0, values = c("tomato", "steelblue")),
        fontWeight = "bold")
}

sci_notation <- function(x) formatC(x, format = "e", digits = 1)
rnd <- function(x) format(round(x, 2), nsmall = 2)

make_tally <- function(p_cutoff){
  TWAS_results <- left_join(read.csv("data/derived/TWAS/TWAS_result_females.csv") %>% 
              select(gene, beta_FE, beta_FL, pval_FE, pval_FL), 
            read.csv("data/derived/TWAS/TWAS_result_males.csv") %>% 
              select(gene, beta_ME, beta_ML, pval_ME, pval_ML), by = "gene") %>%
    mutate(sum_p = pval_FE + pval_FL + pval_ME + pval_ML) %>% 
    arrange(sum_p) %>% select(-sum_p) %>% 
    mutate(
      type = 
        factor(case_when(
          
          (sign(beta_FE) != sign(beta_ME) & pval_FE < p_cutoff & pval_ME < p_cutoff) & 
            (sign(beta_FL) != sign(beta_ML) & pval_FL < p_cutoff & pval_ML < p_cutoff) ~ "Sex antagonistic, both ages",
          
          (sign(beta_FE) == sign(beta_ME) & pval_FE < p_cutoff & pval_ME < p_cutoff) & 
            (sign(beta_FL) == sign(beta_ML) & pval_FL < p_cutoff & pval_ML < p_cutoff) ~ "Sex concordant, both ages",
          
          (sign(beta_FE) != sign(beta_FL) & pval_FE < p_cutoff & pval_FL < p_cutoff) & 
            (sign(beta_ME) != sign(beta_ML) & pval_ME < p_cutoff & pval_ML < p_cutoff) ~ "Age antagonistic, both sexes",
          
          (sign(beta_FE) == sign(beta_FL) & pval_FE < p_cutoff & pval_FL < p_cutoff) &
            (sign(beta_ME) == sign(beta_ML) & pval_ME < p_cutoff & pval_ML < p_cutoff) ~ "Age concordant, both sexes",
          
          sign(beta_FE) != sign(beta_ME) & pval_FE < p_cutoff & pval_ME < p_cutoff ~ "Sex antagonistic, early",
          sign(beta_FL) != sign(beta_ML) & pval_FL < p_cutoff & pval_ML < p_cutoff ~ "Sex antagonistic, late",
          
          sign(beta_FE) == sign(beta_ME) & pval_FE < p_cutoff & pval_ME < p_cutoff ~ "Sex concordant, early",
          sign(beta_FL) == sign(beta_ML) & pval_FL < p_cutoff & pval_ML < p_cutoff ~ "Sex concordant, late",
          
          sign(beta_FE) != sign(beta_FL) & pval_FE < p_cutoff & pval_FL < p_cutoff ~ "Age antagonistic, females",
          sign(beta_ME) != sign(beta_ML) & pval_ME < p_cutoff & pval_ML < p_cutoff ~ "Age antagonistic, males",
          
          sign(beta_FE) == sign(beta_FL) & pval_FE < p_cutoff & pval_FL < p_cutoff ~ "Age concordant, females",
          sign(beta_ME) == sign(beta_ML) & pval_ME < p_cutoff & pval_ML < p_cutoff ~ "Age concordant, males",
          
          pval_FE < p_cutoff ~ "Female early only",
          pval_FL < p_cutoff ~ "Female late only",
          pval_ME < p_cutoff ~ "Male early only",
          pval_ML < p_cutoff ~ "Male late only",
          TRUE ~ "Uncorrelated with fitness"
        ), levels = c("Sex antagonistic, both ages", "Sex concordant, both ages", 
                      "Age antagonistic, both sexes", "Age concordant, both sexes",
                      "Sex antagonistic, early", "Sex concordant, early", "Sex antagonistic, late", "Sex concordant, late",
                      "Age antagonistic, females", "Age antagonistic, males", "Age concordant, females", "Age concordant, males",
                      "Female early only", "Male early only", "Female late only", "Male late only", "Uncorrelated with fitness")))
  
  tally <- table(TWAS_results$type) %>% sort(decreasing = T) %>% enframe()
  names(tally)[2] <- paste("Using p <", p_cutoff)
  tally
}

TWAS_tally_table <- make_tally(0.05) %>% 
  left_join(make_tally(0.01), by = "name") %>% 
  left_join(make_tally(0.001), by = "name") %>% 
  left_join(make_tally(0.0001), by = "name") %>% 
  left_join(make_tally(0.00001), by = "name") %>% 
  rename(`Relationship to fitness` = name)

saveRDS(TWAS_tally_table, "data/derived/TWAS_tally_table.rds")

p_cutoff <- 0.01
sig_TWAS_results <- 
  left_join(read_csv("data/derived/TWAS/TWAS_result_females.csv") %>% 
              select(gene, beta_FE, beta_FL, pval_FE, pval_FL), 
            read_csv("data/derived/TWAS/TWAS_result_males.csv") %>% 
              select(gene, beta_ME, beta_ML, pval_ME, pval_ML), by = "gene") %>%
  mutate(sum_p = pval_FE + pval_FL + pval_ME + pval_ML) %>% arrange(sum_p) %>% 
  select(-sum_p) %>% 
  filter(pval_FE < p_cutoff | pval_FL < p_cutoff | pval_ME < p_cutoff | pval_ML < p_cutoff) 


big_table <- sig_TWAS_results %>% 
  mutate(pval_FE = sci_notation(pval_FE), 
         pval_FL = sci_notation(pval_FL),
         pval_ME = sci_notation(pval_ME),
         pval_ML = sci_notation(pval_ML)) %>%
  left_join(read_csv("data/derived/TWAS/TWAS_results.csv") %>% 
              select(FBID, gene_name, chromosome, male_bias_in_expression, AveExpr), by = c("gene" = "FBID")) %>% 
  select(gene, gene_name, chromosome, male_bias_in_expression, AveExpr,  starts_with("beta"), starts_with("pval")) %>%
  # mutate(across(where(is.numeric), rnd)) %>%
  rename(`Gene name` = gene_name,
         Chromosome = chromosome,
         `Female early effect` = beta_FE,
         `Female late effect` = beta_FL,
         `Male early effect` = beta_ME,
         `Male late effect` = beta_ML,
         `Female early pval` = pval_FE,
         `Female late pval` = pval_FL,
         `Male early pval` = pval_ME,
         `Male late pval` = pval_ML,
         `Male bias in expression (logFC)` = male_bias_in_expression,
         `Average expression level` = AveExpr) 

write_csv(big_table, "data/derived/Supplementary Dataset S1.csv")

big_table <- big_table %>% 
  mutate(across(where(is.numeric), rnd)) 

Table showing significant transcripts from the TWAS

The table shows the 517 transcripts that were associated with one or more of the four phenotypes, with a p-value less than 0.01. The sex difference in expression (column 4) and the average expression level (across both sexes; column 5) were calculated from the DGRP expression data from Huang et al. 2015 (PNAS). Columns 6-9 show the regression coefficients that relate the line mean transcript abundance to the line mean fitnesses (with positive and negative values highlighted in colour), while columns 10-13 give the associated \(p\)-values.

big_table %>% my_data_table()

Table showing the numbers of transcripts related to fitness

The table tallies up the numbers of transcripts showing a particular significant relationship with fitness (rows) for various different \(p\)-value thresholds (columns). For example, ‘Female early only’ counts the number of transcripts whose expression significantly correlated with mean female early life fitness across lines. ‘Age concordant, males’ counts transcripts whose expression correlated with early- and late-life fitness in males, in the same direction, while ‘Age antagonistic, males’ counts transcripts showing a significant, opposite relationships with male early- and male late-life fitness. Similarly, ‘Sex concordant, early’ counts transcripts showing a concordant effect on early life fitness in both males and females, and ‘Sex antagonistic, late’ counts those showing opposing effects on late-life fitness in males and females. The categories are all mutually exclusive, such that transcripts are only counted towards the most specific category that applies to them. Note that this method has low power to detect transcripts that correlate with two or more fitness metrics, because there are two or more opportunities to make a ‘false negative’ error, and the power is low for any given transcripts (this is why we use a range of \(p\)-values, including some low ones).

TWAS_tally_table %>% kable_table() 

Relationship to fitness	Using p < 0.05	Using p < 0.01	Using p < 0.001	Using p < 1e-04	Using p < 1e-05
Uncorrelated with fitness	12143	13769	14222	14280	14285
Female early only	472	123	14	1	0
Female late only	384	112	15	3	1
Age concordant, males	356	59	7	0	0
Male late only	301	76	10	1	0
Age concordant, females	292	44	3	0	0
Male early only	268	89	15	1	0
Sex concordant, early	24	5	0	0	0
Sex concordant, late	18	7	0	0	0
Sex antagonistic, early	12	2	0	0	0
Sex antagonistic, late	8	0	0	0	0
Sex concordant, both ages	6	0	0	0	0
Sex antagonistic, both ages	2	0	0	0	0
Age antagonistic, both sexes	0	0	0	0	0
Age concordant, both sexes	0	0	0	0	0
Age antagonistic, females	0	0	0	0	0
Age antagonistic, males	0	0	0	0	0

Session information

sessionInfo()

R version 4.2.2 (2022-10-31)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur ... 10.16

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] kableExtra_1.3.4  DT_0.27           lubridate_1.9.2   forcats_1.0.0    
 [5] stringr_1.5.0     dplyr_1.1.0       purrr_1.0.1       readr_2.1.4      
 [9] tidyr_1.3.0       tibble_3.1.8      ggplot2_3.4.1     tidyverse_2.0.0  
[13] workflowr_1.7.0.4

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.10       svglite_2.1.1     getPass_0.2-2     ps_1.7.2         
 [5] rprojroot_2.0.3   digest_0.6.31     utf8_1.2.2        R6_2.5.1         
 [9] evaluate_0.20     highr_0.10        httr_1.4.5        pillar_1.8.1     
[13] rlang_1.0.6       rstudioapi_0.14   whisker_0.4.1     callr_3.7.3      
[17] jquerylib_0.1.4   rmarkdown_2.20    webshot_0.5.4     htmlwidgets_1.6.1
[21] bit_4.0.5         munsell_0.5.0     compiler_4.2.2    httpuv_1.6.9     
[25] xfun_0.37         pkgconfig_2.0.3   systemfonts_1.0.4 htmltools_0.5.4  
[29] tidyselect_1.2.0  fansi_1.0.4       viridisLite_0.4.1 crayon_1.5.2     
[33] tzdb_0.3.0        withr_2.5.0       later_1.3.0       grid_4.2.2       
[37] jsonlite_1.8.4    gtable_0.3.1      lifecycle_1.0.3   git2r_0.31.0     
[41] magrittr_2.0.3    scales_1.2.1      vroom_1.6.1       cli_3.6.0        
[45] stringi_1.7.12    cachem_1.0.7      fs_1.6.1          promises_1.2.0.1 
[49] xml2_1.3.3        bslib_0.4.2       ellipsis_0.3.2    generics_0.1.3   
[53] vctrs_0.5.2       tools_4.2.2       bit64_4.0.5       glue_1.6.2       
[57] crosstalk_1.2.0   hms_1.1.2         parallel_4.2.2    processx_3.8.0   
[61] fastmap_1.1.1     yaml_2.3.7        timechange_0.2.0  colorspace_2.1-0 
[65] rvest_1.0.3       knitr_1.42        sass_0.4.5