#' seaMass-Σ
#'
#' Fits the seaMass-Σ Bayesian group-level quantification model to imported data.
#'
#' @include seaMass.R
setClass("seaMass_sigma", contains = "seaMass_group_quants", slots = c(
filepath = "character"
))
#' @describeIn seaMass_sigma-class Runs seaMass-Σ.
#' @param data A \link{data.frame} of input data as returned by \link{import_MaxQuant}, \link{import_OpenSWATH},
#' \link{import_ProteinPilot} or \link{import_ProteomeDiscoverer}, .
#' @param data.design A \link{data.frame} created by \link{new_assay_design} and then customised, which specifies
#' assay names and block design.
#' @param run Run seaMass-Σ now, or just prepare it for later execution on e.g. a HPC cluster?
#' @param control A control object created with \link{sigma_control} specifying control parameters for the model.
#' @param path Name of folder prefix on disk where all intermediate and output data will be stored.
#' @return A \code{seaMass_sigma} object, which allows access to metadata and each block's \link{sigma_block} object to access
#' various results.
#' @import data.table
#' @export seaMass_sigma
seaMass_sigma <- function(
data,
data.design = new_assay_design(data),
path = "fit",
user = Sys.info()[["user"]],
run = TRUE,
control = sigma_control(),
...
) {
# check for finished output
object <- open_sigma(path, quiet = T)
if (!is.null(object)) {
stop(paste0("ERROR: completed seaMass-sigma already found at ", path))
}
### INIT
control@version <- as.character(packageVersion("seaMass"))
cat(paste0("[", Sys.time(), "] seaMass-sigma v", control@version, "\n"))
data.table::setDTthreads(control@nthread)
fst::threads_fst(control@nthread)
# setup seaMass folder structure
if (!grepl("\\.seaMass$", path)) path <- paste0(path, ".seaMass")
if (file.exists(path)) unlink(path, recursive = T)
if (!dir.create(file.path(path, "sigma"), recursive = T))
stop("ERROR: problem creating folder")
# setup seaMass-sigma folder structure
dir.create(file.path(path, "csv"))
path <- file.path(normalizePath(path), "sigma")
# init DT
data.is.data.table <- is.data.table(data)
DT <- setDT(data)
# get design into the format we need
DT.design <- as.data.table(data.design)
if (!is.factor(DT.design[, Assay])) DT.design[, Assay := factor(as.character(Assay), levels = unique(as.character(Assay)))]
if (all(is.na(DT.design[, Run]))) {
DT.design[, Run := NULL]
DT.design[, Run := "1"]
}
if (!is.factor(DT.design[, Run])) DT.design[, Run := factor(as.character(Run), levels = levels(DT[, Run]))]
if (all(is.na(DT.design[, Channel]))) {
DT.design[, Channel := NULL]
DT.design[, Channel := "1"]
}
if (!is.factor(DT.design[, Channel])) DT.design[, Channel := factor(as.character(Channel), levels = levels(DT[, Channel]))]
# extract blocks and save control
block.cols <- colnames(DT.design)[grep("^Block\\.(.*)$", colnames(DT.design))]
if (is.null(attr(DT, "group"))) {
control@group <- c("Group", "Groups")
} else {
control@group <- attr(DT, "group")
}
if (is.null(attr(DT, "component"))) {
control@component <- c("Component", "Components")
} else {
control@component <- attr(DT, "component")
}
if (is.null(attr(DT, "measurement"))) {
control@measurement <- c("Measurement", "Measurements")
} else {
control@measurement <- attr(DT, "measurement")
}
control@user <- user
control@blocks <- sub("^Block\\.(.*)$", "\\1", block.cols)
control@ellipsis <- list(...)
validObject(control)
saveRDS(control, file.path(path, "control.rds"))
# check ref weights
if ("RefWeight" %in% names(DT.design)) {
for (block.col in block.cols) {
if (all(DT.design$RefWeight[DT.design[[block.col]]] == 0)) DT.design$RefWeight[DT.design[[block.col]]] <- 1
}
} else {
DT.design[, RefWeight := 1]
}
# process each block independently
DT.design <- rbindlist(lapply(1:length(control@blocks), function(i) {
cat(paste0("[", Sys.time(), "] preparing block=", control@blocks[i], "...\n"))
# create output directory
blockpath <- file.path(path, control@blocks[i])
dir.create(blockpath)
# design for this block
DT.design.block <- DT.design[as.logical(get(block.cols[i]))]
DT.design.block[, (block.cols) := NULL]
DT.design.block[, Block := factor(control@blocks[i], levels = control@blocks)]
setcolorder(DT.design.block, "Block")
# data for this block
DT.block <- merge(DT, DT.design.block[, .(Run, Channel, Assay)], by = c("Run", "Channel"), sort = F)
setcolorder(DT.block, "Assay")
# if poisson model only integers are allowed, plus apply missingness.threshold
if (control@error.model == "poisson") {
DT.block[, Count0 := round(Count)]
DT.block[, Count0 := ifelse(Count0 <= control@missingness.threshold, NA_real_, Count0)]
} else {
DT.block[, Count0 := ifelse(Count <= control@missingness.threshold, NA_real_, Count)]
}
fst::write.fst(DT.block, file.path(blockpath, "data.fst"))
DT.block[, Count := NULL]
# write Assay index (design)
DT.design.block <- merge(DT.design.block, DT.block[, .(
A.qG = uniqueN(Group[Use & !is.na(Count0)]),
A.uG = uniqueN(Group[Use]),
A.nG = uniqueN(Group),
A.qC = uniqueN(Component[Use & !is.na(Count0)]),
A.uC = uniqueN(Component[Use]),
A.nC = uniqueN(Component),
A.qM = uniqueN(Measurement[Use & !is.na(Count0)]),
A.uM = uniqueN(Measurement[Use]),
A.nM = uniqueN(Measurement),
A.qD = sum(Use & !is.na(Count0)),
A.uD = sum(Use),
A.nD = length(Count0)
), by = .(Run, Channel)], by = c("Run", "Channel"))
setorder(DT.design, Assay, na.last = T)
setcolorder(DT.design, c("Assay", "Run", "Channel", colnames(DT.design)[grep("^Block\\.", colnames(DT.design))], "RefWeight"))
fst::write.fst(DT.design.block, file.path(blockpath, "design.fst"))
# write Group index
DT.groups <- DT.block[, .(
GroupInfo = GroupInfo[1],
G.qC = uniqueN(Component[Use & !is.na(Count0)]),
G.uC = uniqueN(Component[Use]),
G.nC = uniqueN(Component),
G.qM = uniqueN(Measurement[Use & !is.na(Count0)]),
G.uM = uniqueN(Measurement[Use]),
G.nM = uniqueN(Measurement),
G.qD = sum(Use & !is.na(Count0)),
G.uD = sum(Use),
G.nD = length(Count0)
), by = Group]
setorder(DT.groups, -G.qC, -G.uC, -G.nC, -G.qM, -G.uM, -G.nM, -G.qD, -G.uD, -G.nD, Group)
#setorder(DT.groups, as.character(Group))
# use pre-trained regression model to estimate how long each Group will take to process
# Intercept, uC, uM, uC^2, uM^2, uC*uM
a <- c(5.338861e-01, 9.991505e-02, 2.871998e-01, 4.294391e-05, 6.903229e-04, 2.042114e-04)
DT.groups[, pred.time := a[1] + a[2]*G.uC + a[3]*G.uM + a[4]*G.uC*G.uC + a[5]*G.uM*G.uM + a[6]*G.uC*G.uM]
fst::write.fst(DT.groups, file.path(blockpath, "groups.fst"))
# write Component index
DT.components <- DT.block[, .(
C.qM = uniqueN(Measurement[Use & !is.na(Count0)]),
C.uM = uniqueN(Measurement[Use]),
C.nM = uniqueN(Measurement),
C.qD = sum(Use & !is.na(Count0)),
C.uD = sum(Use),
C.nD = length(Count0)
), by = .(Group, Component)]
#setorder(DT.components, -C.qM, -C.uM, -C.nM, -C.qD, -C.uD, -C.nD, Component)
DT.components <- merge(DT.groups[, .(Group)], DT.components, by = "Group", sort = F)
fst::write.fst(DT.components, file.path(blockpath, "components.fst"))
# write Measurement index
DT.measurements <- DT.block[, .(
M.qD = sum(Use & !is.na(Count0)),
M.uD = sum(Use),
M.nD = sum(!is.na(Count0))
), by = .(Group, Component, Measurement)]
#setorder(DT.measurements, -M.qD, -M.uD, -M.nD, Measurement)
DT.measurements <- merge(DT.components[, .(Group, Component)], DT.measurements, by = c("Group", "Component"), sort = F)
fst::write.fst(DT.measurements, file.path(blockpath, "measurements.fst"))
rm(DT.measurements)
# now can get rid of not used
DT.block <- DT.block[Use == T & !is.na(Assay)]
DT.block[, Use := NULL]
# missingness model
if (control@missingness.model == "rm") {
DT.block <- DT.block[!is.na(Count0)]
} else {
# add in NAs to be imputed, removing Measurements where all NA
DT.block <- DT.block[, Use := !all(is.na(Count0)), by = .(Group, Component, Measurement)]
DT.block[, Assay := factor(Assay, levels = unique(DT.design.block$Assay))]
DT.block <- dcast(DT.block[Use == T], Group + Component + Measurement ~ Assay, value.var = "Count0", drop = c(T, F))
DT.block <- melt(DT.block, variable.name = "Assay", value.name = "Count0", id.vars = c("Group", "Component", "Measurement"))
DT.block[, Assay := factor(Assay, levels = levels(DT.design.block$Assay))]
if (substr(control@missingness.model, 1, 8) == "censored") {
DT.block[, Count1 := min(Count0, na.rm = T), by = .(Group, Component, Measurement)]
DT.block[, Count1 := ifelse(is.na(Count0), Count1, Count0)]
if (control@missingness.model == "censored_") DT.block[is.na(Count0), Count0 := ifelse(is.na(Count1), NA_real_, pmin(1.0, Count1))]
if (control@missingness.model == "censored0") DT.block[is.na(Count0), Count0 := Count1 / 2^0]
if (control@missingness.model == "censored1") DT.block[is.na(Count0), Count0 := Count1 / 2^1]
if (control@missingness.model == "censored2") DT.block[is.na(Count0), Count0 := Count1 / 2^2]
if (control@missingness.model == "censored3") DT.block[is.na(Count0), Count0 := Count1 / 2^3]
if (control@missingness.model == "censored4" || control@missingness.model == "censored") DT.block[is.na(Count0), Count0 := Count1 / 2^4]
if (control@missingness.model == "censored5") DT.block[is.na(Count0), Count0 := Count1 / 2^5]
if (control@missingness.model == "censored6") DT.block[is.na(Count0), Count0 := Count1 / 2^6]
if (control@missingness.model == "censored7") DT.block[is.na(Count0), Count0 := Count1 / 2^7]
if (control@missingness.model == "censored8") DT.block[is.na(Count0), Count0 := Count1 / 2^8]
if (control@missingness.model == "censored9") DT.block[is.na(Count0), Count0 := Count1 / 2^9]
if (control@error.model == "poisson") DT.block[, Count0 := round(Count)]
} else {
if (control@missingness.model == "one") DT.block[is.na(Count0), Count0 := 1.0]
if (control@missingness.model == "minimum") {
DT.block[, Count1 := min(Count0, na.rm = T), by = .(Group, Component, Measurement)]
DT.block[is.na(Count0), Count0 := Count1]
DT.block[, Count1 := NULL]
}
}
}
# set ordering for indexing
DT.block <- merge(DT.block, DT.groups[, .(Group, pred.time)], by = "Group", sort = F)
setorder(DT.block, -pred.time, Group, Component, Measurement, Assay)
DT.block[, pred.time := NULL]
# filter DT for Empirical Bayes model
DT0 <- unique(DT.block[, .(Group, Component, Measurement)])
DT0[, nMeasurement := .N, by = .(Group, Component)]
DT0 <- DT0[nMeasurement >= control@measurement.eb.min]
DT0[, nMeasurement := NULL]
DT0.components <- unique(DT0[, .(Group, Component)])
DT0.components[, nComponent := .N, by = Group]
if (control@component.model != "") DT0.components <- DT0.components[nComponent >= control@component.eb.min]
DT0.components[, nComponent := NULL]
DT0 <- merge(DT0, DT0.components, by = c("Group", "Component"), sort = F)
DT0 <- merge(DT.block, DT0, by = c("Group", "Component", "Measurement"), sort = F)
DT0.assays <- unique(DT0[, .(Group, Assay)])
DT0.assays[, nAssay := .N, by = Group]
DT0.assays <- DT0.assays[nAssay >= control@assay.eb.min]
DT0.assays[, nAssay := NULL]
DT0 <- merge(DT0, DT0.assays, by = c("Group", "Assay"), sort = F)
# filter to eb.max
if (control@component.model == "") {
DT0[, I := as.integer(factor(Measurement, levels = unique(Measurement)))]
DT0[, I := I[1], by = Measurement]
} else {
DT0[, I := as.integer(factor(Component, levels = unique(Component)))]
DT0[, I := I[1], by = Component]
}
DT0 <- DT0[I <= control@eb.max]
DT0[, I := NULL]
setcolorder(DT0, c("Group", "Component", "Measurement", "Assay"))
# save random access indices
if (!dir.create(file.path(blockpath, "model0"))) stop()
DT0.index <- DT0[, .(Group = unique(Group), file = factor("input.fst"), from = .I[!duplicated(Group)], to = .I[rev(!duplicated(rev(Group)))])]
fst::write.fst(DT0.index, file.path(blockpath, "model0", "input.index.fst"))
if (!dir.create(file.path(blockpath, "model1"))) stop()
DT.index <- DT.block[, .(Group = unique(Group), file = factor("input.fst"), from = .I[!duplicated(Group)], to = .I[rev(!duplicated(rev(Group)))])]
fst::write.fst(DT.index, file.path(blockpath, "model1", "input.index.fst"))
# convert factors to integers for efficiency for saving
DT0[, Group := as.integer(Group)]
DT0[, Component := as.integer(Component)]
DT0[, Measurement := as.integer(Measurement)]
DT0[, Assay := as.integer(Assay)]
fst::write.fst(DT0, file.path(blockpath, "model0", "input.fst"))
DT.block[, Group := as.integer(Group)]
DT.block[, Component := as.integer(Component)]
DT.block[, Measurement := as.integer(Measurement)]
DT.block[, Assay := as.integer(Assay)]
fst::write.fst(DT.block, file.path(blockpath, "model1", "input.fst"))
return(DT.design.block)
}))
### RUN
object <- new("seaMass_sigma", filepath = path)
prepare_sigma(control@schedule, object)
if (run) {
run(control@schedule, object)
} else {
cat(paste0("[", Sys.time(), "] queued\n"))
}
if (!data.is.data.table) setDF(data)
return(invisible(object))
}
#' @import data.table
#' @export
#' @include generics.R
setMethod("report", "seaMass_sigma", function(object, job.id) {
ctrl <- control(object)
if (ctrl@version != as.character(packageVersion("seaMass")))
stop(paste0("version mismatch - '", filepath(object), "' was prepared with seaMass v", ctrl@version, " but is running on v", packageVersion("seaMass")))
cat(paste0("[", Sys.time(), "] REPORT\n"))
# delete model0
if (!("model0" %in% ctrl@keep)) for (block in blocks(object)) unlink(file.path(filepath(block), "model0", "component.stdevs*"), recursive = T)
if (!("model0" %in% ctrl@keep)) for (block in blocks(object)) unlink(file.path(filepath(block), "model0", "measurement.stdevs*"), recursive = T)
if (!("model0" %in% ctrl@keep)) for (block in blocks(object)) unlink(file.path(filepath(block), "model0", "assay.deviations*"), recursive = T)
# delete model1
if (!("measurement.means" %in% ctrl@keep)) for (block in blocks(object)) unlink(file.path(filepath(block), "model1", "measurement.means*"), recursive = T)
if (!("measurement.stdevs" %in% ctrl@keep)) for (block in blocks(object)) unlink(file.path(filepath(block), "model1", "measurement.stdevs*"), recursive = T)
if (!("component.means" %in% ctrl@keep)) for (block in blocks(object)) unlink(file.path(filepath(block), "model1", "component.means*"), recursive = T)
if (!("component.stdevs" %in% ctrl@keep)) for (block in blocks(object)) unlink(file.path(filepath(block), "model1", "component.stdevs*"), recursive = T)
if (!("assay.deviations" %in% ctrl@keep)) for (block in blocks(object)) unlink(file.path(filepath(block), "model1", "assay.deviations*"), recursive = T)
if (!("component.deviations" %in% ctrl@keep)) for (block in blocks(object)) unlink(file.path(filepath(block), "model1", "component.deviations*"), recursive = T)
if (!("group.means" %in% ctrl@keep)) for (block in blocks(object)) unlink(file.path(filepath(block), "model1", "group.means*"), recursive = T)
if (!("group.quants" %in% ctrl@keep)) for (block in blocks(object)) unlink(file.path(filepath(block), "model1", "group.quants*"), recursive = T)
# assemble report
cat(paste0("[", Sys.time(), "] assembling html report...\n"))
assemble_report(object)
return(invisible(NULL))
})
#' @describeIn seaMass_sigma-class Open a complete \code{seaMass_sigma} run from the supplied \code{path}.
#' @export
open_sigma <- function(
path = "fit.seaMass",
quiet = FALSE,
force = FALSE
) {
path2 <- file.path(ifelse(dir.exists(path), path, paste0(path, ".seaMass")), "sigma")
if (!dir.exists(path2)) {
if (quiet) {
return(NULL)
} else {
stop("'", path, "' does not exist")
}
}
object <- new("seaMass_sigma", filepath = normalizePath(path2))
if (!force && read_completed(filepath(object)) == 0) {
if (quiet) {
return(NULL)
} else {
stop("'", path, "' is not complete")
}
}
return(object)
}
#' @describeIn seaMass_sigma-class Delete the \code{seaMass_sigma} run from disk.
#' @export
#' @include generics.R
setMethod("del", "seaMass_sigma", function(object) {
return(unlink(dirname(filepath(object)), recursive = T))
})
#' @describeIn seaMass_sigma-class Get the \code{seaMass_sigma} object.
#' @export
#' @include generics.R
setMethod("root", "seaMass_sigma", function(object) {
return(object)
})
#' @describeIn seaMass_sigma-class Get name.
#' @export
#' @include generics.R
setMethod("name", "seaMass_sigma", function(object) {
return(sub("\\.seaMass$", "", basename(dirname(filepath(object)))))
})
#' @describeIn seaMass_sigma-class Get path.
#' @export
#' @include generics.R
setMethod("filepath", "seaMass_sigma", function(object) {
return(object@filepath)
})
#' @describeIn seaMass_sigma-class Run.
#' @export
#' @include generics.R
setMethod("run", "seaMass_sigma", function(object) {
run(control(object)@schedule, object)
return(invisible(object))
})
#' @describeIn seaMass_sigma-class Get the \link{sigma_control}.
#' @export
#' @include generics.R
setMethod("control", "seaMass_sigma", function(object) {
return(readRDS(file.path(filepath(object), "control.rds")))
})
#' @describeIn seaMass_sigma-class Get the groups as a \code{data.frame}.
#' @import data.table
#' @export
#' @include generics.R
setMethod("groups", "seaMass_sigma", function(object, summary = FALSE, as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) groups(block, as.data.table = T)))
if (summary) DT <- DT[, .(GroupInfo = GroupInfo[1], G.qC = max(G.qC), G.uC = max(G.uC), G.nC = max(G.nC), G.qM = max(G.qM), G.uM = max(G.uM), G.nM = max(G.nM), G.qD = max(G.qD), G.uD = max(G.uD), G.nD = max(G.nD)), by = Group]
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' @describeIn seaMass_sigma-class Get the n groups with the most quantified components that are present in all blocks \code{data.frame}.
#' @import data.table
#' @export
#' @include generics.R
setMethod("top_groups", "seaMass_sigma", function(object, n = 512) {
DT <- groups(object, as.data.table = T)
DT <- DT[, .(.N, qC = min(G.qC)), by = Group][N == nlevels(DT$Block)]
setorder(DT, -qC)
return(as.character(DT[1:min(n, nrow(DT)), Group]))
})
#' @describeIn seaMass_sigma-class Get the groups as a \code{data.frame}.
#' @import data.table
#' @export
#' @include generics.R
setMethod("components", "seaMass_sigma", function(object, summary = FALSE, as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) components(block, as.data.table = T)))
if (summary) DT <- DT[, .(C.qM = max(C.qM), C.uM = max(C.uM), C.nM = max(C.nM), C.qD = max(C.qD), C.uD = max(C.uD), C.nD = max(C.nD)), by = .(Group, Component)]
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' @describeIn seaMass_sigma-class Get the assay groups as a \code{data.frame}.
#' @import data.table
#' @export
#' @include generics.R
setMethod("assay_groups", "seaMass_sigma", function(object, summary = FALSE, as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) assay_groups(block, as.data.table = T)))
if (summary) DT <- DT[, .(AG.qC = max(AG.qC), AG.uC = max(AG.uC), AG.nC = max(AG.nC), AG.qM = max(AG.qM), AG.uM = max(AG.uM), AG.nM = max(AG.nM), AG.qD = max(AG.qD), AG.uD = max(AG.uD), AG.nD = max(AG.nD)), by = .(Group, Assay)]
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' @describeIn seaMass_sigma-class Get the assay groups as a \code{data.frame}.
#' @import data.table
#' @export
#' @include generics.R
setMethod("assay_components", "seaMass_sigma", function(object, summary = FALSE, as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) assay_components(block, as.data.table = T)))
if (summary) DT <- DT[, .(AC.qM = max(AC.qM), AC.uM = max(AC.uM), AC.nM = max(AC.nM), AC.qD = max(AC.qD), AC.uD = max(AC.uD), AC.nD = max(AC.nD)), by = .(Group, Component, Assay)]
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' @describeIn seaMass_sigma-class Get the groups as a \code{data.frame}.
#' @import data.table
#' @export
#' @include generics.R
setMethod("measurements", "seaMass_sigma", function(object, summary = FALSE, as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) measurements(block, as.data.table = T)))
if (summary) DT <- DT[, .(M.qD = max(M.qD), M.uD = max(M.uD), M.nD = max(M.nD)), by = .(Group, Component, Measurement)]
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' @describeIn seaMass_sigma-class Get the list of \link{sigma_block} objects for the blocks.
#' @export
#' @include generics.R
setMethod("blocks", "seaMass_sigma", function(object) {
blocks <- control(object)@blocks
if (length(blocks) == 0)
stop(paste0("seaMass-sigma output '", sub("\\.seaMass$", "", basename(filepath(object))), "' is missing"))
names(blocks) <- blocks
blocks <- lapply(blocks, function(block) new("sigma_block", filepath = file.path(filepath(object), block)))
return(blocks)
})
#' @describeIn seaMass_sigma-class Open the list of \link{seaMass_delta} objects.
#' @export
#' @include generics.R
setMethod("open_thetas", "seaMass_sigma", function(object, quiet = FALSE, force = FALSE) {
deltas <- lapply(sub("^delta\\.", "", list.files(filepath(object), "^delta\\.*")), function(name) open_delta(object, name, quiet, force))
names(deltas) <- lapply(deltas, function(delta) name(delta))
return(deltas)
})
#' @describeIn seaMass_sigma-class Open the list of \link{seaMass_delta} objects.
#' @export
#' @include generics.R
setMethod("open_deltas", "seaMass_sigma", function(object, quiet = FALSE, force = FALSE) {
deltas <- lapply(sub("^delta\\.", "", list.files(filepath(object), "^delta\\.*")), function(name) open_delta(object, name, quiet, force))
names(deltas) <- lapply(deltas, function(delta) name(delta))
return(deltas)
})
#' @describeIn seaMass_sigma-class Get the study design as a \code{data.frame}.
#' @import data.table
#' @export
#' @include generics.R
setMethod("assay_design", "seaMass_sigma", function(object, as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) assay_design(block, as.data.table = T)))
if (nrow(DT) == 0) return(NULL)
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' @describeIn seaMass_sigma-class Get the model timings as a \link{data.frame}.
#' @import data.table
#' @export
#' @include generics.R
setMethod("timings", "seaMass_sigma", function(object, input = "model1", as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) timings(block, input, as.data.table = T)))
if (nrow(DT) == 0) return(NULL)
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' @import data.table
#' @export
#' @include generics.R
setMethod("read", "seaMass_sigma", function(object, input, type, items = NULL, chains = 1:control(object)@nchain, summary = NULL, summary.func = "robust_normal", as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) read(block, input, type, items, chains, summary, summary.func, as.data.table = T)))
if (is.null(DT)) return(NULL)
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' @describeIn seaMass_sigma-class Get priors as a \link{data.frame}.
#' @import data.table
#' @export
#' @include generics.R
setMethod("priors", "seaMass_sigma", function(object, input = "model1", as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) priors(block, input, as.data.table = T)))
if (nrow(DT) == 0) return(NULL)
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' @describeIn seaMass_sigma-class Get the model assay stdevs as a \link{data.frame}.
#' @import data.table
#' @export
#' @include generics.R
setMethod("assay_stdevs", "seaMass_sigma", function(object, input = "model1", as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) assay_stdevs(block, input, as.data.table = T)))
if (nrow(DT) == 0) return(NULL)
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' @describeIn seaMass_sigma-class Get the model assay deviations as a \link{data.frame}.
#' @import doRNG
#' @import data.table
#' @export
#' @include generics.R
setMethod("assay_deviations", "seaMass_sigma", function(object, assays = NULL, summary = TRUE, input = "model1", chains = 1:control(object)@nchain, as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) assay_deviations(block, assays, summary, input, chains, as.data.table = T)))
if (nrow(DT) == 0) return(NULL)
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' @describeIn seaMass_sigma-class Get the model measurement means as a \link{data.frame}.
#' @import data.table
#' @export
#' @include generics.R
setMethod("measurement_means", "seaMass_sigma", function(object, measurements = NULL, summary = TRUE, input = "model1", chains = 1:control(object)@nchain, as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) measurement_means(block, measurements, summary, input, chains, as.data.table = T)))
if (nrow(DT) == 0) return(NULL)
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' @describeIn seaMass_sigma-class Get the model measurement stdevs as a \link{data.frame}.
#' @import data.table
#' @export
#' @include generics.R
setMethod("measurement_stdevs", "seaMass_sigma", function(object, measurements = NULL, summary = TRUE, input = "model1", chains = 1:control(object)@nchain, as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) measurement_stdevs(block, measurements, summary, input, chains, as.data.table = T)))
if (nrow(DT) == 0) return(NULL)
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' @describeIn seaMass_sigma-class Get the model component means as a \link{data.frame}.
#' @import data.table
#' @export
#' @include generics.R
setMethod("component_means", "seaMass_sigma", function(object, components = NULL, summary = TRUE, input = "model1", chains = 1:control(object)@nchain, as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) component_means(block, components, summary, input, chains, as.data.table = T)))
if (nrow(DT) == 0) return(NULL)
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' @describeIn seaMass_sigma-class Get the model component stdevs as a \link{data.frame}.
#' @import data.table
#' @export
#' @include generics.R
setMethod("component_stdevs", "seaMass_sigma", function(object, components = NULL, summary = TRUE, input = "model1", chains = 1:control(object)@nchain, as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) component_stdevs(block, components, summary, input, chains, as.data.table = T)))
if (nrow(DT) == 0) return(NULL)
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' @describeIn seaMass_sigma-class Gets the model component deviations as a \link{data.frame}.
#' @import data.table
#' @export
#' @include generics.R
setMethod("component_deviations", "seaMass_sigma", function(object, components = NULL, summary = TRUE, input = "model1", chains = 1:control(object)@nchain, as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) component_deviations(block, components, summary, input, chains, as.data.table = T)))
if (nrow(DT) == 0) return(NULL)
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' @describeIn seaMass_sigma-class Get the model group quantifications as a \link{data.frame}.
#' @import data.table
#' @export
#' @include generics.R
setMethod("group_quants", "seaMass_sigma", function(object, groups = NULL, summary = TRUE, input = "model1", chains = 1:control(object)@nchain, as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) group_quants(block, groups, summary, input, chains, as.data.table = T)))
if (nrow(DT) == 0) return(NULL)
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' @describeIn seaMass_sigma-class Get the model group means as a \link{data.frame}.
#' @import data.table
#' @export
#' @include generics.R
setMethod("group_means", "seaMass_sigma", function(object, groups = NULL, summary = TRUE, input = "model1", chains = 1:control(object)@nchain, as.data.table = FALSE) {
DT <- rbindlist(lapply(blocks(object), function(block) group_means(block, groups, summary, input, chains, as.data.table = T)))
if (nrow(DT) == 0) return(NULL)
if (!as.data.table) setDF(DT)
else DT[]
return(DT)
})
#' #' @import data.table
#' #' @export
#' #' @include generics.R
#' setMethod("plot_priors", "seaMass_sigma", function(
#' object,
#' data = list(
#' priors(object, as.data.table = T)[is.na(Assay)][, .(Block, Effect, s, df)],
#' priors(object, as.data.table = T)[is.na(Assay)][, .(Block, Effect, s = s0, df = df0)],
#' rbind(
#' measurement_stdevs(object, input = "model0", summary = T, as.data.table = T)[, .(Block, Effect = "Measurements", value = rinaka(length(s), s, df))],
#' component_stdevs(object, input = "model0", summary = T, as.data.table = T)[, .(Block, Effect = "Components", value = rinaka(length(s), s, df))]
#' )
#' ),
#' horizontal = TRUE,
#' draw_outline = list(TRUE, TRUE, FALSE),
#' draw_quantiles = list(0.5, NULL, NULL),
#' trim = c(0.05, 0.95),
#' colour = list("blue", "black", NULL),
#' fill = list("lightblue", NULL, "grey"),
#' alpha = list(0.5, 0.5, 0.5),
#' facets = "Effect",
#' value.label = "stdev",
#' value.limits = limits_dists(data, trim, c(0, 1), include.zero = T),
#' value.length = 160,
#' variable.sort.cols = NULL,
#' variable.label.cols = "Block",
#' variable.interval = 5,
#' show.legend = TRUE,
#' file = NULL
#' ) {
#' return(plot_dists(object, data, horizontal, draw_outline, draw_quantiles, trim, colour, fill, alpha, facets, value.label, value.limits, value.length, variable.sort.cols, variable.label.cols, variable.interval, show.legend, file))
#' })
#' @import data.table
#' @export
#' @include generics.R
setMethod("plot_assay_stdevs", "seaMass_sigma", function(
object,
data = list(
assay_stdevs(object, as.data.table = T)[, list(Block, Assay, s, df)],
assay_stdevs(object, as.data.table = T)[, list(Block, Assay, s = B.sC, df = B.dfC)],
assay_stdevs(object, as.data.table = T)[, list(Block, Assay, s = B.sM, df = B.dfM)]
),
draw_quantiles = list(0.5, NULL, NULL),
trim = c(0.05, 0.95),
colour = list("A.qM", NULL, NULL),
fill = list(NULL, "darkgreen", "black"),
alpha = list(0.75, 0.2, 0.2),
value.label = "stdev",
value.limits = limits_dists(data, trim, include.zero = T, non.negative = T),
variable.summary.cols = c("Block", "Run", "Channel", "Assay", "RefWeight", "Sample", "Condition", "A.qG", "A.qC", "A.qM", "A.qD"),
variable.label.cols = c("Sample", "Assay", "Block"),
...
) {
return(plot_dists(
object = object,
data = data,
draw_quantiles = draw_quantiles,
trim = trim,
colour = colour,
fill = fill,
alpha = alpha,
value.label = value.label,
value.limits = value.limits,
variable.summary.cols = variable.summary.cols,
variable.label.cols = variable.label.cols,
...
))
})
#' @import data.table
#' @export
#' @include generics.R
setMethod("plot_group_means", "seaMass_sigma", function(
object,
groups = NULL,
summary = TRUE,
colour = "Block",
value.label = "mean",
variable.summary.cols = c("Group", "Block", "G.qC", "G.qM", "G.qD"),
variable.label.cols = "Group",
...
) {
return(plot_dists(
object,
data = group_means(object, groups, summary = summary, as.data.table = T),
colour = colour,
value.label = value.label,
variable.summary.cols = variable.summary.cols,
variable.label.cols = variable.label.cols,
...
))
})
#' @import data.table
#' @export
#' @include generics.R
setMethod("plot_group_quants", "seaMass_sigma", function(
object,
group,
summary = TRUE,
colour = "Condition",
variable.summary.cols = c("Group", "Block", "Run", "Channel", "Assay", "RefWeight", "Sample", "Condition", "AG.qC", "AG.qM", "AG.qD"),
variable.label.cols = c("Sample", "Assay", "Block"),
value.label = "quant",
...
) {
return(plot_dists(
object,
data = group_quants(object, group, summary = summary, as.data.table = T),
colour = colour,
variable.summary.cols = variable.summary.cols,
variable.label.cols = variable.label.cols,
value.label = value.label,
...
))
})
#' @import data.table
#' @export
#' @include generics.R
setMethod("plot_component_means", "seaMass_sigma", function(
object,
group,
summary = TRUE,
colour = "Block",
value.label = "mean",
variable.summary.cols = c("Group", "Component", "Block", "C.qM", "C.qD"),
variable.label.cols = "Component",
...
) {
return(plot_dists(
object,
data = component_means(object, group, summary = summary, as.data.table = T),
colour = colour,
value.label = value.label,
variable.summary.cols = variable.summary.cols,
variable.label.cols = variable.label.cols,
...
))
})
#' @import data.table
#' @export
#' @include generics.R
setMethod("plot_component_stdevs", "seaMass_sigma", function(
object,
group,
summary = TRUE,
colour = "Block",
value.label = "stdev",
variable.summary.cols = c("Group", "Component", "Block", "C.qM", "C.qD"),
variable.label.cols = "Component",
...
) {
return(plot_dists(
object,
data = component_stdevs(object, group, summary = summary, as.data.table = T),
colour = colour,
value.label = value.label,
variable.summary.cols = variable.summary.cols,
variable.label.cols = variable.label.cols,
...
))
})
#' @import data.table
#' @export
#' @include generics.R
setMethod("plot_component_deviations", "seaMass_sigma", function(
object,
group,
summary = TRUE,
colour = "Condition",
value.label = "deviation",
variable.summary.cols = c("Group", "Component", "Block", "Run", "Channel", "Assay", "RefWeight", "Sample", "Condition", "AC.qM", "AC.qD"),
variable.label.cols = c("Component", "Sample", "Assay", "Block"),
...
) {
return(plot_dists(
object,
data = component_deviations(object, group, summary = summary, as.data.table = T),
colour = colour,
value.label = value.label,
variable.summary.cols = variable.summary.cols,
variable.label.cols = variable.label.cols,
...
))
})
#' @import data.table
#' @export
#' @include generics.R
setMethod("plot_measurement_means", "seaMass_sigma", function(
object,
group,
summary = TRUE,
colour = "Block",
value.label = "mean",
variable.summary.cols = c("Group", "Component", "Measurement", "Block", "M.qD"),
variable.label.cols = c("Component", "Measurement"),
...
) {
return(plot_dists(
object,
data = measurement_means(object, group, summary = summary, as.data.table = T),
colour = colour,
value.label = value.label,
variable.summary.cols = variable.summary.cols,
variable.label.cols = variable.label.cols,
...
))
})
#' @import data.table
#' @export
#' @include generics.R
setMethod("plot_measurement_stdevs", "seaMass_sigma", function(
object,
group,
summary = TRUE,
colour = "Block",
value.label = "stdev",
variable.summary.cols = c("Group", "Component", "Measurement", "Block", "M.qD"),
variable.label.cols = c("Component", "Measurement"),
...
) {
return(plot_dists(
object,
data = measurement_stdevs(object, group, summary = summary, as.data.table = T),
colour = colour,
value.label = value.label,
variable.summary.cols = variable.summary.cols,
variable.label.cols = variable.label.cols,
...
))
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.