An example for LDL ~ Predictdb eQTL&sQTL from Liver
XSun
2024-09-22
Last updated: 2024-09-23
Checks: 6 1
Knit directory: multigroup_ctwas_analysis/
This reproducible R Markdown analysis was created with workflowr (version 1.7.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
The R Markdown is untracked by Git. To know which version of the R
Markdown file created these results, you’ll want to first commit it to
the Git repo. If you’re still working on the analysis, you can ignore
this warning. When you’re finished, you can run
wflow_publish to commit the R Markdown file and build the
HTML.
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.
The command set.seed(20231112) 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.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
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 d7cd32f. 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: .Rhistory
Ignored: analysis/figure/
Ignored: results/
Untracked files:
Untracked: analysis/LDL_predictdb_esQTL.Rmd
Unstaged changes:
Modified: analysis/index.Rmd
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.
There are no past versions. Publish this analysis with
wflow_publish() to start tracking its development.
We present a sample cTWAS report based on real data analysis. The analyzed trait is LDL cholesterol, the prediction models are liver gene expression and splicing models trained on GTEx v8 in the PredictDB format.
Analysis settings
Input data
- GWAS Z-scores
The analyzed trait is LDL cholesterol. The summary statistics are
downloaded from https://gwas.mrcieu.ac.uk, using dataset ID:
ukb-d-30780_irnt. The number of SNPs it contains is
13,586,016.
The sample size is
[1] "gwas_n = 343621"- Prediction models
The prediction models used in this analysis are liver gene expression and splicing models, trained on GTEx v8 in the PredictDB format. These models can be downloaded from https://predictdb.org/post/2021/07/21/gtex-v8-models-on-eqtl-and-sqtl/
[1] "The number of eQTLs per gene = 1.5078"[1] "Total genes = 12714"[1] "The number of sQTLs per intron = 1.2151"[1] "Total introns = 29250"- Reference data
The reference data include genomic region definitions and an LD
reference. We use the genomic regions provided by the package and the LD
reference in b38, located at
/project2/mstephens/wcrouse/UKB_LDR_0.1/. Alternatively,
the LD reference can be downloaded from this link:https://uchicago.app.box.com/s/jqocacd2fulskmhoqnasrknbt59x3xkn.
Data processing and harmonization
We map the reference SNPs and LD matrices to regions following the instructions from the cTWAS tutorial.
When processing z-scores, we exclude multi-allelic and
strand-ambiguous variants by setting
drop_multiallelic = TRUE and
drop_strand_ambig = TRUE.
Users can expand the code snippets below to view the exact code used.
## input data
weight_files <- c("/project2/xinhe/shared_data/multigroup_ctwas/weights/expression_models/expression_Liver.db","/project2/xinhe/shared_data/multigroup_ctwas/weights/splicing_models/splicing_Liver.db")
z_snp_file <- "/project2/xinhe/shared_data/multigroup_ctwas/gwas/ctwas_inputs_zsnp/LDL-ukb-d-30780_irnt.z_snp.RDS"
region_file <- system.file("extdata/ldetect", paste0("EUR.", genome_version, ".ldetect.regions.RDS"), package = "ctwas")
region_info <- readRDS(region_file)
genome_version <- "b38"
LD_dir <- "/project2/mstephens/wcrouse/UKB_LDR_0.1/"
## output dir
outputdir <- "/project/xinhe/xsun/multi_group_ctwas/examples/results_predictdb_main/LDL-ukb-d-30780_irnt/"
dir.create(outputdir, showWarnings=F, recursive=T)
## other parameters
ncore <- 5
## process inputs
### Preprocess LD_map & SNP_map
region_metatable <- region_info
region_metatable$LD_file <- file.path(LD_dir, paste0(LD_filestem, ".RDS"))
region_metatable$SNP_file <- file.path(LD_dir, paste0(LD_filestem, ".Rvar"))
res <- create_snp_LD_map(region_metatable)
region_info <- res$region_info
snp_map <- res$snp_map
LD_map <- res$LD_map
### Preprocess GWAS z-scores
z_snp <- readRDS(z_snp_file)
z_snp <- preprocess_z_snp(z_snp = z_snp,
snp_map = snp_map,
drop_multiallelic = TRUE,
drop_strand_ambig = TRUE)
### Preprocess weights
weights_expression1 <- preprocess_weights(weight_file = weight_files[1],
region_info = region_info,
gwas_snp_ids = z_snp$id,
snp_map = snp_map,
LD_map = LD_map,
type = "eQTL",
context = tissue,
weight_format = "PredictDB",
drop_strand_ambig = TRUE,
scale_predictdb_weights = T,
load_predictdb_LD = F, #### F for fusion converted weights
filter_protein_coding_genes = TRUE,
ncore = ncore)
weights_splicing1 <- preprocess_weights(weight_file = weight_files[2],
region_info = region_info,
gwas_snp_ids = z_snp$id,
snp_map = snp_map,
LD_map = LD_map,
type = "sQTL",
context = tissue,
weight_format = "PredictDB",
drop_strand_ambig = TRUE,
scale_predictdb_weights = T, #### F for fusion converted weights
load_predictdb_LD = F,
filter_protein_coding_genes = TRUE,
ncore = ncore)
weights <- c(weights_expression1,weights_splicing1) Running cTWAS analysis
We use the function ctwas_sumstats() to run cTWAS analysis with LD. Details are shown in the tutorial. https://xinhe-lab.github.io/multigroup_ctwas/articles/running_ctwas_analysis.html#running-ctwas-main-function
The arguments are all in defaul settings, more specifically,
- we set
group_prior_var_structure = "shared_type"to allow all groups in one molecular QTL type to share the same variance parameter - we first estimate the number of causal signals (L) for each region
by setting
filter_L = TRUE - The package will not compute the non-SNP PIP by setting
filter_nonSNP_PIP = TRUE - We select regions with non-SNP PIP > 0.5 by setting
min_nonSNP_PIP = 0.5
We use the function ctwas_sumstats() to run the cTWAS
analysis with LD. For more details on this function, refer to the cTWAS
tutorial: https://xinhe-lab.github.io/multigroup_ctwas/articles/running_ctwas_analysis.html#running-ctwas-main-function
All arguments are set to their default values, with the following specific settings:
group_prior_var_structure = "shared_type": Allows all groups within a molecular QTL type to share the same variance parameter.filter_L = TRUE: Estimates the number of causal signals (L) for each region.filter_nonSNP_PIP = TRUE: Prevents the computation of non-SNP PIP values.min_nonSNP_PIP = 0.5: Selects regions where the non-SNP PIP is greater than 0.5.
Users can expand the code snippets below to view the exact code used.
thin <- 0.1
maxSNP <- 20000
ctwas_res <- ctwas_sumstats(z_snp,
weights,
region_info,
LD_map,
snp_map,
thin = thin,
maxSNP = maxSNP,
group_prior_var_structure = "shared_type",
filter_L = TRUE,
filter_nonSNP_PIP = FALSE,
min_nonSNP_PIP = 0.5,
ncore = ncore,
ncore_LD = ncore,
save_cor = TRUE,
cor_dir = paste0(outputdir,"/cor_matrix"),
verbose = T)cTWAS results
ctwas_res is the object contains the outputs of
cTWAS
Parameter estimation
We extract the estimated parameters by
param <- ctwas_res$param
we make plots using the function
make_convergence_plots(param, gwas_n) to see how estimated
parameters converge during the execution of the program:
These plots show the estimated prior inclusion probability, prior effect size variance, enrichment and proportion of variance explained (PVE) over the iterations of parameter estimation. The enrichment is defined as the ratio of the prior inclusion probability of molecular traits over the prior inclusion probability of variants. We generally expect molecular traits to have higher prior inclusion probability than variants. Enrichment values typically range from 20 - 100 for expression traits.
Then, we use summarize_param(param, gwas_n) to obtain
estimated parameters (from the last iteration) and to compute the PVE by
variants and molecular traits.
[1] "The number of genes/introns/SNPs used in the analysis is:"Liver|eQTL Liver|sQTL SNP
8775 18136 7405450 ctwas_parameters$attributable_pve contains the
proportion of heritability mediated by molecular traits and variants, we
visualize it using pie chart.
data <- data.frame(
category = names(ctwas_parameters$attributable_pve),
percentage = ctwas_parameters$attributable_pve
)
# Calculate percentage labels for the chart
data$percentage_label <- paste0(round(data$percentage * 100, 1), "%")
ggplot(data, aes(x = "", y = percentage, fill = category)) +
geom_bar(stat = "identity", width = 1) +
coord_polar("y", start = 0) +
theme_void() + # Remove background and axes
geom_text(aes(label = percentage_label),
position = position_stack(vjust = 0.5), size = 5) +
scale_fill_manual(values = c("#FF9999", "#66B2FF", "#99FF99")) + # Custom colors
labs(fill = "Category") +
ggtitle("Attributable PVE")
Fine-mapping results
We process the fine-mapping results here.
We first add gene annotations to cTWAS results
load("/project2/xinhe/shared_data/multigroup_ctwas/weights/E_S_A_mapping_updated.RData")
colnames(E_S_A_mapping)[1] <- "molecular_id"
finemap_res <- ctwas_res$finemap_res
finemap_res$molecular_id <- get_molecular_ids(finemap_res)
snp_map <- readRDS(paste0(results_dir,trait,".snp_map.RDS"))
finemap_res <- anno_finemap_res(finemap_res,
snp_map = snp_map,
mapping_table = E_S_A_mapping,
add_gene_annot = TRUE,
map_by = "molecular_id",
drop_unmapped = TRUE,
add_position = TRUE,
use_gene_pos = "mid")2024-09-23 16:25:07 INFO::Annotating fine-mapping result ...
2024-09-23 16:25:08 INFO::Map molecular traits to genes
2024-09-23 16:25:11 INFO::Split PIPs for molecular traits mapped to multiple genes
2024-09-23 16:25:16 INFO::Add gene positions
2024-09-23 16:25:16 INFO::Add SNP positionsDT::datatable(finemap_res[finemap_res$susie_pip > 0.8,],caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','The annotated fine-mapping results, ones with susie pip > 0.8 are shown'),options = list(pageLength = 5) )For all genes analyzed, we compare the z-scores and fine-mapping PIPs
ggplot(data = finemap_res[finemap_res$type!="SNP",], aes(x = abs(z), y = susie_pip)) +
geom_point() +
labs(x = "abs(z-scores)", y = "PIPs") +
theme_minimal()
Next, we compute gene PIPs across different types of molecular traits
combined_pip_by_type <- combine_gene_pips(finemap_res,
group_by = "gene_name",
by = "type",
method = "combine_cs",
filter_cs = TRUE)2024-09-23 16:25:30 INFO::Limit gene results to credible setsDT::datatable(combined_pip_by_type,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','Gene PIPs, only genes within credible sets are shown'),options = list(pageLength = 5) )Visualization
- The top one shows -log10(p-value) of the association of variants (from LDL GWAS) and molecular traits (from the package computed z-scores) with the phenotype
- The next track shows the PIPs of variants and molecular traits. By
default, we limit PIP results to credible sets in the PIP track
(
filter_cs = TRUE) - The next track shows the QTLs of the focal gene.
- The bottom is the gene track.
sessionInfo()R version 4.2.0 (2022-04-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)
Matrix products: default
BLAS/LAPACK: /software/openblas-0.3.13-el7-x86_64/lib/libopenblas_haswellp-r0.3.13.so
locale:
[1] C
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] ggplot2_3.5.1 EnsDb.Hsapiens.v86_2.99.0
[3] ensembldb_2.20.2 AnnotationFilter_1.20.0
[5] GenomicFeatures_1.48.3 AnnotationDbi_1.58.0
[7] Biobase_2.56.0 GenomicRanges_1.48.0
[9] GenomeInfoDb_1.39.9 IRanges_2.30.0
[11] S4Vectors_0.34.0 BiocGenerics_0.42.0
[13] ctwas_0.4.14
loaded via a namespace (and not attached):
[1] colorspace_2.0-3 rjson_0.2.21
[3] ellipsis_0.3.2 rprojroot_2.0.3
[5] XVector_0.36.0 locuszoomr_0.2.1
[7] fs_1.5.2 rstudioapi_0.13
[9] farver_2.1.0 DT_0.22
[11] ggrepel_0.9.1 bit64_4.0.5
[13] fansi_1.0.3 xml2_1.3.3
[15] codetools_0.2-18 logging_0.10-108
[17] cachem_1.0.6 knitr_1.39
[19] jsonlite_1.8.0 workflowr_1.7.0
[21] Rsamtools_2.12.0 dbplyr_2.1.1
[23] png_0.1-7 readr_2.1.2
[25] compiler_4.2.0 httr_1.4.3
[27] assertthat_0.2.1 Matrix_1.5-3
[29] fastmap_1.1.0 lazyeval_0.2.2
[31] cli_3.6.1 later_1.3.0
[33] htmltools_0.5.2 prettyunits_1.1.1
[35] tools_4.2.0 gtable_0.3.0
[37] glue_1.6.2 GenomeInfoDbData_1.2.8
[39] dplyr_1.1.4 rappdirs_0.3.3
[41] Rcpp_1.0.12 jquerylib_0.1.4
[43] vctrs_0.6.5 Biostrings_2.64.0
[45] rtracklayer_1.56.0 crosstalk_1.2.0
[47] xfun_0.41 stringr_1.5.1
[49] lifecycle_1.0.4 irlba_2.3.5
[51] restfulr_0.0.14 XML_3.99-0.14
[53] zlibbioc_1.42.0 zoo_1.8-10
[55] scales_1.3.0 gggrid_0.2-0
[57] hms_1.1.1 promises_1.2.0.1
[59] MatrixGenerics_1.8.0 ProtGenerics_1.28.0
[61] parallel_4.2.0 SummarizedExperiment_1.26.1
[63] LDlinkR_1.2.3 yaml_2.3.5
[65] curl_4.3.2 memoise_2.0.1
[67] sass_0.4.1 biomaRt_2.54.1
[69] stringi_1.7.6 RSQLite_2.3.1
[71] highr_0.9 BiocIO_1.6.0
[73] filelock_1.0.2 BiocParallel_1.30.3
[75] rlang_1.1.2 pkgconfig_2.0.3
[77] matrixStats_0.62.0 bitops_1.0-7
[79] evaluate_0.15 lattice_0.20-45
[81] purrr_1.0.2 labeling_0.4.2
[83] GenomicAlignments_1.32.0 htmlwidgets_1.5.4
[85] cowplot_1.1.1 bit_4.0.4
[87] tidyselect_1.2.0 magrittr_2.0.3
[89] R6_2.5.1 generics_0.1.2
[91] DelayedArray_0.22.0 DBI_1.2.2
[93] withr_2.5.0 pgenlibr_0.3.3
[95] pillar_1.9.0 KEGGREST_1.36.3
[97] RCurl_1.98-1.7 mixsqp_0.3-43
[99] tibble_3.2.1 crayon_1.5.1
[101] utf8_1.2.2 BiocFileCache_2.4.0
[103] plotly_4.10.0 tzdb_0.4.0
[105] rmarkdown_2.25 progress_1.2.2
[107] grid_4.2.0 data.table_1.14.2
[109] blob_1.2.3 git2r_0.30.1
[111] digest_0.6.29 tidyr_1.3.0
[113] httpuv_1.6.5 munsell_0.5.0
[115] viridisLite_0.4.0 bslib_0.3.1