Nothing
library(testthat)
acontext("geom line")
data(intreg, package = "animint2")
min.logratio <- min(intreg$signals$logratio)-0.2
max.logratio <- max(intreg$signals$logratio)
intreg$breaks$min.logratio <- min.logratio
intreg$breaks$max.logratio <- max.logratio
signal.colors <- c(estimate="#0adb0a", latent="#0098ef")
breakpoint.colors <- c("1breakpoint"="#ff7d7d", "0breakpoints"='#f6f4bf')
models.by.signal <- with(intreg, split(segments, segments$signal))
signals.by.signal <- with(intreg, split(signals, signals$signal))
model.selection.list <- list()
sig.labels.list <- list()
sig.seg.names.list <- list()
sig.names.list <- list()
for(signal.name in names(models.by.signal)){
signal <- signals.by.signal[[signal.name]]
signal.models <- models.by.signal[[signal.name]]
models.by.segments <- split(signal.models, signal.models$segments)
sig.seg.by.segments <- list()
for(segments.str in names(models.by.segments)){
segments <- as.numeric(segments.str)
model <- models.by.segments[[segments.str]]
model$rss <- NA
model$data <- NA
for(segment.i in 1:nrow(model)){
one.segment <- model[segment.i, ]
seg.data <-
subset(signal,
one.segment$first.base < base &
base < one.segment$last.base)
model$data[segment.i] <- nrow(seg.data)
residual.vec <- seg.data$logratio - one.segment$mean
model$rss[segment.i] <- sum(residual.vec * residual.vec)
}
stopifnot(sum(model$data) == nrow(signal))
model.stats <-
data.frame(signal=signal.name,
segments,
error=sum(model$rss),
data=sum(model$data))
model.selection.list[[paste(signal.name, segments.str)]] <- model.stats
if(segments.str == "1"){
sig.labels.list[[signal.name]] <- model.stats
}
sig.seg.by.segments[[segments.str]] <-
data.frame(signal=signal.name,
segments, min.logratio, max.logratio,
base=(min(signal$base)+max(signal$base))/2)
}
sig.seg.tall <- do.call(rbind, sig.seg.by.segments)
sig.seg.names.list[[signal.name]] <- sig.seg.tall
sig.names.list[[signal.name]] <- sig.seg.tall[1,]
}
model.selection.sorted <- do.call(rbind, model.selection.list)
set.seed(1)
model.selection <-
model.selection.sorted[sample(1:nrow(model.selection.sorted)),]
sig.seg.names <- do.call(rbind, sig.seg.names.list)
sig.names <- do.call(rbind, sig.names.list)
sig.labels <- do.call(rbind, sig.labels.list)
## Plot segments rather than penalty.
mmir.selection <-
list(error=ggplot()+
ggtitle("Select profile and number of segments")+
make_tallrect(model.selection, "segments",
colour=signal.colors[["estimate"]])+
theme_bw()+
theme_animint(width=600)+
theme(panel.margin=grid::unit(0, "lines"))+
facet_grid(. ~ geom)+
geom_text(aes(0, error, label=signal),
clickSelects="signal",
data=sig.labels, hjust=1)+
scale_x_continuous("segments", breaks=c(1, 5, 10, 20),
limits=c(-2, 20))+
xlab("squared error")+
geom_line(aes(segments, error,
group=signal),
data=data.frame(model.selection, geom="line"),
clickSelects="signal",
alpha=0.6, size=8)+
geom_path(aes(segments, error,
group=signal),
data=data.frame(model.selection, geom="path"),
clickSelects="signal",
alpha=0.6, size=8),
signal=ggplot()+
theme_bw()+
theme_animint(width=800)+
scale_x_continuous("position on chromosome (mega base pairs)",
breaks=c(100,200))+
scale_fill_manual(values=breakpoint.colors,guide="none")+
geom_blank(aes(first.base/1e6, logratio+2/8), data=intreg$ann)+
ggtitle("Copy number profile and maximum likelihood segmentation")+
ylab("logratio")+
geom_point(aes(base/1e6, logratio),
data=intreg$sig,
showSelected="signal")+
geom_segment(aes(first.base/1e6, mean, xend=last.base/1e6, yend=mean),
data=intreg$seg, colour=signal.colors[["estimate"]],
showSelected=c("signal", "segments"))+
geom_segment(aes(base/1e6, min.logratio,
xend=base/1e6, yend=max.logratio),
colour=signal.colors[["estimate"]],
showSelected=c("signal", "segments"),
linetype="dashed",
data=intreg$breaks)+
geom_text(aes(base/1e6, max.logratio, label=paste("signal", signal)),
data=sig.names,
showSelected="signal")+
geom_text(aes(base/1e6, min.logratio,
label=sprintf("%d segment%s", segments,
ifelse(segments==1, "", "s"))),
data=sig.seg.names,
showSelected=c("signal", "segments"),
color=signal.colors[["estimate"]]),
first=list(signal="4.2", segments=4))
tdir <- tempfile()
animint2dir(mmir.selection, tdir, open.browser=FALSE)
all.increasing <- function(num.vec){
stopifnot(is.numeric(num.vec))
all(0 < diff(num.vec))
}
expected.list <-
list(geom3_line_error=TRUE,
geom4_path_error=FALSE)
result.list <- list()
for(g.class in names(expected.list)){
expected <- expected.list[[g.class]]
tsv.path <- Sys.glob(file.path(tdir, paste0(g.class, "*")))
g.data <- read.table(tsv.path, header=TRUE, comment.char = "")
tsv.by.signal <- split(g.data, g.data$clickSelects)
for(signal.name in names(tsv.by.signal)){
one.signal <- tsv.by.signal[[signal.name]]
computed <- all.increasing(one.signal$x)
result.list[[paste(g.class, signal.name)]] <-
data.frame(g.class, signal.name, computed, expected)
}
}
result <- do.call(rbind, result.list)
test_that("line sorts tsv data by x value, path does not", {
with(result, {
expect_identical(computed, expected)
})
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.