Spectre - advanced spatial analysis

---

<- Back to Spectre home page

Introduction

---

Here we provide instructions for analysing spatial data after cell segmentation as part of Spectre’s advanced segmentation and spatial analysis workflow. This workflow is split into 3x scripts, which can be run after performing cell segmentation:

1. Extract cellular data from ROIs and masks
2. Cellular analysis
3. Spatial analysis

Setup

---

If you go to https://github.com/ImmuneDynamics/Spectre you can download the repository, including the analysis workflow scripts.

You can find the advanced spatial analysis workflow scripts in this folder.

You will need a folder called ‘data’ containing a folder called ‘ROIs’. Each ROI should be a folder with a unique name, and the TIFF files for that ROI should be in the corresponding folder.

Within the ‘data’ folder, you will need a folder called ‘masks’. Each of your mask types should be stored here – with the ROI name, and then some pattern that can separate which type of mask it is.

You will also need a ‘metadata’ folder containing a CSV file called ‘sample.metadata’. In this file you should include a column called ‘ROI’ with the name of each ROI, and then as many other columns as you like containing metadata for each ROI. We recommend including ‘Sample’ (e.g. if multiple ROIs are taken from one sample, you can note this here), ‘Group’ (the experiment group each ROI belongs to), and ‘Batch’ (the sample preparation/acquisition batch for the ROI – this could be some kind of reference, or a date).

Script 1. Extract cellular data

---

Load libraries

…

###################################################################################
### Spectre: spatial 1 - add masks and extract cellular data
###################################################################################

First, load the Spectre package and associated packages.

    ### Load libraries

    library('Spectre')

    Spectre::package.check(type = 'spatial')
    Spectre::package.load(type = 'spatial')

Next, set directories.

    ### Set PrimaryDirectory

    dirname(rstudioapi::getActiveDocumentContext()$path)
    setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
    getwd()
    PrimaryDirectory <- getwd()
    PrimaryDirectory

    ### Set InputDirectory (ROI TIFFs)

    setwd(PrimaryDirectory)
    setwd("data/ROIs/")
    InputDirectory <- getwd()
    InputDirectory

    ### Set MaskDirectory (ROI mask TIFFs)

    setwd(PrimaryDirectory)
    setwd("data/masks")
    MaskDirectory <- getwd()
    MaskDirectory

    ### Create output directory

    setwd(PrimaryDirectory)
    dir.create("Output 1 - add masks")
    setwd("Output 1 - add masks")
    OutputDirectory <- getwd()
    OutputDirectory

Check ROIs and TIFFs

###################################################################################
### Check ROIs and TIFFs
###################################################################################

### Initialise the spatial data object with channel TIFF files

    setwd(InputDirectory)

    rois <- list.dirs(full.names = FALSE, recursive = FALSE)
    as.matrix(rois)

##      [,1]
## [1,] "ROI002"
## [2,] "ROI004"
## [3,] "ROI006"
## [4,] "ROI008"
## [5,] "ROI010"
## [6,] "ROI012"

    ### Check channel names

    tiff.list <- list()

    for(i in rois){
      setwd(InputDirectory)
      setwd(i)
      tiff.list[[i]] <- list.files(getwd())
    }

    t(as.data.frame(tiff.list))

##        [,1]               [,2]              [,3]             [,4]
## ROI002 "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff"
## ROI004 "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff"
## ROI006 "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff"
## ROI008 "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff"
## ROI010 "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff"
## ROI012 "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff"
##        [,5]              [,6]              [,7]              [,8]
## ROI002 "CD45_Sm152.tiff" "CD8a_Dy162.tiff" "DNA1_Ir191.tiff" "DNA3_Ir193.tiff"
## ROI004 "CD45_Sm152.tiff" "CD8a_Dy162.tiff" "DNA1_Ir191.tiff" "DNA3_Ir193.tiff"
## ROI006 "CD45_Sm152.tiff" "CD8a_Dy162.tiff" "DNA1_Ir191.tiff" "DNA3_Ir193.tiff"
## ROI008 "CD45_Sm152.tiff" "CD8a_Dy162.tiff" "DNA1_Ir191.tiff" "DNA3_Ir193.tiff"
## ROI010 "CD45_Sm152.tiff" "CD8a_Dy162.tiff" "DNA1_Ir191.tiff" "DNA3_Ir193.tiff"
## ROI012 "CD45_Sm152.tiff" "CD8a_Dy162.tiff" "DNA1_Ir191.tiff" "DNA3_Ir193.tiff"

Read in TIFF files

###################################################################################
### Read in TIFF files and create spatial objects
################################################################################### 

    ### Read in ROI channel TIFFs

    setwd(InputDirectory)
    spatial.dat <- read.spatial.files(dir = InputDirectory)

As the files are read in, you will see feedback like this:

## Reading TIFFs from:/Users/thomasa/OneDrive - The University of Sydney (Staff)/Library/Github (public)/Spectre/workflows/Spatial workflows/Spatial - advanced/data/ROIs

##   ...with extent correction

##   ...with y-axis orientation flipping

## Reading ROI: ROI002

##   -- reading single band in TIFF:CD11b_Sm149.tiff

##   -- reading single band in TIFF:CD20_Dy161.tiff

##   -- reading single band in TIFF:CD3_Er170.tiff

##   -- reading single band in TIFF:CD4_Gd156.tiff

##   -- reading single band in TIFF:CD45_Sm152.tiff

##   -- reading single band in TIFF:CD8a_Dy162.tiff

##   -- reading single band in TIFF:DNA1_Ir191.tiff

##   -- reading single band in TIFF:DNA3_Ir193.tiff

    ### Check results

    str(spatial.dat, 3)

## List of 6
##  $ ROI002:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  : list()
##   .. ..@ DATA   : list()
##  $ ROI004:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  : list()
##   .. ..@ DATA   : list()
##  $ ROI006:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  : list()
##   .. ..@ DATA   : list()
##  $ ROI008:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  : list()
##   .. ..@ DATA   : list()
##  $ ROI010:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  : list()
##   .. ..@ DATA   : list()
##  $ ROI012:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  : list()
##   .. ..@ DATA   : list()

    spatial.dat[[1]]@RASTERS

## class      : RasterStack
## dimensions : 501, 500, 250500, 8  (nrow, ncol, ncell, nlayers)
## resolution : 1, 1  (x, y)
## extent     : 0, 500, 0, 501  (xmin, xmax, ymin, ymax)
## crs        : NA
## names      : CD11b_Sm149, CD20_Dy161, CD3_Er170, CD4_Gd156, CD45_Sm152, CD8a_Dy162, DNA1_Ir191, DNA3_Ir193
## min values :           0,          0,         0,         0,          0,          0,          0,          0
## max values :          33,        779,        41,        40,        734,        109,       1670,       3007

Read in masks files

###################################################################################
### Read in masks files
###################################################################################

    ### Define cell mask extension for different mask types

        setwd(MaskDirectory)

        all.masks <- list.files(pattern = '.tif')
        as.matrix(all.masks)

##       [,1]
##  [1,] "ROI002_ilastik_s2_Object Identities.tif"
##  [2,] "ROI002_ilastik_s2_Object Predictions.tif"
##  [3,] "ROI002_ilastik_s2_Simple Segmentation.tif"
##  [4,] "ROI004_ilastik_s2_Object Identities.tif"
##  [5,] "ROI004_ilastik_s2_Object Predictions.tif"
##  [6,] "ROI004_ilastik_s2_Simple Segmentation.tif"
##  [7,] "ROI006_ilastik_s2_Object Identities.tif"
##  [8,] "ROI006_ilastik_s2_Object Predictions.tif"
##  [9,] "ROI006_ilastik_s2_Simple Segmentation.tif"
## [10,] "ROI008_ilastik_s2_Object Identities.tif"
## [11,] "ROI008_ilastik_s2_Object Predictions.tif"
## [12,] "ROI008_ilastik_s2_Simple Segmentation.tif"
## [13,] "ROI010_ilastik_s2_Object Identities.tif"
## [14,] "ROI010_ilastik_s2_Object Predictions.tif"
## [15,] "ROI010_ilastik_s2_Simple Segmentation.tif"
## [16,] "ROI012_ilastik_s2_Object Identities.tif"
## [17,] "ROI012_ilastik_s2_Object Predictions.tif"
## [18,] "ROI012_ilastik_s2_Simple Segmentation.tif"

        mask.types <- list('cell.mask' = '_ilastik_s2_Object Identities.tif',
                           'cell.type' = '_ilastik_s2_Object Predictions.tif',
                           'regions' = '_ilastik_s2_Simple Segmentation.tif')

        mask.types

## $cell.mask
## [1] "_ilastik_s2_Object Identities.tif"
##
## $cell.type
## [1] "_ilastik_s2_Object Predictions.tif"
##
## $regions
## [1] "_ilastik_s2_Simple Segmentation.tif"

    ### Read in masks

        for(i in names(mask.types)){
              spatial.dat <- do.add.masks(dat = spatial.dat,
                                          mask.dir = MaskDirectory,
                                          mask.pattern = mask.types[[i]],
                                          mask.label = i)
        }

        str(spatial.dat, 3)

## List of 6
##  $ ROI002:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  :List of 3
##   .. ..@ DATA   : list()
##  $ ROI004:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  :List of 3
##   .. ..@ DATA   : list()
##  $ ROI006:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  :List of 3
##   .. ..@ DATA   : list()
##  $ ROI008:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  :List of 3
##   .. ..@ DATA   : list()
##  $ ROI010:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  :List of 3
##   .. ..@ DATA   : list()
##  $ ROI012:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  :List of 3
##   .. ..@ DATA   : list()

        str(spatial.dat[[1]]@MASKS, 3)

## List of 3
##  $ cell.mask:List of 1
##   ..$ maskraster:Formal class 'RasterLayer' [package "raster"] with 12 slots
##  $ cell.type:List of 1
##   ..$ maskraster:Formal class 'RasterLayer' [package "raster"] with 12 slots
##  $ regions  :List of 1
##   ..$ maskraster:Formal class 'RasterLayer' [package "raster"] with 12 slots

Rename rasters (if required)

###################################################################################
### Rename rasters (if required)
###################################################################################

Here you can check the consistency of the raster names, and if required, adjust those names.

    ### Check channel names

        channel.names <- list()

        for(i in names(spatial.dat)){
          channel.names[[i]] <- names(spatial.dat[[i]]@RASTERS)
        }

        t(as.data.frame(channel.names))

##        [,1]          [,2]         [,3]        [,4]        [,5]
## ROI002 "CD11b_Sm149" "CD20_Dy161" "CD3_Er170" "CD4_Gd156" "CD45_Sm152"
## ROI004 "CD11b_Sm149" "CD20_Dy161" "CD3_Er170" "CD4_Gd156" "CD45_Sm152"
## ROI006 "CD11b_Sm149" "CD20_Dy161" "CD3_Er170" "CD4_Gd156" "CD45_Sm152"
## ROI008 "CD11b_Sm149" "CD20_Dy161" "CD3_Er170" "CD4_Gd156" "CD45_Sm152"
## ROI010 "CD11b_Sm149" "CD20_Dy161" "CD3_Er170" "CD4_Gd156" "CD45_Sm152"
## ROI012 "CD11b_Sm149" "CD20_Dy161" "CD3_Er170" "CD4_Gd156" "CD45_Sm152"
##        [,6]         [,7]         [,8]
## ROI002 "CD8a_Dy162" "DNA1_Ir191" "DNA3_Ir193"
## ROI004 "CD8a_Dy162" "DNA1_Ir191" "DNA3_Ir193"
## ROI006 "CD8a_Dy162" "DNA1_Ir191" "DNA3_Ir193"
## ROI008 "CD8a_Dy162" "DNA1_Ir191" "DNA3_Ir193"
## ROI010 "CD8a_Dy162" "DNA1_Ir191" "DNA3_Ir193"
## ROI012 "CD8a_Dy162" "DNA1_Ir191" "DNA3_Ir193"

 ### List of corrections (first entry is the 'correct' one)

        # corrections <- list(c('CD4','Cd4'),
        #                     c('CD8','CD8a')
        #                     )

    ### Replace the 'incorrect' names

        # for(i in names(spatial.dat)){
        #   # i <- names(spatial.dat)[[1]]
        #
        #   for(a in c(1:length(corrections))){
        #     # a <- 1
        #
        #     trg <- which(names(spatial.dat[[i]]@RASTERS) == corrections[[a]][2])
        #     if(length(trg) != 0){
        #       names(spatial.dat[[i]]@RASTERS)[trg] <- corrections[[a]][1]
        #     }
        #   }
        # }

    ### Check channel names

        # channel.names <- list()
        #
        # for(i in names(spatial.dat)){
        #   channel.names[[i]] <- names(spatial.dat[[i]]@RASTERS)
        # }
        #
        # t(as.data.frame(channel.names))      

Generate outlines

###################################################################################
### Generate polygons and outlines
###################################################################################

    ### Generate polygons and outlines

        for(i in names(mask.types)){
          spatial.dat <- do.create.outlines(dat = spatial.dat, mask.name = i)
        }

    ### Checks

        str(spatial.dat, 3)

## List of 6
##  $ ROI002:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  :List of 3
##   .. ..@ DATA   : list()
##  $ ROI004:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  :List of 3
##   .. ..@ DATA   : list()
##  $ ROI006:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  :List of 3
##   .. ..@ DATA   : list()
##  $ ROI008:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  :List of 3
##   .. ..@ DATA   : list()
##  $ ROI010:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  :List of 3
##   .. ..@ DATA   : list()
##  $ ROI012:Formal class 'spatial' [package "Spectre"] with 3 slots
##   .. ..@ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
##   .. ..@ MASKS  :List of 3
##   .. ..@ DATA   : list()

        str(spatial.dat[[1]]@MASKS, 2)

## List of 3
##  $ cell.mask:List of 4
##   ..$ maskraster:Formal class 'RasterLayer' [package "raster"] with 12 slots
##   ..$ polygons  :Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..$ outlines  :'data.frame':   102730 obs. of  7 variables:
##   ..$ centroids :Formal class 'SpatialPoints' [package "sp"] with 3 slots
##  $ cell.type:List of 4
##   ..$ maskraster:Formal class 'RasterLayer' [package "raster"] with 12 slots
##   ..$ polygons  :Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..$ outlines  :'data.frame':   23699 obs. of  7 variables:
##   ..$ centroids :Formal class 'SpatialPoints' [package "sp"] with 3 slots
##  $ regions  :List of 4
##   ..$ maskraster:Formal class 'RasterLayer' [package "raster"] with 12 slots
##   ..$ polygons  :Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..$ outlines  :'data.frame':   8452 obs. of  7 variables:
##   ..$ centroids :Formal class 'SpatialPoints' [package "sp"] with 3 slots

Mask QC plots

###################################################################################
### Mask QC plots
###################################################################################      

Here you can choose a ‘base’ raster to use for the spatial plot, and which mask you would like to plot.

    ### Mask plot setup

        setwd(OutputDirectory)
        dir.create('Plots - cell masks')
        setwd('Plots - cell masks')

        as.matrix(names(spatial.dat[[1]]@RASTERS))

##      [,1]
## [1,] "CD11b_Sm149"
## [2,] "CD20_Dy161"
## [3,] "CD3_Er170"
## [4,] "CD4_Gd156"
## [5,] "CD45_Sm152"
## [6,] "CD8a_Dy162"
## [7,] "DNA1_Ir191"
## [8,] "DNA3_Ir193"

        base <- 'DNA1_Ir191'
        base

## [1] "DNA1_Ir191"

        as.matrix(names(spatial.dat[[1]]@MASKS))

##      [,1]
## [1,] "cell.mask"
## [2,] "cell.type"
## [3,] "regions"

        mask <- 'cell.mask'
        mask    

## [1] "cell.mask"

    ### Create plots

        for(i in names(spatial.dat)){
          make.spatial.plot(dat = spatial.dat,
                            image.roi = i,
                            image.channel = base,
                            mask.outlines = mask)
        }

You will see plots that look like this for each ROI, allowing you to assess the suitability of the mask.

    ### Create plots

        for(i in names(spatial.dat)[1]){
          make.spatial.plot(dat = spatial.dat,
                            image.roi = i,
                            image.channel = base,
                            mask.outlines = mask)
        }

Calculate cellular data

###################################################################################
### Calculate cellular data and plot
###################################################################################   

    ### Calculate cellular data for each cell mask (this step may take some time)

        spatial.dat <- do.extract(spatial.dat, 'cell.mask')
        str(spatial.dat, 3)

        spatial.dat[[1]]@DATA

    ### Factor plot setup

        setwd(OutputDirectory)
        dir.create('Plots - factors')
        setwd('Plots - factors')

        plot.rois <- names(spatial.dat)[c(1:2)]
        plot.rois

## [1] "ROI002" "ROI004"

        as.matrix(names(spatial.dat[[1]]@RASTERS))

##      [,1]
## [1,] "CD11b_Sm149"
## [2,] "CD20_Dy161"
## [3,] "CD3_Er170"
## [4,] "CD4_Gd156"
## [5,] "CD45_Sm152"
## [6,] "CD8a_Dy162"
## [7,] "DNA1_Ir191"
## [8,] "DNA3_Ir193"

        base <- 'DNA1_Ir191'
        base

## [1] "DNA1_Ir191"

        plot.factors <- names(spatial.dat[[1]]@MASKS)[-which('cell.mask' == names(spatial.dat[[1]]@MASKS))]
        plot.factors

## [1] "cell.type" "regions"

        plot.exp <- names(spatial.dat[[1]]@RASTERS)
        plot.exp

## [1] "CD11b_Sm149" "CD20_Dy161"  "CD3_Er170"   "CD4_Gd156"   "CD45_Sm152"
## [6] "CD8a_Dy162"  "DNA1_Ir191"  "DNA3_Ir193"

    ### Make factor plots

        for(i in plot.rois){

          setwd(OutputDirectory)
          setwd('Plots - factors')
          dir.create(i)
          setwd(i)

          for(a in plot.factors){
            make.spatial.plot(dat = spatial.dat,
                              image.roi = i,
                              image.channel = base,
                              mask.outlines = mask,
                              cell.dat = 'CellData',
                              cell.col = a,
                              cell.col.type = 'factor',
                              title = paste0(a, ' - ', i))
          }
        }

    ### Make exp plots

        for(i in plot.rois){

          setwd(OutputDirectory)
          setwd('Plots - factors')
          dir.create(i)
          setwd(i)

          for(a in plot.exp){
            make.spatial.plot(dat = spatial.dat,
                              image.roi = i,
                              image.channel = base,
                              mask.outlines = mask,
                              cell.dat = 'CellData',
                              cell.col = a,
                              title = paste0(a, ' - ', i))
          }
        }

Extract cellular data and annotate

###################################################################################
### Extract cellular data and annotate
###################################################################################  

    ### Extract cellular data

        cell.dat <- do.pull.data(spatial.dat, 'CellData')
        cell.dat

##           ROI   ID          x          y Area CD11b_Sm149 CD20_Dy161 CD3_Er170
##     1: ROI002    1 147.666667 500.166667    6  0.33333334 10.1666670 0.3333333
##     2: ROI002    2  97.000000 500.000000   16  0.37500000  1.3750000 0.5000000
##     3: ROI002    3 470.815789 499.710526   19  0.57894737  0.1052632 0.2631579
##     4: ROI002    4 310.500000 499.833333   12  0.08333334 60.5833321 0.2500000
##     5: ROI002    5 245.166667 499.833333   12  0.33333334 28.2500000 2.0000000
##    ---
## 25348: ROI012 2907  45.569767   2.081395   43  0.34883720  3.6744187 0.5581396
## 25349: ROI012 2908  36.675676   4.500000   74  0.27027026  0.8918919 1.9459460
## 25350: ROI012 2909  27.035714   2.446429   56  0.32142857  0.5357143 0.3571429
## 25351: ROI012 2910  16.928571   2.128571   35  0.45714286  0.9142857 1.3428571
## 25352: ROI012 2911   5.986486  10.918919  222  0.31081080  0.7567568 0.2027027
##        CD4_Gd156 CD45_Sm152 CD8a_Dy162 DNA1_Ir191 DNA3_Ir193 cell.type regions
##     1: 0.0000000  5.6666665  0.5000000  54.666668   98.50000     65534   65535
##     2: 1.7500000  1.8125000  0.8125000   9.750000   15.62500     65534   65534
##     3: 0.7368421  0.2631579  0.6842105   8.578947   16.31579     65534   65534
##     4: 0.1666667  7.0833335  1.5833334  50.833332  104.41666     65534   65535
##     5: 0.0000000 18.6666660  4.2500000  95.916664  171.58333     65534   65535
##    ---
## 25348: 0.5581396  2.0465117  0.9767442  52.651161   99.34884     65534   65534
## 25349: 0.4324324  7.4594593  3.1216216  61.621620  115.50000     65534   65534
## 25350: 1.1964285  1.3214285  1.0892857  95.482140  173.32143     65534   65534
## 25351: 1.4571428  3.2285714  0.6000000  65.428574  122.00000     65534   65534
## 25352: 0.8693694  0.5855856  0.9729730   7.554054   14.17117     65535   65533

    ### Annotate cell types

        cell.type.annot <- list('Cell' = c(65534),
                                'Non-cell' = c(65535)
                                )

        cell.type.annot <- do.list.switch(cell.type.annot)
        names(cell.type.annot) <- c('Values', 'Annotated cell type')

        cell.dat <- do.add.cols(cell.dat, 'cell.type', cell.type.annot, 'Values')

## Step 1/3. Mapping data

## Step 2/3. Merging data

## Step 3/3. Returning data

        cell.dat

##           ROI   ID          x          y Area CD11b_Sm149 CD20_Dy161 CD3_Er170
##     1: ROI002    1 147.666667 500.166667    6  0.33333334 10.1666670 0.3333333
##     2: ROI002    2  97.000000 500.000000   16  0.37500000  1.3750000 0.5000000
##     3: ROI002    3 470.815789 499.710526   19  0.57894737  0.1052632 0.2631579
##     4: ROI002    4 310.500000 499.833333   12  0.08333334 60.5833321 0.2500000
##     5: ROI002    5 245.166667 499.833333   12  0.33333334 28.2500000 2.0000000
##    ---
## 25348: ROI012 2907  45.569767   2.081395   43  0.34883720  3.6744187 0.5581396
## 25349: ROI012 2908  36.675676   4.500000   74  0.27027026  0.8918919 1.9459460
## 25350: ROI012 2909  27.035714   2.446429   56  0.32142857  0.5357143 0.3571429
## 25351: ROI012 2910  16.928571   2.128571   35  0.45714286  0.9142857 1.3428571
## 25352: ROI012 2911   5.986486  10.918919  222  0.31081080  0.7567568 0.2027027
##        CD4_Gd156 CD45_Sm152 CD8a_Dy162 DNA1_Ir191 DNA3_Ir193 cell.type regions
##     1: 0.0000000  5.6666665  0.5000000  54.666668   98.50000     65534   65535
##     2: 1.7500000  1.8125000  0.8125000   9.750000   15.62500     65534   65534
##     3: 0.7368421  0.2631579  0.6842105   8.578947   16.31579     65534   65534
##     4: 0.1666667  7.0833335  1.5833334  50.833332  104.41666     65534   65535
##     5: 0.0000000 18.6666660  4.2500000  95.916664  171.58333     65534   65535
##    ---
## 25348: 0.5581396  2.0465117  0.9767442  52.651161   99.34884     65534   65534
## 25349: 0.4324324  7.4594593  3.1216216  61.621620  115.50000     65534   65534
## 25350: 1.1964285  1.3214285  1.0892857  95.482140  173.32143     65534   65534
## 25351: 1.4571428  3.2285714  0.6000000  65.428574  122.00000     65534   65534
## 25352: 0.8693694  0.5855856  0.9729730   7.554054   14.17117     65535   65533
##        Annotated cell type
##     1:                Cell
##     2:                Cell
##     3:                Cell
##     4:                Cell
##     5:                Cell
##    ---
## 25348:                Cell
## 25349:                Cell
## 25350:                Cell
## 25351:                Cell
## 25352:            Non-cell

    ### Annotate regions

        region.annot <- list('Background' = c(65533),
                                'Red pulp' = c(65534),
                                'White pup' = c(65535)
                                )

        region.annot <- do.list.switch(region.annot)
        names(region.annot) <- c('Values', 'Annotated region')

        cell.dat <- do.add.cols(cell.dat, 'regions', region.annot, 'Values')

## Step 1/3. Mapping data

## Step 2/3. Merging data

## Step 3/3. Returning data

        cell.dat

##           ROI   ID          x          y Area CD11b_Sm149 CD20_Dy161 CD3_Er170
##     1: ROI002    1 147.666667 500.166667    6  0.33333334 10.1666670 0.3333333
##     2: ROI002    2  97.000000 500.000000   16  0.37500000  1.3750000 0.5000000
##     3: ROI002    3 470.815789 499.710526   19  0.57894737  0.1052632 0.2631579
##     4: ROI002    4 310.500000 499.833333   12  0.08333334 60.5833321 0.2500000
##     5: ROI002    5 245.166667 499.833333   12  0.33333334 28.2500000 2.0000000
##    ---
## 25348: ROI012 2907  45.569767   2.081395   43  0.34883720  3.6744187 0.5581396
## 25349: ROI012 2908  36.675676   4.500000   74  0.27027026  0.8918919 1.9459460
## 25350: ROI012 2909  27.035714   2.446429   56  0.32142857  0.5357143 0.3571429
## 25351: ROI012 2910  16.928571   2.128571   35  0.45714286  0.9142857 1.3428571
## 25352: ROI012 2911   5.986486  10.918919  222  0.31081080  0.7567568 0.2027027
##        CD4_Gd156 CD45_Sm152 CD8a_Dy162 DNA1_Ir191 DNA3_Ir193 cell.type regions
##     1: 0.0000000  5.6666665  0.5000000  54.666668   98.50000     65534   65535
##     2: 1.7500000  1.8125000  0.8125000   9.750000   15.62500     65534   65534
##     3: 0.7368421  0.2631579  0.6842105   8.578947   16.31579     65534   65534
##     4: 0.1666667  7.0833335  1.5833334  50.833332  104.41666     65534   65535
##     5: 0.0000000 18.6666660  4.2500000  95.916664  171.58333     65534   65535
##    ---
## 25348: 0.5581396  2.0465117  0.9767442  52.651161   99.34884     65534   65534
## 25349: 0.4324324  7.4594593  3.1216216  61.621620  115.50000     65534   65534
## 25350: 1.1964285  1.3214285  1.0892857  95.482140  173.32143     65534   65534
## 25351: 1.4571428  3.2285714  0.6000000  65.428574  122.00000     65534   65534
## 25352: 0.8693694  0.5855856  0.9729730   7.554054   14.17117     65535   65533
##        Annotated cell type Annotated region
##     1:                Cell        White pup
##     2:                Cell         Red pulp
##     3:                Cell         Red pulp
##     4:                Cell        White pup
##     5:                Cell        White pup
##    ---
## 25348:                Cell         Red pulp
## 25349:                Cell         Red pulp
## 25350:                Cell         Red pulp
## 25351:                Cell         Red pulp
## 25352:            Non-cell       Background

Area calculation

###################################################################################
### Area calculation
###################################################################################       

        area.totals <- do.calculate.area(spatial.dat)
        area.totals

##       ROI    Total
## 1: ROI002 500.4998
## 2: ROI004 500.4998
## 3: ROI006 500.0000
## 4: ROI008 500.4998
## 5: ROI010 500.4998
## 6: ROI012 500.4998

        area.table <- do.calculate.area(spatial.dat, region = 'regions')
        area.table

##       ROI    65533    65534     65535
## 1: ROI002 105.8489 369.4577 320.61971
## 2: ROI004 113.2652 128.5535 470.26057
## 3: ROI006 126.2062 236.1123 422.28308
## 4: ROI008 184.9595 403.7499 230.81594
## 5: ROI010 157.2514 370.8005 297.11782
## 6: ROI012 266.4001 422.2606  35.02856

        for(i in region.annot[['Values']]){
          # i <- 65533

          nm <- region.annot[region.annot[['Values']] == i,'Values'][[1]]
          new.nm <- region.annot[region.annot[['Values']] == i,'Annotated region'][[1]]

          names(area.table)[which(names(area.table) == i)] <- new.nm

        }

        area.table

##       ROI Background Red pulp White pup
## 1: ROI002   105.8489 369.4577 320.61971
## 2: ROI004   113.2652 128.5535 470.26057
## 3: ROI006   126.2062 236.1123 422.28308
## 4: ROI008   184.9595 403.7499 230.81594
## 5: ROI010   157.2514 370.8005 297.11782
## 6: ROI012   266.4001 422.2606  35.02856

Save data

###################################################################################
### Save data
###################################################################################  

    ### Output QS and CSV file

        setwd(OutputDirectory)
        dir.create('Data')
        setwd('Data')

        qsave(spatial.dat, "spatial.dat.qs")
        fwrite(cell.dat, 'cell.dat.csv')
        fwrite(area.totals, 'area.totals.csv')
        fwrite(area.table, 'area.table.csv')

    ### Pull cellular data and write FCS file from each ROI independently

        setwd(OutputDirectory)
        dir.create('FCS files')
        setwd('FCS files')

        for(i in names(spatial.dat)){

          ## Extract data and setup cols

              tmp <- list()
              tmp[[i]] <- spatial.dat[[i]]

              cell.dat <- do.pull.data(tmp, 'CellData')
              cell.dat <- do.asinh(cell.dat, names(spatial.dat[[i]]@RASTERS), cofactor = 1)

          ### Invert y axis

              all.neg <- function(test) -1*abs(test)

              y_invert <- cell.dat[['y']]
              y_invert <- all.neg(y_invert)
              cell.dat[['y_invert']] <- y_invert

          ### Write FCS files

              write.files(cell.dat, i, write.csv = FALSE, write.fcs = TRUE)
              rm(cell.dat)
              rm(i)
        }

Script 2. Cellular analysis

---

Coming soon!

Script 3. Spatial analysis

---

Coming soon!