Last updated: 2021-03-04

Checks: 2 0

Knit directory: fitnessGWAS/

This reproducible R Markdown analysis was created with workflowr (version 1.6.2). 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.

Repository version: c606d3d

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 c606d3d. 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:    analysis/correlations_SNP_effects_cache/
    Ignored:    analysis/plot_models_variant_effects_cache/
    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:  Rplots.pdf
    Untracked:  analysis/GO_KEGG_enrichment.Rmd
    Untracked:  analysis/TWAS_tables.Rmd
    Untracked:  code/GO_and_KEGG_gsea.R
    Untracked:  code/eQTL_analysis.Rmd
    Untracked:  code/make_sites_files_for_ARGweaver.R
    Untracked:  code/plotting_results.Rmd
    Untracked:  code/run_argweaver.R
    Untracked:  code/run_mashr.R
    Untracked:  data/argweaver/
    Untracked:  data/input/TF_binding_sites.csv
    Untracked:  figures/GWAS_mixture_proportions.rds
    Untracked:  figures/GWAS_stats_figure.pdf
    Untracked:  figures/TWAS_mixture_proportions.rds
    Untracked:  figures/TWAS_stats_figure.pdf
    Untracked:  figures/antagonism_ratios.pdf
    Untracked:  figures/composite_mixture_figure.pdf
    Untracked:  figures/eff_size_histos.pdf
    Untracked:  figures/stats_figure.pdf
    Untracked:  manuscript/

Unstaged changes:
    Modified:   .gitignore
    Modified:   README.md
    Modified:   analysis/GWAS_tables.Rmd
    Modified:   analysis/TWAS.Rmd
    Modified:   analysis/_site.yml
    Deleted:    analysis/about.Rmd
    Deleted:    analysis/eQTL_analysis.Rmd
    Modified:   analysis/gwas_adaptive_shrinkage.Rmd
    Modified:   analysis/license.Rmd
    Modified:   analysis/make_annotation_database.Rmd
    Modified:   analysis/plot_line_means.Rmd
    Modified:   analysis/plotting_results.Rmd
    Deleted:    code/gwas_adaptive_shrinkage.R
    Modified:   data/input/female_fitness.csv
    Modified:   data/input/female_fitness_CLEANED.csv
    Modified:   data/input/male_fitness.csv
    Modified:   data/input/male_fitness_CLEANED.csv
    Modified:   figures/figure1.eps
    Deleted:    figures/figure1.png
    Modified:   figures/figure2.eps

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/index.Rmd) and HTML (docs/index.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
html 836a780 lukeholman 2021-03-04 Build site.
html 359ff37 lukeholman 2021-03-04 Build site.
html 5506c4b lukeholman 2021-03-04 Build site.
Rmd 390a393 lukeholman 2021-03-04 big first commit 2021
Rmd 8d54ea5 Luke Holman 2018-12-23 Initial commit
html 8d54ea5 Luke Holman 2018-12-23 Initial commit

Running the analysis

1. Setting up a database to hold variant/gene annotations and GWAS results

This script creates a SQLite3 database holding two tables: one with annotations for each variant (provided by the Mackay lab), and one with annotations for each gene (from annotation hub). In step 3, we also add the GWAS results to this database, allowing memory-efficient handling of the results. Using this database also avoids having to load Bioconductor packages that conflict with dplyr::select(), etc.

2. Estimating line mean fitness using Bayesian models

This script uses Bayesian mixed models implemented in brms to estimate the line means for our four fitness traits.

3. Calculating quantitative genetic parameters

Using linear mixed models implemented in the R package sommer to estimate the heritability of, and genetic covariance among, our four fitness traits.

4. Running the GWAS

The script first performs quality control and imputation on the dataset of SNPs and indels for the DGRP. Second, it runs association tests on our four phenotypes using the software command-line software GEMMA. Third, it uses PLINK to identify groups of loci that are in complete linkage disequilibrium (‘SNP clumps’).

5. Applying adaptive shrinkage to the GWAS results

5a. Applying adaptive shrinkage to the GWAS results

Here, we use the R package mashr to perform multivariate adaptive shrinkage on the results of the GWAS. We ran mashr using both possible modes: canonical, and data-driven. The canonical approach is useful for testing hypotheses about the relative abundance of different types of fitness-affecting SNPs, e.g. sexually antagonistic, sexually concordant, sex-limited, or null. The data-driven approach adjusts the GWAS estimates most accurately, as it estimates the patterns of covariances in the SNPs’ effects on the 4 traits from the data; this is useful for deriving shrinked estimates of all the SNPs’ effect sizes.

5b. Plots illustrating how adaptive shrinkage affected the results

This script generates some plots of the SNP effect sizes under the two forms of shrinkage (canonical and data-driven) to check that mashr is performing as expected.

6. Running the GWAS, identifying eQTLs, and applying mashr to the TWAS results

This script uses the transcriptomic data from Huang et al. 2015 PNAS on the DGRP. In this script, we

  • Calculate the average expression and sex bias in expression for each transcript in Huang et al.
  • Identify loci that associated with fitness in our GWAS, and which were also identified as eQTLs in Huang et al.
  • Perform ‘TWAS’, i.e. searching for transcripts whose expression correlates with our fitness measures across lines
  • Adjust the results of the TWAS linear models using mashr in canonical and data-driven modes.

Results of the analysis

5. Plot and model the variant effect sizes

Here, we plot the correlations in the estimated effects of each variant (or group of variants in perfect linkage disequilibrium with one another). For the GWAS results, we run linear models to test whether variants differ in their effect on fitness according to chromosome, site class, and minor allele frequency. For the TWAS results, we run linear models to test whether transcripts differ in their effect on fitness according to chromosome, average expression level, and the sex bias in expression.

6. GO and KEGG enrichment analyses

Text here

Heritability - Correlated with selection experienced by the gene? Angatonism? - Correlated with sex-specificity in expression?

eQTLs - Are there eQTLs with matching, independent, and opposite-sex effects? YES TO ALL. NEED TO RE-DO eQTL ANALYSIS - Are there more overlaps between the eQTLs and fitness QTLs than expected by chance? Correlation in p-vals?

Predictions - eQTLs that affect transciption the same way in both sexes will have SA fitness effects for transcripts under SA selection (or concordant for concordant) We need: - Effect size for the eQTL on the transcript, in both sexes (re-run eQTL) - Effect size for the eQTL on fitness, in both sexes (done) - Effect size for the transcript on fitness (done) - “Allele 1 made the transcript increase in both sexes, good for males and bad for females” - eQTLs that only affect transciption in one sex will have a more sex-specific fitness effect as well (not concordant, not antag) - eQTLs where the + and - alleles are swapped between sexes will be antagonistic or concordant, depending on sex-specific selection on that transcript - Transcripts will show the twin peaks relationship between sex bias in expression and SA selection - Genes targetted by dsx will be extra antagonistic