# beadplexr with CBA and MACSPlex assays

#### Ulrik Stervbo

#### 2018-03-06

# Introduction

This document demonstrates the potential usefulness of beadplexr with data generated with the CBA assay system from BD, and the MACSPlex assay system from Miltenyi Biotec. The package has been tested with assay data from the two systems. Unfortunately, I cannot share the data - for demonstration purposes I have simulated data for the two systems. The simulated data is deliberately noisy, and do not reflect the true quality of the respective assays. The sole purpose of the simulated data is to mimic the two dimensional bead identification.

# Load and view the data

    library(beadplexr)
    library(dplyr)

    ##
    ## Attaching package: 'dplyr'

    ## The following objects are masked from 'package:stats':
    ##
    ##     filter, lag

    ## The following objects are masked from 'package:base':
    ##
    ##     intersect, setdiff, setequal, union

    library(ggplot2)
    data(simplex)

## The MACSPlex-like data

The data from the CBA system is similar to the data from the MACSPlex system: there is a single population in the forward-side scatter and 10 populations separated by brightness of FITC and PE. The concentration of the analytes are given in the APC channel.

    mplex_data <- simplex[["mplex"]]

    mplex_data %>%
    facs_plot(.x = "FSC", .y = "SSC", .type = "hex")

    mplex_data %>%
    facs_plot(.x = "FITC", .y = "PE", .type = "hex")

## The CBA-like data

There is a single bead population in the forward-side scatter and up to 30 bead populations separated by the brightness of APC and APC-Cy7. The concentration of the analytes are given in the PS channel.

    cba_data <- simplex[["cba"]]

    cba_data %>%
    facs_plot(.x = "FSC", .y = "SSC", .type = "hex")

    cba_data %>%
    facs_plot(.x = "APC", .y = "APC-Cy7", .type = "hex")

# Analyte identification

Since there is just a single population in the forward-side scatter we can jump straight to identification of the analytes. It is of course possible to decrease spurious analyte signals by focusing on true bead events in the forward-side scatter, and it might be necessary in a true experiment.

## Analyte identification

### MACSPlex-like

We use the `identify_analyte()` to annotate the correct analytes in the MACSPlex-like data. The function requires the parameter `.analyte_id`, which – in a true experiment – correspond to some identifier given by Miltenyi Biotec. Here we just assign a number.

The function also allows for removing events that are far from the centre of the analyte clusters. This is done by the parameter `.trim`. Here we set the value to `0.1` which means that 10% of the events are removed from the clusters. Excluded events are given the analyte ID `NA` and can easily be excluded from further processing.

    mplex_analyte <- mplex_data %>%
      identify_analyte(.parameter = c("FITC", "PE"), .analyte_id = as.character(c(1:10)))

    mplex_analyte %>%
      facs_plot(.x = "FITC", .y = "PE", .beads = "analyte")

    mplex_analyte <- mplex_data %>%
      identify_analyte(.parameter = c("FITC", "PE"), .analyte_id = as.character(c(1:10)), .trim = 0.1)

    mplex_analyte %>%
      facs_plot(.x = "FITC", .y = "PE", .beads = "analyte")

We can now look at the intensity of the APC channel, to find the relative amount of analyte bound to each bead. In this case it makes very little sense, but we’ll do it for the purpose of completeness

    mplex_analyte %>%
      filter(!is.na(analyte)) %>%
      mutate(analyte = factor(analyte, levels = c(1:10))) %>%
      ggplot() +
      aes(x = APC, fill = analyte) +
      geom_density()+
      facet_wrap(~ analyte)+
      theme(legend.position = "none")

The trimming can also be done on the identified analytes. Here we remove 5% of the events.

    mplex_analyte %>%
      # filter(!is.na(analyte)) %>%
      group_by(analyte) %>%
      trim_population(.parameter = "APC", .column_name = "analyte", .trim = 0.05) %>%
      ungroup() %>%
      mutate(analyte = factor(analyte, levels = c(1:10))) %>%
      ggplot() +
      aes(x = APC, fill = analyte) +
      geom_density() +
      facet_wrap(~ analyte) +
      theme(legend.position = "none")

### CBA-like

For the CBA assay system, the approach is quite similar to the MACSPlex:

    cba_analyte <- cba_data %>%
      identify_analyte(.parameter = c("APC", "APC-Cy7"), .analyte_id = as.character(c(1:30)), .trim = 0.1)

    cba_analyte %>%
      facs_plot(.x = "APC", .y = "APC-Cy7", .beads = "analyte")

## ID assignment

The analyte ID is assigned to each cluster, based on the order of the centre of the cluster. For analytes where the bead is identified by two fluorescent parameters it could be a little tricky, because the centres are first sorted by the first parameter given, and then the second. The order of parameters and analyte IDs matter a lot! I urge you to double check that the beads are identified as expected. Once you have everything sorted out, it should remain stable for your cytometer.

For the simulated MACSPlex-like data the IDs are assigned as below.

    mplex_analyte %>%
      filter(!is.na(analyte)) %>%
      group_by(analyte) %>%
      summarise(`FITC mean` = mean(FITC), `PE mean` = mean(PE)) %>%
      arrange(`FITC mean`, `PE mean`) %>%
      knitr::kable()

analyte | FITC mean | PE mean
--- | --- | ---
1 | 5.253017 | 5.343757
2 | 5.342780 | 12.427284
3 | 6.408917 | 7.362167
4 | 7.709976 | 10.075022
5 | 8.494268 | 14.271564
6 | 8.827282 | 7.689722
7 | 12.166583 | 6.347947
8 | 12.292569 | 9.710517
9 | 14.361117 | 9.792471
10 | 14.880174 | 14.453440

## Noise reduction

For the CBA-like data, some clusters are still a bit too wide, while the trimming of `0.1` made others look quite fine. The simulated data contains a bit of noise hat might interfere the the cluster identification. This noise might not be present in the a true experiment, but is included to demonstrate the removal of lonely, noisy events by the `despeckel()` functionality of beadplexr.

    cba_analyte <- cba_data %>%
      despeckle(.parameter = c("APC", "APC-Cy7"), .neighbours = 2) %>%
      identify_analyte(.parameter = c("APC", "APC-Cy7"), .analyte_id = as.character(c(1:30)), .trim = 0.01)

    ## Loading required namespace: igraph

    cba_analyte %>%
      facs_plot(.x = "APC", .y = "APC-Cy7", .beads = "analyte")

In this particular example, the `despeckle()` is a little too rough on the population in the lower left corner, and it is probably more prudent to accept a little noise on one analyte than decimation of another.

# Session info

    ## R version 3.3.1 (2016-06-21)
    ## Platform: x86_64-pc-linux-gnu (64-bit)
    ## Running under: Ubuntu 15.10
    ##
    ## locale:
    ##  [1] LC_CTYPE=de_DE.UTF-8       LC_NUMERIC=C
    ##  [3] LC_TIME=de_DE.UTF-8        LC_COLLATE=de_DE.UTF-8
    ##  [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=de_DE.UTF-8
    ##  [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C
    ##  [9] LC_ADDRESS=C               LC_TELEPHONE=C
    ## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C
    ##
    ## attached base packages:
    ## [1] stats     graphics  grDevices utils     datasets  methods   base
    ##
    ## other attached packages:
    ## [1] bindrcpp_0.2    hexbin_1.27.1   ggplot2_2.2.1   dplyr_0.7.4
    ## [5] beadplexr_0.1.0
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] nlme_3.1-131        matrixStats_0.52.2  pbkrtest_0.4-7
    ##  [4] devtools_1.13.0     rprojroot_1.3-2     prabclus_2.2-6
    ##  [7] hunspell_2.4        tools_3.3.1         backports_1.1.2
    ## [10] R6_2.2.2            lazyeval_0.2.0      BiocGenerics_0.20.0
    ## [13] mgcv_1.8-17         colorspace_1.3-2    trimcluster_0.1-2
    ## [16] nnet_7.3-12         raster_2.5-8        withr_2.1.1
    ## [19] sp_1.2-4            curl_3.1            git2r_0.18.0
    ## [22] graph_1.52.0        quantreg_5.33       Biobase_2.34.0
    ## [25] SparseM_1.77        xml2_1.1.1          desc_1.1.1
    ## [28] sandwich_2.3-4      labeling_0.3        diptest_0.75-7
    ## [31] scales_0.4.1        DEoptimR_1.0-8      mvtnorm_1.0-6
    ## [34] robustbase_0.92-7   commonmark_1.2      stringr_1.2.0
    ## [37] digest_0.6.12       minqa_1.2.4         rmarkdown_1.5
    ## [40] rrcov_1.4-3         pkgconfig_2.0.1     htmltools_0.3.6
    ## [43] lme4_1.1-13         plotrix_3.6-5       highr_0.6
    ## [46] rlang_0.2.0         rstudioapi_0.6      flowCore_1.40.6
    ## [49] bindr_0.1           zoo_1.8-0           mclust_5.2.3
    ## [52] gtools_3.5.0        car_2.1-4           magrittr_1.5
    ## [55] modeltools_0.2-21   Matrix_1.2-8        Rcpp_0.12.10
    ## [58] munsell_0.4.3       stringi_1.1.5       multcomp_1.4-6
    ## [61] yaml_2.1.14         MASS_7.3-45         flexmix_2.3-14
    ## [64] plyr_1.8.4          grid_3.3.1          parallel_3.3.1
    ## [67] crayon_1.3.4        lattice_0.20-34     splines_3.3.1
    ## [70] knitr_1.15.1        pillar_1.2.1        igraph_1.0.1
    ## [73] fpc_2.1-10          corpcor_1.6.9       codetools_0.2-15
    ## [76] stats4_3.3.1        glue_1.2.0          drc_3.0-1
    ## [79] packrat_0.4.8-1     evaluate_0.10       nloptr_1.0.4
    ## [82] testthat_1.0.2      MatrixModels_0.4-1  gtable_0.2.0
    ## [85] purrr_0.2.2         tidyr_0.6.2         kernlab_0.9-25
    ## [88] assertthat_0.2.0    roxygen2_6.0.1      class_7.3-14
    ## [91] survival_2.41-2     pcaPP_1.9-61        tibble_1.4.2
    ## [94] memoise_1.1.0       rversions_1.0.3     cluster_2.0.6
    ## [97] TH.data_1.0-8