Commit 5d4ff407 authored by Kunst, Jonathan's avatar Kunst, Jonathan
Browse files

Added parts for plotly visualisation

Added script to work through plotly book where needed
parent 65e3b4e2
# plotly ----
data("diamonds", package = 'ggplot2')
diamonds
plot_ly(diamonds, x = ~cut)
plot_ly(diamonds, x = ~cut, y = ~clarity)
plot_ly(diamonds, x = ~cut, y = ~clarity, colors = "Accent")
# plotly uses purely functional programming style
# input -> function modifies -> modified output
# layout expect a plot_ly object and modifies it according to function call
layout(plot_ly(diamonds, x = ~cut),
title = 'beautiful histogram')
# above code is less readable than
diamonds %>%
plot_ly(x = ~cut) %>%
layout(title = 'beautiful histogram')
# we can control layout even more by using add_* (i.e. add_histogram)
diamonds %>%
plot_ly() %>%
add_histogram(x = ~cut)
# or with bars where you need to calculate the statistics beforehand
diamonds %>%
count(cut) %>%
plot_ly() %>%
add_bars( x = ~cut, y = ~n)
# non-statistical layers (i.e. bars) are faster and more responsive but less flexible client-side
# another example of:
# - globally assigning x to cut
# - modify data after histogram is plotted
# - add another layer of text from modified summarised data
# note: make sure you want to display this information on the same axes
diamonds %>%
plot_ly(x = ~cut) %>%
add_histogram() %>%
group_by(cut) %>%
summarise(n = n()) %>%
add_text(
text = ~scales::comma(n), y = ~n,
textposition = "top middle",
cliponaxis = FALSE
)
# underlying plot_ly() is the json code which is captured in a list in plotly_build()
plot <- plot_ly(diamonds, x = ~cut, color = ~clarity, colors = "Accent")
# the underlying json
plotly_json(plot)
# the json figure contains data (traces) and layout
# a trace defines mapping from data and visuals
# every trace has a type and determines attributes
# using plotly_build we can debug a plot
b_plot <- plotly_build(plot)
length(b_plot$x$data)
# extract name of each trace (8 traces for 8 clarities) appearing in the figure definition x in data
map_chr(b_plot$x$data, "name")
length(unique(diamonds$clarity))
# the colors of the clarity are not directly coded in the json list
# instead plotly_build() designates colors with marker.color
m <- lm(log(price) ~ log(carat), data = diamonds)
diamonds <- modelr::add_residuals(diamonds, m)
diamonds <- modelr::add_predictions(diamonds, m)
head(diamonds)
p <- plot_ly(mpg, x = ~cty, y = ~hwy, alpha = 0.3)
subplot(
add_markers(p, symbol = ~cyl, name = "A single trace"),
add_markers(p, symbol = ~factor(cyl), color = I("black"))
)
# pass any number of plotly objects to subplot()
p1 <- plot_ly(economics, x = ~date, y = ~uempmed)
p2 <- plot_ly(economics, x = ~date, y = ~unemploy)
subplot(p1, p2, p1, p2, nrows = 2, margin = 0.05)
#' # anchor multiple traces on the same legend entry
p1 <- add_lines(p1, color = I("black"), name = "1st", legendgroup = "1st")
p2 <- add_lines(p2, color = I("red"), name = "2nd", legendgroup = "2nd")
subplot(
p1, style(p1, showlegend = FALSE),
p2, style(p2, showlegend = FALSE),
nrows = 2, margin = 0.05
)
\ No newline at end of file
library(polyqtlR)
library(plotly)
library(tidyverse)
library(scales)
library(trelliscopejs)
data("phased_maplist.4x", "SNP_dosages.4x", "Phenotypes_4x")
# IBD estimation
......@@ -10,12 +14,6 @@ IBD_4x <- estimate_IBD(phased_maplist = phased_maplist.4x,
bivalent_decoding = TRUE,
ncores = 2)
names(IBD_4x)
IBD_4x[['LG1']]
test <- 1:length(IBD_4x)
which(names(IBD_4x) == 'LG2')
# haplo visualisation
vis_hap <- visualiseHaplo(IBD_list = IBD_4x,
display_by = "name",
......@@ -46,6 +44,9 @@ plotQTL(LOD_data = qtl_LODs.4x,
col = "darkgreen")
# linear plot
# is made with standard plotting function, padding for splitting linkage groups
# e.g. length previous LG + padding = 110 = 0 linkage group 1
# See whether LGs can be selected with layout in plot_ly object
plotLinearQTL(LOD_data = qtl_LODs.4x,
col = "dodgerblue")
......@@ -62,7 +63,7 @@ plotLinearQTL(LOD_data = qtl_LODs.4x,
col = "dodgerblue")
# co-factor analysis
# Manual co-factor analysis ----
findPeak(qtl_LODs.4x, linkage_group = 1)
qtl_LODs.4x_cofactor <- QTLscan(IBD_list = IBD_4x,
......@@ -73,100 +74,91 @@ qtl_LODs.4x_cofactor <- QTLscan(IBD_list = IBD_4x,
cofactor_df = data.frame("LG" = 1,
"cM" = 12.3),
perm_test = FALSE,
ncores = 2)#nc is the number of cores, defined earlier
# plotly ----
data("diamonds", package = 'ggplot2')
diamonds
plot_ly(diamonds, x = ~cut)
plot_ly(diamonds, x = ~cut, y = ~clarity)
plot_ly(diamonds, x = ~cut, y = ~clarity, colors = "Accent")
# plotly uses purely functional programming style
# input -> function modifies -> modified output
# layout expect a plot_ly object and modifies it according to function call
layout(plot_ly(diamonds, x = ~cut),
title = 'beautiful histogram')
# above code is less readable than
diamonds %>%
plot_ly(x = ~cut) %>%
layout(title = 'beautiful histogram')
# we can control layout even more by using add_* (i.e. add_histogram)
diamonds %>%
plot_ly() %>%
add_histogram(x = ~cut)
# or with bars where you need to calculate the statistics beforehand
diamonds %>%
count(cut) %>%
plot_ly() %>%
add_bars( x = ~cut, y = ~n)
# non-statistical layers (i.e. bars) are faster and more responsive but less flexible client-side
# another example of:
# - globally assigning x to cut
# - modify data after histogram is plotted
# - add another layer of text from modified summarised data
# note: make sure you want to display this information on the same axes
diamonds %>%
plot_ly(x = ~cut) %>%
add_histogram() %>%
group_by(cut) %>%
summarise(n = n()) %>%
add_text(
text = ~scales::comma(n), y = ~n,
textposition = "top middle",
cliponaxis = FALSE
)
# underlying plot_ly() is the json code which is captured in a list in plotly_build()
plot <- plot_ly(diamonds, x = ~cut, color = ~clarity, colors = "Accent")
# the underlying json
plotly_json(plot)
ncores = 2)
# the json figure contains data (traces) and layout
# a trace defines mapping from data and visuals
# every trace has a type and determines attributes
new_pheno <- qtl_LODs.4x_cofactor$Residuals
head(new_pheno)
# using plotly_build we can debug a plot
b_plot <- plotly_build(plot)
length(b_plot$x$data)
# extract BLUES from new_pheno
blues <- BLUE(data = new_pheno,
model = Pheno~geno,
random = ~1|Block,
genotype.ID = "geno")
# extract name of each trace (8 traces for 8 clarities) appearing in the figure definition x in data
map_chr(b_plot$x$data, "name")
length(unique(diamonds$clarity))
# the colors of the clarity are not directly coded in the json list
# instead plotly_build() designates colors with marker.color
m <- lm(log(price) ~ log(carat), data = diamonds)
diamonds <- modelr::add_residuals(diamonds, m)
diamonds <- modelr::add_predictions(diamonds, m)
head(diamonds)
# qtl LOD plot two models example ----
lodsmod1 <- qtl_LODs.4x$QTL.res %>%
dplyr::mutate(model = 'model1')
# cofactor analysis with new_pheno blues
qtl_LODs.4x_cofactor <- QTLscan(IBD_list = IBD_4x,
Phenotype.df = blues,
genotype.ID = "geno",
trait.ID = "blue",
perm_test = TRUE,
ncores = 4)
plotLinearQTL(LOD_data = qtl_LODs.4x_cofactor,
col = "red")
# automatic cofactor analysis ----
# first extract blues to reduce computation time by simplifying data structure
# also by using the blues with co-factors analysis, the major QTL is "removed"
blues <- BLUE(data = Phenotypes_4x,
model = pheno~geno,
random = ~1|year,
genotype.ID = "geno")
# run check_cofactors for automatic cofactor analysis
QTLmodel <- check_cofactors(IBD_list = IBD_4x,
Phenotype.df = blues,
genotype.ID = "geno",
trait.ID = "blue",
LOD_data = qtl_LODs.4x,
ncores = 4)
QTLmodel
# handy to directly run PVE for seeing PVE per model/submodel
PVE(IBD_list = IBD_4x,
Phenotype.df = blues,
genotype.ID = "geno",
trait.ID = "blue",
QTL_df = QTLmodel)
# comparison of mulitple analyses in one plot
par(mfrow = c(1, 1))
plotLinearQTL_list(LOD_data.ls = list(qtl_LODs.4x,
qtl_LODs.4x_cofactor),
plot_type = "lines")
# further QTL analysis
visualiseQTLeffects(IBD_list = IBD_4x,
Phenotype.df = blues,
genotype.ID = "geno",
trait.ID = "blue",
linkage_group = 1,
LOD_data = qtl_LODs.4x)
# plotly tests ----
modellist <- list()
modellist[['nonfactor']] <- qtl_LODs.4x$QTL.res %>%
left_join(bind_rows(phased_maplist.4x), by = 'position')
lodsmod2 <- lodsmod1 %>%
mutate(model = 'model2',
LOD = jitter(LOD, 1, 0.25),
.keep = 'unused')
modellist[['co-factor']] <- qtl_LODs.4x_cofactor$QTL.res %>%
left_join(select(bind_rows(phased_maplist.4x), !marker), by = 'position')
lodsmod <- rbind(lodsmod1, lodsmod2) %>%
left_join(bind_rows(phased_maplist.4x), by = 'position')
p <- ggplot(lodsmod, aes(x = position, y = LOD, color = model, linetype = model)) +
p <- ggplot(bind_rows(modellist, .id = 'model'), aes(x = position,
y = LOD,
color = model,
linetype = model)) +
geom_point(aes(
text = paste(paste('LOD:', round(LOD, 2)),
paste('Marker:', marker),
sep = '\n'
))) +
geom_line() +
facet_grid(~chromosome,
switch = 'x',
space = 'free_x') +
facet_wrap(~chromosome) +
theme(
panel.grid = element_blank(),
panel.background = element_blank(),
......@@ -174,40 +166,19 @@ p <- ggplot(lodsmod, aes(x = position, y = LOD, color = model, linetype = model)
strip.background = element_blank(),
strip.placement = 'outside'
)
p
gp <- ggplotly(p) %>%
gp <- ggplotly(p, tooltip = 'text') %>%
highlight('plotly_selected') %>%
layout(title = list(text = 'Chromosome', xanchor = 'center'))
gp
# qtl LOD plot firstexample ----
dfs <- qtl_LODs.4x$QTL.res %>%
nest_by(chromosome)
lapply(seq_along(dfs$data), function(i){
plot_ly(dfs$data[[i]], x = ~position, y = ~LOD, type = 'scatter', mode = 'lines') %>%
layout(yaxis = list(title = 'Logarithm of Odds (LOD)', showgrid = FALSE),
xaxis = list(title = str(i), showgrid = FALSE))
}) %>%
subplot(., shareX = FALSE, shareY = TRUE, titleX = TRUE)
gp
# qtl LOD plot group_map----
palet <- park_palette('Hawaii')
panel <- . %>%
plot_ly(x = ~position, y = ~LOD) %>%
add_lines() %>%
add_annotations(
text = ~unique(chromosome),
x = 0.5,
y = 0,
yref = "paper",
xref = "paper",
yanchor = "bottom",
xanchor = 'center',
yshift = -35,
showarrow = FALSE,
font = list(size = 15)
) %>%
add_lines(color = ~model, legendgroup = ~model) %>%
layout(
showlegend = TRUE,
shapes = list(
......@@ -227,15 +198,34 @@ panel <- . %>%
yaxis = list(showgrid = FALSE)
)
qtl_LODs.4x$QTL.res %>%
mutate(chromosome = paste('Linkage Group', chromosome)) %>%
group_by(chromosome) %>%
do(p = panel(.)) %>%
subplot(nrows = 1, shareX = FALSE, shareY = TRUE)
bind_rows(modellist, .id = 'model' ) %>%
group_by(chromosome) %>%
group_map(~panel(.)) %>%
subplot(nrows = 1, shareX = TRUE, shareY=TRUE) %>%
highlight('plotly_selected')
plot_ly(qtl_LODs.4x$QTL.res,
x = ~position,
y = ~LOD)
# trelliscope example ----
p <- ggplot(bind_rows(modellist, .id = 'model'), aes(x = position,
y = LOD,
color = model,
linetype = model)) +
geom_point(aes(
text = paste(paste('LOD:', round(LOD, 2)),
paste('Marker:', marker),
sep = '\n'
))) +
geom_line() +
facet_trelliscope(~chromosome,
ncol = length(unique(bind_rows(modellist, .id = 'model')$chromosome)),
as_plotly = TRUE) +
theme(
panel.grid = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
strip.background = element_blank(),
strip.placement = 'outside'
)
p
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment