Summary of analysis plan

• predict whole blood expression
• check how well the prediction works with GEUVADIS expression data
• run association between predicted expression and a simulated phenotype
• calculate association between expression levels and coronary artery disease risk using s-predixcan
• fine-map the coronary artery disease gwas results using torus (need some preformatting)
• calculate colocalization probability using fastenloc
• run transcriptome-wide mendelian randomization in one locus of interest

library(tidyverse)

── Attaching packages ──────────────────────────────────────────────── tidyverse 1.3.0 ──

✓ ggplot2 3.3.0     ✓ purrr   0.3.4
✓ tibble  3.0.0     ✓ dplyr   0.8.5
✓ tidyr   1.0.2     ✓ stringr 1.4.0
✓ readr   1.3.1     ✓ forcats 0.5.0

── Conflicts ─────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

predict expression

Visual summary of predixcan runs

• Remember you need to copy and paste this code chunk into the terminal to run it. Also make sure you activated the imlabtools environment which has all the necessary python modules.

• Make sure all the paths and file names are correct. This run should take about one minute.

<!-- PRE="/home/student/QGT-Columbia-HKI" -->
<!-- DATA=$PRE/data -->
<!-- MODEL=$PRE/models -->
<!-- METAXCAN=$PRE/repos/MetaXcan-master/software -->
<!-- RESULTS=$PRE/results -->

printf "Predict expression\n\n"
python3 $METAXCAN/Predict.py \
--model_db_path $PRE/models/gtex_v8_en/en_Whole_Blood.db \
--vcf_genotypes $DATA/predixcan/genotype/filtered.vcf.gz \
--vcf_mode genotyped \
--variant_mapping $DATA/predixcan/gtex_v8_eur_filtered_maf0.01_monoallelic_variants.txt.gz id rsid \
--on_the_fly_mapping METADATA "chr{}_{}_{}_{}_b38" \
--prediction_output $RESULTS/predixcan/Whole_Blood__predict.txt \
--prediction_summary_output $RESULTS/predixcan/Whole_Blood__summary.txt \
--verbosity 9 \
--throw

assess prediction performance (optional)

predicted_expression = read_tsv(glue::glue("{results.dir}/predixcan/Whole_Blood__predict.txt"))
dim(predicted_expression)
head(predicted_expression[,1:5])
prediction_summary = read_tsv(glue::glue("{results.dir}/predixcan/Whole_Blood__summary.txt"))
dim(prediction_summary)
head(prediction_summary)

## merge with GEUVADIS expression data

## calculate spearman correlation

## select a few genes and plot predicted vs observed expression

run association with a phenotype

PHENO="sim.infinitesimal_pve0.1"
printf "association\n\n"
python3 $METAXCAN/PrediXcanAssociation.py \
--expression_file $RESULTS/predixcan/Whole_Blood__predict.txt \
--input_phenos_file $DATA/predixcan/phenotype/$PHENO.txt \
--input_phenos_column pheno \
--output $RESULTS/predixcan/$PHENO/Whole_Blood__association.txt \
--verbosity 9 \
--throw

PHENO="sim.spike_n_slab_0.01_pve0.1"
printf "association\n\n"
python3 $METAXCAN/PrediXcanAssociation.py \
--expression_file $RESULTS/predixcan/Whole_Blood__predict.txt \
--input_phenos_file $DATA/predixcan/phenotype/$PHENO.txt \
--input_phenos_column pheno \
--output $RESULTS/predixcan/$PHENO/Whole_Blood__association.txt \
--verbosity 9 \
--throw


PHENO="sim.spike_n_slab_0.1_pve0.05"
printf "association\n\n"
python3 $METAXCAN/PrediXcanAssociation.py \
--expression_file $RESULTS/predixcan/Whole_Blood__predict.txt \
--input_phenos_file $DATA/predixcan/phenotype/$PHENO.txt \
--input_phenos_column pheno \
--output $RESULTS/predixcan/$PHENO/Whole_Blood__association.txt \
--verbosity 9 \
--throw

read association results

## read association results
PHENO="sim.spike_n_slab_0.01_pve0.1"


predixcan_association = read_tsv(glue::glue("{results.dir}/predixcan/{PHENO}/Whole_Blood__association.txt"))
## take a look at the results
dim(predixcan_association)
predixcan_association %>% arrange(pvalue) %>% head
predixcan_association %>% arrange(pvalue) %>% ggplot(aes(pvalue)) + geom_histogram(bins=20)
## compare distribution against the null (uniform)
qq(predixcan_association$pvalue)

truebetas = read_tsv(glue::glue("{DATA}/predixcan/phenotype/gene-effects/{PHENO}.txt"))
predixcan_association %>% inner_join(truebetas,by=c("gene"="gene_id")) %>% ggplot(aes(effect,effect_size))+geom_point()+geom_abline()

Exercise

-[ ] show top genes -[ ] compare with true effect sizes

library(qqman)

For example usage please run: vignette('qqman')

Citation appreciated but not required:

Turner, S.D. qqman: an R package for visualizing GWAS results using Q-Q and manhattan plots. biorXiv DOI: 10.1101/005165 (2014).

——-

run s-predixcan

export PRE="/home/student/QGT-Columbia-HKI"
export DATA=$PRE/data/s-predixcan
export MODEL=$PRE/models
export METAXCAN=$PRE/repos/MetaXcan-master/software
export RESULTS=$PRE/results

python $METAXCAN/SPrediXcan.py \
--gwas_file  $DATA/imputed_CARDIoGRAM_C4D_CAD_ADDITIVE.txt.gz \
--snp_column panel_variant_id --effect_allele_column effect_allele --non_effect_allele_column non_effect_allele --zscore_column zscore \
--model_db_path $MODEL/gtex_v8_mashr/mashr_Whole_Blood.db \
--covariance $MODEL/gtex_v8_mashr/mashr_Whole_Blood.txt.gz \
--keep_non_rsid --additional_output --model_db_snp_key varID \
--throw \
--output_file $RESULTS/spredixcan/eqtl/CARDIoGRAM_C4D_CAD_ADDITIVE__PM__Whole_Blood.csv

plot and interpret s-predixcan results

spredixcan_association = read_csv(glue::glue("{results.dir}/spredixcan/eqtl/CARDIoGRAM_C4D_CAD_ADDITIVE__PM__Whole_Blood.csv"))
dim(spredixcan_association)
spredixcan_association %>% arrange(pvalue) %>% head
spredixcan_association %>% arrange(pvalue) %>% ggplot(aes(pvalue)) + geom_histogram(bins=20)

SORT1, considered to be a causal gene for LDL cholesterol and as a consequence of coronary artery disease, is not found here. Why? (tissue)

run multixcan (optional)

export MODEL=$PRE/models
export DATA=$PRE/data/s-predixcan

python $METAXCAN/SMulTiXcan.py \
--models_folder $MODEL/gtex_v8_mashr \
--models_name_pattern "mashr_(.*).db" \
--snp_covariance $MODEL/gtex_v8_expression_mashr_snp_covariance.txt.gz \
--metaxcan_folder $RESULTS/spredixcan/eqtl/ \
--metaxcan_filter "CARDIoGRAM_C4D_CAD_ADDITIVE__PM__(.*).csv" \
--metaxcan_file_name_parse_pattern "(.*)__PM__(.*).csv" \
--gwas_file $DATA/imputed_CARDIoGRAM_C4D_CAD_ADDITIVE.txt.gz \
--snp_column panel_variant_id --effect_allele_column effect_allele --non_effect_allele_column non_effect_allele --zscore_column zscore --keep_non_rsid --model_db_snp_key varID \
--cutoff_condition_number 30 \
--verbosity 7 \
--throw \
--output $RESULTS/smultixcan/eqtl/CARDIoGRAM_C4D_CAD_ADDITIVE_smultixcan.txt

fine-map GWAS results

We will run torus due to time limitation but ideally we would like to run a method that allows multiple causal variants per locus.

#torus -d Height.torus.zval.gz --load_zval -dump_pip Height.gwas.pip
#gzip Height.gwas.pip
TORUSOFT=torus

$TORUSOFT -d $PRE/data/fastenloc/Height.torus.zval.gz --load_zval -dump_pip $PRE/data/fastenloc/Height.gwas.pip
cd $PRE/data/fastenloc
gzip Height.gwas.pip
cd $PRE 

calculate colocalization with fastENLOC

## check out tutorial https://github.com/xqwen/fastenloc/tree/master/tutorial

export eqtl_annotation_gzipped=$PRE/data/fastenloc/FASTENLOC-gtex_v8.eqtl_annot.vcf.gz
export gwas_data_gzipped=$PRE/data/fastenloc/Height.gwas.pip.gz
export TISSUE=Whole_Blood
export FASTENLOCSOFT=fastenloc

mkdir $RESULTS/fastenloc/
cd $RESULTS/fastenloc/
$FASTENLOCSOFT -eqtl $eqtl_annotation_gzipped -gwas $gwas_data_gzipped -t $TISSUE 

#[-total_variants total_snp] [-thread n] [-prefix prefix_name] [-s shrinkage]

analyze results

## optional - compare with s-predixcan results

-[] prepare

run SMR (optional)

run TWMR (for a locus)

TWMR

export TWMR=$PRE/repos/TWMR-master
export OUTPUT=$PRE/results
GENE=ENSG00000002919

cd $TWMR

R < $TWMR/MR.R --no-save $GENE

cd $PRE

## output: /home/student/QGT-Columbia-HKI/repos/TWMR-master/ENSG00000002919.alpha

Session information

sessionInfo()

R version 3.6.3 (2020-02-29)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Catalina 10.15.2

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/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] qqman_0.1.4     forcats_0.5.0   stringr_1.4.0   dplyr_0.8.5    
 [5] purrr_0.3.4     readr_1.3.1     tidyr_1.0.2     tibble_3.0.0   
 [9] ggplot2_3.3.0   tidyverse_1.3.0

loaded via a namespace (and not attached):
 [1] tidyselect_1.0.0 xfun_0.13        haven_2.2.0      lattice_0.20-41 
 [5] colorspace_1.4-1 vctrs_0.2.4      generics_0.0.2   htmltools_0.4.0 
 [9] yaml_2.2.1       rlang_0.4.5      later_1.0.0      pillar_1.4.3    
[13] withr_2.1.2      glue_1.4.1       DBI_1.1.0        dbplyr_1.4.2    
[17] calibrate_1.7.5  readxl_1.3.1     modelr_0.1.6     lifecycle_0.2.0 
[21] cellranger_1.1.0 munsell_0.5.0    gtable_0.3.0     workflowr_1.6.2 
[25] rvest_0.3.5      evaluate_0.14    knitr_1.28       httpuv_1.5.2    
[29] fansi_0.4.1      broom_0.5.5      Rcpp_1.0.4.6     promises_1.1.0  
[33] backports_1.1.6  scales_1.1.0     jsonlite_1.6.1   fs_1.4.1        
[37] hms_0.5.3        digest_0.6.25    stringi_1.4.6    rprojroot_1.3-2 
[41] grid_3.6.3       cli_2.0.2        tools_3.6.3      magrittr_1.5    
[45] crayon_1.3.4     whisker_0.4      pkgconfig_2.0.3  MASS_7.3-51.5   
[49] ellipsis_0.3.0   xml2_1.3.1       reprex_0.3.0     lubridate_1.7.8 
[53] rstudioapi_0.11  assertthat_0.2.1 rmarkdown_2.1    httr_1.4.1      
[57] R6_2.4.1         nlme_3.1-147     git2r_0.26.1     compiler_3.6.3