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
[https://zenodo.org/record/6043600](https://zenodo.org/record/6043600).
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/3.14/OSCA.intro/the-singlecellexperiment-class.html).

The [SpatialExperiment](https://bioconductor.org/packages/release/bioc/html/SingleCellExperiment.html)
class [@Righelli2021] 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` function, 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 in the current setting 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 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
"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 are ordered by
isotope mass and therefore match the channel order of the multi-channel images
(see Section \@ref(read-images)).

```{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)`.

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 excel file:

```{r add-indication}
library(openxlsx)
library(stringr)
meta <- read.xlsx("data/sample_metadata.xlsx")

spe$patient_id <- as.vector(str_extract_all(spe$sample_id, "Patient[1-4]", simplify = TRUE))
spe$ROI <- as.vector(str_extract_all(spe$sample_id, "00[1-8]", simplify = TRUE))
spe$indication <- meta$Indication[match(spe$patient_id, meta$Sample.ID)]

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.

```{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` package 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(dittoColors(reps = 1)[seq_along(unique(spe$sample_id))], 
                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 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 save 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))
patient_id <- str_extract_all(names(images), "Patient[1-4]", simplify = TRUE)
indication <- meta$Indication[match(patient_id, meta$Sample.ID)] 

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
```

## 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(names(colData(spe)), c("sample_id", "ObjectNumber", "area", "major_axis_length", "minor_axis_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, 
#    listData = list(sample_id = c("Patient1_001", "Patient1_001", 
#    "Patient1_001", "Patient1_001", "Patient1_001"), ObjectNumber = c(1, 
#    2, 3, 4, 5), area = c(22, 8, 28, 23, 27), major_axis_length = c(12.6491106406735, 
#    4.47213595499958, 10.6701384606448, 8.87999016101037, 10.7326126996496
#    ), minor_axis_length = c(2, 2, 3.52327281843438, 3.21095561912589, 
#    3.13466326780999), eccentricity = c(0.987420882906575, 0.894427190999916, 
#    0.943911217824785, 0.932335463337873, 0.956397255130736), 
#        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")), elementType = "ANY", elementMetadata = NULL, 
#     metadata = list()))
# expect_equal(names(rowData(spe)), c("channel", "name", "keep", "ilastik", "deepcell", "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, 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), 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)), elementType = "ANY", elementMetadata = NULL, metadata = list()))
# 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.5, 11.6130212111907, 1.82404176213525, 3.22545072165402, 
# 0.809814534404061, 0.375, 3.12646059691906, 0.283107548952103, 
# 0.83465687930584, 0.427393049001694, 0.34352376631328, 4.82839978592736, 
# 0.929368317127228, 1.01492113300732, 0.431800333516938, 0.340718041295591, 
# 6.75667715072632, 0.974456113317738, 6.96184600436169, 1.14356865571893, 
# 0.33489340322989, 6.11668086934973, 5.71861606174045, 0.594917981712906, 
# 0.39980254791401))
# expect_equal(as.vector(assay(spe, "exprs")[1:5, 1:5]), c(0.481211825059603, 3.14712277710762, 1.36205720976067, 1.88742752156615, 
# 0.740313787257858, 0.366724604230137, 1.85769590688433, 0.279455938469543, 
# 0.759502582492501, 0.415347442233046, 0.337102761865711, 2.27821722251645, 
# 0.830537454818704, 0.891885162165194, 0.419396869687053, 0.3344480991532, 
# 2.60911009071464, 0.86319550027886, 2.63871050928924, 0.979339314203556, 
# 0.328929815775784, 2.51078282058994, 2.44443246736646, 0.564462221647734, 
# 0.389851983870733))
# expect_equal(colPairNames(spe), "neighborhood")
# expect_equal(colPair(spe)[1:5,], new("SelfHits", from = c(1L, 1L, 2L, 2L, 2L), to = c(41L, 69L, 
# 18L, 58L, 85L), nLnode = 43289L, nRnode = 43289L, elementMetadata = NULL, 
#     metadata = list()))
# expect_equal(spatialCoords(spe)[1:5,], structure(c(517, 231.5, 271.464285714286, 287.347826086957, 410.62962962963, 
# 0.5, 0.5, 1, 0.956521739130435, 0.888888888888889), .Dim = c(5L, 
# 2L), .Dimnames = list(c("1", "2", "3", "4", "5"), 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, 
#     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.V1 = 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(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.V1 = 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(52L, 52L, 52L, 52L, 52L, 52L, 52L, 52L, 52L, 52L, 
# 52L, 52L, 52L, 52L, 52L, 52L, 52L, 52L, 52L, 52L, 52L, 52L, 52L, 
# 52L, 52L), .Dim = c(5L, 5L)))
```

## Session Info

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