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)

suppressPackageStartupMessages(library(dplyr))
library(formattable)
library(phenoptr)

export_path = params$export_path
export_data = params$export_data

# Quantiles to show on the plots
quantiles = sort(params$quantiles)

This report shows component levels for samples from

r export_path

comp_files = list.files(export_path, pattern='component_data.tif', full.names=TRUE) %>% 
  purrr::set_names(., basename(.) %>% stringr::str_remove('_component_data.tif'))

if (length(comp_files) == 0)
  stop('No component data files found at "', export_path, '"')

Pixel intensity by image

These plots show the number of pixels at each intensity for each component of each image. Intensity is shown on a log scale.

# Initial data acquisition and creation of the first round of plots
# all_data is a list of lists of plot, data, quants, clipping
all_data = purrr::imap(comp_files, function(path, name) {
                      message('Reading ', basename(path))
                      phenoptrReports:::get_component_level_data(path, name, 
                      quantiles=quantiles)
})

# Save some data
all_quants = purrr::map_dfr(all_data, 'quants')
readr::write_csv(all_quants, file.path(export_path, 'all_quants.csv'))

all_clipping = purrr::map_dfr(all_data, 'clipping')
readr::write_csv(all_clipping, file.path(export_path, 'all_clipping.csv'))

all_components = purrr::map_dfr(all_data, 'data')
density_fig_height =  1 + 0.6 * n_distinct(all_components$Fluor)

# Output density plots by file
message('Creating histograms by source')
purrr::iwalk(all_data, function(d, name) {
  message(name)
  print(phenoptrReports:::component_ridge_plot(d$data, d$quants, d$clipping, name))
  cat('\n\n  \n----\n  \n')
})

component_fig_height = 1 + 0.5 * n_distinct(all_components$Source)

Pixel intensity by component

message('Creating histograms by component')
fluors = sort(unique(all_components$Fluor))
num_fluors = length(fluors)
palette = phenoptrReports::discrete_colors(num_fluors)

purrr::walk2(fluors, palette, 
   function(fluor, fill) {
     message(fluor)
     print(phenoptrReports:::fluor_ridge_plot(
       all_components %>% filter(Fluor==fluor),
       all_quants %>% filter(Fluor==fluor), 
       all_clipping %>% filter(Fluor==fluor),
       fluor, fill))
    cat('\n\n  \n----\n  \n')
})

corr_plot_size = num_fluors/2

Correlation between markers

These heatmaps show the Pearson's correlation between per-pixel expression values for each pair of markers. This may help identify crosstalk between markers.

Clustered

The first set of heatmaps group correlated components together

message('Creating clustered correlation heatmaps')
purrr::iwalk(all_data, function(d, name) {
  message(name)
  corr = d$data %>% 
    dplyr::select(-Source) %>% 
    tidyr::pivot_wider(names_from='Fluor', values_from='value', values_fn=list) %>% 
    tidyr::unnest(cols=everything()) %>% 
    stats::cor()

  p = pheatmap::pheatmap(corr, color = rev(RColorBrewer::brewer.pal(n = 7, name = 'RdYlBu')), 
           breaks=seq(-1, 1, length.out=8),
           treeheight_row=0,
           clustering_distance_rows='correlation',
           clustering_distance_cols='correlation',
           main=name)
  print(p)
  cat('\n\n  \n----\n  \n')
})

Natural order

The second set of heatmaps shows components in the order in which they appear in the source files.

message('Creating ordered correlation heatmaps')
purrr::iwalk(all_data, function(d, name) {
  message(name)
  corr = d$data %>% 
    dplyr::select(-Source) %>% 
    tidyr::pivot_wider(names_from='Fluor', values_from='value', values_fn=list) %>% 
    tidyr::unnest(cols=everything()) %>% 
    stats::cor()

  p = pheatmap::pheatmap(corr, color = rev(RColorBrewer::brewer.pal(n = 7, name = 'RdYlBu')), 
           breaks=seq(-1, 1, length.out=8),
           cluster_rows=FALSE,
           cluster_cols=FALSE,
           main=name)
  print(p)
  cat('\n\n  \n----\n  \n')
})

Pair plots

Pair plots show the expression of each marker against each other marker. They may help identify crosstalk between markers.

Pair plots are shown only for images with ten or fewer components.

message('Creating pair plots')
purrr::iwalk(all_data, function(d, name) {
  # Unstack and subsample
  message(name)
  d = d$data %>% select(-Source) %>% 
    filter(!Fluor %in% c('Sample AF', 'AF', 'Autofluorescence')) %>% 
    unstack(value ~ Fluor) %>% 
    sample_n(2000)

  if (export_data) {
    # Write pairs plot data to a file
    pairs_name = paste0(name, '_pairs_data.txt')
    pairs_path = file.path(export_path, pairs_name)
    readr::write_tsv(d, pairs_path)
  }

  # Only show pairs plots if the number of fluors is <= 10, they are not
  # helpful with many parts
  if (ncol(d) <= 10) {
    axis_max = max(d)
    pairs(d, pch='.', cex=3, col=scales::alpha('black', 0.2),
                main=name, xlim=c(0, axis_max), ylim=c(0, axis_max))
    cat('\n\n  \n----\n  \n')
  }
})

# If there are at least two quantiles, show a signal-to-noise table
# with the top and bottom quantiles.
if (length(quantiles) >= 2) {
  cat('\n\n### Signal to noise\n\n')
  cat(stringr::str_glue(
  'This table shows the ratio of the ',
  '{scales::percent(quantiles[length(quantiles)], accuracy=0.1)}ile signal ',
  'level to the {scales::percent(quantiles[1])}ile signal level.\n\n'))

  # What columns of all_quants do we want?
  qlow_col = 3
  qhi_col = qlow_col + length(quantiles) - 1
  all_quants %>% dplyr::mutate(SN=round(.[[qhi_col]]/.[[qlow_col]], 1)) %>% 
    dplyr::select(Fluor, Source, SN) %>% 
    tidyr::spread(Fluor, SN) %>% 
    knitr::kable()
}

Guidance

Pixel intensity by image: The top 99.9^th^ percentile pixel is usually reflective of 'positive' signal. If your multiplex panel is well-balanced, the 99.9^th^ percentile values should be similar between channels.
Pixel intensity by component: Allows you to quickly assess to see if there are any outliers or large-scale differences in staining intensity across the image set.
Pair plots: Tight correlations may indicate crosstalk between channels.

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