alicemacqueen/snpdiver
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Edits to mash evaluation functions to improve plotting

Browse Source
This commit is contained in:
Alice MacQueen
2021-04-07 15:53:59 -05:00
parent 04fe4f1281
commit 284dd5ef14
3 changed files with 209 additions and 15 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
NAMESPACE
View File

@@ -79,6 +79,7 @@ importFrom(stats,ppoints)
importFrom(stats,predict) importFrom(stats,predict)
importFrom(stats,qbeta) importFrom(stats,qbeta)
importFrom(stats,uniroot) importFrom(stats,uniroot)
importFrom(stringr,str_length)
importFrom(stringr,str_replace) importFrom(stringr,str_replace)
importFrom(stringr,str_replace_all) importFrom(stringr,str_replace_all)
importFrom(stringr,str_sub) importFrom(stringr,str_sub)

152
R/mash-evaluation.R
View File

@@ -539,7 +539,7 @@ reorder_cormat <- function(cormat){
#' @importFrom tidyr separate #' @importFrom tidyr separate
#' @import ggplot2 #' @import ggplot2
#' @importFrom purrr as_vector #' @importFrom purrr as_vector
#' @importFrom stringr str_replace str_replace_all str_sub #' @importFrom stringr str_replace str_replace_all str_sub str_length
#' #'
#' @export #' @export
mash_plot_marker_effect <- function(m, snp = snp, n = NA, i = NA, marker = TRUE, mash_plot_marker_effect <- function(m, snp = snp, n = NA, i = NA, marker = TRUE,
@@ -576,7 +576,28 @@ mash_plot_marker_effect <- function(m, snp = snp, n = NA, i = NA, marker = TRUE,
x = .data$Conditions), width = 0.3) + x = .data$Conditions), width = 0.3) +
geom_hline(yintercept = 0, lty = 2) + geom_hline(yintercept = 0, lty = 2) +
labs(x = "Conditions", y = "Effect Size") + labs(x = "Conditions", y = "Effect Size") +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) theme_classic() +
theme(axis.text.x = element_text(size = 10, angle = 45, vjust = 1, hjust = 1),
panel.background = element_rect(fill=NA),
legend.position = "none",
axis.title = element_text(size = 10),
axis.text.y = element_text(size = 10),
axis.line.x = element_line(size = 0.35, colour = 'grey50'),
axis.line.y = element_line(size = 0.35, colour = 'grey50'),
axis.ticks = element_line(size = 0.25, colour = 'grey50'),
legend.justification = c(1, 0.75),
legend.key.size = unit(0.35, 'cm'),
legend.title = element_blank(),
legend.text = element_text(size = 9),
legend.text.align = 0, legend.background = element_blank(),
plot.subtitle = element_text(size = 10, vjust = 0),
strip.background = element_blank(),
strip.text = element_text(hjust = 0.5, size = 10, vjust = 0),
strip.placement = 'outside', panel.spacing.x = unit(-0., 'cm'))
if(str_length(effectplot$Conditions[1]) > 10){
ggobject <- ggobject +
theme(plot.margin = unit(c(0.2, 0.2, 0.5, 1.5), 'cm'),)
}
if(marker == TRUE){ if(marker == TRUE){
ggobject <- ggobject + ggtitle(label = marker_name) ggobject <- ggobject + ggtitle(label = marker_name)
} }
@@ -633,7 +654,25 @@ mash_plot_covar <- function(m, saveoutput = FALSE, suffix = ""){
levels = (df$`Covariance Matrix`)) levels = (df$`Covariance Matrix`))
ggobject <- ggplot(df) + ggobject <- ggplot(df) +
geom_bar(aes(x = .data$Mass, y = .data$`Covariance Matrix`), geom_bar(aes(x = .data$Mass, y = .data$`Covariance Matrix`),
stat = "identity") stat = "identity") +
theme_classic() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.background = element_rect(fill=NA),
axis.title = element_text(size = 10),
axis.text = element_text(size = 10),
axis.line.x = element_line(size = 0.35, colour = 'grey50'),
axis.line.y = element_line(size = 0.35, colour = 'grey50'),
axis.ticks = element_line(size = 0.25, colour = 'grey50'),
legend.justification = c(1, 0.75),
legend.key.size = unit(0.35, 'cm'),
legend.title = element_blank(),
legend.text = element_text(size = 9),
legend.text.align = 0, legend.background = element_blank(),
plot.subtitle = element_text(size = 10, vjust = 0),
strip.background = element_blank(),
strip.text = element_text(hjust = 0.5, size = 10 ,vjust = 0),
strip.placement = 'outside', panel.spacing.x = unit(-0.1, 'cm'))
if(saveoutput == TRUE){ if(saveoutput == TRUE){
if(!(str_sub(suffix, end = 1) %in% c("", "_"))){ if(!(str_sub(suffix, end = 1) %in% c("", "_"))){
@@ -699,12 +738,16 @@ mash_plot_Ulist <- function(m, range = NA, saveoutput = FALSE, suffix = "",
if(ncol(U1) == length(get_colnames(m))){ if(ncol(U1) == length(get_colnames(m))){
colnames(U1) <- str_replace(get_colnames(m), "(Stand_)??Bhat_?", "") %>% colnames(U1) <- str_replace(get_colnames(m), "(Stand_)??Bhat_?", "") %>%
str_replace("-mean$", "") str_replace("-mean$", "")
name_order <- str_replace(get_colnames(m), "(Stand_)??Bhat_?", "") %>%
str_replace("-mean$", "")
U1 <- as_tibble(U1, .name_repair = "unique") %>% U1 <- as_tibble(U1, .name_repair = "unique") %>%
mutate(rowU = str_replace(get_colnames(m), "(Stand_)??Bhat_?", ""), mutate(rowU = str_replace(get_colnames(m), "(Stand_)??Bhat_?", ""),
rowU = str_replace(.data$rowU, "-mean$", "")) %>% rowU = str_replace(.data$rowU, "-mean$", "")) %>%
pivot_longer(cols = -.data$rowU, names_to = "colU", pivot_longer(cols = -.data$rowU, names_to = "colU",
values_to = "covar") %>% values_to = "covar") %>%
filter(!is.na(.data$covar)) filter(!is.na(.data$covar))
U1$rowU <- factor(U1$rowU, levels = name_order)
U1$colU <- factor(U1$colU, levels = name_order)
} else { } else {
U1 <- as_tibble(U1, rownames = "rowU", .name_repair = "unique") %>% U1 <- as_tibble(U1, rownames = "rowU", .name_repair = "unique") %>%
pivot_longer(cols = -.data$rowU, names_to = "colU", pivot_longer(cols = -.data$rowU, names_to = "colU",
@@ -736,8 +779,24 @@ mash_plot_Ulist <- function(m, range = NA, saveoutput = FALSE, suffix = "",
"#5DC863FF", "#AADC32FF", "#FDE725FF"), "#5DC863FF", "#AADC32FF", "#FDE725FF"),
limits = c(lim_lower, lim_upper)) + limits = c(lim_lower, lim_upper)) +
#geom_text(aes(label = round(.data$covar, 1)), color = "darkgrey") + #geom_text(aes(label = round(.data$covar, 1)), color = "darkgrey") +
theme(legend.position = "right", theme_classic() +
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) + theme(panel.background = element_rect(fill=NA),
legend.position = "right",
axis.title = element_text(size = 10),
axis.text.y = element_text(size = 10),
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
axis.line.x = element_line(size = 0.35, colour = 'grey50'),
axis.line.y = element_line(size = 0.35, colour = 'grey50'),
axis.ticks = element_line(size = 0.25, colour = 'grey50'),
legend.justification = c(1, 0.75),
legend.key.size = unit(0.35, 'cm'),
legend.title = element_blank(),
legend.text = element_text(size = 9),
legend.text.align = 0, legend.background = element_blank(),
plot.subtitle = element_text(size = 10, vjust = 0),
strip.background = element_blank(),
strip.text = element_text(hjust = 0.5, size = 10 ,vjust = 0),
strip.placement = 'outside', panel.spacing.x = unit(-0.1, 'cm')) +
xlab("") + ylab("") + ggtitle(names(Ulist)[u]) xlab("") + ylab("") + ggtitle(names(Ulist)[u])
} }
if (labels == TRUE) { if (labels == TRUE) {
@@ -833,13 +892,25 @@ mash_plot_manhattan_by_condition <- function(m, snp, cond = NA,
facet_wrap(~ .data$CHR, nrow = 1, scales = "free_x", strip.position = "bottom") + facet_wrap(~ .data$CHR, nrow = 1, scales = "free_x", strip.position = "bottom") +
scale_color_viridis_c(option = "B", end = 0.95) + scale_color_viridis_c(option = "B", end = 0.95) +
scale_fill_viridis_c(option = "B", end = 0.95) + scale_fill_viridis_c(option = "B", end = 0.95) +
theme(strip.text = element_text(size = 8), theme_classic() +
axis.text.y = element_text(size = 8), theme(axis.text.x = element_blank(),
axis.text.x = element_blank(),
legend.text = element_text(size = 8),
axis.ticks.x = element_blank(), axis.ticks.x = element_blank(),
panel.background = element_rect(fill=NA), panel.background = element_rect(fill=NA),
legend.position = "right") + legend.position = "right",
axis.title = element_text(size = 10),
axis.text = element_text(size = 10),
axis.line.x = element_line(size = 0.35, colour = 'grey50'),
axis.line.y = element_line(size = 0.35, colour = 'grey50'),
axis.ticks = element_line(size = 0.25, colour = 'grey50'),
legend.justification = c(1, 0.75),
legend.key.size = unit(0.35, 'cm'),
legend.title = element_blank(),
legend.text = element_text(size = 9),
legend.text.align = 0, legend.background = element_blank(),
plot.subtitle = element_text(size = 10, vjust = 0),
strip.background = element_blank(),
strip.text = element_text(hjust = 0.5, size = 10 ,vjust = 0),
strip.placement = 'outside', panel.spacing.x = unit(-0.05, 'cm')) +
labs(x = "Chromosome", y = log10BF) + labs(x = "Chromosome", y = log10BF) +
scale_x_continuous(expand = c(0.32, 0.32)) + scale_x_continuous(expand = c(0.32, 0.32)) +
scale_shape_manual(values = rep(c(21, 22, 23), 9), guide = FALSE) scale_shape_manual(values = rep(c(21, 22, 23), 9), guide = FALSE)
@@ -949,14 +1020,50 @@ mash_plot_pairwise_sharing <- function(m = NULL, effectRDS = NULL,
label_size = label_size, hjust = hjust, vjust = vjust, label_size = label_size, hjust = hjust, vjust = vjust,
layout.exp = layout.exp, min_size = min_size, layout.exp = layout.exp, min_size = min_size,
max_size = max_size) + max_size = max_size) +
scale_color_viridis_c(option = option) scale_color_viridis_c(option = option) +
theme_classic() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.background = element_rect(fill=NA),
axis.title = element_text(size = 10),
axis.text = element_text(size = 10),
axis.line.x = element_line(size = 0.35, colour = 'grey50'),
axis.line.y = element_line(size = 0.35, colour = 'grey50'),
axis.ticks = element_line(size = 0.25, colour = 'grey50'),
legend.justification = c(1, 0.75),
legend.key.size = unit(0.35, 'cm'),
legend.title = element_blank(),
legend.text = element_text(size = 9),
legend.text.align = 0, legend.background = element_blank(),
plot.subtitle = element_text(size = 10, vjust = 0),
strip.background = element_blank(),
strip.text = element_text(hjust = 0.5, size = 10 ,vjust = 0),
strip.placement = 'outside', panel.spacing.x = unit(-0.1, 'cm'))
} else { } else {
corrplot <- ggcorr(data = NULL, cor_matrix = shared_effects, geom = geom, corrplot <- ggcorr(data = NULL, cor_matrix = shared_effects, geom = geom,
label = label, label_alpha = label_alpha, label = label, label_alpha = label_alpha,
label_size = label_size, hjust = hjust, vjust = vjust, label_size = label_size, hjust = hjust, vjust = vjust,
layout.exp = layout.exp, min_size = min_size, layout.exp = layout.exp, min_size = min_size,
max_size = max_size) + max_size = max_size) +
scale_color_viridis_c(option = option) scale_color_viridis_c(option = option) +
theme_classic() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.background = element_rect(fill=NA),
axis.title = element_text(size = 10),
axis.text = element_text(size = 10),
axis.line.x = element_line(size = 0.35, colour = 'grey50'),
axis.line.y = element_line(size = 0.35, colour = 'grey50'),
axis.ticks = element_line(size = 0.25, colour = 'grey50'),
legend.justification = c(1, 0.75),
legend.key.size = unit(0.35, 'cm'),
legend.title = element_blank(),
legend.text = element_text(size = 9),
legend.text.align = 0, legend.background = element_blank(),
plot.subtitle = element_text(size = 10, vjust = 0),
strip.background = element_blank(),
strip.text = element_text(hjust = 0.5, size = 10 ,vjust = 0),
strip.placement = 'outside', panel.spacing.x = unit(-0.1, 'cm'))
} }
if(saveoutput == TRUE){ if(saveoutput == TRUE){
@@ -1019,7 +1126,26 @@ mash_plot_sig_by_condition <- function(m, conditions = NA, saveoutput = FALSE,
geom_point() + geom_point() +
geom_hline(yintercept = 0, lty = 2) + geom_hline(yintercept = 0, lty = 2) +
xlab(label = "Number of Conditions") + xlab(label = "Number of Conditions") +
ylab(label = "Number of Significant SNPs") ylab(label = "Number of Significant SNPs")+
theme_classic() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.background = element_rect(fill=NA),
legend.position = "none",
axis.title = element_text(size = 10),
axis.text = element_text(size = 10),
axis.line.x = element_line(size = 0.35, colour = 'grey50'),
axis.line.y = element_line(size = 0.35, colour = 'grey50'),
axis.ticks = element_line(size = 0.25, colour = 'grey50'),
legend.justification = c(1, 0.75),
legend.key.size = unit(0.35, 'cm'),
legend.title = element_blank(),
legend.text = element_text(size = 9),
legend.text.align = 0, legend.background = element_blank(),
plot.subtitle = element_text(size = 10, vjust = 0),
strip.background = element_blank(),
strip.text = element_text(hjust = 0.5, size = 10 ,vjust = 0),
strip.placement = 'outside', panel.spacing.x = unit(-0.1, 'cm'))
if(saveoutput == TRUE){ if(saveoutput == TRUE){
if(!(str_sub(suffix, end = 1) %in% c("", "_"))){ if(!(str_sub(suffix, end = 1) %in% c("", "_"))){

71
man/dive_phe2mash.Rd
View File

@@ -1,5 +1,6 @@
% Generated by roxygen2: do not edit by hand % Generated by roxygen2: do not edit by hand
% Please edit documentation in R/wrapper.R % Please edit documentation in R/dive_effects2mash.R, R/dive_phe2effects.R,
% R/wrapper.R
\name{dive_phe2mash} \name{dive_phe2mash}
\alias{dive_phe2mash} \alias{dive_phe2mash}
\title{Wrapper to run mash given a phenotype data frame} \title{Wrapper to run mash given a phenotype data frame}
@@ -14,7 +15,47 @@ dive_phe2mash(
min.phe = 200, min.phe = 200,
save.plots = TRUE, save.plots = TRUE,
thr.r2 = 0.2, thr.r2 = 0.2,
thr.m = c("sum", "max"), thr.m = c("max", "sum"),
num.strong = 1000,
num.random = NA,
scale.phe = TRUE,
roll.size = 50,
U.ed = NA,
U.hyp = NA,
verbose = TRUE
)
dive_phe2mash(
df,
snp,
type = "linear",
svd = NULL,
suffix = "",
outputdir = ".",
min.phe = 200,
save.plots = TRUE,
thr.r2 = 0.2,
thr.m = c("max", "sum"),
num.strong = 1000,
num.random = NA,
scale.phe = TRUE,
roll.size = 50,
U.ed = NA,
U.hyp = NA,
verbose = TRUE
)
dive_phe2mash(
df,
snp,
type = "linear",
svd = NULL,
suffix = "",
outputdir = ".",
min.phe = 200,
save.plots = TRUE,
thr.r2 = 0.2,
thr.m = c("max", "sum"),
num.strong = 1000, num.strong = 1000,
num.random = NA, num.random = NA,
scale.phe = TRUE, scale.phe = TRUE,
@@ -76,12 +117,38 @@ matrices must have dimensions that match the number of phenotypes where
univariate GWAS ran successfully.} univariate GWAS ran successfully.}
\item{verbose}{Output some information on the iterations? Default is \code{TRUE}.} \item{verbose}{Output some information on the iterations? Default is \code{TRUE}.}
\item{effects}{fbm created using 'dive_phe2effects' or 'dive_phe2mash'.
Saved under the name "gwas_effects_{suffix}.rds" and can be loaded into
R using the bigstatsr function "big_attach".}
} }
\value{ \value{
A mash object made up of all phenotypes where univariate GWAS ran
successfully.
A mash object made up of all phenotypes where univariate GWAS ran
successfully.
A mash object made up of all phenotypes where univariate GWAS ran A mash object made up of all phenotypes where univariate GWAS ran
successfully. successfully.
} }
\description{ \description{
Though step-by-step GWAS, preparation of mash inputs, and mash
allows you the most flexibility and opportunities to check your results
for errors, once those sanity checks are complete, this function allows
you to go from a phenotype data.frame of a few phenotypes you want to
compare to a mash result. Some exception handling has been built into
this function, but the user should stay cautious and skeptical of any
results that seem 'too good to be true'.
This function allows
you to go from a phenotype data.frame of a few phenotypes you want to
compare to filebacked matrix of univariate GWAS effects, standard errors,
and -log10pvalues. This output object can be used in "dive_effects2mash"
function. Some exception handling has been built into
this function, but the user should stay cautious and skeptical of any
results that seem 'too good to be true'.
Though step-by-step GWAS, preparation of mash inputs, and mash Though step-by-step GWAS, preparation of mash inputs, and mash
allows you the most flexibility and opportunities to check your results allows you the most flexibility and opportunities to check your results
for errors, once those sanity checks are complete, this function allows for errors, once those sanity checks are complete, this function allows