Introduction to cytominer

Load data

First, load the data, which is stored in a database backend created using https://github.com/cytomining/cytominer-database.

fixture <-
  system.file("extdata", "fixture_intensities_shapes.sqlite",
              package = "cytominer")

db <- DBI::dbConnect(RSQLite::SQLite(), fixture)

Then load associated metadata if it exists, and copy it to the backend so that we can use it later.

ext_metadata <-
  readr::read_csv(system.file("extdata", "metadata.csv",
                              package = "cytominer")) %>%
  dplyr::rename(g_well = Well)

## Parsed with column specification:
## cols(
##   Well = col_character(),
##   Type = col_character()
## )

ext_metadata <- dplyr::copy_to(db, ext_metadata)

Next, select a measurement table that you want to work. Here, we will pick intensities but we can easily extend to using multiple or all measurement classes if needed by creating new views.

intensities <-
  dplyr::tbl(src = db, "view_intensities") %>%
  dplyr::compute()

For this example, lets filter the data down to a few wells.

measurements <-
  intensities %>%
  dplyr::filter(g_well %in% c("A01", "A02", "A10", "A11"))

How many rows does this table have?

measurements %>%
  dplyr::tally() %>%
  knitr::kable()

That actually 9 times the number of cells in this experiment (n = 40): each compartment and each channel gets its own row. Here, we have 3 compartments (cell, nucleus, cytoplasm) and 3 channels (CellMask, Hoechst, Alexa568). So that’s 3 x 3 x 40 = 360.

Next, do some setup stuff that we will need later

qualities <- c("q_debris")

groupings <-
  c("g_plate",
    "g_well",
    "g_image",
    "g_pattern",
    "g_channel")

variables <-
  colnames(measurements) %>%
  stringr::str_subset("^m_")

measurements %<>%
  dplyr::select(dplyr::one_of(c(groupings, qualities, variables)))

Clean

Let’s remove cells that come from images that were marked as having debris

debris_removed <-
  measurements %>% dplyr::filter(q_debris == 0)

Then, remove cells where all the measurements are NA’s - this may happen if the identified cell mask was too small to measure any of the features.

na_rows_removed <-
  cytominer::drop_na_rows(
    population = debris_removed,
    variables = variables
  ) %>%
  dplyr::compute()

Normalize

We need to normalize the data so that

• features are on the same scale

• plate-to-plate variation is reduced

The default for doing this is standardization. Here, we take all the cells from control wells in the experiment (this is where the external metadata gets used) and compute normalizations parameters from that (in this case, just the mean and s.d.) and then apply it to the whole dataset (i.e. the population)

normalized <-
  cytominer::normalize(
    population = na_rows_removed %>% 
      dplyr::collect(),
    variables = variables,
    strata =  c("g_plate", "g_pattern", "g_channel"),
    sample =
      na_rows_removed %>%
      dplyr::inner_join(
        ext_metadata %>% 
          dplyr::filter(Type == "ctrl") %>% 
          dplyr::select(g_well) 
      ) %>% dplyr::collect()
  )

## Joining, by = "g_well"

normalized %<>% dplyr::collect()

In some cases, we may have features that have no variance at all (e.g. Euler number). If these features have not already been removed by this stage, the standardization step will results in all values for that feature being NA ( because s.d. = 0). Lets remove them:

First, count how many cells have NA values per feature:

na_frequency <-
  cytominer::count_na_rows(
    population = normalized,
    variables = variables)

na_frequency %>%
  tidyr::gather(feature, na_count) %>%
  knitr::kable()

As it turns out, no feature has NA in this example. But lets run this cleaning operation anyway (no features will be dropped)

cleaned <-
  cytominer::variable_select(
    population = normalized,
    variables = variables,
    operation = "drop_na_columns"
)

variables <-
  colnames(cleaned) %>%
  stringr::str_subset("^m_")

Transform

Tranform the data so that assumptions we may later make about the data distribution are satisfied (e.g. Gaussianity). The default here is generalized_log.

transformed <-
  cytominer::transform(
    population = cleaned,
    variables = variables
  )

Select features

Finally, we typically perform feature selection on the data. Feature selection is an expensive operation, so we usually want to train the feature selection model on a sample of the dataset. Here, we choose to aggregate the data instead of sampling it (i.e. collapse it to per-well aggregates)

aggregated <-
  cytominer::aggregate(
    population = transformed,
    variables = variables,
    strata = groupings
  ) %>%
  dplyr::collect()

variables <-
  colnames(aggregated) %>%
  stringr::str_subset("^m_")

… and then apply feature selection on the per-cell data. Here correlation_threshold - a method that reduces the redundancy of features - is used.

selected <-
  cytominer::variable_select(
    population = transformed,
    variables = variables,
    sample = aggregated,
    operation = "correlation_threshold"
  ) %>%
  dplyr::collect()

And now lets take a glimpse at the data!

selected %>%
  dplyr::glimpse()

## Rows: 315
## Columns: 10
## $ g_plate                                 <chr> "110000106828", "110000106828…
## $ g_well                                  <chr> "A01", "A01", "A01", "A01", "…
## $ g_image                                 <chr> "5a822385194aceb347a07f27285d…
## $ g_pattern                               <chr> "Cells", "Cells", "Cells", "C…
## $ g_channel                               <chr> "Alexa568", "Alexa568", "Alex…
## $ q_debris                                <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ m_intensities_integrated                <dbl> -1.48831626, -0.06761634, -0.…
## $ m_intensities_maximum                   <dbl> -1.25783096, -1.51483166, -1.…
## $ m_intensities_median_absolute_deviation <dbl> -1.4847843, -0.5372476, -1.10…
## $ m_intensities_minimum                   <dbl> -2.3552231, -1.9025552, -0.10…

  DBI::dbDisconnect(db)