#' Fit a 4-parameter log-logistic model for each analyte in
#' a MagPix kit
#'
#' @md
#' @param d Output from `assayr2::()`
#' @param wells_to_omit If not NA, the standard wells to remove.
#' @return list with 4 elements: data, curves, models, and omitted data
#' @examples
#' \dontrun{
#' file <- system.file("extdata", "Magpix_example.csv", package = "assayr2")
#' d <- read_magpix(file)
#' sc <- fitStandardCurves(d)
#' }
#' @importFrom dplyr group_by filter inner_join mutate select_if
#' @importFrom tidyr nest unnest
#' @importFrom purrr map
#' @importFrom drc drm
#' @export
mp_fit_standards <- function(d, wells_to_omit=NA) {
model_data <- tibble::as.tibble(d[[2]]) %>%
dplyr::group_by(target) %>%
tidyr::nest()
obs_standard <- dplyr::filter(d[[1]], grepl("Standard", content))
# Remove cases with 0 intensity
to_remove <- dplyr::filter(obs_standard, is.infinite(log10_intensity))
if (nrow(to_remove) > 0) {
msg <- paste0("Removed ", nrow(to_remove), " standard observations for 0 signal")
warning(msg)
obs_standard <- dplyr::filter(obs_standard, !(is.infinite(log10_intensity)))
}
# Remove cases specified by user
if (!is.na(wells_to_omit)) {
obs_standard %<>% dplyr::filter(!well %in% wells_to_omit)
}
standard_data <- dplyr::inner_join(d[[2]], obs_standard) %>%
tibble::as.tibble() %>%
dplyr::group_by(target) %>%
tidyr::nest() %>%
dplyr::mutate(model = purrr::map(data,
~ drc::drm(log10_intensity ~ concentration, data = ., fct = drc::LL.4())),
curves = purrr::map(model,
~ get_curves(.)))
models <- standard_data$model
names(models) <- standard_data$target
# returns a 4 element list
list(
data = tidyr::unnest(standard_data, data) %>%
dplyr::select_if(~ !is.list(.)),
curves = tidyr::unnest(standard_data, curves) %>%
dplyr::select_if(~ !is.list(.)),
models = models,
omitted_standards = to_remove
)
}
#' Predict analyte concentrations from observed intensities
#'
#' @md
#' @param d Output from `assayr2::readMagPix()`
#' @param sc Output from `assayr2::fitStandardCurves()`
#' @return A `data.frame` of samples processed with predicted concentrations, QC flags and target ID.
#' @examples
#' file <- system.file("extdata", "Magpix_example.csv", package = "assayr2")
#' magpix <- read_magpix(file)
#' sc <- mp_fit_standards(magpix)
#' samples <- mp_predict_concentrations(magpix, sc)
#' @importFrom dplyr filter group_by summarise inner_join mutate select_if
#' @importFrom tidyr nest unnest
#' @importFrom tibble tibble
#' @export
mp_predict_concentrations <- function(d, sc) {
# flag samples below LLoQ
bounds <- dplyr::filter(d[[1]],
grepl("Standard",
content), intensity > 0) %>%
dplyr::group_by(target) %>%
dplyr::summarise(min_standard_intensity = min(intensity),
max_standard_intensity = max(intensity))
d[[1]] %<>% dplyr::inner_join(bounds) %>%
dplyr::mutate(below_lloq = intensity < min_standard_intensity,
above_uloq = intensity > max_standard_intensity) %>%
dplyr::mutate(flag = ifelse(below_lloq, "Below Range", "In Range")) %>%
dplyr::mutate(flag = ifelse(above_uloq, "Above Range", flag))
tmp <- tibble::tibble(sc$models)
dplyr::filter(d[[1]], grepl("Unknown", content)) %>%
as.data.frame() %>%
dplyr::group_by(target) %>%
tidyr::nest() %>%
dplyr::inner_join( tibble::tibble(target = names(sc$models), model = sc$models) ) %>%
dplyr::mutate(data = purrr::map2(data, model,
~ dplyr::mutate(.x,
pg_ml = inverse_predict(.y, .x$log10_intensity),
lloq_concentration = inverse_predict(.y, log10(.x$min_standard_intensity))))) %>%
tidyr::unnest(data) %>%
dplyr::select_if(~ !is.list(.))
}
#' High level wrapper to process raw MagPix files.
#'
#' @param infiles Character vector denoting file paths for raw data files
#' @param plate_ids Character vector of plate identifiers (since sometimes, it is not always in order)
#' @param concentrations Optional argument, list where keys are analyte identifiers and values
#' are numeric vectors indicating standard concentrations used.
#' @param wells_to_omit If not NA, a list or equal length to plate_ids with the
#' standard wells to remove. Elements are passed to `mp_fit_standards(wells_to_omit=...)`
#' individually, so an element of `NA` will result in no wells being dropped
#' from the corresponding plate.
#' @param value Name of the absorbance value to use ('Median' or 'Net MFI').
#' @return MagPix object
#' @examples
#' \dontrun{
#' file <- system.file("extdata", "Magpix_example.csv", package = "assayr2")
#' infiles <- c(file, file)
#' d <- mp_express(infiles, 1:2)
#' # If the line above fails, then we might need to pass an extra argument
#' # for concentrations we can see in the template file
#' d <- mp_express(files, 1:2, list('BMP-10' = c(3250, 1083, 361, 120, 40, 13, 4)))
#' }
#' @importFrom assertthat assert_that
#' @importFrom evaluate evaluate
#' @importFrom purrr map map_lgl map2 map2_df
#' @importFrom dplyr mutate
#' @importFrom ggplot2 ggplot geom_line geom_rug scale_x_log10 scale_y_log10 annotation_logticks facet_wrap labs theme element_text geom_text geom_point
#' @export
mp_express <- function(infiles, plate_ids, concentrations, wells_to_omit=NA,
value='Net MFI') {
if (missing(concentrations)) {
concentrations <- list()
}
# unique plate identifiers should be present for each infile
assertthat::assert_that(length(infiles) == length(plate_ids),
msg = 'Length of infiles should be the same as the number of plate identifers')
# See if the standard data in the file is correct
a <- evaluate::evaluate("temp <- purrr::map(infiles, ~ assayr::readMagPix(.))")
missing_standards <- purrr::map(a, ~ .$message) %>%
{ Filter(function(x) (!(is.null(x))), .)}
# If the standard data are incorrect and the user has not provided
# an alternative vector, then print the problematic analytes to console
if (any(purrr::map_lgl(missing_standards,
~ grepl('Invalid standard concentrations', .))) &&
length(concentrations) == 0) {
for (i in 1:length(missing_standards)) {
if (grepl('Invalid standard concentrations', missing_standards[i])) {
cat(paste0(missing_standards[[i]], '\n'))
}
}
stop("Raw standard data are missing/incorrect. Pass in concentrations as a named list,
where the keys are analytes and elements are numeric vectors.")
}
# Read the raw data if the data are correct and/or the user
# has provided vectors for the problematic analytes
d <- purrr::map(infiles, ~ read_magpix(., value=value))
# Change concentrations to those that user provides
if (length(concentrations) > 0) {
for (j in 1:length(concentrations)) {
this_analyte <- names(concentrations)[j]
for (i in 1:length(d)) {
d[[i]][[2]]$concentration[which(d[[i]][[2]]$target == this_analyte)] <-
concentrations[[this_analyte]]
}
}
}
# Fit standard curves
sc <- tryCatch({purrr::map2(d, wells_to_omit,
~ mp_fit_standards(.x, wells_to_omit = .y))},
error = function(e) {NA})
if (is.na(sc)) {
stop("Could not fit standard curves. Try passing in a concentration list for the problematic analyte(s).")
}
# Predict concentration from intensity and add plate ids
sample_data <- purrr::map2(d, sc,
~ mp_predict_concentrations(.x, .y)) %>%
purrr::map2_df(plate_ids,
~ dplyr::mutate(.x, plate = .y))
# Collapse standard curve data and data.frame with curves to plot
sc_data <- purrr::map2_df(sc, plate_ids,
~ .x$data %>% dplyr::mutate(plate = .y))
sc_curves <- purrr::map2_df(sc, plate_ids,
~ .x$curves %>% dplyr::mutate(plate = .y))
list(sample = sample_data, sc_data = sc_data, curves = sc_curves)
}
.standard_curve_plot <- function(sample_data, standard_data, standard_curves,
geom, title, subtitle) {
p <- ggplot2::ggplot(standard_data, ggplot2::aes(x=concentration, y=intensity)) +
ggplot2::geom_line(data=standard_curves, aes(y=10^log10_intensity)) +
ggplot2::geom_rug(data=sample_data, ggplot2::aes(x=pg_ml),
color='red') +
ggplot2::scale_x_log10() +
ggplot2::scale_y_log10() +
ggplot2::annotation_logticks() +
theme_assayr() +
ggplot2::facet_wrap(~ plate) +
ggplot2::labs(x="Standard Concentration (pg/ml)",
y="Intensity", title=title,subtitle=subtitle) +
ggplot2::theme(axis.text.x = element_text(size=rel(1.5)),
axis.text.y = ggplot2::element_text(size=rel(1.5)),
axis.title.y = ggplot2::element_text(size=rel(1.5)),
axis.title.x = ggplot2::element_text(size=rel(1.5)))
if (geom == 'point') {
p + ggplot2::geom_point()
} else {
p + ggplot2::geom_text(aes(label=paste(plate, well)))
}
}
#' Generate a list of standard curve plots, where each element is a plot for one analyte
#' @param magpix MagPix object generated by magPixExpress.
#' @param title Character vector for plot title.
#' @param geom Control the geom used. If "point" will be `geom_point()`, if anything
#' else will be `geom_text()`.
#' @param chunk_size Number of analytes that should be displayed in each plot
#' @importFrom dplyr filter
#' @importFrom ggplot2 facet_wrap
#' @export
mp_plot <- function(magpix, title, chunk_size, geom) {
sanalytes <- sort(unique(magpix[[1]]$target))
if (chunk_size > 1) {
sanalytes <- split(sanalytes, ceiling(seq_along(sanalytes) / chunk_size))
}
plots <- vector("list", length(sanalytes))
k <- 1
for (analyte in sanalytes) {
these_sample_data <- dplyr::filter(magpix[[1]], target %in% analyte)
these_std_data <- dplyr::filter(magpix[[2]], target %in% analyte)
these_std_curves <- dplyr::filter(magpix[[3]], target %in% analyte)
plots[[k]] <- .standard_curve_plot(these_sample_data, these_std_data, these_std_curves,
geom=geom, title = title, subtitle = analyte)
if (is.list(sanalytes)) {
plots[[k]] <- plots[[k]] + ggplot2::facet_wrap(target ~ plate) +
labs(subtitle=NULL)
}
k <- k + 1
}
plots
}
#' Get character vector of analytes measured in
#' MagPix assay
#'
#' @param x A data frame of MagPix data
#' @return A character vector of detected analytes
#' @export
analytes <- function(x) {
sort(unique(x$sample$target))
}
#' Add metadata to a MagPix object
#'
#' @md
#' @param x Magpix object
#' @param templates data.frame with columns found in device layout table
#' in addition to : `time`, `sample_id`, `plate`, and `well`.
#' `time` represents the number of days which proteins are allowed to accumulate
#' `sample_id` is an intermediate mapping from device to well
#' `plate` is an identifier for the assay plate
#' `well` is a string denoting a unique position on the assay plate
#' @return MagPix object with metadata added to the `sample` attribute
#' and a new attribute `media_blanks`, which is a data.frame of the media
#' blanks run on the plates. `time` is standardized to read as 'Day x' for
#' consistency.
#' @importFrom assertthat assert_that
#' @importFrom dplyr filter left_join
#' @export
mp_add_metadata <- function(x, templates) {
# Normalize data
assertthat::assert_that(
'time' %in% colnames(templates),
msg = '`time` must be a column in `templates`.')
templates$time %<>% gsub("day", "Day", .) %>% trimws
templates$time %<>% factor(levels = c("Day 7", "Day 10"))
assertthat::assert_that(
all(c('plate', 'well') %in% colnames(templates)),
msg = '`templates` is missing columns `plate` and `well` to join on.')
x$sample %<>% dplyr::left_join(templates)
if ('device' %in% colnames(templates)) {
x$media_blanks <- dplyr::filter(x$sample, device %in% LETTERS)
assertthat::assert_that(nrow(x$media_blanks) > 0,
msg = 'No media blanks found in templates.')
x$sample %<>% dplyr::filter(!(device %in% LETTERS))
} else {
x$media_blanks <- NA
}
x
}
#' Print method for MagPix objects
#' @param x A data frame of MagPix data
#' @export
mp_print <- function(x) {
cat(paste0(length(analytes(x)), ' unique analytes.\n'))
cat(paste0(length(unique(x$sample$plate)), ' plates.\n'))
if ('time' %in% colnames(x$sample)) {
cat(paste0(length(unique(x$sample$time)), ' unique time points.'))
} else {
cat("No metadata added. See ?addMetadata()\n")
}
}
#' A color pallete for flagging MagPix samples
#'
#' Uses `ggsci::pal_d3()`
#' @md
#' @importFrom ggsci pal_d3
#' @export
pal_magpix_flag <- function() {
x <- ggsci::pal_d3()(4)
x <- x[c(1, 2, 4)]
names(x) <- c('In Range', 'Below Range', 'Above Range')
x
}
#' Generate one or multiple dotplots for a single analyte in a MagPix assay
#'
#' @md
#' @param x Magpix object
#' @param mapping ggplot2 mapping created with `aes()`, must include x, y, and color aesthetics
#' @param analyte String denoting analyte to be plotted
#' @param facets formula for faceting
#' @param geom can be 'point' or 'text'
#' @param color Needs to be flag if color is mapped to flag
#' @param x_text_angle Angle of text on the x axis
#' @param title Title of the plot
#' @param label_pos Where should the labels go?
#' @param data_geom_size Passed onto `geom_point(size = ...)`
#' @param label_geom_size Passed onto `geom_label(size = ...)`
#' @param log_y_axis Boolean, should the y-axis be logged via `+ scale_y_log10()`?
#' @param chunk_size Integer, the number of analytes plotted in each plot.
#' @param text_label Column in `x` for mapping to `aes(label= ...)`.
#' @name mp_dotplot
#' @importFrom dplyr filter mutate
#' @importFrom ggplot2 ggplot scale_color_manual theme scale_shape_manual scale_shape_discrete stat_summary labs geom_line geom_point aes position_dodge
#' @importFrom ggsci scale_color_d3
#' @export
mp_dotplot <- function(x, mapping, analyte, facets, geom='point',
color='flag', x_text_angle=0,
title='', label_pos='top',
data_geom_size=3.5, label_geom_size=3.5,
log_y_axis=F, text_label='device') {
d <- dplyr::filter(x$sample, target %in% analyte)
# Make paste of device and plate for plotting technical replicates
d %<>% dplyr::mutate(device_plate = paste0(device, plate))
# Color plots by QC flags
if (color == 'flag') {
p <- ggplot2::ggplot(d, mapping) +
ggplot2::scale_color_manual(name=NULL, values=pal_magpix_flag()) +
ggplot2::theme(legend.position = 'top')
} else if ('color' %in% names(mapping)) { # Other color scheme
print('color in aes')
p <- ggplot2::ggplot(d, mapping) +
ggsci::scale_color_d3(name=NULL, 'category20')
} else {
p <- ggplot2::ggplot(d, mapping)
}
if (x_text_angle > 0) {
p <- p + ggplot2::theme(axis.text.x = ggplot2::element_text(angle=x_text_angle, vjust=1, hjust=1))
}
# Add custom shape pal
if ('shape' %in% names(mapping)) {
if (mapping['shape'] == 'flag') {
p <- p + ggplot2::scale_shape_manual(name=NULL,
values = c(`Above Range` = 2,
`In Range` = 1,
`Below Range` = 3))
} else {
p <- p + ggplot2::scale_shape_discrete(name=NULL)
}
}
# Facet by user-specified formula
if (!(missing(facets))) {
p <- p + facets
}
if (length(analyte) == 1) {
# Add label based on length of `analyte`
p <- p + ggplot2::labs(y= paste0('[', unique(d$target), '] (pg/ml)'),
x=NULL, subtitle=unique(d$target), title=title)
} else {
p <- p + ggplot2::labs(y='[Analyte] (pg/ml)', x=NULL,
title=title)
}
# Add geom for data
if (geom == 'point') {
p <- p +
ggplot2::geom_point(ggplot2::aes(group=device_plate), size=data_geom_size,
position = ggplot2::position_dodge(width=0.5)) +
ggplot2::geom_line(ggplot2::aes(group=device_plate),
position = ggplot2::position_dodge(width=0.5))
} else if (geom == 'text') {
if (text_label == 'device') {
p <- p +
ggplot2::stat_summary(fun.y = mean, ggplot2::aes(group=device_plate, label=device),
geom='text', position = ggplot2::position_dodge(width=0.3))
}
else if (text_label == 'plate_well') {
p <- p +
ggplot2::geom_text(
ggplot2::aes(group=device_plate, label=paste(plate, "-", well)),
position = ggplot2::position_dodge(width=0.3))
}
else {
stop("if geom == 'text', then text_label must be 'device' or 'plate_well'")
}
}
# Add confidence intervals and p-values if they have been computed
if ('stats' %in% attributes(x)$names) {
# only draw stats if they exist and at least one CI exists
s <- dplyr::filter(x$stats, target %in% analyte)
v <- dplyr::filter(x$vehicle_means, target %in% analyte)
if (label_pos == 'top') {
y_val <- Inf
vjust_val <- 1.2
} else {
y_val <- -Inf
vjust_val <- -0.4
}
if (nrow(v) > 0) {
p <- p +
ggplot2::geom_crossbar(data=s, ggplot2::aes(ymin=lwr, y=est, ymax=upr,
color=NULL)) +
ggplot2::geom_hline(data=v, ggplot2::aes(yintercept=est, color=NULL),
linetype=3) +
ggplot2::geom_label(data=s, ggplot2::aes(label = label, y=y_val,
color=NULL),
alpha=0.6, vjust=vjust_val,
size = label_geom_size) +
ggplot2::scale_y_continuous(expand = c(a=0.05, b=0.05, c=0.2, d=0.05))
}
}
# Log y-axis
if (log_y_axis) {
p <- p +
ggplot2::scale_y_log10(expand = c(a=0.05, b=0.05, c=0.2, d=0.05)) +
ggplot2::annotation_logticks(sides='l')
}
p
}
#'@rdname mp_dotplot
#'@export
mp_dotplots <- function(x, mapping, facets, geom='point', color='flag',
x_text_angle=0, title='', label_pos='top', chunk_size=1,
data_geom_size = 3.5, label_geom_size = 3.5,
log_y_axis=F, text_label='device') {
plots <- list()
if (chunk_size == 1) {
for (i in 1:length(analytes(x))) {
plots[[i]] <- mp_dotplot(x, mapping, analytes(x)[i],
facets, geom, color,
x_text_angle, title, label_pos,
data_geom_size, label_geom_size,
log_y_axis=F, text_label)
}
} else {
chunks <- split(analytes(x), ceiling(seq_along(analytes(x))/chunk_size))
for (i in 1:length(chunks)) {
plots[[i]] <- mp_dotplot(x, mapping, chunks[[i]],
facets, geom, color,
x_text_angle, title, label_pos,
data_geom_size, label_geom_size,
log_y_axis, text_label)
}
}
plots
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.