Key links:

Page contents:

← back to Spectre Home Page

Introduction

Introduction

Spectre's simple spatial analysis workflow is intended to provide a straightforward method of spatial analysis for high-dimensional imaging analysis. Critically, this workflow does not include spatially-classified cell phenotypes or tissue regions, but largely relies on the labelling of cell types by annotating clusters, and then assessing changes in these cells across images/ROIs. For more comprehensive spatial analysis approaches, including regional composition etc, see our advanced spatial analysis workflow.

Overview

This workflow is broken up into three stages, managed by three R scripts:

  1. Add masks
  2. Cellular analysis (clustering, dimensionality reduction)
  3. Spatial analysis

Following stage 1 (add masks), FCS files are also generated that can be used in FlowJo for analysis if desired.


Setup

For this workflow, you will need to use Spectre v0.5.0 or higher. For installation and update instructions, as well as introductory tutorials, see this page. Alternatively, you can check out our spatial analysis protocols in FlowJo.


Script 1: add masks and extract cellular data


Load packages and set directories

x

###################################################################################
### Spectre: spatial 1 - add masks and extract cellular data
###################################################################################
    ### Load libraries
    
        library('Spectre')
        Spectre::package.check(type = 'spatial')
        Spectre::package.load(type = 'spatial')
    ### Set PrimaryDirectory
    
        dirname(rstudioapi::getActiveDocumentContext()$path)            # Finds the directory where this script is located
        setwd(dirname(rstudioapi::getActiveDocumentContext()$path))     # Sets the working directory to where the script is located
        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 TIFFs

x

###################################################################################
### 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"
        tiff.list <- list()
    ### Check channel names
    
        for(i in rois){
          setwd(InputDirectory)
          setwd(i)
          tiff.list[[i]] <- list.files(getwd())
        }
        t(as.data.frame(tiff.list))
       [,1]               [,2]                 [,3]                  [,4]               [,5]              [,6]             [,7]             [,8]              [,9]              [,10]             [,11]                  [,12]            
ROI002 "140Ce_Ce140.tiff" "80ArAr_ArAr80.tiff" "alphaSMA_Pr141.tiff" "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff" "CD45_Sm152.tiff" "CD68_Tb159.tiff" "CD8a_Dy162.tiff" "CollagenI_Tm169.tiff" "DNA1_Ir191.tiff"
ROI004 "140Ce_Ce140.tiff" "80ArAr_ArAr80.tiff" "alphaSMA_Pr141.tiff" "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff" "CD45_Sm152.tiff" "CD68_Tb159.tiff" "CD8a_Dy162.tiff" "CollagenI_Tm169.tiff" "DNA1_Ir191.tiff"
ROI006 "140Ce_Ce140.tiff" "80ArAr_ArAr80.tiff" "alphaSMA_Pr141.tiff" "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff" "CD45_Sm152.tiff" "CD68_Tb159.tiff" "CD8a_Dy162.tiff" "CollagenI_Tm169.tiff" "DNA1_Ir191.tiff"
ROI008 "140Ce_Ce140.tiff" "80ArAr_ArAr80.tiff" "alphaSMA_Pr141.tiff" "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff" "CD45_Sm152.tiff" "CD68_Tb159.tiff" "CD8a_Dy162.tiff" "CollagenI_Tm169.tiff" "DNA1_Ir191.tiff"
ROI010 "140Ce_Ce140.tiff" "80ArAr_ArAr80.tiff" "alphaSMA_Pr141.tiff" "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff" "CD45_Sm152.tiff" "CD68_Tb159.tiff" "CD8a_Dy162.tiff" "CollagenI_Tm169.tiff" "DNA1_Ir191.tiff"
ROI012 "140Ce_Ce140.tiff" "80ArAr_ArAr80.tiff" "alphaSMA_Pr141.tiff" "CD11b_Sm149.tiff" "CD20_Dy161.tiff" "CD3_Er170.tiff" "CD4_Gd156.tiff" "CD45_Sm152.tiff" "CD68_Tb159.tiff" "CD8a_Dy162.tiff" "CollagenI_Tm169.tiff" "DNA1_Ir191.tiff"
       [,13]             [,14]                  [,15]           
ROI002 "DNA3_Ir193.tiff" "HistoneH3_Yb176.tiff" "VIM_Nd143.tiff"
ROI004 "DNA3_Ir193.tiff" "HistoneH3_Yb176.tiff" "VIM_Nd143.tiff"
ROI006 "DNA3_Ir193.tiff" "HistoneH3_Yb176.tiff" "VIM_Nd143.tiff"
ROI008 "DNA3_Ir193.tiff" "HistoneH3_Yb176.tiff" "VIM_Nd143.tiff"
ROI010 "DNA3_Ir193.tiff" "HistoneH3_Yb176.tiff" "VIM_Nd143.tiff"
ROI012 "DNA3_Ir193.tiff" "HistoneH3_Yb176.tiff" "VIM_Nd143.tiff"


Import channel TIFFs

x

###################################################################################
### Read in TIFF files and create spatial objects
###################################################################################   
    ### Read in ROI channel TIFFs
    
        setwd(InputDirectory)
        spatial.dat <- read.spatial.files(rois = rois, roi.loc = getwd())
    ### Check results
    
        str(spatial.dat, 3)
List of 6
 $ ROI002:List of 1
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
 $ ROI004:List of 1
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
 $ ROI006:List of 1
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
 $ ROI008:List of 1
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
 $ ROI010:List of 1
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
 $ ROI012:List of 1
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
        spatial.dat[[1]]$RASTERS
class      : RasterStack 
dimensions : 501, 500, 250500, 15  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : 0, 500, 0, 501  (xmin, xmax, ymin, ymax)
crs        : NA 
names      : X140Ce_Ce140f, X80ArAr_ArAr80f, alphaSMA_Pr141f, CD11b_Sm149f, CD20_Dy161f, CD3_Er170f, CD4_Gd156f, CD45_Sm152f, CD68_Tb159f, CD8a_Dy162f, CollagenI_Tm169f, DNA1_Ir191f, DNA3_Ir193f, HistoneH3_Yb176f, VIM_Nd143f 
min values :             0,            3399,               0,            0,           0,          0,          0,           0,           0,           0,                0,           0,           0,                0,          0 
max values :           443,            5297,            1231,           33,         779,         41,         40,         734,        1223,         109,              183,        1670,        3007,             3167,       7764


Import masks

x

###################################################################################
### 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_Cell_mask.tiff"     
 [2,] "ROI002_Cytoplasm_mask.tiff"
 [3,] "ROI002_Nuclei_mask.tiff"   
 [4,] "ROI004_Cell_mask.tiff"     
 [5,] "ROI004_Cytoplasm_mask.tiff"
 [6,] "ROI004_Nuclei_mask.tiff"   
 [7,] "ROI006_Cell_mask.tiff"     
 [8,] "ROI006_Cytoplasm_mask.tiff"
 [9,] "ROI006_Nuclei_mask.tiff"   
[10,] "ROI008_Cell_mask.tiff"     
[11,] "ROI008_Cytoplasm_mask.tiff"
[12,] "ROI008_Nuclei_mask.tiff"   
[13,] "ROI010_Cell_mask.tiff"     
[14,] "ROI010_Cytoplasm_mask.tiff"
[15,] "ROI010_Nuclei_mask.tiff"   
[16,] "ROI012_Cell_mask.tiff"     
[17,] "ROI012_Cytoplasm_mask.tiff"
[18,] "ROI012_Nuclei_mask.tiff" 
    ### Import CELL masks
        
        as.matrix(all.masks)
        cell.mask.ext <- '_Cell_mask.tiff'
        
        cell.masks <- list.files(pattern = cell.mask.ext)
        as.matrix(cell.masks)
     [,1]                   
[1,] "ROI002_Cell_mask.tiff"
[2,] "ROI004_Cell_mask.tiff"
[3,] "ROI006_Cell_mask.tiff"
[4,] "ROI008_Cell_mask.tiff"
[5,] "ROI010_Cell_mask.tiff"
[6,] "ROI012_Cell_mask.tiff"
        spatial.dat <- do.add.masks(spatial.dat = spatial.dat,
                                    mask.loc = getwd(),
                                    masks = cell.masks,
                                    mask.label = 'cell.mask',
                                    mask.ext = cell.mask.ext)
Adding masks to cell.mask
  -- processing ROI002
     ... correcting extent
     ... flipping Y
  -- processing ROI004
     ... correcting extent
     ... flipping Y
  -- processing ROI006
     ... correcting extent
     ... flipping Y
  -- processing ROI008
     ... correcting extent
     ... flipping Y
  -- processing ROI010
     ... correcting extent
     ... flipping Y
  -- processing ROI012
     ... correcting extent
     ... flipping Y
Returning spatial data object with added masks
    ### Review data
        
        str(spatial.dat, 3)
List of 6
 $ ROI002:List of 2
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 1
 $ ROI004:List of 2
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 1
 $ ROI006:List of 2
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 1
 $ ROI008:List of 2
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 1
 $ ROI010:List of 2
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 1
 $ ROI012:List of 2
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 1
        spatial.dat[[1]]$MASKS
$cell.mask
$cell.mask$maskraster
class      : RasterLayer 
dimensions : 501, 500, 250500  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
extent     : 0, 500, 0, 501  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : cell.mask 
values     : 0, 65535  (min, max)


Create polygons for cell

x

###################################################################################
### Generate polygons and outlines
###################################################################################
    ### Generate polygons and outlines
    
        spatial.dat <- do.create.outlines(spatial.dat, 'cell.mask')
Creating polygons, outlines, and centroids using 'stars' method.
Processing masks for ROI ROI002
although coordinates are longitude/latitude, st_union assumes that they are planar
although coordinates are longitude/latitude, st_union assumes that they are planar
although coordinates are longitude/latitude, st_union assumes that they are planar

...

... polygons complete
Regions defined for each Polygons
... outlines complete
... centroids complete
Returning spatial data
    ### Checks
    
        str(spatial.dat, 3)
List of 6
 $ ROI002:List of 2
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 4
 $ ROI004:List of 2
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 4
 $ ROI006:List of 2
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 4
 $ ROI008:List of 2
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 4
 $ ROI010:List of 2
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 4
 $ ROI012:List of 2
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 4
        spatial.dat[[1]]$MASKS
$cell.mask
$cell.mask$maskraster
class      : RasterLayer 
dimensions : 501, 500, 250500  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
extent     : 0, 500, 0, 501  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : cell.mask 
values     : 0, 65535  (min, max)


$cell.mask$polygons
class       : SpatialPolygonsDataFrame 
features    : 2750 
extent      : 0, 500, 0, 501  (xmin, xmax, ymin, ymax)
crs         : NA 
variables   : 1
names       : cell.mask 
min values  :         0 
max values  :     65535 

$cell.mask$outlines
    long lat order  hole piece id group
1     19 103     1 FALSE     1  1   1.1
2     19 104     2 FALSE     1  1   1.1
3     20 104     3 FALSE     1  1   1.1


Mask QC plots

x

###################################################################################
### Mask QC plots
###################################################################################      
    ### Mask plots
        
        setwd(OutputDirectory)
        dir.create('Plots - cell masks')
        setwd('Plots - cell masks')
		as.matrix(names(spatial.dat[[1]]$RASTERS))

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


Calculate 'per cell' data from masks

x

###################################################################################
### Extract 'per cell' data
################################################################################### 
    ### Extract cellular data for each cell mask
        
        spatial.dat <- do.extract(spatial.dat, 'cell.mask')
Processing ROI002
... X140Ce_Ce140f
|==============================================================================================================| 100%
... X80ArAr_ArAr80f
|==============================================================================================================| 100%
... alphaSMA_Pr141f
|==============================================================================================================| 100%

...
        str(spatial.dat, 3)
List of 6
 $ ROI002:List of 3
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 4
  ..$ DATA   :List of 1
  .. ..$ CellData:Classes ‘data.table’ and 'data.frame':	2750 obs. of  19 variables:
  .. .. ..- attr(*, ".internal.selfref")=<externalptr> 
 $ ROI004:List of 3
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 4
  ..$ DATA   :List of 1
  .. ..$ CellData:Classes ‘data.table’ and 'data.frame':	4205 obs. of  19 variables:
  .. .. ..- attr(*, ".internal.selfref")=<externalptr> 
 $ ROI006:List of 3
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 4
  ..$ DATA   :List of 1
  .. ..$ CellData:Classes ‘data.table’ and 'data.frame':	4200 obs. of  19 variables:
  .. .. ..- attr(*, ".internal.selfref")=<externalptr> 
 $ ROI008:List of 3
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 4
  ..$ DATA   :List of 1
  .. ..$ CellData:Classes ‘data.table’ and 'data.frame':	1685 obs. of  19 variables:
  .. .. ..- attr(*, ".internal.selfref")=<externalptr> 
 $ ROI010:List of 3
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 4
  ..$ DATA   :List of 1
  .. ..$ CellData:Classes ‘data.table’ and 'data.frame':	1301 obs. of  19 variables:
  .. .. ..- attr(*, ".internal.selfref")=<externalptr> 
 $ ROI012:List of 3
  ..$ RASTERS:Formal class 'RasterStack' [package "raster"] with 11 slots
  ..$ MASKS  :List of 1
  .. ..$ cell.mask:List of 4
  ..$ DATA   :List of 1
  .. ..$ CellData:Classes ‘data.table’ and 'data.frame':	1427 obs. of  19 variables:
  .. .. ..- attr(*, ".internal.selfref")=<externalptr>


Save spatial data object

x

###################################################################################
### Save spatial data object as qs file
###################################################################################         
    ### Output
        
        setwd(OutputDirectory)
        dir.create('QS file')
        setwd('QS file')
        
        qsave(spatial.dat, "spatial.dat.qs")    


Write FCS files

x

###################################################################################
### Write FCS files
###################################################################################       
    ### Pull cellular data
        
        cell.dat <- do.pull.data(spatial.dat, 'CellData')
        cell.dat
         ROI   ID         x          y   Area X140Ce_Ce140f X80ArAr_ArAr80f alphaSMA_Pr141f CD11b_Sm149f CD20_Dy161f	...
    1: ROI002    1 276.24140 218.403834 100982    0.11077221        4291.782        9.292290    0.3623715   6.8142638	...
    2: ROI002    2 222.54000   3.200000     50    0.06000000        4143.240        8.840000    0.2800000   1.9000000 	... 
    3: ROI002    3 247.44643   3.660714     56    0.03571429        4150.036        3.392857    0.3928571   0.4464286 	... 
    4: ROI002    4 333.27941   4.250000     68    0.08823530        4184.838        5.720588    0.3529412   1.2205882 	... 
    5: ROI002    5 388.37755   4.785714     49    0.08163265        4123.775        2.102041    0.1428571   0.5714286 	... 
   ---                                                                                                               
15564: ROI012 1423 216.41892 494.986486     37    0.05405406        4016.324       64.027023    0.2972973   0.2972973	... 
15565: ROI012 1424  59.82099 495.919753     81    0.16049382        4092.963        5.938272    0.3333333   1.0000000	... 
15566: ROI012 1425 432.20732 495.731707     82    0.08536585        4042.378        9.524390    0.6829268   0.8780488	... 
15567: ROI012 1426 397.91463 495.548780     41    0.00000000        4051.805        8.341463    0.3170732   0.5121951	... 
15568: ROI012 1427 420.92105 496.578947     38    0.13157895        4051.763        4.315790    0.2105263   1.1578947 	...
        as.matrix(names(cell.dat))
      [,1]              
 [1,] "ROI"             
 [2,] "ID"              
 [3,] "x"               
 [4,] "y"               
 [5,] "Area"            
 [6,] "X140Ce_Ce140f"   
 [7,] "X80ArAr_ArAr80f" 
 [8,] "alphaSMA_Pr141f" 
 [9,] "CD11b_Sm149f"    
[10,] "CD20_Dy161f"     
[11,] "CD3_Er170f"      
[12,] "CD4_Gd156f"      
[13,] "CD45_Sm152f"     
[14,] "CD68_Tb159f"     
[15,] "CD8a_Dy162f"     
[16,] "CollagenI_Tm169f"
[17,] "DNA1_Ir191f"     
[18,] "DNA3_Ir193f"     
[19,] "HistoneH3_Yb176f"
[20,] "VIM_Nd143f" 
        cellular.cols <- names(cell.dat)[c(6:20)]         
		as.matrix(cellular.cols)
      [,1]              
 [1,] "X140Ce_Ce140f"   
 [2,] "X80ArAr_ArAr80f" 
 [3,] "alphaSMA_Pr141f" 
 [4,] "CD11b_Sm149f"    
 [5,] "CD20_Dy161f"     
 [6,] "CD3_Er170f"      
 [7,] "CD4_Gd156f"      
 [8,] "CD45_Sm152f"     
 [9,] "CD68_Tb159f"     
[10,] "CD8a_Dy162f"     
[11,] "CollagenI_Tm169f"
[12,] "DNA1_Ir191f"     
[13,] "DNA3_Ir193f"     
[14,] "HistoneH3_Yb176f"
[15,] "VIM_Nd143f"
    ### Arcsinh transformation
    
        cell.dat <- do.asinh(cell.dat, cellular.cols, 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
        
        cell.dat
    ### Write FCS files  
        
        setwd(OutputDirectory)
        dir.create('FCS files')
        setwd('FCS files')
        
        write.files(cell.dat, 'spatial_all', write.csv = FALSE, write.fcs = TRUE)
        write.files(cell.dat, 'spatial', divide.by = 'ROI', write.csv = FALSE, write.fcs = TRUE)


Further plots

x

###################################################################################
### Further plots
###################################################################################       
    ### Plot settings
        
        as.matrix(names(spatial.dat))
        
        plot.rois <- c("ROI002", "ROI004")
        plot.rois
[1] "ROI002" "ROI004"
        as.matrix(names(spatial.dat[[1]]$RASTERS))
        exp.plot <- names(spatial.dat[[1]]$RASTERS)[1:15]
        exp.plot
 [1] "X140Ce_Ce140f"    "X80ArAr_ArAr80f"  "alphaSMA_Pr141f"  "CD11b_Sm149f"     "CD20_Dy161f"      "CD3_Er170f"      
 [7] "CD4_Gd156f"       "CD45_Sm152f"      "CD68_Tb159f"      "CD8a_Dy162f"      "CollagenI_Tm169f" "DNA1_Ir191f"     
[13] "DNA3_Ir193f"      "HistoneH3_Yb176f" "VIM_Nd143f"     
        as.matrix(names(spatial.dat[[1]]$RASTERS))
        factor.plot <- names(spatial.dat[[1]]$RASTERS)[0]
        factor.plot
character(0)
    ### Expression data point plots
        
        setwd(OutputDirectory)
        dir.create('Plots - expression')
        setwd('Plots - expression')
        
        for(i in plot.rois){
          for(a in exp.plot){
            make.spatial.plot(spatial.dat = spatial.dat, 
                              image.roi = i, 
                              image.channel = a, 
                              mask.outlines = 'cell.mask')
          }
        }

    ### Factor data point plots
        
        setwd(OutputDirectory)
        dir.create('Plots - data point expression')
        setwd('Plots - data point expression')
        
        for(i in plot.rois){
          for(a in exp.plot){
            make.spatial.plot(spatial.dat = spatial.dat, 
                              image.roi = i, 
                              image.channel = a, 
                              mask.outlines = 'cell.mask', 
                              cell.dat = 'CellData', 
                              cell.col = a)
          }
          
          for(a in factor.plot){
            make.spatial.plot(spatial.dat = spatial.dat, 
                              image.roi = i, 
                              image.channel = a, 
                              mask.outlines = 'cell.mask', 
                              cell.dat = 'CellData', 
                              cell.col = a, 
                              cell.col.type = 'factor')
          }
        }








Script 2: Cellular analysis

x

###################################################################################
### Spectre: spatial 2 - cellular analysis
###################################################################################
    ### Load libraries

        library('Spectre')

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

        dirname(rstudioapi::getActiveDocumentContext()$path)            # Finds the directory where this script is located
        setwd(dirname(rstudioapi::getActiveDocumentContext()$path))     # Sets the working directory to where the script is located
        getwd()
        PrimaryDirectory <- getwd()
        PrimaryDirectory
    ### Set InputDirectory

        setwd(PrimaryDirectory)
        setwd("Output 1 - add masks")
        InputDirectory <- getwd()
        InputDirectory
    ### Set metadata directory

        setwd(PrimaryDirectory)
        setwd("metadata")
        MetaDirectory <- getwd()
        MetaDirectory
    ### Create output directory

        setwd(PrimaryDirectory)
        dir.create("Output 2 - cellular analysis")
        setwd("Output 2 - cellular analysis")
        OutputDirectory <- getwd()
        OutputDirectory


###################################################################################
### Read in spatial data object
###################################################################################
    ### Read in spatial data object

        setwd(InputDirectory)
        setwd("QS file/")
        
        list.files(getwd(), '.qs')
        
        spatial.dat <- qread('spatial.dat.qs')
        str(spatial.dat, 3)
        
        spatial.dat[[1]]$DATA$CellData
    ### Metadata
        
        setwd(MetaDirectory)
        
        sample.meta <- fread("ROIs and samples.csv")
        sample.meta


###################################################################################
### Extract cellular data from images using masks and add annotations
###################################################################################
    ### Pull cellular data into a separate data.table

        cell.dat <- do.pull.data(spatial.dat, 'CellData')
        cell.dat
    ### Add sample metadata

        cell.dat <- do.add.cols(cell.dat, base.col = 'ROI', add.dat = sample.meta, add.by = 'ROI')
        cell.dat


###################################################################################
### Filtering
###################################################################################
    ### Filter
        
        #
        #
        #


###################################################################################
### Setup columns
###################################################################################
    ### Select 'cellular' columns and perform arcsinh transformation

        as.matrix(names(cell.dat))

        cellular.cols <- names(cell.dat)[c(6:20)]
        as.matrix(cellular.cols)
    ### Define key columns
        
        as.matrix(names(cell.dat))
        
        roi.col <- 'ROI'
        sample.col <- 'Sample'
        group.col <- 'Group'
        batch.col <- 'Batch'
    ### Asinh transformation

        cell.dat <- do.asinh(cell.dat, cellular.cols, cofactor = 1)
        cell.dat

        cellular.cols <- paste0(cellular.cols, "_asinh")
        as.matrix(cellular.cols)
    ### Re-scaling

        cell.dat <- do.rescale(cell.dat, cellular.cols)
        cell.dat

        cellular.cols <- paste0(cellular.cols, "_rescaled")
        as.matrix(cellular.cols)
    ### Define clustering cols

        as.matrix(names(cell.dat))

        cluster.cols <- names(cell.dat)[c(42:48)]
        as.matrix(cluster.cols)


###################################################################################
### Clustering and dimensionality reduction
###################################################################################
    setwd(OutputDirectory)
    ### Run FlowSOM

        cell.dat <- run.flowsom(cell.dat, cluster.cols, meta.k = 15)
        cell.dat
    ### Run DR

        cell.dat <- run.fitsne(cell.dat, cluster.cols, perplexity = 200)
        cell.dat
    ### Individual plots

        make.colour.plot(cell.dat, 'FItSNE_X', 'FItSNE_Y', roi.col, 'factor')
        make.colour.plot(cell.dat, 'FItSNE_X', 'FItSNE_Y', sample.col, 'factor')
        make.colour.plot(cell.dat, 'FItSNE_X', 'FItSNE_Y', group.col, 'factor')
        make.colour.plot(cell.dat, 'FItSNE_X', 'FItSNE_Y', batch.col, 'factor')
    
        make.colour.plot(cell.dat, 'FItSNE_X', 'FItSNE_Y', 'FlowSOM_metacluster', 'factor', add.label = TRUE)
    ### Multi plots

        make.multi.plot(cell.dat, 'FItSNE_X', 'FItSNE_Y', plot.by = cellular.cols, figure.title = 'Cellular cols')
        make.multi.plot(cell.dat, 'FItSNE_X', 'FItSNE_Y', plot.by = cluster.cols, col.type = 'factor', figure.title = 'Clustering cols')

        make.multi.plot(cell.dat, 'FItSNE_X', 'FItSNE_Y', 'FlowSOM_metacluster', 'ROI', col.type = 'factor', figure.title = 'FlowSOM_metacluster by ROI')
    ### Expression heatmap

        exp <- do.aggregate(cell.dat, cellular.cols, 'FlowSOM_metacluster', 'mean')
        make.pheatmap(exp, 'FlowSOM_metacluster', cellular.cols)
        rm(exp)


###################################################################################
### Annotate clusters
###################################################################################
    ### Cluster annotations
        
        cluster.annots <- list('CD4 T cells' = c(1,5),
                               'CD8 T cells' = c(6,2),
                               'B cells' = c(7,12,8,9,13,14,15,10,11))
        
        cluster.annots <- do.list.switch(cluster.annots)
        names(cluster.annots) <- c('Values', 'Annotated metacluster')
        cluster.annots
    ### Add annotations, fill in NAs
        
        cell.dat <- do.add.cols(cell.dat, 'FlowSOM_metacluster', cluster.annots, 'Values')
        cell.dat[['Annotated metacluster']][is.na(cell.dat[['Annotated metacluster']])] <- 'Other'
        cell.dat
    ### Extra plots
        
        make.colour.plot(cell.dat, 'FItSNE_X', 'FItSNE_Y', 'Annotated metacluster', 'factor', add.label = TRUE)
        make.multi.plot(cell.dat, 'FItSNE_X', 'FItSNE_Y', 'Annotated metacluster', 'ROI', col.type = 'factor', figure.title = 'Annotated metacluster by ROI')


###################################################################################
### Save data
###################################################################################
setwd(OutputDirectory)
    ### data.table
    
        fwrite(cell.dat, 'cell.dat.csv')
    ### FCS files
        
        setwd(OutputDirectory)
        dir.create('FCS')
        setwd('FCS')
        
        write.files(cell.dat, file.prefix = 'IMC_All', write.csv = FALSE, write.fcs = TRUE)
        write.files(cell.dat, file.prefix = 'IMC_', divide.by = sample.col, write.csv = FALSE, write.fcs = TRUE)
        


###################################################################################
### Spatial plotting of cluster data
###################################################################################
    ### Clusters
        
        setwd(OutputDirectory)
        dir.create('Clusters')
        setwd('Clusters')
        
            for(i in unique(cell.dat$ROI)){
                temp <- cell.dat[cell.dat[['ROI']] == i,]
                make.spatial.plot(spatial.dat, i, 'DNA1_Ir191', mask.outlines = 'cell.mask', temp, cell.col = 'FlowSOM_metacluster', cell.col.type = 'factor')
            }
    ### Cell type
        
        setwd(OutputDirectory)
        dir.create('Annotated')
        setwd('Annotated')
        
            for(i in unique(cell.dat$ROI)){
                temp <- cell.dat[cell.dat[['ROI']] == i,]
                make.spatial.plot(spatial.dat, i, 'DNA1_Ir191', mask.outlines = 'cell.mask', temp, cell.col = 'Annotated metacluster', cell.col.type = 'factor')
            }




Script 3: spatial analysis

x




  • No labels