Ganymede genetic report for: IL22
2025-01-03
knitr::opts_chunk$set(fig.width = 12, fig.height = 8, fig.retina = 2, error = TRUE)
# Load libs
library(callisto)
library(kableExtra)
library(knitr)
library(DT)
# Set logging to errors only
futile.logger::flog.threshold(ERROR)
# Include this here to suppress a warning
DT::datatable(data = tibble(var = 1:3))
# Connect to Ganymede
dbc <- db_connect()
# QC that the input_gene is not the default value
if(params$input_gene == "DOUBLE_CLICK_HERE_TO_CHANGE"){
log_fatal_stop("The parameter `input_gene` is still the default value (DOUBLE_CLICK_HERE_TO_CHANGE). Please try again to pass a valid gene_name to the notebook.")
}
# QC that the input_gene is in our database
if(nrow(find_gene_names(name_like = params$input_gene)) == 0){
log_fatal_stop("The gene input: ({params$input_gene}) isn't found in our database, please check the input.")
}
tags_2_prioritize <-
tibble(tags1 = c("qtls_all", "qtls_imm", "qtls_eec", "qtls_gtex", "qtls_pqtl", "qtls_pgx"),
tags1_prio = 1:6)callisto_version()## [1] "Package: callisto|Version: 0.8.1|Built: R 3.6.1; ; 2025-01-03 13:49:49 UTC; unix|Packaged: 2025-01-03 13:49:48 UTC; ubuntu"cat(glue("Report version: {params$auto_rpt_version}"))## Report version: 2.01 Available instruments
These genetic instruments are potential hypotheses for how a gene can be linked to disease. This section is to help us understand how many “shots on goal” we have for a given gene. These instruments don’t represent robust variant to gene mapping.
1.1 Number of instruments
# Print output summary
tribble(~data_type, ~num_instruments,
"Coding variants", n_coding_vars,
"QTL tissue/cell types", n_tissue_qtls,
"Nearby GWAS", n_nearby_gwas) %>%
DT::datatable(rownames = FALSE, filter = "top")1.2 List of nearby GWAS
Nearby GWAS are genome wide significant top hits within +/- 1 MB of IL22.
if(n_nearby_gwas > 0){
nearby_gwas %>%
DT::datatable(rownames = FALSE, filter = "top", options = list(pageLength = 15, scrollX = TRUE))
} else {
cat("No nearby GWAS signals to list.")
}2 Regional architecture
This is a high level illustration of the variant to gene (V2G) mapping at a genomic region across all genes and phenotypes. More detailed follow up interogation of the v2g will come next.
if(!exists("all_v2g_data") | all_v2g_input_rows == 0){
coding_eval_lgl <- FALSE
coloc_eval_lgl <- FALSE
near_eval_lgl <- FALSE
cat(glue("Ganymede doesn't have any variant to gene mapping for: {params$input_gene}"))
} else {
# Use for future eval checks
coding_eval_lgl <-
all_v2g_data %>%
filter(gene_name == params$input_gene &
v2g_type == "coding" &
v2g_value == 1) %>%
nrow() > 0
coloc_eval_lgl <-
all_v2g_data %>%
filter(gene_name == params$input_gene &
v2g_type == "coloc" &
v2g_value == 1) %>%
nrow() > 0
near_eval_lgl <-
all_v2g_data %>%
filter(gene_name == params$input_gene &
v2g_type == "nearest" &
v2g_value == 1) %>%
nrow() > 0
}2.1 IL22 & competing hypothesis only
if(!exists("all_v2g_data")){
cat(glue("Ganymede doesn't have any variant to gene mapping for: {params$input_gene}"))
} else if(params$input_gene %in% all_v2g_data$gene_name){
gather_region_arch.plot(all_v2g_data) %>%
gather_region_arch.zoom(., gene_name = params$input_gene)
} else {
cat(glue("Ganymede doesn't have any variant to gene mapping for: {params$input_gene}"))
}
2.2 All nearby V2G
# Plot the data
gather_region_arch.plot(all_v2g_data)
3 Coding evidence
Coding variant to gene mapping occurs when there is a high or medium consequence coding variant is within the credible set and has a posterior probability of causality >= 1%. Coding variants impact on the gene(s) should be confirmed using external resources.
3.1 Phenome-wide coding evidence
Illustrates the coding variant(s) (right y-axis) that are in the credible set and >= 1% probability in any disease or biomarker credible set available in Ganymede (y axis-left). Foreach variant and phenotype the odd’s ratio (x-axis), the variant impact (shape) are illustrated. In addition, the ratio of the variant’s probability relative to the credible set max probability is illustrated by the color to illustrate the liklihood that the variant is driving genetic signal.
3.1.1 Latest disease & biomarkers
if(coding_eval_lgl == TRUE){
if(is.data.frame(gwas_coding_results) & !is_empty(gwas_coding_results)){
instrument_gene_coding.plot(gwas_coding_results)
}
} else {
cat(glue("No coding variant to gene mapping for: {params$input_gene}."))
}## No coding variant to gene mapping for: IL22.3.1.2 cis-QTLs
if(!is_empty(cis_coding) & nrow(cis_coding) > 0){
instrument_gene_coding.plot(cis_coding)
} else {
cat(glue("No coding cis-QTL variant to gene mapping for: {params$input_gene}."))
}## Error in is_empty(cis_coding): object 'cis_coding' not found3.1.3 trans-QTLs
if(!is_empty(trans_coding) & nrow(trans_coding) > 0){
instrument_gene_coding.plot(trans_coding)
} else {
cat(glue("No coding trans-QTL variant to gene mapping for: {params$input_gene}."))
}## Error in is_empty(trans_coding): object 'trans_coding' not found3.2 Phenotypes
# TODO: Account for conditional signals in future
if(coding_eval_lgl == FALSE){
cat(glue("No coding variant to gene mapping for: {params$input_gene}."))
} else {
results <-
map(seq_along(gwas_coding_results$analysis_id), # .x
~{knitr::knit_child(text =
c('### `r gwas_coding_results$analysis_label[.x]` {.tabset .tabset-pills}',
'#### Regionplot',
'```{r}',
'rp_data <-
regionplot.data(analysis_id = gwas_coding_results$analysis_id[.x],
chrom = gwas_coding_results$chrom[.x],
locus = gwas_coding_results$locus[.x])
regionplot.plot(rp_data)',
'```',
'',
'#### Credible Set',
'```{r}',
'rp_data %>%
pluck("cred_set", 1) %>%
filter(cs == TRUE & pp >= 0.01) %>%
arrange(desc(pp)) %>%
DT::datatable(rownames = FALSE,
filter = "top",
options = list(pageLength = 15, scrollX = TRUE))',
'```',
''),
envir = environment(),
quiet = TRUE)
})
cat(unlist(results), sep = '\n')
}No coding variant to gene mapping for: IL22.
3.3 Variants
if(coding_eval_lgl == FALSE){
cat(glue("No coding variant(s) to follow-up for: {params$input_gene}."))
} else {
phewas_vars <-
gwas_coding_results %>%
distinct(chrom, pos, ref, alt) %>%
mutate(var_label = glue("chr{chrom}:{pos}-{ref}-{alt}"))
results <-
pmap(list(phewas_vars$chrom, # ..1
phewas_vars$pos, # ..2
phewas_vars$ref, # ..3
phewas_vars$alt, # ..4
phewas_vars$var_label), # ..5
~{
knitr::knit_child(text =
c('### `r ..5` {.tabset .tabset-pills}',
'#### Credible Phewas {.tabset .tabset-pills}',
'##### Disease & biomarkers',
'```{r}',
glue('credible_phewas(chrom = {..1}, pos = {..2}, ref = "{..3}", alt = "{..4}",
ignore_xqtl = TRUE)'),
'```',
'',
'##### All QTLs',
'```{r}',
glue('credible_phewas(chrom = {..1}, pos = {..2}, ref = "{..3}", alt = "{..4}",
tag = "qtls_all", ignore_xqtl = FALSE, pc_genes_only = TRUE)'),
'```',
'',
'#### Phewas {.tabset .tabset-pills}',
'##### Disease & biomarkers',
'```{r}',
glue('phewas(chrom = {..1}, pos = {..2}, ref = "{..3}", alt = "{..4}",
ignore_xqtl = TRUE)'),
'```',
'',
'##### All QTLs',
'```{r}',
glue('phewas(chrom = {..1}, pos = {..2}, ref = "{..3}", alt = "{..4}",
tag = "qtls_all", ignore_xqtl = FALSE, pc_genes_only = FALSE)'),
'```',
'',
'#### External links',
glue('- [Gnomad](https://gnomad.broadinstitute.org/variant/{..1}-{..2}-{..3}-{..4}?dataset=gnomad_r3)'),
glue('- [Finngen PheWeb](https://results.finngen.fi/variant/{..1}-{..2}-{..3}-{..4})'),
glue('- [Open Targets](https://genetics.opentargets.org/variant/{..1}_{..2}_{..3}_{..4})'),
glue('- [ClinVar](https://www.ncbi.nlm.nih.gov/clinvar/variation/3396/?oq={..1}[Chromosome]+AND+{..2}[Base+Position+for+Assembly+GRCh38])'),
''),
envir = environment(),
quiet = TRUE)
})
cat(unlist(results), sep = '\n')
}No coding variant(s) to follow-up for: IL22.
4 Coloc evidence
Evidence linking IL22 xQTLs with disease/biomarker phenotypes through colocalization. Colocalization provides the probability that a variant in the GWAS and xQTL is shared and therefore that the gene is driving disease. In addition, it provides the direction of effect - where associations can be either positive (higher disease risk ~ higher gene expression) or inversely related (higher disease risk ~ lower gene expression).
aba_coloc_query_results <-
safely_exec(instrument_gene_coloc.query,
gene_name = params$input_gene,
tag = c("gwas_latest", "rivas", "neale_cell_quant"),
pc_genes_only = TRUE)4.1 Phenome-wide colocs
These plots illustrate all the colocalization data in the region near IL22, including other genes, with colocalizations from all phenotypes and all tissues/cell types available. These plots only highlight positive results for illustrative purposes. All possible combinations of signals have been calculated and tested.
4.1.1 All tissue/cell type xQTLs
if(coloc_eval_lgl == TRUE & exists("aba_coloc_query_results")){
aba_coloc_query_results %>%
instrument_gene_coloc.wrangle(proximity = 'cis') %>%
instrument_gene_coloc.plot()
} else {
cat(glue("No coloc results in any tissue/cell type for: {params$input_gene}"))
}
4.1.2 Immune cells
if(exists("aba_coloc_imm_tag_res")){
if(!is_empty(aba_coloc_imm_tag_res)){
if(nrow(aba_coloc_imm_tag_res) > 0 ){
instrument_gene_coloc.plot(aba_coloc_imm_tag_res)
} else {
cat(glue("No coloc results for IMM related tissues/cell types."))
}
} else {
cat(glue("No coloc results for IMM related tissues/cell types."))
}
} else {
cat(glue("No coloc results for IMM related tissues/cell types."))
}## No coloc results for IMM related tissues/cell types.4.1.3 cis-pQTLs
if(exists("aba_coloc_pqtl_cis_tag_res")){
if(!is_empty(aba_coloc_pqtl_cis_tag_res)){
if(nrow(aba_coloc_pqtl_cis_tag_res) > 0){
instrument_gene_coloc.plot(aba_coloc_pqtl_cis_tag_res)
} else {
cat(glue("No colocs with cis-pQTLs phenotypes for: {params$input_gene}."))
}
} else {
cat(glue("No colocs with cis-pQTLs phenotypes for: {params$input_gene}."))
}
} else {
cat(glue("No colocs with cis-pQTLs phenotypes for: {params$input_gene}."))
}
4.1.4 GTEx
if(exists("aba_coloc_gtex_tag_res")){
if(!is_empty(aba_coloc_gtex_tag_res)){
if(nrow(aba_coloc_gtex_tag_res) > 0){
instrument_gene_coloc.plot(aba_coloc_gtex_tag_res)
} else {
cat(glue("No GTEx colocs for: {params$input_gene}."))
}
} else {
cat(glue("No GTEx colocs for: {params$input_gene}."))
}
} else {
cat(glue("No GTEx colocs for: {params$input_gene}."))
}## No GTEx colocs for: IL22.4.1.5 EEC
if(exists("aba_coloc_eec_tag_res")){
if(!is_empty(aba_coloc_eec_tag_res)){
if(nrow(aba_coloc_eec_tag_res) > 0 ){
instrument_gene_coloc.plot(aba_coloc_eec_tag_res)
} else {
cat(glue("No EEC colocs for: {params$input_gene}."))
}
} else {
cat(glue("No EEC colocs for: {params$input_gene}."))
}
} else {
cat(glue("No EEC colocs for: {params$input_gene}."))
}## No EEC colocs for: IL22.4.1.6 Biomarkers
if(exists("aba_coloc_rivas_tag_res")){
if(!is_empty(aba_coloc_rivas_tag_res)){
if(nrow(aba_coloc_rivas_tag_res) > 0 ){
instrument_gene_coloc.plot(aba_coloc_rivas_tag_res)
} else {
cat(glue("No colocs with biomarkers phenotypes for: {params$input_gene}."))
}
} else {
cat(glue("No colocs with biomarkers phenotypes for: {params$input_gene}."))
}
} else {
cat(glue("No colocs with biomarkers phenotypes for: {params$input_gene}."))
}## No colocs with biomarkers phenotypes for: IL22.4.1.7 trans-pQTLs
if(exists("aba_coloc_pqtl_trans_tag_res")){
if(!is_empty(aba_coloc_pqtl_trans_tag_res)){
if(nrow(aba_coloc_pqtl_trans_tag_res) > 0 ){
instrument_gene_coloc.plot(aba_coloc_pqtl_trans_tag_res)
} else {
cat(glue("No colocs with trans-pQTLs phenotypes for: {params$input_gene}."))
}
} else {
cat(glue("No colocs with trans-pQTLs phenotypes for: {params$input_gene}."))
}
} else {
cat(glue("No colocs with trans-pQTLs phenotypes for: {params$input_gene}."))
}## No colocs with trans-pQTLs phenotypes for: IL22.4.1.8 all-pQTLs
if(exists("aba_coloc_pqtl_all_tag_res")){
if(!is_empty(aba_coloc_pqtl_all_tag_res)){
if(nrow(aba_coloc_pqtl_all_tag_res) > 0 ){
instrument_gene_coloc.plot(aba_coloc_pqtl_all_tag_res)
} else {
cat(glue("No colocs with all-pQTLs phenotypes for: {params$input_gene}."))
}
} else {
cat(glue("No colocs with all-pQTLs phenotypes for: {params$input_gene}."))
}
} else {
cat(glue("No colocs with all-pQTLs phenotypes for: {params$input_gene}."))
}
4.2 Locus
These plots illustrate the colocalization results in more depth by analyzing just 1 GWAS locus with all genes and tissue/cell types provided. While all genes/cell types possible have been calculated and tested, a gene/tissue pair is only shown if at least 1 result is >= 50%. Select a phenotype (below) to investigate the locus.
coloc_loci_knitr_fxn <- function(gwas_aid, gwas_chrom, gwas_locus, gwas_label, tags_list){
tag4loop <-
tags_list %>%
left_join(tags_2_prioritize, by = "tags1") %>%
arrange(tags1_prio) %>%
pull(tags1)
ret <-
c(glue('### {{gwas_label}} {.tabset .tabset-pills}', .open = "{{", .close = "}}"),
'```{r}',
glue('mc_query_res <- coloc_multi.query(',
'analysis_id2 = "{{gwas_aid}}", ',
'chrom = "{{gwas_chrom}}", ',
'locus2 = "{{gwas_locus}}")',
.open = "{{",
.close = "}}"),
'```',
'',
'Select a tag below to determine which xQTLs to view.\n',
map(tag4loop, # .x
~{c(glue('#### {{.x}}', .open = "{{", .close = "}}"),
'```{r}',
glue('coloc_multi.wrangle(',
'.data = mc_query_res, ',
'tag = "{{.x}}", ',
'pc_genes_only = TRUE) %>% ',
'coloc_multi.plot()',
.open = "{{",
.close = "}}"
),
'```',
'')})
) %>%
unlist()
return(ret)
}if(coloc_eval_lgl == TRUE){
input_gene_gwas_aids_with_colocs <-
# Select coloc data for input gene only
aba_coloc_query_results %>%
instrument_gene_coloc.wrangle() %>%
filter(gene_name1 == params$input_gene & h4 >= 0.80) %>%
# Annotate xQTLs with their respective tags
mutate(tags1_list = map(analysis_id1, ~find_tags(analysis_id = .))) %>%
unnest(tags1_list) %>%
rename(tags1 = tag) %>%
semi_join(tags_2_prioritize, by = "tags1") %>%
# Select the list of GWAS and the xQTL tags
distinct(analysis_id2, chrom, locus2, analysis_label2,
data_type2, n_cases2, n_cohort2, tags1) %>%
distinct(analysis_id2, chrom, locus2, analysis_label2, tags1) %>%
mutate(analysis_id2 = as.character(analysis_id2)) %>%
group_by(analysis_id2, chrom, locus2) %>%
nest(tags1_list = tags1) %>%
ungroup() %>%
arrange(analysis_id2)
input_knitr_txt <-
pmap(list(input_gene_gwas_aids_with_colocs$analysis_id2, # ..1
input_gene_gwas_aids_with_colocs$chrom, # ..2
input_gene_gwas_aids_with_colocs$locus2, # ..3
input_gene_gwas_aids_with_colocs$analysis_label2, # ..4
input_gene_gwas_aids_with_colocs$tags1_list), # ..5
~coloc_loci_knitr_fxn(gwas_aid = ..1,
gwas_chrom = ..2,
gwas_locus = ..3,
gwas_label = ..4,
tags_list = ..5)) %>%
unlist()
results <- knitr::knit_child(text = input_knitr_txt, envir = environment(), quiet = TRUE)
cat(unlist(results), sep = '\n')
} else {
cat(glue("No coloc results for: {params$input_gene}."))
}4.2.1 FinnGen R10 - Atopic dermatitis
mc_query_res <- coloc_multi.query(analysis_id2 = "finngen_r10_ad", chrom = "12", locus2 = "2")Select a tag below to determine which xQTLs to view.
4.2.1.1 qtls_all
coloc_multi.wrangle(.data = mc_query_res, tag = "qtls_all", pc_genes_only = TRUE) %>% coloc_multi.plot()
4.2.1.2 qtls_pqtl
coloc_multi.wrangle(.data = mc_query_res, tag = "qtls_pqtl", pc_genes_only = TRUE) %>% coloc_multi.plot()
4.2.2 FinnGen R10 - Dermatitis & eczema
mc_query_res <- coloc_multi.query(analysis_id2 = "finngen_r10_eczema", chrom = "12", locus2 = "1")Select a tag below to determine which xQTLs to view.
4.2.2.1 qtls_all
coloc_multi.wrangle(.data = mc_query_res, tag = "qtls_all", pc_genes_only = TRUE) %>% coloc_multi.plot()
4.2.2.2 qtls_pqtl
coloc_multi.wrangle(.data = mc_query_res, tag = "qtls_pqtl", pc_genes_only = TRUE) %>% coloc_multi.plot()
4.2.3 FinnGen R10 Meta - Atopic dermatitis
mc_query_res <- coloc_multi.query(analysis_id2 = "finngen_r10_meta_ad", chrom = "12", locus2 = "2")Select a tag below to determine which xQTLs to view.
4.2.3.1 qtls_all
coloc_multi.wrangle(.data = mc_query_res, tag = "qtls_all", pc_genes_only = TRUE) %>% coloc_multi.plot()
4.2.3.2 qtls_pqtl
coloc_multi.wrangle(.data = mc_query_res, tag = "qtls_pqtl", pc_genes_only = TRUE) %>% coloc_multi.plot()
4.2.4 FinnGen R9 - Atopic dermatitis
mc_query_res <- coloc_multi.query(analysis_id2 = "finngen_r9_ad", chrom = "12", locus2 = "2")Select a tag below to determine which xQTLs to view.
4.2.4.1 qtls_all
coloc_multi.wrangle(.data = mc_query_res, tag = "qtls_all", pc_genes_only = TRUE) %>% coloc_multi.plot()
4.2.4.2 qtls_pqtl
coloc_multi.wrangle(.data = mc_query_res, tag = "qtls_pqtl", pc_genes_only = TRUE) %>% coloc_multi.plot()
4.2.5 Galileo ’22 meta - Atopic dermatitis
mc_query_res <- coloc_multi.query(analysis_id2 = "galileo_2022_meta_ad", chrom = "12", locus2 = "2")Select a tag below to determine which xQTLs to view.
4.2.5.1 qtls_all
coloc_multi.wrangle(.data = mc_query_res, tag = "qtls_all", pc_genes_only = TRUE) %>% coloc_multi.plot()
4.2.5.2 qtls_pqtl
coloc_multi.wrangle(.data = mc_query_res, tag = "qtls_pqtl", pc_genes_only = TRUE) %>% coloc_multi.plot()
4.2.6 Galileo ’23 meta - Atopic dermatitis
mc_query_res <- coloc_multi.query(analysis_id2 = "galileo_2023_meta_ad", chrom = "12", locus2 = "2")Select a tag below to determine which xQTLs to view.
4.2.6.1 qtls_all
coloc_multi.wrangle(.data = mc_query_res, tag = "qtls_all", pc_genes_only = TRUE) %>% coloc_multi.plot()
4.2.6.2 qtls_pqtl
coloc_multi.wrangle(.data = mc_query_res, tag = "qtls_pqtl", pc_genes_only = TRUE) %>% coloc_multi.plot()
4.2.7 Paternoster ’15 - Atopic dermatitis
mc_query_res <- coloc_multi.query(analysis_id2 = "paternoster2015_ad", chrom = "12", locus2 = "1")Select a tag below to determine which xQTLs to view.
4.2.7.1 qtls_all
coloc_multi.wrangle(.data = mc_query_res, tag = "qtls_all", pc_genes_only = TRUE) %>% coloc_multi.plot()
4.2.7.2 qtls_pqtl
coloc_multi.wrangle(.data = mc_query_res, tag = "qtls_pqtl", pc_genes_only = TRUE) %>% coloc_multi.plot()
4.3 Pairwise
These plots illustrate the colocalization of 1 GWAS locus with the tissue/cell type that had the highest colocalization probability for IL22. For the GWAS and the xQTL the regionplots illustrate which variants are in each respective credible set with red color. The green circles around the variants highlight credible set variants that are shared betwen the two sets. These plots allow us to visually confirm the colocalization results.
coloc_pairwise_knitr_fxn <- function(chrom, analysis_id1, analysis_id2, locus1, locus2, entity1){
gwas_label <- describe_analysis_ids(analysis_id = analysis_id2) %>% pull(analysis_label)
ret <-
c(glue('### {gwas_label}'),
'```{r}',
glue('coloc_pair(analysis_id1 = "{analysis_id1}", analysis_id2 = "{analysis_id2}", ',
'chrom = {chrom}, entity1 = "{entity1}", ',
'locus1 = {locus1}, locus2 = {locus2})'),
'```',
'')
return(ret)
}if(coloc_eval_lgl == TRUE){
pairwise_plot_input <-
aba_coloc_query_results %>%
instrument_gene_coloc.wrangle() %>%
# Get our gene of interest
filter(gene_name1 == params$input_gene & h4 >= 0.80) %>%
# Find highest H4 tissue for each GWAS w/ coloc
group_by(analysis_id2) %>%
mutate(gwas_max_h4 = max(h4)) %>%
ungroup() %>%
filter(h4 == gwas_max_h4) %>%
# Select just 1 tissue in case of exact ties
group_by(analysis_id2) %>%
mutate(row_n = ntile(n = n())) %>%
ungroup() %>%
filter(row_n == 1) %>%
# Arrange by GWAS sample size
select(chrom, analysis_id2, locus2, analysis_id1, locus1, entity1) %>%
mutate(analysis_id2 = as.character(analysis_id2)) %>%
arrange(analysis_id2)
input_knitr_txt <-
pmap(dots_list(
pairwise_plot_input$analysis_id1, # ..1
pairwise_plot_input$analysis_id2, # ..2
pairwise_plot_input$chrom, # ..3
pairwise_plot_input$locus1, # ..4
pairwise_plot_input$locus2, # ..5
pairwise_plot_input$entity1, # ..6
.named = TRUE),
~coloc_pairwise_knitr_fxn(
analysis_id1 = ..1,
analysis_id2 = ..2,
chrom = ..3,
locus1 = ..4,
locus2 = ..5,
entity1 = ..6)) %>%
unlist()
results <- knitr::knit_child(text = input_knitr_txt, envir = environment(), quiet = TRUE) #
cat(unlist(results), sep = '\n')
} else {
cat(glue("No coloc results for: {params$input_gene}."))
}4.3.1 FinnGen R10 - Atopic dermatitis
coloc_pair(analysis_id1 = "ppp_3072_pqtl_plasma", analysis_id2 = "finngen_r10_ad", chrom = 12, entity1 = "OID31407", locus1 = 1, locus2 = 2)
4.3.2 FinnGen R10 - Dermatitis & eczema
coloc_pair(analysis_id1 = "ppp_3072_pqtl_plasma", analysis_id2 = "finngen_r10_eczema", chrom = 12, entity1 = "OID31407", locus1 = 1, locus2 = 1)
4.3.3 FinnGen R10 Meta - Atopic dermatitis
coloc_pair(analysis_id1 = "ppp_3072_pqtl_plasma", analysis_id2 = "finngen_r10_meta_ad", chrom = 12, entity1 = "OID31407", locus1 = 1, locus2 = 2)
4.3.4 FinnGen R9 - Atopic dermatitis
coloc_pair(analysis_id1 = "ppp_3072_pqtl_plasma", analysis_id2 = "finngen_r9_ad", chrom = 12, entity1 = "OID31407", locus1 = 1, locus2 = 2)
4.3.5 Galileo ’22 meta - Atopic dermatitis
coloc_pair(analysis_id1 = "ppp_3072_pqtl_plasma", analysis_id2 = "galileo_2022_meta_ad", chrom = 12, entity1 = "OID31407", locus1 = 1, locus2 = 2)
4.3.6 Galileo ’23 meta - Atopic dermatitis
coloc_pair(analysis_id1 = "ppp_3072_pqtl_plasma", analysis_id2 = "galileo_2023_meta_ad", chrom = 12, entity1 = "OID31407", locus1 = 1, locus2 = 2)
4.3.7 Paternoster ’15 - Atopic dermatitis
coloc_pair(analysis_id1 = "ppp_3072_pqtl_plasma", analysis_id2 = "paternoster2015_ad", chrom = 12, entity1 = "OID31407", locus1 = 1, locus2 = 1)
4.4 Genome-wide pQTL analysis
These plots illustrate direction of effect between the pQTL & GWAS data. Each dot illustrates the beta for the lead cred set variant of the colocalized GWAS & pQTL. The bars illustrate the summary stat standard error around the beta. Variants with a blue box surround them are in cis to the target gene. GWAS betas have been adjusted so that they are relative to the pQTL beta being positive. Only GWAS with >= 1 cis coloc are plotted.
coloc_gw_pqtl_knitr_fxn <- function(analysis_id1, analysis_id2, gene_name, analysis_label2){
ret <-
c(glue('### {analysis_label2}'),
'```{r}',
glue('coloc_gw_pqtl(analysis_id1 = "{analysis_id1}", analysis_id2 = "{analysis_id2}", ',
'gene_name = "{params$input_gene}")'),
'```',
'')
return(ret)
}if(coloc_eval_lgl == TRUE & exists("aba_coloc_query_results")){
pqtl_colocs <-
aba_coloc_query_results %>%
filter(data_type1 == 'pqtl' & h4 >= 0.80 & gene_name1 == params$input_gene)
if(nrow(pqtl_colocs) == 0){
cat(glue("No pQTL coloc results for: {params$input_gene}."))
} else {
input_knitr_txt <-
pmap(dots_list(
pqtl_colocs$analysis_id1, # ..1
pqtl_colocs$analysis_id2, # ..2
pqtl_colocs$analysis_label2, # ..3
.named = TRUE),
~coloc_gw_pqtl_knitr_fxn(
analysis_id1 = ..1,
analysis_id2 = ..2,
gene_name = params$input_gene,
analysis_label2 = ..3)) %>%
unlist()
results <- knitr::knit_child(text = input_knitr_txt, envir = environment(), quiet = TRUE) #
cat(unlist(results), sep = '\n')
}
} else {
cat(glue("No coloc results for: {params$input_gene}."))
}4.4.1 FinnGen R10 Meta - Atopic dermatitis
coloc_gw_pqtl(analysis_id1 = "ppp_3072_pqtl_plasma", analysis_id2 = "finngen_r10_meta_ad", gene_name = "IL22")
4.4.2 FinnGen R10 - Dermatitis & eczema
coloc_gw_pqtl(analysis_id1 = "ppp_3072_pqtl_plasma", analysis_id2 = "finngen_r10_eczema", gene_name = "IL22")
4.4.3 FinnGen R10 - Atopic dermatitis
coloc_gw_pqtl(analysis_id1 = "ppp_3072_pqtl_plasma", analysis_id2 = "finngen_r10_ad", gene_name = "IL22")
4.4.4 FinnGen R9 - Atopic dermatitis
coloc_gw_pqtl(analysis_id1 = "ppp_3072_pqtl_plasma", analysis_id2 = "finngen_r9_ad", gene_name = "IL22")
4.4.5 Paternoster ’15 - Atopic dermatitis
coloc_gw_pqtl(analysis_id1 = "ppp_3072_pqtl_plasma", analysis_id2 = "paternoster2015_ad", gene_name = "IL22")
4.4.6 Galileo ’22 meta - Atopic dermatitis
coloc_gw_pqtl(analysis_id1 = "ppp_3072_pqtl_plasma", analysis_id2 = "galileo_2022_meta_ad", gene_name = "IL22")
4.4.7 Galileo ’23 meta - Atopic dermatitis
coloc_gw_pqtl(analysis_id1 = "ppp_3072_pqtl_plasma", analysis_id2 = "galileo_2023_meta_ad", gene_name = "IL22")
5 Nearest gene
These data illustrate the nearest gene variant to gene mapping results.
5.1 Regionplots
# TODO: Update for signals
# TODO: Add to callisto?
nearest_plot_fxn <- function(gwas_aid, gwas_chrom, gwas_locus){
# Query other nearby genes
gene_data <-
sw_dbGetQuery(
conn = dbc,
statement = glue("SELECT
index_analysis_id, index_chrom, index_locus, index_pos, index_mlog_p,
near_gene_name, near_gene_tss, near_gene_pos_start, near_gene_pos_end,
near_tss_dist, near_gene_dist, near_tss_rank, near_gene_rank
FROM nearest_pc_gene
WHERE
index_entity IS NULL AND
index_analysis_id = '{gwas_aid}' AND
index_chrom = '{gwas_chrom}' AND
index_locus = '{gwas_locus}' AND
(near_tss_rank <= 3 OR near_gene_rank <= 3)")) %>%
mutate(near_gene_tes = if_else(near_gene_tss == near_gene_pos_start,
near_gene_pos_end,
near_gene_pos_start),
near_type = case_when(near_tss_rank == near_gene_rank ~ "tss+gene",
near_tss_rank <= near_gene_rank ~ "tss",
TRUE ~ "gene"),
near_best_rank = if_else(str_detect(near_type, "tss"), near_tss_rank, near_gene_rank),
nearest_dist = if_else(near_tss_dist < near_gene_dist, near_tss_dist, near_gene_dist),
nearest_dist_rank = ntile(nearest_dist, n()),
plot_x = case_when(near_tss_rank == near_gene_rank ~ index_pos,
near_tss_rank <= near_gene_rank ~ near_gene_tss,
TRUE ~ near_gene_tes),
plot_y =
max(index_mlog_p) -
(((max(index_mlog_p)-pval2mlog(5e-8))/n_distinct(nearest_dist_rank))*(nearest_dist_rank-1)),
label_txt = glue("{near_gene_name}\nRank: {near_best_rank}\nMethod: {near_type}\nDist: {nearest_dist}"))
rp_fig <-
regionplot.data(analysis_id = gwas_aid, chrom = gwas_chrom, locus = gwas_locus) %>%
regionplot.plot(., append_xaxis_genes_labels = FALSE)
near_fig <-
rp_fig +
# Index variant vertical line
geom_segment(data = gene_data %>% distinct(index_pos, index_mlog_p),
mapping = aes(x = index_pos, xend = index_pos,
y = 0, yend = index_mlog_p),
alpha = 0.8,
color = "red",
linetype = 2) +
# Horizontal top lines for genes
geom_segment(data = gene_data,
mapping = aes(x = index_pos, xend = plot_x,
y = plot_y, yend = plot_y,
color = as.character(nearest_dist_rank)),
alpha = 0.8,
show.legend = FALSE) +
# Vertical lines for genes
geom_segment(data = gene_data,
mapping = aes(x = plot_x, xend = plot_x,
y = 0, yend = plot_y,
color = as.character(nearest_dist_rank)),
alpha = 0.8,
show.legend = FALSE) +
# Adjust colors manually to color blind friendly colors.
## Ranks will be 1-4 will be normal, maybe 5 rare, 6 should be *very* rare.
scale_discrete_manual(
aesthetics = c("color", "segment.color"),
values = c(
"1" = rgb(086, 180, 233, maxColorValue = 255), # sky blue
"2" = rgb(230, 159, 000, maxColorValue = 255), # orange (yellow orange)
"3" = rgb(204, 121, 167, maxColorValue = 255), # "reddish-purple" (pink)
"4" = rgb(000, 158, 115, maxColorValue = 255), # blueish green
"5" = rgb(213, 094, 000, maxColorValue = 255), # Vermillion (reddish orange)
"6" = rgb(000, 114, 178, maxColorValue = 255) # blue
))
# Upstream labels
if(gene_data %>% filter(plot_x <= index_pos) %>% nrow() > 0){
near_fig <-
near_fig +
geom_label_repel(data = gene_data %>% filter(plot_x <= index_pos),
mapping = aes(x = plot_x, y = plot_y,
label = label_txt,
segment.color = as.character(nearest_dist_rank)),
box.padding = 0.75,
min.segment.length = 0,
xlim = c(NA,
gene_data %>%
filter(plot_x <= index_pos) %>%
summarize(max_x = max(plot_x, na.rm = TRUE)) %>%
pull()),
show.legend = FALSE
)
}
# Downstream labels
if(gene_data %>% filter(plot_x > index_pos) %>% nrow() > 0){
near_fig <-
near_fig +
geom_label_repel(data = gene_data %>% filter(plot_x > index_pos),
mapping = aes(x = plot_x, y = plot_y,
label = label_txt,
segment.color = as.character(nearest_dist_rank)),
box.padding = 0.75,
min.segment.length = 0,
xlim = c(gene_data %>%
filter(plot_x > index_pos) %>%
summarize(min_x = min(plot_x, na.rm = TRUE)) %>%
pull(),
NA),
show.legend = FALSE
)
}
# Append genes below x-axis
near_fig <- callisto::append_locus_gene_plot(.data = near_fig, style = "contrast", pc_genes_only = TRUE)
# Return figure
return(near_fig)
}if(near_eval_lgl == TRUE){
nearest_data <- instrument_gene_nearest(gene_name = params$input_gene, ignore_xqtl = TRUE)
# Sort GWAS by sample size
nearest_data <-
nearest_data %>%
left_join(describe_analysis_ids(analysis_id = nearest_data %>%
distinct(index_analysis_id) %>%
pull(index_analysis_id)),
by = c("index_analysis_id" = "analysis_id")) %>%
arrange(desc(data_type), desc(n_cases), desc(n_cohort)) %>%
select(near_gene_name:n_genes_in_locus, analysis_label)
nearest_gwas_ids <-
nearest_data %>%
distinct(index_analysis_id, index_chrom, index_locus, analysis_label) %>%
arrange(index_analysis_id)
results <-
pmap(list(nearest_gwas_ids$index_analysis_id, # ..1
nearest_gwas_ids$index_chrom, # ..2
nearest_gwas_ids$index_locus, # ..3
nearest_gwas_ids$analysis_label), # ..4
~{knitr::knit_child(text =
c(glue('### {..4}'),
'```{r}',
glue('nearest_plot_fxn(gwas_aid = "{..1}",
gwas_chrom = "{..2}",
gwas_locus = {..3})'),
'```',
''
),
envir = environment(),
quiet = TRUE)
})
cat(unlist(results), sep = '\n')
} else {
cat(glue("No nearest gene evidence for: {params$input_gene}."))
}5.1.1 FinnGen R10 - Atopic dermatitis
nearest_plot_fxn(gwas_aid = "finngen_r10_ad",
gwas_chrom = "12",
gwas_locus = 2)
5.1.2 FinnGen R10 - Dermatitis & eczema
nearest_plot_fxn(gwas_aid = "finngen_r10_eczema",
gwas_chrom = "12",
gwas_locus = 1)
5.1.3 FinnGen R10 Meta - Atopic dermatitis
nearest_plot_fxn(gwas_aid = "finngen_r10_meta_ad",
gwas_chrom = "12",
gwas_locus = 2)
5.1.4 FinnGen R9 - Atopic dermatitis
nearest_plot_fxn(gwas_aid = "finngen_r9_ad",
gwas_chrom = "12",
gwas_locus = 2)
5.1.5 Galileo ’22 meta - Atopic dermatitis
nearest_plot_fxn(gwas_aid = "galileo_2022_meta_ad",
gwas_chrom = "12",
gwas_locus = 2)
5.1.6 Galileo ’23 meta - Atopic dermatitis
nearest_plot_fxn(gwas_aid = "galileo_2023_meta_ad",
gwas_chrom = "12",
gwas_locus = 2)
5.1.7 Neale ’18 - White blood cell
nearest_plot_fxn(gwas_aid = "neale2018_wbc_cnt",
gwas_chrom = "12",
gwas_locus = 13)
5.1.8 Paternoster ’15 - Atopic dermatitis
nearest_plot_fxn(gwas_aid = "paternoster2015_ad",
gwas_chrom = "12",
gwas_locus = 1)
5.2 Tabular results
if(near_eval_lgl == TRUE){
# Create pretty table for html
nearest_data %>%
DT::datatable(rownames = FALSE, filter = "top", options = list(pageLength = 5, scrollX = TRUE))
} else {
cat(glue("No nearest gene evidence for: {params$input_gene}."))
}6 External links
These are some useful links to external data to further investigate other genetics information about IL22.
input_ensg <-
tibble(entities = gene2entity(gene_name = params$input_gene)) %>%
filter(str_detect(entities, "ENSG\\d+$")) %>%
pull(entities)
ot_url <-
if_else(is_empty(input_ensg),
"https://genetics.opentargets.org/",
"https://genetics.opentargets.org/gene/`r input_ensg`")
genebasss_url <-
if_else(is_empty(input_ensg),
"https://app.genebass.org/",
"https://app.genebass.org/gene/`r input_ensg`?burdenSet=pLoF&phewasOpts=1&resultLayout=full")- Open Targets
- OMIM
- FinnGen
- Genebass tool & Publication
- Reminder, gene based thresholds:
- SKAT-O = 2.5x10^-8
- Burden = 6.7x10^-7
- Single var = 8.0x10^-9
- Reminder, gene based thresholds: