predict expression

Visual summary of predixcan runs

conda activate imlabtools

export PRE="/home/student/QGT-Columbia-HKI"
export DATA=$PRE/data/predixcan
export MODEL=$PRE/models
export METAXCAN=$PRE/repos/MetaXcan-master/software
export 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/genotype/ALL.chr22.shapeit2_integrated_snvindels_v2a_27022019.GRCh38.phased.vcf.gz \
--vcf_mode genotyped \
--variant_mapping $DATA/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 phenotype

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

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

read results

predixcan_association = read_tsv(glue::glue("{results.dir}/predixcan/random_pheno/Whole_Blood__association.txt"))
dim(predixcan_association)
predixcan_association %>% arrange(pvalue) %>% head
predixcan_association %>% arrange(pvalue) %>% ggplot(aes(pvalue)) + geom_histogram(bins=20)

Exercise

-[ ] Run association with another phenotype in $PRE/predixcan/data/phenotype/ALL.chr22.shapeit2_integrated_snvindels_v2a_27022019.GRCh38_x_en_Whole_Blood.simulated_phenotype.spike_n_slab_0.5_x_pve0.6.txt

——-

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.2 (2019-12-12)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6

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] forcats_0.4.0   stringr_1.4.0   dplyr_0.8.3     purrr_0.3.3    
 [5] readr_1.3.1     tidyr_1.0.2     tibble_2.1.3    ggplot2_3.2.1  
 [9] tidyverse_1.3.0 workflowr_1.6.0

loaded via a namespace (and not attached):
 [1] tidyselect_1.0.0 xfun_0.12        haven_2.2.0      lattice_0.20-38 
 [5] colorspace_1.4-1 vctrs_0.2.2      generics_0.0.2   htmltools_0.4.0 
 [9] yaml_2.2.0       rlang_0.4.3      later_1.0.0      pillar_1.4.3    
[13] withr_2.1.2      glue_1.3.1       DBI_1.1.0        dbplyr_1.4.2    
[17] modelr_0.1.5     readxl_1.3.1     lifecycle_0.1.0  munsell_0.5.0   
[21] gtable_0.3.0     cellranger_1.1.0 rvest_0.3.5      evaluate_0.14   
[25] knitr_1.27       httpuv_1.5.2     fansi_0.4.1      broom_0.5.4     
[29] Rcpp_1.0.3       promises_1.1.0   backports_1.1.5  scales_1.1.0    
[33] jsonlite_1.6     fs_1.3.1         hms_0.5.3        digest_0.6.23   
[37] stringi_1.4.5    grid_3.6.2       rprojroot_1.3-2  cli_2.0.1       
[41] tools_3.6.2      magrittr_1.5     lazyeval_0.2.2   crayon_1.3.4    
[45] whisker_0.4      pkgconfig_2.0.3  xml2_1.2.2       reprex_0.3.0    
[49] lubridate_1.7.4  rstudioapi_0.10  assertthat_0.2.1 rmarkdown_2.1   
[53] httr_1.4.1       R6_2.4.1         nlme_3.1-143     git2r_0.26.1    
[57] compiler_3.6.2