Skip to content

Commit

Permalink
add make_convergence_plots() to make convergence plots and simplify s…
Browse files Browse the repository at this point in the history 
…ummarize_param()
  • Loading branch information
@kevinlkx
kevinlkx committed Jun 17, 2024
1 parent 2833698 commit f5a7a46
Show file tree
Hide file tree
Showing 8 changed files with 163 additions and 149 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ Type: Package
Title: Adjusting for Genetic Confounders in Transcriptome-Wide Association
Studies Improves Discovery of Risk Genes of Complex Traits
Date: 2024-06-14
Version: 0.2.40
Version: 0.2.41
Authors@R: c(person("Siming","Zhao",role="aut",email="siming.zhao06@gmail.com"),
person("Wesley","Crouse",role="aut"),
person("Sheng","Qian",role="aut",email="shengqian@uchicago.edu"),
Expand Down
2 changes: 1 addition & 1 deletion NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ export(get_gene_annot_from_ens_db)
export(get_problematic_genes)
export(get_region_cor)
export(load_weights)
export(make_convergence_plots)
export(make_locusplot)
export(map_snp_info_regions)
export(merge_finemap_regions)
Expand Down Expand Up @@ -77,7 +78,6 @@ importFrom(ggplot2,labs)
importFrom(ggplot2,margin)
importFrom(ggplot2,scale_alpha_manual)
importFrom(ggplot2,scale_color_manual)
importFrom(ggplot2,scale_colour_discrete)
importFrom(ggplot2,scale_fill_manual)
importFrom(ggplot2,scale_shape_manual)
importFrom(ggplot2,scale_size_manual)
Expand Down
138 changes: 137 additions & 1 deletion R/ctwas_plots.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -238,7 +238,7 @@ make_locusplot <- function(finemap_res,

# PIP panel
p_pip <- ggplot(loc$data, aes(x=.data$pos/1e6, y=.data$susie_pip, shape=.data$type,
size=.data$object_type, alpha=.data$object_type)) +
size=.data$object_type, alpha=.data$object_type)) +
geom_point(aes(color=.data$r2_levels, fill=.data$r2_levels)) +
geom_text_repel(aes(label=.data$label), size=label.text.size, color="black") +
scale_shape_manual(values = point.shapes) +
Expand Down Expand Up @@ -317,3 +317,139 @@ make_locusplot <- function(finemap_res,
cowplot::plot_grid(p_pvalue, p_pip, p_qtl, p_genes, ncol = 1,
rel_heights = panel.heights, align = "v", axis="tblr")
}

#' @title Make convergence plots for the estimated parameters
#'
#' @param param a list of cTWAS parameter estimation result from \code{est_param}
#'
#' @param gwas_n the sample size of the GWAS summary statistics
#'
#' @param title.size font size of the plot title
#'
#' @param legend.size font size of the plot legend title
#'
#' @param colors colors for different groups, passed as the
#' \code{values} input to \code{\link[ggplot2]{scale_color_manual}}.
#' If fewer colors than "fits" are given, the colors are recycled.
#'
#' @param ncol number of columns in the output plot
#'
#' @importFrom ggplot2 ggplot
#' @importFrom ggplot2 aes
#' @importFrom ggplot2 geom_point
#' @importFrom ggplot2 geom_line
#' @importFrom ggplot2 xlab ylab ggtitle
#' @importFrom ggplot2 expand_limits
#' @importFrom ggplot2 theme
#' @importFrom ggplot2 theme_bw
#' @importFrom ggplot2 element_text
#' @importFrom ggplot2 guides
#' @importFrom ggplot2 guide_legend
#' @importFrom ggplot2 scale_color_manual
#' @importFrom cowplot plot_grid theme_cowplot
#' @importFrom rlang .data
#'
#' @export
make_convergence_plots <- function(param,
gwas_n,
title.size = 10,
legend.size = 8,
colors = c("#E69F00","#56B4E9","#009E73","#F0E442",
"#0072B2","#D55E00","#CC79A7", "#999999"),
ncol = 2){

# estimated group prior (all iterations)
group_prior_iters <- param$group_prior_iters
# estimated group prior variance (all iterations)
group_prior_var_iters <- param$group_prior_var_iters

# group size
group_size <- param$group_size[rownames(group_prior_iters)]

# estimated group PVE (all iterations)
group_pve_iters <- group_prior_var_iters*group_prior_iters*group_size/gwas_n
# estimated enrichment of genes (all iterations)
enrichment_iters <- t(sapply(rownames(group_prior_iters)[rownames(group_prior_iters)!="SNP"], function(x){
group_prior_iters[rownames(group_prior_iters)==x,]/group_prior_iters[rownames(group_prior_iters)=="SNP"]}))

# make convergence plots

# inclusion plot
df <- data.frame(niter = rep(1:ncol(group_prior_iters), nrow(group_prior_iters)),
value = unlist(lapply(1:nrow(group_prior_iters), function(x){group_prior_iters[x,]})),
group = rep(rownames(group_prior_iters), each=ncol(group_prior_iters)))
factor_levels <- c(setdiff(rownames(group_prior_iters), "SNP"), "SNP")
df$group <- factor(df$group, levels = factor_levels)

p_pi <- ggplot(df, aes(x=.data$niter, y=.data$value, group=.data$group, color=.data$group)) +
geom_line() +
geom_point() +
scale_color_manual(values = colors) +
xlab("Iteration") + ylab(bquote(pi)) +
ggtitle("Proportion Causal") +
theme_cowplot() +
theme(plot.title=element_text(size=title.size)) +
expand_limits(y=0) +
guides(color = guide_legend(title = "Group")) +
theme(legend.title = element_text(size=legend.size, face="bold"),
legend.text = element_text(size=legend.size))

# effect size plot
df <- data.frame(niter = rep(1:ncol(group_prior_var_iters), nrow(group_prior_var_iters)),
value = unlist(lapply(1:nrow(group_prior_var_iters), function(x){group_prior_var_iters[x,]})),
group = rep(rownames(group_prior_var_iters), each=ncol(group_prior_var_iters)))
df$group <- factor(df$group, levels = factor_levels)

p_sigma2 <- ggplot(df, aes(x=.data$niter, y=.data$value, group=.data$group, color=.data$group)) +
geom_line() +
geom_point() +
scale_color_manual(values = colors) +
xlab("Iteration") + ylab(bquote(sigma^2)) +
ggtitle("Effect Size") +
theme_cowplot() +
theme(plot.title=element_text(size=title.size)) +
expand_limits(y=0) +
guides(color = guide_legend(title = "Group")) +
theme(legend.title = element_text(size=legend.size, face="bold"),
legend.text = element_text(size=legend.size))

# PVE plot
df <- data.frame(niter = rep(1:ncol(group_pve_iters), nrow(group_pve_iters)),
value = unlist(lapply(1:nrow(group_pve_iters), function(x){group_pve_iters[x,]})),
group = rep(rownames(group_pve_iters), each=ncol(group_pve_iters)))
df$group <- factor(df$group, levels = factor_levels)

p_pve <- ggplot(df, aes(x=.data$niter, y=.data$value, group=.data$group, color=.data$group)) +
geom_line() +
geom_point() +
scale_color_manual(values = colors) +
xlab("Iteration") + ylab(bquote(h[G]^2)) +
ggtitle("PVE") +
theme_cowplot() +
theme(plot.title=element_text(size=title.size)) +
expand_limits(y=0) +
guides(color = guide_legend(title = "Group")) +
theme(legend.title = element_text(size=legend.size, face="bold"),
legend.text = element_text(size=legend.size))

# enrichment plot
df <- data.frame(niter = rep(1:ncol(enrichment_iters), nrow(enrichment_iters)),
value = unlist(lapply(1:nrow(enrichment_iters), function(x){enrichment_iters[x,]})),
group = rep(rownames(enrichment_iters), each=ncol(enrichment_iters)))
df$group <- factor(df$group, levels = factor_levels[factor_levels!="SNP"])

p_enrich <- ggplot(df, aes(x=.data$niter, y=.data$value, group=.data$group, color=.data$group)) +
geom_line() +
geom_point() +
scale_color_manual(values = colors) +
xlab("Iteration") + ylab(bquote(pi[G]/pi[S])) +
ggtitle("Enrichment") +
theme_cowplot() +
theme(plot.title=element_text(size=title.size)) +
expand_limits(y=0) +
guides(color = guide_legend(title = "Group")) +
theme(legend.title = element_text(size=legend.size, face="bold"),
legend.text = element_text(size=legend.size))

cowplot::plot_grid(p_pi, p_sigma2, p_enrich, p_pve, ncol = ncol)
}
132 changes: 12 additions & 120 deletions R/ctwas_summarize_parameters.R
View file
Original file line number Diff line number Diff line change
@@ -1,43 +1,14 @@
#' @title Summarizes and plots cTWAS parameter estimates
#' @title Summarizes estimated parameters
#'
#' @param param a list of cTWAS parameter estimation result from \code{est_param}
#' @param param a list of parameter estimation result from \code{est_param}
#'
#' @param gwas_n the sample size of the GWAS summary statistics
#'
#' @param plot if TRUE, return a plot of the estimated parameters
#'
#' @param title.size font size of the plot title
#'
#' @param legend.size font size of the plot legend title
#'
#' @param title.font font of the plot title
#'
#' @importFrom logging loginfo
#' @importFrom ggplot2 ggplot
#' @importFrom ggplot2 aes
#' @importFrom ggplot2 geom_point
#' @importFrom ggplot2 geom_line
#' @importFrom ggplot2 xlab ylab ggtitle
#' @importFrom ggplot2 expand_limits
#' @importFrom ggplot2 theme
#' @importFrom ggplot2 theme_bw
#' @importFrom ggplot2 element_text
#' @importFrom ggplot2 guides
#' @importFrom ggplot2 guide_legend
#' @importFrom ggplot2 scale_colour_discrete
#' @importFrom cowplot plot_grid theme_cowplot
#' @importFrom rlang .data
#' @return a list of summarized parameters
#'
#' @export
#'
summarize_param <- function(param,
gwas_n,
plot = TRUE,
title.size = 10,
legend.size = 8,
title.font = c("bold", "plain", "italic", "bold.italic")){

title.font <- match.arg(title.font)
summarize_param <- function(param, gwas_n){

group_prior <- param$group_prior
# estimated group prior (all iterations)
Expand All @@ -60,92 +31,13 @@ summarize_param <- function(param,
group_prior_iters[rownames(group_prior_iters)==x,]/group_prior_iters[rownames(group_prior_iters)=="SNP"]}))
enrichment <- enrichment_iters[,ncol(enrichment_iters)]

outlist <- list(group_size = group_size,
group_prior = group_prior,
group_prior_var = group_prior_var,
enrichment = enrichment,
group_pve = group_pve,
total_pve = sum(group_pve),
attributable_pve = group_pve/sum(group_pve))

if (plot){

# inclusion plot
df <- data.frame(niter = rep(1:ncol(group_prior_iters), nrow(group_prior_iters)),
value = unlist(lapply(1:nrow(group_prior_iters), function(x){group_prior_iters[x,]})),
group = rep(rownames(group_prior_iters), each=ncol(group_prior_iters)))
factor_levels <- c(setdiff(rownames(group_prior_iters), "SNP"), "SNP")
df$group <- factor(df$group, levels = factor_levels)

p_pi <- ggplot(df, aes(x=.data$niter, y=.data$value, group=.data$group, color=.data$group)) +
geom_line() +
geom_point() +
xlab("Iteration") + ylab(bquote(pi)) +
ggtitle("Proportion Causal") +
theme_cowplot() +
theme(plot.title=element_text(size=title.size)) +
expand_limits(y=0) +
guides(color = guide_legend(title = "Group")) +
theme(legend.title = element_text(size=legend.size, face=title.font),
legend.text = element_text(size=legend.size))

# effect size plot
df <- data.frame(niter = rep(1:ncol(group_prior_var_iters), nrow(group_prior_var_iters)),
value = unlist(lapply(1:nrow(group_prior_var_iters), function(x){group_prior_var_iters[x,]})),
group = rep(rownames(group_prior_var_iters), each=ncol(group_prior_var_iters)))
df$group <- factor(df$group, levels = factor_levels)

p_sigma2 <- ggplot(df, aes(x=.data$niter, y=.data$value, group=.data$group, color=.data$group)) +
geom_line() +
geom_point() +
xlab("Iteration") + ylab(bquote(sigma^2)) +
ggtitle("Effect Size") +
theme_cowplot() +
theme(plot.title=element_text(size=title.size)) +
expand_limits(y=0) +
guides(color = guide_legend(title = "Group")) +
theme(legend.title = element_text(size=legend.size, face=title.font),
legend.text = element_text(size=legend.size))

# PVE plot
df <- data.frame(niter = rep(1:ncol(group_pve_iters), nrow(group_pve_iters)),
value = unlist(lapply(1:nrow(group_pve_iters), function(x){group_pve_iters[x,]})),
group = rep(rownames(group_pve_iters), each=ncol(group_pve_iters)))
df$group <- factor(df$group, levels = factor_levels)

p_pve <- ggplot(df, aes(x=.data$niter, y=.data$value, group=.data$group, color=.data$group)) +
geom_line() +
geom_point() +
xlab("Iteration") + ylab(bquote(h[G]^2)) +
ggtitle("PVE") +
theme_cowplot() +
theme(plot.title=element_text(size=title.size)) +
expand_limits(y=0) +
guides(color = guide_legend(title = "Group")) +
theme(legend.title = element_text(size=legend.size, face=title.font),
legend.text = element_text(size=legend.size))

# enrichment plot
df <- data.frame(niter = rep(1:ncol(enrichment_iters), nrow(enrichment_iters)),
value = unlist(lapply(1:nrow(enrichment_iters), function(x){enrichment_iters[x,]})),
group = rep(rownames(enrichment_iters), each=ncol(enrichment_iters)))
df$group <- factor(df$group, levels = factor_levels[factor_levels!="SNP"])

p_enrich <- ggplot(df, aes(x=.data$niter, y=.data$value, group=.data$group, color=.data$group)) +
geom_line() +
geom_point() +
xlab("Iteration") + ylab(bquote(pi[G]/pi[S])) +
ggtitle("Enrichment") +
theme_cowplot() +
theme(plot.title=element_text(size=title.size)) +
expand_limits(y=0) +
guides(color = guide_legend(title = "Group")) +
theme(legend.title = element_text(size=legend.size, face=title.font),
legend.text = element_text(size=legend.size)) +
scale_colour_discrete(drop = FALSE)

outlist$convergence_plot <- plot_grid(p_pi, p_sigma2, p_enrich, p_pve)
}
res <- list(group_size = group_size,
group_prior = group_prior,
group_prior_var = group_prior_var,
enrichment = enrichment,
group_pve = group_pve,
total_pve = sum(group_pve),
attributable_pve = group_pve/sum(group_pve))

return(outlist)
return(res)
}
26 changes: 7 additions & 19 deletions man/summarize_param.Rd
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 1 addition & 1 deletion tests/testthat/test-ctwas_summarize_parameters.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ test_that("summarize_param works", {
param <- ctwas_res$param
precomputed_ctwas_parameters <- readRDS(system.file("extdata/sample_data", "LDL_example.ctwas_parameters_no_plots.RDS", package = "ctwas"))

ctwas_parameters <- summarize_param(param, gwas_n, plot = FALSE)
ctwas_parameters <- summarize_param(param, gwas_n)
expect_equal(ctwas_parameters, precomputed_ctwas_parameters)

})
5 changes: 2 additions & 3 deletions vignettes/simple_ctwas_tutorial.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -245,10 +245,9 @@ ctwas_parameters$total_pve
ctwas_parameters$attributable_pve
```


Plot convergence
Make convergence plots of estimated parameters
```{r convergence_plot, fig.width = 7, fig.height=5}
ctwas_parameters$convergence_plot
make_convergence_plots(param, gwas_n)
```


Expand Down
5 changes: 2 additions & 3 deletions vignettes/summarizing_ctwas_results.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,13 +92,12 @@ ctwas_parameters$attributable_pve
```


Plot convergence
Make convergence plots of estimated parameters
```{r convergence_plot, fig.width = 7, fig.height=5}
ctwas_parameters$convergence_plot
make_convergence_plots(param, gwas_n)
```



## Summarizing cTWAS results


Expand Down

0 comments on commit f5a7a46

Please sign in to comment.