In akoyabio/phenoptrReports: Create reports using Phenoptics data

knitr::opts_chunk$set(echo=FALSE,fig.width=9, fig.height=6, 
                      comment=NA, warning=FALSE, message=FALSE)

# The default plot hook for vignettes creates a markdown image 'tag'
# which breaks in pandoc if there are any spaces in the image path.
# The standard HTML plot hook inserts an <img> tag which renders correctly.
knitr::knit_hooks$set(plot = knitr:::hook_plot_html)

# Need that pipe, don't want all of magrittr
suppressPackageStartupMessages(library(dplyr))
library(ggplot2)

Unmixing quality report for singly-stained samples from

r params$export_path

n_bright = 3200 # Number of bright pixels to find
skip_fraction = 0.0001

export_path = params$export_path
comp_files = list.files(export_path, pattern='component_data.tif', full.names=TRUE)

# Put DAPI first
dapi_files = stringr::str_detect(basename(comp_files), 'DAPI')
comp_files = c(comp_files[dapi_files], comp_files[!dapi_files])
comp_names = basename(comp_files)
primary_fluors = purrr::map_chr(comp_names, phenoptrReports:::file_to_fluor)
all_data = purrr::map2(comp_files, primary_fluors, 
                       phenoptrReports:::process_singleplex_image, 
                       n_bright=n_bright, skip_fraction=skip_fraction) %>% 
    purrr::set_names(comp_names)

# If the images do not all have the same components, give a warning
# Get all component names (and a few extra names)
image_components = purrr::map(all_data, 'signals') %>% 
  purrr::map(names)
all_component_names = unlist(image_components) %>% unique()

# Check for any missing
missing = purrr::map_lgl(image_components, 
                         ~length(setdiff(all_component_names, .x))>0)
if (any(missing)) {
  cat('\n\n### Notice - missing components! {.warning}\n\n')
  cat('Not all supplied component files have the same components.
      Missing components cause NA values in the tables below.\n\n')

  purrr::walk2(names(image_components), image_components,
    function(name, components) {
      # Remove the extra names
      components = head(components, -2)
      components = paste(components, collapse=', ')
      cat(stringr::str_glue('- <strong>{name}</strong>: {components}\n\n'))
    })

}

Crosstalk by image

This table shows crosstalk from each primary fluor to the other components in the same image. This is a measure of unmixing error.

message('Creating report')

crosstalk = all_data %>%
  purrr::map('sn') %>% dplyr::bind_rows() %>% 
  dplyr::select(Primary, Source, dplyr::matches('DAPI'), dplyr::everything())

to_display = crosstalk %>% 
  dplyr::select(-Primary, -dplyr::matches('Autofluorescence')) %>% 
  dplyr::mutate_at(-1, formattable::percent, digits=1)

# Put DAPI first
to_display = to_display %>% 
  dplyr::select(Source, dplyr::matches('DAPI'), dplyr::everything())

my_gradient = function(colors) {
  grad = colorRamp(colors, space='Lab')
  function(x) {
    color = grad(x)
    colnames(color)= c('red', 'green', 'blue')
    formattable::csscolor(round(t(color), 0), format='rgb')
  }
}

# Overall maximum
limit = max({x<-unlist(to_display[,-1]); x[x<1]})

# Colors for the green-yellow-red/orange ramp
# Excel colors from Carla
# green = '#00B050'
# yellow = '#FFEB84'
# red_orange = '#ED7D31'

# These are a bit more saturated
green = '#08b151'
yellow = 'gold'
red_orange = 'red'

my_color = function(x) {
  cutpt = 0.015
  green_ramp = my_gradient(c(green, yellow))
  yellow_ramp = my_gradient(c(yellow, red_orange))
  ifelse(x==1, formattable::csscolor('white'), # Don't show 100%
         ifelse(x<=cutpt, green_ramp(x/cutpt),
                yellow_ramp((x-cutpt)/(limit-cutpt))))
}

my_font_weight= function(x) ifelse(x>0.005 & x<1, 'bold', 'normal')

# Make a color_tile variation
my_color_tile = function (...) 
{
    formattable::formatter("span", 
                           style = function(x) 
                             formattable::style(display = "block", 
        padding = "3px", `border-radius` = "0px",
        `background-color` = my_color(x), 
        color = ifelse(x==1, 'white', 'black'),
        `font-weight`=my_font_weight(x)))
}

padded_cell = function (...) 
{
    formattable::formatter("span", 
                           style = function(x) 
                             formattable::style(display = "block", 
        padding = "3px"))
}

formats = purrr::map(names(to_display)[-1], my_color_tile) %>% 
  purrr::set_names(names(to_display)[-1]) %>% 
  c(list(Source=padded_cell('Source')))

formattable::formattable(to_display, formats, 
                         align=c('l', rep('r', ncol(to_display)-1)),
            table.attr = "class=\"table table-condensed fill-cell\"")

Crosstalk by component

This table estimates crosstalk from other fluors to the primary fluor as a percent of the primary fluor. This estimates the effect of unmixing error on other signals.

signals = all_data %>%
  purrr::map('signals') %>% dplyr::bind_rows() %>% 
  dplyr::select(Primary, Source, dplyr::matches('DAPI'), dplyr::everything())

to_display = signals %>% 
  dplyr::select(-Primary, -dplyr::matches('Autofluorescence')) %>% 
  dplyr::mutate_at(-1, ~formattable::percent(.x/max(.x), digits=1))

# Put DAPI first
to_display = to_display %>% 
  dplyr::select(Source, dplyr::matches('DAPI'), dplyr::everything())

# We need to know the smallest relative value in own fluor
min_own_value = to_display %>% dplyr::mutate(Fluor=primary_fluors) %>% 
  dplyr::select(-Source) %>% 
  tidyr::gather('Component', 'Value', -Fluor) %>% 
  dplyr::filter(stringr::str_detect(Component, Fluor)) %>% 
  dplyr::pull(Value) %>% 
  min()

limit = max({x<-unlist(to_display[,-1]); x[x<min_own_value]})

formats = purrr::map(names(to_display)[-1], my_color_tile) %>% 
  purrr::set_names(names(to_display)[-1]) %>% 
  c(list(Source=padded_cell('Source')))

formattable::formattable(to_display, formats, 
                         align=c('l', rep('r', ncol(to_display)-1)),
            table.attr = "class=\"table table-condensed fill-cell\"")

Guidance

Crosstalk values for a given component	Typical performance of component for different tasks
Crosstalk values for a given component	Expression level assessment	Cell lineage/phenotype classification
All values in a column are < 2%	Good	Good
All values in a column are < 5% Some values are > 2%	Marginal, depending on dynamic range of expression	Good
All values in a column are < 10% Some values are > 5%	Poor	Ok, may require more extensive classifier training
Any crosstalk values > 10%	Poor	Marginal, see Opal Assay Development Guide for advice

This table suggests guidelines to help interpret the Crosstalk by Component table. Components used to assess expression levels have stricter requirements than components used for phenotyping.

Pixel intensity in primary component

This plot shows the intensity distribution of pixels of the primary fluor in each image. Intensity is shown on a log scale.

density_fig_height = 0.5 + 0.5 * length(comp_files)

# Make density plots using random sampled data
all_primaries = tibble::tibble(
  Source=factor(comp_names, levels=comp_names), 
  Fluor=factor(primary_fluors, levels=unique(primary_fluors)),
  data=purrr::map(all_data, 
                  ~tibble::tibble(primary=as.numeric(.x$primary)) %>%
                     dplyr::sample_frac(size=0.1))) %>% 
  tidyr::unnest(cols=c(data))

palette = discrete_colors(dplyr::n_distinct(all_primaries$Fluor))

if (require(ggridges)){
  all_primaries %>% 
    dplyr::filter(primary > 0) %>% 
    ggplot(aes(primary, forcats::fct_rev(Source), fill=Fluor)) +
    geom_density_ridges(color='gray50', na.rm=TRUE) +
    labs(x='Pixel intensity', y='') +
    scale_fill_manual(values=palette) +
    guides(color='none', fill='none') +
    theme_minimal() +
    theme(axis.text.y=element_text(face='bold', hjust=0, vjust=0)) +
    scale_x_log10(expand = c(0.01, 0), limits=c(0.1, NA)) +
    scale_y_discrete(expand = c(0.01, 0))
} else {
  # ggplot version
  ggplot(all_primaries, aes(primary, color=Fluor, fill=Fluor)) +
           geom_density(na.rm=TRUE) +
    facet_wrap(~Source, ncol=1, strip.position='left', scales='free_y') +
    labs(x='Pixel intensity', y='') +
    scale_x_log10(limits=c(0.1, NA)) +
    scale_color_manual(values=palette) +
    scale_fill_manual(values=palette) +
    guides(color='none', fill='none') +
    theme_minimal() +
    theme(strip.text.y=element_text(face='bold', angle=180, hjust=0),
          axis.ticks.y=element_blank(),
          axis.text.y=element_blank(),
          panel.grid.major.y=element_blank(),
          panel.grid.minor=element_blank(),)
}

Signal strength

This table shows signal strength for each component, measured at the brightest pixels of the primary fluor. This is the raw data underlying the crosstalk tables.

to_display = signals %>% 
  dplyr::select(-Primary, -dplyr::matches('Autofluorescence')) %>% 
  dplyr::mutate_at(-1, round, digits=2)

# Put DAPI first
to_display = to_display %>% 
  dplyr::select(Source, dplyr::matches('DAPI'), dplyr::everything())

knitr::kable(to_display)

What is this?

Image analysis

For each image, the brightest r scales::comma(n_bright) pixels in the primary component are found (ignoring the r scales::percent(skip_fraction) very brightest).

For these pixels, the mean signal level in each component is computed. The values are shown in the Signal Strength table.

Crosstalk by image

This table shows the relative signal level of each fluor in the components where the fluor should not express.

The ratio of non-signal component strength to signal component strength within a single image is reported as a percent. This is a measure of unmixing error.

Each column in the report shows crosstalk into the component named at the top of the column. Crosstalk comes from the fluors represented by the files listed in the first column.

Crosstalk by component

This table estimates the amount of crosstalk from non-signal fluors to each component.

For each component, it shows the crosstalk from other fluors as a percent of the mean signal for the expressing component. This estimates the effect of unmixing error on other signals.

If multiple samples are available for the expressing component, the largest mean is used.

Pixel intensity

These plots show the relative numbers of pixels at each pixel intensity. The Y-axis varies between plots; the X-axis is the log base 10 of pixel intensity.

Signal strength

The signal strength table shows the mean expression of each fluor in the bright pixels. This is the raw data underlying the crosstalk tables. Crosstalk by image normalizes this data by row; crosstalk by component normalizes by column.

{height=50px style='border:none;'}