Spectre - simple discovery workflow

---

<- Back to Spectre home page

Introduction

---

Overview

Here we provide a worked example of a ‘simple’ discovery analysis workflow, where the entire process (data prep, clustering, dimensionality reduction, cluster annotation, plotting, summary data, and statistical analysis) is contained within a single script. The ‘simple’ workflow is most suitable for fast analysis of small datasets. For larger or more complex datasets, or datasets with multiple batches, we recommend the general discovery workflow, where the data preparation, batch alignment, clustering/dimensionality reduction, and quantitative analysis are separated into separate scripts. The demo dataset used for this worked example are cells extracted from mock- or virally-infected mouse brains, measured by flow cytometry.

Strategy

The ‘simple’ and ‘general’ discovery workflows are designed to facilitate the analysis of large and complex cytometry datasets using the Spectre R package. We’ve tested up to 30 million cells in a single analysis session so far. The workflow is designed to get around the cell number limitations of tSNE/UMAP. The analysis starts with clustering with FlowSOM – which is fast and scales well to large datasets. The clustered data is then downsampled, and dimensionality reduction is performed with tSNE/UMAP. This allows for visualisation of the data, and the clusters present in the dataset. Once the possible cell types in the datasets have been explored, the clusters can be labelled with the appropriate cellular identities. Finally, we can use the clusters/populations to generate summary statistics (expression levels, frequencies, total counts etc), which allows us to create graphs and heatmaps, facilitating statistical analysis.

Multiple samples

To analyse multiple samples, all the files must be imported into the one analysis session. This allows cells from each session to be clustered and analysed together, and allows us to examine the differential expression of markers, or the differences in cell proportions, between experimental groups.

Batch alignment

The ‘simple’ discovery workflow does not include any batch alignment steps. If batch correction needs to be applied, we recommend using the general discovery workflow.

---

---

Before you start

---

If you haven’t installed Spectre, please visit our Spectre installation page. If you are new to R and Spectre, we recommend trying out the R/RStudio and Spectre tutorials available on our getting st arted page, to familiarise yourself with R/RStudio first.

Citation and methods

---

Citation

If you use Spectre in your work, please consider citing Ashhurst TM, Marsh-Wakefield F, Putri GH et al. (2022). Cytometry Part A 101 (3), 237-253. To continue providing open-source tools such as Spectre, it helps us if we can demonstrate that our efforts are contributing to analysis efforts in the community. Please also consider citing the authors of the individual packages or tools (e.g. CytoNorm, FlowSOM, tSNE, UMAP, etc) that are critical elements of your analysis work. We have provided some generic text that you can use for your methods section with each protocol and on the ‘about’ page.

Sample methods blurb

Here is a sample methods blurb for this workflow. You may need to adapt this text to reflect any changes made in your analysis.

Computational analysis of data was performed using the Spectre R package (Ashhurst et al., 2022), with instructions and source code provided at https://github.com/ImmuneDynamics/spectre. Samples were initially prepared in FlowJo, and the population of interest was exported as raw value CSV files. Arcsinh transformation was performed on the data in R using a co-factor of 15 to redistribute the data on a linear scale and compress low end values near zero. The dataset was then merged into a single data.table, with keywords denoting the sample, group, and other factors added to each row (cell). The FlowSOM algorithm (Van Gassen et al., 2015) was then run on the merged dataset to cluster the data, where every cell is assigned to a specific cluster and metacluster. Subsequently, the data was downsampled and analysed by the dimensionality reduction algorithm Uniform Manifold Approximation and Projection (UMAP) (McInnes, Healy, Melville, 2018) for cellular visualisation.

Setup

---

Directories

Create a master folder with a meaningful name. Then inside that folder, insert the following:

• One folder called ‘data’ – this will contain your data CSV or FCS files
• One folder called ‘metadata’ – this will contain a CSV containg your sample metadata
• One folder called ‘Spectre simple discovery’ or similar – place this analysis script there

Example:

  CNS analysis
    /data
        -- Contains data files, one CSV or FCS per sample
    /metadata
        -- Contains a CSV containing sample metadata (group, batch, etc)
    /Spectre simple discovery
        -- Spectre simple.discovery.R

Analysis script

You can download the simple discovery script from this link – place this inside the Spectre simple discovery folder.

Data files

If you would to use the demo data as a test run for the simple discovery workflow, nothing to do at this step. Simply follow the relevant instructions further down this page to download the demo data (under 2. Import and prep data).

If you would like to use your own data, add your data and metadata files:

• Place the sample CSV or FCS files in the data folder you created above.
• Place the metadata CSV files in the metadata folder you created above.

Please see this page for detailed instructions on exporting data for Spectre and setting up a metadata file.

1. Load packages and set directories

---

#######################################################################################################
#### 1. Load packages, and set working directory
#######################################################################################################

Running library(Spectre) will load the Spectre package (also known as a ‘library’). We can then use package.check() to see if the standard dependency packages are installed, and package.load() to load those packages.

    ### Load libraries

        library(Spectre)
        Spectre::package.check()    # Check that all required packages are installed
        Spectre::package.load()     # Load required packages

Here we can set our ‘primary’ directory, which is going to be the location where the R script is saved. This path will be stored as PrimaryDirectory.

Note: if you aren’t sure how to navigate directories in R, check out our brief introduction to R tutorial.

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

We can then set our input directory, which will be the ‘data’ folder where we placed our files during setup. To do this we ask R to start at the location of the PrimaryDirectory, go up one level .., and then find the data folder.

    ### Set 'input' directory
        
        setwd(PrimaryDirectory)
        dir.create('../data', showWarnings = FALSE)
        setwd("../data/")
        InputDirectory <- getwd()
        setwd(PrimaryDirectory)

We need to set the location of the ‘metadata’ folder. This is where we can store a CSV file that contains any relevant metadata that we want to embed in our samples. In this example, it is located in a sub-folder called ‘metadata’.

    ### Set 'metadata' directory
        
        setwd(PrimaryDirectory)
        dir.create('../metadata', showWarnings = FALSE)
        setwd("../metadata/")
        MetaDirectory <- getwd()
        setwd(PrimaryDirectory)

We need to create a folder where our output data can go once our analysis is finished. In this example we will call this ‘Output_Spectre’.

    ### Create output directory
        
        setwd(PrimaryDirectory)
        dir.create("Output_Spectre", showWarnings = FALSE)
        setwd("Output_Spectre")
        OutputDirectory <- getwd()
        setwd(PrimaryDirectory)

2. Import and prep data

---

#######################################################################################################
#### 2. Import and prep data
#######################################################################################################

Demo data

If you need the demo dataset, uncomment the following code in the analysis script (select all, CMD+SHIFT+C) and run to download. If you are using your own datasets, then skip this step.

This code will download the demo dataset files and metadata file, and place them in the data and metadata folders respectively.

# setwd(PrimaryDirectory)
# setwd("../")
# getwd()
# download.file(url = "https://github.com/ImmuneDynamics/data/blob/main/msCNS.zip?raw=TRUE", destfile = 'msCNS.zip', mode = 'wb')
# unzip(zipfile = 'msCNS.zip')
# for(i in list.files('msCNS/data', full.names = TRUE)){
#   file.rename(from = i,  to = gsub('msCNS/', '', i))
# }
# for(i in list.files('msCNS/metadata', full.names = TRUE)){
#   file.rename(from = i,  to = gsub('msCNS/', '', i))
# }
# unlink(c('msCNS/', 'msCNS.zip', '__MACOSX'), recursive = TRUE)

Import data

To begin, we will change our working directory to ‘InputDirectory’ and list all the CSV files in that directory – these should be the sample CSV files. We can then read in all of our samples (in this example, one CSV file per sample) into a list called ‘data.list’. Spectre uses the data.table framework to store data, which reads, writes, and performs operations on data very quickly.

    ### Import data

        setwd(InputDirectory)
        list.files(InputDirectory, ".csv")

##  [1] "CNS_Mock_01.csv"   "CNS_Mock_02.csv"   "CNS_Mock_03.csv"  
##  [4] "CNS_Mock_04.csv"   "CNS_Mock_05.csv"   "CNS_Mock_06.csv"  
##  [7] "CNS_WNV_D7_01.csv" "CNS_WNV_D7_02.csv" "CNS_WNV_D7_03.csv"
## [10] "CNS_WNV_D7_04.csv" "CNS_WNV_D7_05.csv" "CNS_WNV_D7_06.csv"

        data.list <- Spectre::read.files(file.loc = InputDirectory,
                                         file.type = ".csv",
                                         do.embed.file.names = TRUE)

By default, the read.files() function will generate some other variables, which you can review, by running the do.list.summary() function.

The ‘name.table’ variable is a table of all the column names for all of your samples (one row per sample, one column per column name). If all of the column names are matching, then this table should be a repeating pattern. If it has been jumbled, then some of your samples have columns that don’t appear in other samples. The ‘ncol.check’ and ‘nrow.check’ are simple tables indicating the number or columns and rows in each sample.

    ### Check the data

        check <- do.list.summary(data.list)

        check$name.table # Review column names and their subsequent values

##      X1  X2   X3   X4    X5   X6   X7   X8  X9      X10    X11
## 1  NK11 CD3 CD45 Ly6G CD11b B220 CD8a Ly6C CD4 FileName FileNo
## 2  NK11 CD3 CD45 Ly6G CD11b B220 CD8a Ly6C CD4 FileName FileNo
## 3  NK11 CD3 CD45 Ly6G CD11b B220 CD8a Ly6C CD4 FileName FileNo
## 4  NK11 CD3 CD45 Ly6G CD11b B220 CD8a Ly6C CD4 FileName FileNo
## 5  NK11 CD3 CD45 Ly6G CD11b B220 CD8a Ly6C CD4 FileName FileNo
## 6  NK11 CD3 CD45 Ly6G CD11b B220 CD8a Ly6C CD4 FileName FileNo
## 7  NK11 CD3 CD45 Ly6G CD11b B220 CD8a Ly6C CD4 FileName FileNo
## 8  NK11 CD3 CD45 Ly6G CD11b B220 CD8a Ly6C CD4 FileName FileNo
## 9  NK11 CD3 CD45 Ly6G CD11b B220 CD8a Ly6C CD4 FileName FileNo
## 10 NK11 CD3 CD45 Ly6G CD11b B220 CD8a Ly6C CD4 FileName FileNo
## 11 NK11 CD3 CD45 Ly6G CD11b B220 CD8a Ly6C CD4 FileName FileNo
## 12 NK11 CD3 CD45 Ly6G CD11b B220 CD8a Ly6C CD4 FileName FileNo

        check$ncol.check # Review number of columns (features, markers) in each sample

##       [,1]
##  [1,] 11  
##  [2,] 11  
##  [3,] 11  
##  [4,] 11  
##  [5,] 11  
##  [6,] 11  
##  [7,] 11  
##  [8,] 11  
##  [9,] 11  
## [10,] 11  
## [11,] 11  
## [12,] 11

        check$nrow.check # Review number of rows (cells) in each sample

##       [,1] 
##  [1,] 9937 
##  [2,] 15415
##  [3,] 14246
##  [4,] 17044
##  [5,] 5459 
##  [6,] 4891 
##  [7,] 17950
##  [8,] 16233
##  [9,] 15999
## [10,] 17131
## [11,] 15926
## [12,] 18773

You can review the first 6 rows of the first sample in your data using the following:

        data.list[[1]]

##            NK11        CD3     CD45       Ly6G     CD11b      B220      CD8a
##           <num>      <num>    <num>      <num>     <num>     <num>     <num>
##    1:   42.3719  40.098700  6885.08  -344.7830 14787.300  -40.2399  83.71750
##    2:   42.9586 119.014000  1780.29  -429.6650  5665.730   86.6673  34.72190
##    3:   59.2366 206.238000 10248.30 -1603.8400 19894.300  427.8310 285.88000
##    4:  364.9480  -0.233878  3740.04  -815.9800  9509.430  182.4200 333.60500
##    5:  440.2470  40.035200  9191.38    40.5055  5745.820 -211.6940 149.22000
##   ---                                                                       
## 9933:   11.2126  36.951600  2515.82  -647.4930  6172.070  221.9380 266.90000
## 9934:  239.9700 440.217000  7247.28 -1449.7200 15355.400  809.3040 456.59900
## 9935: -134.9650 111.350000  2472.85    81.5975  9657.160 -113.1320   3.79607
## 9936:   86.3333  28.286900  5745.27 -1284.0800 18303.100  353.5290 262.96300
## 9937:   10.1467 122.255000  1971.69  -215.7660   727.708  506.8580 113.14400
##            Ly6C      CD4    FileName FileNo
##           <num>    <num>      <char>  <int>
##    1:  958.7000  711.072 CNS_Mock_01      1
##    2:  448.2590  307.272 CNS_Mock_01      1
##    3: 1008.8300  707.094 CNS_Mock_01      1
##    4:  440.0710  249.784 CNS_Mock_01      1
##    5:   87.4815  867.570 CNS_Mock_01      1
##   ---                                      
## 9933:  141.4200  708.348 CNS_Mock_01      1
## 9934: 2093.6900 2119.270 CNS_Mock_01      1
## 9935: -114.1510  110.743 CNS_Mock_01      1
## 9936:  745.8080  537.750 CNS_Mock_01      1
## 9937:  244.2210 2334.800 CNS_Mock_01      1

Merge data.tables

Once the metadata has been added, we can then merge the data into a single data.table using do.merge.files(). By default, columns with matching names will be aligned in the new table, and any columns that are present in some samples, but not others, will be added and filled with ‘NA’ for any samples that didn’t have that column initially. Once the data has been merged, we can review the data:

    ### Merge data

        cell.dat <- Spectre::do.merge.files(dat = data.list)
        cell.dat

##              NK11        CD3     CD45       Ly6G    CD11b      B220      CD8a
##             <num>      <num>    <num>      <num>    <num>     <num>     <num>
##      1:   42.3719  40.098700  6885.08  -344.7830 14787.30  -40.2399   83.7175
##      2:   42.9586 119.014000  1780.29  -429.6650  5665.73   86.6673   34.7219
##      3:   59.2366 206.238000 10248.30 -1603.8400 19894.30  427.8310  285.8800
##      4:  364.9480  -0.233878  3740.04  -815.9800  9509.43  182.4200  333.6050
##      5:  440.2470  40.035200  9191.38    40.5055  5745.82 -211.6940  149.2200
##     ---                                                                      
## 169000:  910.8890  72.856100 31466.20  -316.5570 28467.80   -7.7972 -271.8040
## 169001:  -10.2642  64.188700 45188.00  -540.5140 22734.00  202.4110 -936.4920
## 169002: -184.2910  -9.445650 11842.60   -97.9383 17237.00  123.4760 -219.9320
## 169003:  248.3860 229.986000 32288.20  -681.1630 19255.80 -656.0540 -201.5880
## 169004:  738.9810  95.470300 46185.10 -1004.6000 22957.80 -661.6280   72.3356
##               Ly6C       CD4      FileName FileNo
##              <num>     <num>        <char>  <int>
##      1:   958.7000  711.0720   CNS_Mock_01      1
##      2:   448.2590  307.2720   CNS_Mock_01      1
##      3:  1008.8300  707.0940   CNS_Mock_01      1
##      4:   440.0710  249.7840   CNS_Mock_01      1
##      5:    87.4815  867.5700   CNS_Mock_01      1
##     ---                                          
## 169000: 12023.7000 1103.0500 CNS_WNV_D7_06     12
## 169001:  4188.3300  315.9400 CNS_WNV_D7_06     12
## 169002:  8923.4000 -453.4640 CNS_WNV_D7_06     12
## 169003: 10365.7000   61.6765 CNS_WNV_D7_06     12
## 169004:  9704.4700  -31.8532 CNS_WNV_D7_06     12

Read in metadata

    ### Read in metadata  
       
        setwd(MetaDirectory)
        
        meta.dat <- fread("sample.details.csv")
        meta.dat

##              Filename     Sample  Group  Batch Cells per sample
##                <char>     <char> <char> <char>            <num>
##  1:   CNS_Mock_01.csv 01_Mock_01   Mock      A           420000
##  2:   CNS_Mock_02.csv 02_Mock_02   Mock      B           240000
##  3:   CNS_Mock_03.csv 03_Mock_03   Mock      B           256000
##  4:   CNS_Mock_04.csv 04_Mock_04   Mock      A           252000
##  5:   CNS_Mock_05.csv 05_Mock_05   Mock      A           345000
##  6:   CNS_Mock_06.csv 06_Mock_06   Mock      B           702000
##  7: CNS_WNV_D7_01.csv  07_WNV_01    WNV      A          5070000
##  8: CNS_WNV_D7_02.csv  08_WNV_02    WNV      B          2940000
##  9: CNS_WNV_D7_03.csv  09_WNV_03    WNV      A          2120000
## 10: CNS_WNV_D7_04.csv  10_WNV_04    WNV      A          4320000
## 11: CNS_WNV_D7_05.csv  11_WNV_05    WNV      B          4080000
## 12: CNS_WNV_D7_06.csv  12_WNV_06    WNV      A          1830000

3. Data transformation

---

#######################################################################################################
#### 3. Data transformation
#######################################################################################################

Before we perform clustering etc, we need to meaningfully transform the data. For more information on why this is necessary, please see this page.

Note: If you have imported CSV (channel value) files exported from FlowJo, then no data transformations are required, and you can skip all of the arcsinh transformation steps and proceed straight to adding the metadata. More information on the FCS, CSV scale, and CSV channel value file types can be found here.

    setwd(OutputDirectory)
    dir.create("Output 1 - transformed plots")
    setwd("Output 1 - transformed plots")

First, check the column names of the dataset.

    ### Arcsinh transformation

        as.matrix(names(cell.dat))

##       [,1]      
##  [1,] "NK11"    
##  [2,] "CD3"     
##  [3,] "CD45"    
##  [4,] "Ly6G"    
##  [5,] "CD11b"   
##  [6,] "B220"    
##  [7,] "CD8a"    
##  [8,] "Ly6C"    
##  [9,] "CD4"     
## [10,] "FileName"
## [11,] "FileNo"

The columns we want to apply arcsinh transformation to are the cellular columns – column 1 to column 9. We can specify those columns using the code below.

    ### Arcsinh transformation

        as.matrix(names(cell.dat))

##       [,1]      
##  [1,] "NK11"    
##  [2,] "CD3"     
##  [3,] "CD45"    
##  [4,] "Ly6G"    
##  [5,] "CD11b"   
##  [6,] "B220"    
##  [7,] "CD8a"    
##  [8,] "Ly6C"    
##  [9,] "CD4"     
## [10,] "FileName"
## [11,] "FileNo"

        to.asinh <- names(cell.dat)[c(1:9)]
        to.asinh

## [1] "NK11"  "CD3"   "CD45"  "Ly6G"  "CD11b" "B220"  "CD8a"  "Ly6C"  "CD4"

Define the cofactor we will use for transformation. As a general recommendation, we suggest using cofactor = 15 for CyTOF data, and cofactor between 100 and 1000 for flow data (we suggest 500 as a starting point). Here is a quick comparison figure showing how different co-factors compare to bi-exponential transformations performed on an LSR-II. For more detailed information on this choice, and for approaches where different cofactors for different columns might be required, see this page.