04-read_data.Rmd

# Read in the data {#read-data}

This section describes how to read in single-cell data and images into `R`
**after** image processing and segmentation (see Section \@ref(processing)).

To highlight examples for IMC data analysis, we provide already processed data at
[10.5281/zenodo.6043599](https://zenodo.org/record/6043599).
This data has already been downloaded in Section \@ref(download-data) and can 
be accessed in the folder `data`.

We use the [imcRtools](https://github.com/BodenmillerGroup/imcRtools) package to
read in single-cell data extracted using the `steinbock` framework or the IMC
Segmentation Pipeline. Both image processing approaches also generate
multi-channel images and segmentation masks that can be read into `R` using the
[cytomapper](https://github.com/BodenmillerGroup/cytomapper) package.

```{r, message=FALSE}
library(imcRtools)
library(cytomapper)
```

## Read in single-cell information

For single-cell data analysis in `R` the
[SingleCellExperiment](https://bioconductor.org/packages/release/bioc/html/SingleCellExperiment.html)
[@Amezquita2019] data container is commonly used within the Bioconductor
framework. It allows standardized access to (i) expression data, (ii) cellular
metadata (e.g., cell type), (iii) feature metadata (e.g., marker name) and (iv)
experiment-wide metadata. For an in-depth introduction to the `SingleCellExperiment`
container, please refer to the [SingleCellExperiment class](https://bioconductor.org/books/release/OSCA.intro/the-singlecellexperiment-class.html).

The [SpatialExperiment](https://bioconductor.org/packages/release/bioc/html/SingleCellExperiment.html)
class [@Righelli2022] is an extension of the `SingleCellExperiment` class. It
was developed to store spatial data in addition to single-cell data and an
extended introduction is accessible
[here](https://bioconductor.org/packages/release/bioc/vignettes/SpatialExperiment/inst/doc/SpatialExperiment.html).

To read in single-cell data generated by the `steinbock` framework or the IMC
Segmentation Pipeline, the `imcRtools` package provides the `read_steinbock` and
`read_cpout` functions, respectively. By default, the data is read into a
`SpatialExperiment` object; however, data can be read in as a
`SingleCellExperiment` object by setting `return_as = "sce"`. All functions
presented in this book are applicable to both data containers.

### steinbock generated data

The downloaded example data (Section \@ref(download-data)) processed with the [steinbock](https://github.com/BodenmillerGroup/steinbock) framework can be read in with the `read_steinbock` function provided by `imcRtools`. For more information, please refer to
`?read_steinbock`.

```{r read-steinbock}
spe <- read_steinbock("data/steinbock/")
spe
```

By default, single-cell data is read in as `SpatialExperiment` object. 
The summarized pixel intensities per channel and cell (here mean intensity) are
stored in the `counts` slot. Columns represent cells and rows represent channels.

```{r counts}
counts(spe)[1:5,1:5]
```

Metadata associated to individual cells are stored in the `colData` slot. After
initial image processing, these metadata include the numeric identifier (`ObjectNumber`),
the area, and morphological features of each cell. In addition, `sample_id` stores
the image name from which each cell was extracted and the width and height of the
corresponding images are stored.

```{r colData}
head(colData(spe))
```

The main difference between the `SpatialExperiment` and the
`SingleCellExperiment` data container is the way spatial
locations of all cells are stored. For the `SingleCellExperiment` container, the
locations are stored in the `colData` slot while the `SpatialExperiment`
container stores them in the `spatialCoords` slot:

```{r spatialCoords}
head(spatialCoords(spe))
```

The _spatial object graphs_ generated by steinbock (see Section
\@ref(feature-extraction) are read into a `colPair` slot with the name
`neighborhood` of the `SpatialExperiment` (or `SingleCellExperiment`) object.
Cell-cell interactions (cells in close spatial proximity) are represented as
"edge list" (stored as `SelfHits` object). Here, the left side represents the
column indices of the `SpatialExperiment` object of the "from" cells and the
right side represents the column indices of the "to" cells. For visualization of
the _spatial object graphs_, please refer to Section \@ref(spatial-viz).

```{r colPair}
colPair(spe, "neighborhood")
```

Finally, metadata regarding the channels are stored in the `rowData` slot. This
information is extracted from the `panel.csv` file.

Channels have the same order as the rows in the `panel.csv` file for which the
keep column is set to 1, and match the order of channels in the multi-channel
images (see Section \@ref(read-images)). For the example data, channels are
ordered by isotope mass.

```{r rowData}
head(rowData(spe))
```

### IMC Segmentation Pipeline generated data

The [IMC Segmentation Pipeline](https://github.com/BodenmillerGroup/ImcSegmentationPipeline) offers an
alternative approach to multiplexed image processing and segmentation. The
default pipeline is also available via `steinbock`. The IMC Segmentation
Pipeline is based on [Ilastik](https://www.ilastik.org/) pixel classification
and image segmentation using [CellProfiler](https://cellprofiler.org/). We recommend
to become familiar with the pipeline as it allows flexible extension to more
complicated image analysis and segmentation tasks. For standard image analysis
and segmentation, `steinbock` is the preferred choice. Please refer to 
the [documentation](https://bodenmillergroup.github.io/ImcSegmentationPipeline/)
to get an overview on the pipeline.

All relevant [output](https://bodenmillergroup.github.io/ImcSegmentationPipeline/output.html)
storing single-cell data is contained in the `cpout` folder. 
For reading in the single-cell measurement, the `imcRtools` package offers the
`read_cpout` function:

```{r read-cpout}
spe2 <- read_cpout("data/ImcSegmentationPipeline/analysis/cpout/")
rownames(spe2) <- rowData(spe2)$Clean_Target
spe2
```

Similar to the `steinbock` output, cell morphological features and image level
metadata are stored in the `colData(spe2)` slot, the interaction information
is contained in `colPair(spe2, type = "neighborhood")` and the mean intensity 
per channel and cell is stored in `counts(spe2)`.

### Reading custom files

When not using `steinbock` or the `ImcSegmentationPipeline`, the single-cell
information has to be read in from custom files. We now demonstrate how
to generate a `SpatialExperiment` object from single-cell data contained
in individual files. As an example, we use files generated by `CellProfiler`
as part of the `ImcSegmentationPipeline`.

First we will read in the single-cell features stored in a CSV file:

```{r read-csv, message=FALSE}
library(readr)

cur_features <- read_csv("data/ImcSegmentationPipeline/analysis/cpout/cell.csv")

dim(cur_features)

head(colnames(cur_features))
```

This file contains a large number of single-cell features including the cell
identifier (`ObjectNumber`), the image identifier (`ImageNumber`), morphological 
features (`AreaShape_*`), the cells' locations (`Location_Center_*`) and the 
mean pixel intensity per cell and per channel (`Intensity_MeanIntensity_FullStack_*`).

Now, we split the features into intensity features, cell-specific metadata and
the physical location of the cells:

```{r split-features}
counts <- cur_features[,grepl("Intensity_MeanIntensity_FullStack", 
                                  colnames(cur_features))]

meta <- cur_features[,c("ImageNumber", "ObjectNumber", "AreaShape_Area",
                            "AreaShape_Eccentricity", "AreaShape_MeanRadius")]

coords <- cur_features[,c("Location_Center_X", "Location_Center_Y")]
```

`CellProfiler` writes out the mean pixel intensities after scaling them 
bit a scaling factor which is bit encoding-specific. The images to which
the IMC Segmentation Pipeline was applied were saved with 16-bit encoding.
This means for the example data, the mean pixel intensities need to 
be scaled by a factor of `2 ^ 16 - 1 = 65535`.

```{r scale-features}
counts <- counts * 65535
```

In addition, `CellProfiler` does not order the channel numerically but rather
as a character; `1, 10, 2, 3, ...` rather than `1, 2, 3, ...`. Therefore we
will need to reorder the channels.

```{r reorder-channels, message=FALSE}
library(stringr)
cur_ch <- str_split(colnames(counts), "_", simplify = TRUE)[,4]
cur_ch <- sub("c", "", cur_ch)

counts <- counts[,order(as.numeric(cur_ch))]
```

From these features we can now construct the `SpatialExperiment` object.

```{r construct-spe}
spe3 <- SpatialExperiment(assays = list(counts = t(counts)),
                          colData = meta, 
                          sample_id = as.character(meta$ImageNumber),
                          spatialCoords = as.matrix(coords))
```

Next, we can store the spatial cell graph generated by `CellProfiler` in the
`colPairs` slot of the object. Spatial cell graphs are usually stored as edge
list in form of a CSV file. The `colPairs` slot requires a `SelfHits` entry
storing an edge list where numeric entries represent the index of the `from` and
`to` cell in the `SpatialExperiment` object. To generate such an edge list, we
need to match the cell IDs contained in the CSV against the cell IDs in the
`SpatialExperiment` object.

```{r create-colPairs, message=FALSE}
cur_pairs <- read_csv("data/ImcSegmentationPipeline/analysis/cpout/Object relationships.csv")

cur_from <- paste(cur_pairs$`First Image Number`, cur_pairs$`First Object Number`)
cur_to <- paste(cur_pairs$`Second Image Number`, cur_pairs$`Second Object Number`)

edgelist <- SelfHits(from = match(cur_from, 
                                  paste(spe3$ImageNumber, spe3$ObjectNumber)),
                     to = match(cur_to, 
                                  paste(spe3$ImageNumber, spe3$ObjectNumber)),
                     nnode = ncol(spe3))

colPair(spe3, "neighborhood") <- edgelist
```

For further downstream analysis, we will use the `steinbock` results.

## Single-cell processing {#cell-processing}

After reading in the single-cell data, few further processing steps need to be
taken.

**Add additional metadata**

We can set the `colnames` of the object to generate unique identifiers per cell:

```{r set-colnames}
colnames(spe) <- paste0(spe$sample_id, "_", spe$ObjectNumber)
```

It is also often the case that sample-specific metadata are available externally.
For the current data, we need to link the cancer type (also referred to as "Indication")
to each sample. This metadata is available as external CSV file:

```{r add-indication, message=FALSE}
library(tidyverse)

# Read patient metadata
meta <- read_csv("data/sample_metadata.csv")

# Extract patient id and ROI id from sample name
spe$patient_id <- str_extract(spe$sample_id, "Patient[1-4]")
spe$ROI <- str_extract(spe$sample_id, "00[1-8]")

# Store cancer type in SPE object
spe$indication <- meta$Indication[match(spe$patient_id, meta$`Sample ID`)]

unique(spe$patient_id)
unique(spe$ROI)
unique(spe$indication)
```

The selected patients were diagnosed with different cancer types:

* SCCHN - head and neck cancer  
* BCC - breast cancer  
* NSCLC - lung cancer  
* CRC - colorectal cancer

**Transform counts**

The distribution of expression counts across cells is often observed to be
skewed towards the right side meaning lots of cells display low counts and few
cells have high counts. To avoid analysis biases from these high-expressing
cells, the expression counts are commonly transformed or clipped.

Here, we perform counts transformation using an inverse hyperbolic sine
function. This transformation is commonly applied to [flow cytometry
data](https://support.cytobank.org/hc/en-us/articles/206148057-About-the-Arcsinh-transform).
The `cofactor` here defines the expression range on which no scaling is
performed. While the `cofactor` for CyTOF data is often set to `5`, IMC data
usually display much lower counts. We therefore apply a `cofactor` of `1`.

However, other transformations such as `log(counts(spe) + 0.01)` should be
tested when analysing IMC data.

```{r transform-counts, message=FALSE}
library(dittoSeq)
dittoRidgePlot(spe, var = "CD3", group.by = "patient_id", assay = "counts") +
    ggtitle("CD3 - before transformation")
assay(spe, "exprs") <- asinh(counts(spe)/1)
dittoRidgePlot(spe, var = "CD3", group.by = "patient_id", assay = "exprs") +
    ggtitle("CD3 - after transformation")
```

**Define interesting channels**

For downstream analysis such as visualization, dimensionality reduction and
clustering, only a subset of markers should be used. As convenience, we can
store an additional entry in the `rowData` slot that specifies the markers of
interest. Here, we deselect the nuclear markers, which were primarily used for
cell segmentation, and keep all other biological targets. However, more informed
marker selection should be performed to exclude lowly expressed marker or
markers with low signal-to-noise ratio. 

```{r select-features}
rowData(spe)$use_channel <- !grepl("DNA|Histone", rownames(spe))
```

**Define color schemes**

We will define color schemes for different metadata entries of the data and
conveniently store them in the `metadata` slot of the `SpatialExperiment` which
will be helpful for downstream data visualizations. We will use colors from the
`RColorBrewer` and `dittoSeq` packages but any other coloring package will
suffice.

```{r define-colors, message=FALSE}
library(RColorBrewer)
color_vectors <- list()

ROI <- setNames(brewer.pal(length(unique(spe$ROI)), name = "BrBG"), 
                unique(spe$ROI))
patient_id <- setNames(brewer.pal(length(unique(spe$patient_id)), name = "Set1"), 
                unique(spe$patient_id))
sample_id <- setNames(c(brewer.pal(6, "YlOrRd")[3:5],
                        brewer.pal(6, "PuBu")[3:6],
                        brewer.pal(6, "YlGn")[3:5],
                        brewer.pal(6, "BuPu")[3:6]),
                unique(spe$sample_id))
indication <- setNames(brewer.pal(length(unique(spe$indication)), name = "Set2"), 
                unique(spe$indication))

color_vectors$ROI <- ROI
color_vectors$patient_id <- patient_id
color_vectors$sample_id <- sample_id
color_vectors$indication <- indication

metadata(spe)$color_vectors <- color_vectors
```

## Read in images {#read-images}

The `cytomapper` package allows multi-channel image handling and visualization
within the Bioconductor framework. The most common data format for multi-channel
images or segmentation masks is the TIFF file format, which is used by `steinbock`
and the `IMC segementation pipeline` to save images. 

Here, we will read in multi-channel images and segmentation masks into a
[CytoImageList](https://www.bioconductor.org/packages/release/bioc/vignettes/cytomapper/inst/doc/cytomapper.html#5_The_CytoImageList_object)
data container. It allows storing multiple multi-channel images and requires
matched channels across all images within the object.

The `loadImages` function is used to read in processed multi-channel images and
their corresponding segmentation masks. **Of note**: the multi-channel images
generated by `steinbock` are saved as 32-bit images while the segmentation masks
are saved as 16-bit images. To correctly scale pixel values of the segmentation
masks when reading them in, we will need to set `as.is = TRUE`.

```{r read-images}
images <- loadImages("data/steinbock/img/")
masks <- loadImages("data/steinbock/masks_deepcell/", as.is = TRUE)
```

In the case of multi-channel images, it is beneficial to set the `channelNames`
for easy visualization. Using the `steinbock` framework, the channel order of
the single-cell data matches the channel order of the multi-channel images.
However, it is recommended to make sure that the channel order is identical
between the single-cell data and the images.

```{r set-channelNames}
channelNames(images) <- rownames(spe)
images
```

For visualization shown in Section \@ref(image-visualization) we will need to
add additional metadata to the `elementMetadata` slot of the `CytoImageList`
objects. This slot is easily accessible using the `mcols` function.

Here, we will store the matched `sample_id`, `patient_id` and `indication`
information within the `elementMetadata` slot of the multi-channel images and
segmentation masks objects. It is crucial that the order of the images in 
both `CytoImageList` objects is the same.

```{r add-metadata}
all.equal(names(images), names(masks))

# Extract patient id from image name
patient_id <- str_extract(names(images), "Patient[1-4]")

# Retrieve cancer type per patient from metadata file
indication <- meta$Indication[match(patient_id, meta$`Sample ID`)] 

# Store patient and image level information in elementMetadata
mcols(images) <- mcols(masks) <- DataFrame(sample_id = names(images),
                                           patient_id = patient_id,
                                           indication = indication)
```

## Generate single-cell data from images

An alternative way of generating a `SingleCellExperiment` object directly 
from the multi-channel images and segmentation masks is supported by the 
[measureObjects](https://bodenmillergroup.github.io/cytomapper/reference/measureObjects.html)
function of the `cytomapper` package. For each cell present in the `masks`
object, the function computes the mean pixel intensity per channel as well as
morphological features (area, radius, major axis length, eccentricity) and the 
location of cells:

```{r measureObjects, message=FALSE}
cytomapper_sce <- measureObjects(masks, image = images, img_id = "sample_id")

cytomapper_sce
```

## Accessing publicly available IMC datasets

The [imcdatasets](https://github.com/BodenmillerGroup/imcdatasets)
R/Bioconductor package provides a number of publicly available IMC datasets. For
a complete introduction to the package, please refer to the
[documentation](https://bioconductor.org/packages/release/data/experiment/vignettes/imcdatasets/inst/doc/imcdatasets.html).
Here, we can read in example data of [@Damond2019] taken from patients diagnosed
with Type I Diabetes. The example here consists of a `CytoImageList` object of
100 images, a `CytoImageList` object of 100 segmentation masks and a
`SingleCellExperiment` object containing 252059 cells. **Of note:** downloading the
images takes quite some time and uses 8GB of memory.

```{r access-damond-data, message=FALSE, eval=FALSE}
library(imcdatasets)

pancreasImages <- Damond_2019_Pancreas(data_type = "images")
pancreasMasks <- Damond_2019_Pancreas(data_type = "masks")
pancreasSCE <- Damond_2019_Pancreas(data_type = "sce")
```

## Save objects

Finally, the generated data objects can be saved for further downstream 
processing and analysis.

```{r save-objects-read-data}
saveRDS(spe, "data/spe.rds")
saveRDS(images, "data/images.rds")
saveRDS(masks, "data/masks.rds")
```

```{r integration-test-read-single-cell, include = FALSE}
library(testthat)

expect_equal(dim(spe), c(40L, 47859L))

expect_equal(names(colData(spe)), c("sample_id", "ObjectNumber", "area", "axis_major_length", "axis_minor_length", 
"eccentricity", "width_px", "height_px", "patient_id", "ROI", 
"indication"))

expect_equal(colData(spe)[1:5, ], new("DFrame", rownames = c("Patient1_001_1", "Patient1_001_2", 
"Patient1_001_3", "Patient1_001_4", "Patient1_001_5"), nrows = 5L, 
    elementType = "ANY", elementMetadata = NULL, metadata = list(), 
    listData = list(sample_id = c("Patient1_001", "Patient1_001", 
    "Patient1_001", "Patient1_001", "Patient1_001"), ObjectNumber = c(1, 
    2, 3, 4, 5), area = c(12, 24, 17, 24, 22), axis_major_length = c(7.40623431863439, 
    16.4800397865237, 9.85084915362206, 8.08290376865476, 8.79366551809981
    ), axis_minor_length = c(1.89529355560732, 1.96283801650549, 
    1.98581719060974, 3.91578004149024, 3.11653157690275), eccentricity = c(0.966701893374626, 
    0.99288178876792, 0.979470319045067, 0.874817765279706, 0.935091484545154
    ), width_px = c(600, 600, 600, 600, 600), height_px = c(600, 
    600, 600, 600, 600), patient_id = c("Patient1", "Patient1", 
    "Patient1", "Patient1", "Patient1"), ROI = c("001", "001", 
    "001", "001", "001"), indication = c("SCCHN", "SCCHN", "SCCHN", 
    "SCCHN", "SCCHN"))))

 expect_equal(names(rowData(spe)), c("channel", "name", "keep", "ilastik", "deepcell", "cellpose", "Tube.Number", 
 "Target", "Antibody.Clone", "Stock.Concentration", "Final.Concentration...Dilution", 
 "uL.to.add", "use_channel"))
 
 expect_equal(rowData(spe)[1:5,], new("DFrame", rownames = c("MPO", "HistoneH3", "SMA", "CD16", 
"CD38"), nrows = 5L, elementType = "ANY", elementMetadata = NULL, 
    metadata = list(), listData = list(channel = c("Y89", "In113", 
    "In115", "Pr141", "Nd142"), name = c("MPO", "HistoneH3", 
    "SMA", "CD16", "CD38"), keep = c(1, 1, 1, 1, 1), ilastik = c(NA, 
    1, NA, NA, NA), deepcell = c(NA, 1, NA, NA, NA), cellpose = c(NA, 
    NA, NA, NA, NA), Tube.Number = c(2101, 2113, 1914, 2079, 
    2095), Target = c("Myeloperoxidase MPO", "Histone H3", "SMA", 
    "CD16", "CD38"), Antibody.Clone = c("Polyclonal MPO", "D1H2", 
    "1A4", "EPR16784", "EPR4106"), Stock.Concentration = c(500, 
    500, 500, 500, 500), Final.Concentration...Dilution = c("4 ug/mL", 
    "1 ug/mL", "0.25 ug/mL", "5 ug/mL", "2.5 ug/mL"), uL.to.add = c("0.8", 
    "0.2", "0.05", "1", "0.5"), use_channel = c(TRUE, FALSE, 
    TRUE, TRUE, TRUE))))
 
 expect_equal(rownames(spe), c("MPO", "HistoneH3", "SMA", "CD16", "CD38", "HLADR", "CD27", 
 "CD15", "CD45RA", "CD163", "B2M", "CD20", "CD68", "Ido1", "CD3", 
 "LAG3 / LAG33", "CD11c", "PD1", "PDGFRb", "CD7", "GrzB", "PDL1", 
 "TCF7", "CD45RO", "FOXP3", "ICOS", "CD8a", "CarbonicAnhydrase", 
 "CD33", "Ki67", "VISTA", "CD40", "CD4", "CD14", "Ecad", "CD303", 
 "CD206", "cleavedPARP", "DNA1", "DNA2"))
 
 expect_equal(colnames(spe)[100:150], c("Patient1_001_100", "Patient1_001_101", "Patient1_001_102", 
 "Patient1_001_103", "Patient1_001_104", "Patient1_001_105", "Patient1_001_106", 
 "Patient1_001_107", "Patient1_001_108", "Patient1_001_109", "Patient1_001_110", 
 "Patient1_001_111", "Patient1_001_112", "Patient1_001_113", "Patient1_001_114", 
 "Patient1_001_115", "Patient1_001_116", "Patient1_001_117", "Patient1_001_118", 
 "Patient1_001_119", "Patient1_001_120", "Patient1_001_121", "Patient1_001_122", 
 "Patient1_001_123", "Patient1_001_124", "Patient1_001_125", "Patient1_001_126", 
 "Patient1_001_127", "Patient1_001_128", "Patient1_001_129", "Patient1_001_130", 
 "Patient1_001_131", "Patient1_001_132", "Patient1_001_133", "Patient1_001_134", 
 "Patient1_001_135", "Patient1_001_136", "Patient1_001_137", "Patient1_001_138", 
 "Patient1_001_139", "Patient1_001_140", "Patient1_001_141", "Patient1_001_142", 
 "Patient1_001_143", "Patient1_001_144", "Patient1_001_145", "Patient1_001_146", 
 "Patient1_001_147", "Patient1_001_148", "Patient1_001_149", "Patient1_001_150"
 ))
 
expect_equal(assayNames(spe), c("counts", "exprs"))

expect_equal(as.vector(counts(spe)[1:5, 1:5]), c(0.575106422106425, 3.12730821967125, 0.260093857844671, 2.03477470080058, 
0.253013650576274, 0.416666666666667, 11.3597883482774, 1.67203828195731, 
2.58805358409882, 0.682666927576065, 0.497549358536215, 2.38414396257962, 
0.15351896426257, 2.2943073721493, 1.19029789812425, 0.890154048800468, 
7.71296111742655, 1.19394763807456, 15.6290833055973, 2.12606023748716, 
0.181818181818182, 1.45127152854746, 0.298670253970406, 0.60842195966027, 
0.29177929054607))

expect_equal(as.vector(assay(spe, "exprs")[1:5, 1:5]), c(0.547361972030777, 1.85795410049082, 0.257247185209247, 1.4590799040046, 
0.250389092048768, 0.405465108108164, 3.12515866853812, 1.28655606552992, 
1.67944602661282, 0.638408955377395, 0.479018831081476, 1.6033222122082, 
0.152922247068558, 1.56800637877875, 1.00974715085955, 0.801530560849479, 
2.74022551398308, 1.01209273247221, 3.44330257094523, 1.49863067518784, 
0.180831044069351, 1.16742600509862, 0.294399149051551, 0.576033240430445, 
0.287790189739289))

expect_equal(colPairNames(spe), "neighborhood")

expect_equal(colPair(spe)[1:5,], new("SelfHits", from = c(1L, 1L, 2L, 2L, 2L), to = c(27L, 55L, 
10L, 44L, 81L), nLnode = 47859L, nRnode = 47859L, elementMetadata = NULL, 
    metadata = list()))

expect_equal(spatialCoords(spe)[1:5,], structure(c(468.583333333333, 515.833333333333, 587.235294117647, 
192.25, 231.772727272727, 0.416666666666667, 0.416666666666667, 
0.470588235294118, 1.25, 0.909090909090909), dim = c(5L, 2L), dimnames = list(NULL, c("Pos_X", "Pos_Y"))))
```

```{r integration-test-read-images, include = FALSE}
expect_equal(names(images)[1:5], c("Patient1_001", "Patient1_002", "Patient1_003", "Patient2_001",
"Patient2_002"))

expect_s4_class(images, "CytoImageList")

expect_equal(mcols(images), new("DFrame", rownames = c("Patient1_001", "Patient1_002", "Patient1_003", 
"Patient2_001", "Patient2_002", "Patient2_003", "Patient2_004", 
"Patient3_001", "Patient3_002", "Patient3_003", "Patient4_005", 
"Patient4_006", "Patient4_007", "Patient4_008"), nrows = 14L, 
    elementType = "ANY", elementMetadata = NULL, metadata = list(), 
    listData = list(sample_id = c("Patient1_001", "Patient1_002", 
    "Patient1_003", "Patient2_001", "Patient2_002", "Patient2_003", 
    "Patient2_004", "Patient3_001", "Patient3_002", "Patient3_003", 
    "Patient4_005", "Patient4_006", "Patient4_007", "Patient4_008"
    ), patient_id = c("Patient1", "Patient1", "Patient1", "Patient2", 
    "Patient2", "Patient2", "Patient2", "Patient3", "Patient3", 
    "Patient3", "Patient4", "Patient4", "Patient4", "Patient4"
    ), indication = c("SCCHN", "SCCHN", "SCCHN", "BCC", "BCC", 
    "BCC", "BCC", "NSCLC", "NSCLC", "NSCLC", "CRC", "CRC", "CRC", 
    "CRC"))))

expect_equal(as.array(images$Patient1_001)[1:5,1:5,1:2],
             structure(c(0, 0, 0, 0, 1.76359105110168, 0, 0, 0, 1.68342781066895,
1, 0, 3.068843126297, 0, 0, 0, 0, 0, 0, 1.68342781066895, 0,
0, 1, 0, 0, 0, 21.8074226379395, 36.0178718566895, 21.4657039642334,
11.2600688934326, 2.51463913917542, 1.4867354631424, 6.26006889343262,
5.14433002471924, 3, 1, 1.4009622335434, 6.20082473754883, 8.63312721252441,
2.20082473754883, 1.00068736076355, 1.14364266395569, 13.5759449005127,
6.88769769668579, 5.66034412384033, 0, 2, 14.6624059677124, 22.2669429779053,
7.00206184387207, 9.37512016296387), .Dim = c(5L, 5L, 2L), .Dimnames = list(
    NULL, NULL, c("MPO", "HistoneH3"))))
```

```{r integration-test-read-masks, include = FALSE}
expect_equal(names(masks)[1:5], c("Patient1_001", "Patient1_002", "Patient1_003", "Patient2_001",
"Patient2_002"))

expect_s4_class(masks, "CytoImageList")

expect_equal(mcols(masks), new("DFrame", rownames = c("Patient1_001", "Patient1_002", "Patient1_003",
"Patient2_001", "Patient2_002", "Patient2_003", "Patient2_004",
"Patient3_001", "Patient3_002", "Patient3_003", "Patient4_005",
"Patient4_006", "Patient4_007", "Patient4_008"), nrows = 14L,
    listData = list(sample_id = c("Patient1_001", "Patient1_002",
    "Patient1_003", "Patient2_001", "Patient2_002", "Patient2_003",
    "Patient2_004", "Patient3_001", "Patient3_002", "Patient3_003",
    "Patient4_005", "Patient4_006", "Patient4_007", "Patient4_008"
    ), patient_id = c("Patient1", "Patient1", "Patient1",
    "Patient2", "Patient2", "Patient2", "Patient2", "Patient3",
    "Patient3", "Patient3", "Patient4", "Patient4", "Patient4",
    "Patient4"), indication = c("SCCHN", "SCCHN", "SCCHN", "BCC",
    "BCC", "BCC", "BCC", "NSCLC", "NSCLC", "NSCLC", "CRC", "CRC",
    "CRC", "CRC")), elementType = "ANY", elementMetadata = NULL,
    metadata = list()))

expect_equal(as.array(masks$Patient1_001)[1:5,1:5], structure(c(58L, 58L, 58L, 58L, 58L, 58L, 58L, 58L, 58L, 58L, 
58L, 58L, 58L, 58L, 58L, 58L, 58L, 58L, 58L, 58L, 58L, 58L, 58L, 
58L, 58L), dim = c(5L, 5L)))
```

```{r integration-test-imcsegmentation-pipeline, include = FALSE}
expect_equal(dim(spe2), c(40L, 43796L))

expect_equal(names(colData(spe2)), c("sample_id", "ObjectNumber", "AreaShape_Area", "Neighbors_NumberOfNeighbors_8", 
"AreaShape_Eccentricity", "AreaShape_MajorAxisLength", "AreaShape_MinorAxisLength", 
"AreaShape_MeanRadius", "ImageNumber", "Metadata_acname", "Metadata_acid", 
"Metadata_description"))

expect_equal(colData(spe2)[1:5, ], new("DFrame", rownames = NULL, nrows = 5L, elementType = "ANY", 
    elementMetadata = NULL, metadata = list(), listData = list(
        sample_id = c("1", "1", "1", "1", "1"), ObjectNumber = c(1, 
        2, 3, 4, 5), AreaShape_Area = c(61, 40, 27, 48, 34), 
        Neighbors_NumberOfNeighbors_8 = c(6, 8, 7, 7, 7), AreaShape_Eccentricity = c(0.808013494907822, 
        0.801208429202, 0.52994193892238, 0.839066042071405, 
        0.774470900280141), AreaShape_MajorAxisLength = c(12.0610890660951, 
        9.34348318770232, 6.39364541933779, 11.3851699913821, 
        9.72619134201331), AreaShape_MinorAxisLength = c(7.10595920371522, 
        5.59100365955202, 5.4220282329127, 6.19386917672975, 
        6.15288138869827), AreaShape_MeanRadius = c(1.56945391923939, 
        1.53106601717798, 1.40344156396267, 1.47571815899011, 
        1.24432752321663), ImageNumber = c(1, 1, 1, 1, 1), Metadata_acname = c("Patient1", 
        "Patient1", "Patient1", "Patient1", "Patient1"), Metadata_acid = c(1, 
        1, 1, 1, 1), Metadata_description = c("pos1_1", "pos1_1", 
        "pos1_1", "pos1_1", "pos1_1"))))

 expect_equal(names(rowData(spe2)), c("Tube.Number", "Metal.Tag", "Target", "Clean_Target", "Antibody.Clone", 
"Stock.Concentration", "Final.Concentration...Dilution", "uL.to.add", 
"full", "ilastik", "deepcell"))
 
 expect_equal(rowData(spe2)[1:5,], new("DFrame", rownames = c("MPO", "HistoneH3", "SMA", "CD16", 
"CD38"), nrows = 5L, elementType = "ANY", elementMetadata = NULL, 
    metadata = list(), listData = list(Tube.Number = c(2101, 
    2113, 1914, 2079, 2095), Metal.Tag = c("Y89", "In113", "In115", 
    "Pr141", "Nd142"), Target = c("Myeloperoxidase MPO", "Histone H3", 
    "SMA", "CD16", "CD38"), Clean_Target = c("MPO", "HistoneH3", 
    "SMA", "CD16", "CD38"), Antibody.Clone = c("Polyclonal MPO", 
    "D1H2", "1A4", "EPR16784", "EPR4106"), Stock.Concentration = c(500, 
    500, 500, 500, 500), Final.Concentration...Dilution = c("4 ug/mL", 
    "1 ug/mL", "0.25 ug/mL", "5 ug/mL", "2.5 ug/mL"), uL.to.add = c("0.8", 
    "0.2", "0.05", "1", "0.5"), full = c(1, 1, 1, 1, 1), ilastik = c(0, 
    1, 0, 0, 0), deepcell = c(NA, 1, NA, NA, NA))))
 
 expect_equal(rownames(spe2), c("MPO", "HistoneH3", "SMA", "CD16", "CD38", "HLADR", "CD27", 
 "CD15", "CD45RA", "CD163", "B2M", "CD20", "CD68", "Ido1", "CD3", 
 "LAG3", "CD11c", "PD1", "PDGFRb", "CD7", "GrzB", "PDL1", 
 "TCF7", "CD45RO", "FOXP3", "ICOS", "CD8a", "CarbonicAnhydrase", 
 "CD33", "Ki67", "VISTA", "CD40", "CD4", "CD14", "Ecad", "CD303", 
 "CD206", "cleavedPARP", "DNA1", "DNA2"))

expect_equal(as.vector(counts(spe2)[1:5, 1:5]), c(0.344262295001812, 4.77049180216797, 0.557377049050553, 0.770491803099294, 
0.393442622859214, 0.374999999912689, 6.74999999842776, 0.724999999831198, 
0.224999999947613, 0.899999999790452, 0.185185185142068, 12.8518518488591, 
1.03703703679558, 0.555555555426205, 0.999999999767169, 0.31249999992724, 
14.5624999966091, 1.27083333303744, 0.791666666482343, 0.666666666511446, 
0.323529411689378, 11.8529411737103, 1.08823529386427, 6.8235294101755, 
0.794117646873929))

expect_equal(colPairNames(spe2), "neighborhood")

expect_equal(colPair(spe2)[1:5,], new("SelfHits", from = c(1L, 1L, 1L, 1L, 1L), to = c(2L, 3L, 
41L, 51L, 79L), nLnode = 43796L, nRnode = 43796L, elementMetadata = NULL, 
    metadata = list()))

expect_equal(spatialCoords(spe2)[1:5,], structure(c(10.2295081967213, 18.7, 25.4074074074074, 32.2708333333333, 
45.0294117647059, 2.62295081967213, 2, 1.81481481481481, 2.22916666666667, 
1.5), dim = c(5L, 2L), dimnames = list(NULL, c("Pos_X", "Pos_Y"))))
```

```{r integration-test-measureObjects, include = FALSE}
expect_equal(counts(spe), counts(cytomapper_sce), check.attributes = FALSE)

expect_equal(spe$area, cytomapper_sce$s.area)
expect_equal(spe$axis_major_length, cytomapper_sce$m.majoraxis)
expect_equal(spe$eccentricity, cytomapper_sce$m.eccentricity)
expect_equal(as.vector(spatialCoords(spe)[,1]), cytomapper_sce$m.cx - 1)
expect_equal(as.vector(spatialCoords(spe)[,2]), cytomapper_sce$m.cy - 1)
```

## Session Info

<details>
   <summary>SessionInfo</summary>
   
```{r, echo = FALSE}
sessionInfo()
```
</details>