For the YouTube livestream recording, see here
For screen-shot slides, click here
Welcome to Week 13 of the course. Over the past two weeks, we have seen how to isolate normalized fluorescent signatures from our unmixing controls. We have also explored ways to make this process more feasible (and less painful) to implement when working with all the fluorophores in our panel.
The focus of today’s walkthrough will be to explore the various ways by which we can compare these normalized fluorescent signatures to each other, and try to infer how our proposed fluorophore panels will perform when we factor in regions of signature overlap that are present between the fluorophores. This in turn can lead to improved panel designs (as we try to reduce the extent of the overlap to improve resolution) in future panels.
By the end of this primary walkthrough (and the bonus content), you will also hopefully have a better understanding about what goes into the various panel metrics displayed by various commercial vendors.
As always, lets start by loading to our local environment any R packages that we are likely to need today.
library(dplyr)
Attaching package: 'dplyr'The following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, unionlibrary(purrr)
library(stringr)
library(ggplot2)
library(Luciernaga)
library(tidyr)
library(tibble)We will start off by first working with the signature matrices for the bead and cell single-color controls from the SDY3080 dataset that we completed retrieving as part of the bonus Signature Matrices walk-through. These were saved as .csv files that can be found in this week’s data folder.
#StorageLocation <- file.path("course", "13_SpectralSimilarities", "data")
StorageLocation <- file.path("data")
#OutputLocation <- file.path("course", "13_SpectralSimilarities", "outputs")
OutputLocation <- file.path("outputs")Lets go ahead and load the signature matrix that we derived from the bead single-color controls.
BeadsPath <- list.files(StorageLocation, pattern="BeadInitialSignature", full.names=TRUE)
BeadData <- read.csv(BeadsPath, check.names=FALSE)
BeadData <- BeadData |> filter(!str_detect(Fluorophore, "Unstained"))
head(BeadData, 3) name
1 DTR_2023_ILT_00-Reference Group-DR_CCR4 BUV615 (Beads).1237096.fcs
2 DTR_2023_ILT_00-Reference Group-DR_CCR6 BV786 (Beads).1237097.fcs
3 DTR_2023_ILT_00-Reference Group-DR_CCR7 BV650 (Beads).1237098.fcs
Fluorophore Antigen Type Detector Negative UV1-A UV2-A UV3-A UV4-A
1 BUV615 CCR4 Beads UV10 BeadUnstained 0.007 0.031 0.025 0.021
2 BV786 CCR6 Beads V15 BeadUnstained 0.000 0.000 0.002 0.001
3 BV650 CCR7 Beads V11 BeadUnstained 0.000 0.000 0.003 0.002
UV5-A UV6-A UV7-A UV8-A UV9-A UV10-A UV11-A UV12-A UV13-A UV14-A UV15-A
1 0.021 0.023 0.016 0.007 0.207 1.000 0.397 0.211 0.121 0.107 0.066
2 0.002 0.002 0.001 0.001 0.001 0.000 0.000 0.000 0.002 0.042 0.155
3 0.002 0.002 0.002 0.001 0.001 0.032 0.169 0.078 0.056 0.045 0.026
UV16-A V1-A V2-A V3-A V4-A V5-A V6-A V7-A V8-A V9-A V10-A V11-A
1 0.035 0.000 0.001 0.001 0.001 0.002 0.002 0.003 0.018 0.103 0.281 0.094
2 0.126 0.053 0.054 0.051 0.022 0.010 0.006 0.006 0.003 0.002 0.003 0.003
3 0.015 0.064 0.067 0.062 0.026 0.011 0.007 0.006 0.005 0.024 0.300 1.000
V12-A V13-A V14-A V15-A V16-A B1-A B2-A B3-A B4-A B5-A B6-A B7-A B8-A
1 0.049 0.028 0.019 0.013 0.005 0 0 0.001 0.005 0.032 0.069 0.031 0.024
2 0.002 0.016 0.238 1.000 0.551 0 0 0.000 0.000 0.000 0.000 0.000 0.000
3 0.444 0.320 0.194 0.122 0.047 0 0 0.000 0.000 0.001 0.003 0.012 0.008
B9-A B10-A B11-A B12-A B13-A B14-A YG1-A YG2-A YG3-A YG4-A YG5-A YG6-A YG7-A
1 0.021 0.009 0.006 0.005 0.003 0.002 0.089 0.471 0.999 0.449 0.339 0.254 0.222
2 0.000 0.000 0.002 0.009 0.020 0.023 0.000 0.000 0.000 0.000 0.000 0.000 0.001
3 0.007 0.004 0.002 0.002 0.001 0.001 0.001 0.003 0.026 0.209 0.133 0.101 0.110
YG8-A YG9-A YG10-A R1-A R2-A R3-A R4-A R5-A R6-A R7-A R8-A
1 0.085 0.057 0.025 0.007 0.008 0.007 0.004 0.003 0.002 0.001 0.001
2 0.006 0.025 0.021 0.000 0.000 0.000 0.001 0.007 0.033 0.090 0.050
3 0.039 0.025 0.012 0.124 0.122 0.103 0.068 0.046 0.026 0.023 0.010And repeat the same process for the signature matrix retrieved from the cell-single color controls.
CellsPath <- list.files(StorageLocation, pattern="CellInitialSignature", full.names=TRUE)
CellsData <- read.csv(CellsPath, check.names=FALSE)
CellsData <- CellsData |> filter(!str_detect(Fluorophore, "Unstained"))
head(CellsData, 3) name
1 DTR_2023_ILT_01-Reference Group-DR_CCR4 BUV615 (Cells).1235678.fcs
2 DTR_2023_ILT_01-Reference Group-DR_CCR6 BV786 (Cells).1235679.fcs
3 DTR_2023_ILT_01-Reference Group-DR_CCR7 BV650 (Cells).1235680.fcs
Fluorophore Antigen Type Detector Negative UV1-A UV2-A UV3-A UV4-A
1 BUV615 CCR4 Cells UV10 PBMC_Unstained 0.010 0.044 0.039 0.038
2 BV786 CCR6 Cells V15 PBMC_Unstained 0.002 0.002 0.005 0.003
3 BV650 CCR7 Cells V11 PBMC_Unstained 0.000 0.001 0.004 0.003
UV5-A UV6-A UV7-A UV8-A UV9-A UV10-A UV11-A UV12-A UV13-A UV14-A UV15-A
1 0.041 0.057 0.080 0.052 0.268 1.000 0.406 0.223 0.128 0.120 0.076
2 0.004 0.004 0.003 0.002 0.003 0.003 0.002 0.001 0.003 0.044 0.154
3 0.004 0.004 0.005 0.006 0.009 0.039 0.166 0.081 0.058 0.048 0.029
UV16-A V1-A V2-A V3-A V4-A V5-A V6-A V7-A V8-A V9-A V10-A V11-A
1 0.041 0.019 0.067 0.123 0.130 0.164 0.144 0.186 0.153 0.196 0.380 0.158
2 0.127 0.064 0.074 0.078 0.048 0.043 0.036 0.045 0.032 0.023 0.028 0.016
3 0.019 0.065 0.074 0.075 0.041 0.030 0.027 0.039 0.037 0.053 0.326 1.000
V12-A V13-A V14-A V15-A V16-A B1-A B2-A B3-A B4-A B5-A B6-A B7-A B8-A
1 0.086 0.062 0.046 0.040 0.017 0.022 0.031 0.050 0.047 0.071 0.101 0.056 0.044
2 0.009 0.023 0.247 1.000 0.555 0.004 0.008 0.014 0.014 0.014 0.013 0.010 0.008
3 0.445 0.323 0.201 0.131 0.054 0.008 0.015 0.026 0.020 0.017 0.017 0.025 0.018
B9-A B10-A B11-A B12-A B13-A B14-A YG1-A YG2-A YG3-A YG4-A YG5-A YG6-A YG7-A
1 0.040 0.021 0.015 0.014 0.012 0.011 0.102 0.432 0.894 0.430 0.326 0.253 0.223
2 0.008 0.005 0.005 0.011 0.020 0.023 0.007 0.013 0.021 0.016 0.011 0.010 0.011
3 0.017 0.010 0.007 0.006 0.005 0.005 0.008 0.013 0.040 0.209 0.136 0.105 0.116
YG8-A YG9-A YG10-A R1-A R2-A R3-A R4-A R5-A R6-A R7-A R8-A
1 0.085 0.059 0.026 0.011 0.012 0.009 0.005 0.010 0.004 0.005 0.002
2 0.010 0.027 0.022 0.002 0.003 0.003 0.003 0.008 0.030 0.079 0.044
3 0.043 0.029 0.015 0.117 0.114 0.097 0.066 0.044 0.027 0.024 0.012Additionally, lets retrieve the reference signatures for our respective fluorophores as they are stored in Luciernaga() for easier comparison of the references vs. the signatures we ourselves derived.
OurFluorophores <- CellsData |> pull(Fluorophore)
Data <- Luciernaga:::InstrumentReferences(64) |> filter(Fluorophore %in% OurFluorophores)
ReferenceData <- Data |> select(-Instrument) |>
tidyr::pivot_wider(names_from="Detector", values_from="AdjustedY")
head(ReferenceData, 3)# A tibble: 3 × 65
Fluorophore UV1 UV2 UV3 UV4 UV5 UV6 UV7 UV8 UV9
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Alexa Fluor 6… 0 0 0 0 0 0 0 0 0
2 APC 0 2.01e-4 2.01e-4 2.98e-4 4.01e-4 6.02e-4 0 0 0
3 APC-Fire 750 0 0 0 0 0 0 0 0 0
# ℹ 55 more variables: UV10 <dbl>, UV11 <dbl>, UV12 <dbl>, UV13 <dbl>,
# UV14 <dbl>, UV15 <dbl>, UV16 <dbl>, V1 <dbl>, V2 <dbl>, V3 <dbl>, V4 <dbl>,
# V5 <dbl>, V6 <dbl>, V7 <dbl>, V8 <dbl>, V9 <dbl>, V10 <dbl>, V11 <dbl>,
# V12 <dbl>, V13 <dbl>, V14 <dbl>, V15 <dbl>, V16 <dbl>, B1 <dbl>, B2 <dbl>,
# B3 <dbl>, B4 <dbl>, B5 <dbl>, B6 <dbl>, B7 <dbl>, B8 <dbl>, B9 <dbl>,
# B10 <dbl>, B11 <dbl>, B12 <dbl>, B13 <dbl>, B14 <dbl>, YG1 <dbl>,
# YG2 <dbl>, YG3 <dbl>, YG4 <dbl>, YG5 <dbl>, YG6 <dbl>, YG7 <dbl>, …Now that we have our 3 signature matrices, to make our lives a little bit easier, lets rearrange the rows so that we see the fluorophores displayed in order of peak detector according to laser and wavelength. We can first modify the tidying code we had to retrieve the Luciernaga references to instead retrieve the peak detectors for each fluorophore.
# Determine Peak Detectors for rearranging later
PeakDetectors <- Data |> group_by(Fluorophore) |>
arrange(desc(AdjustedY)) |> slice(1) |> ungroup() |>
select(!c(Instrument, AdjustedY))
PeakDetectors# A tibble: 29 × 2
Fluorophore Detector
<chr> <chr>
1 APC R1
2 APC-Fire 750 R7
3 APC-Fire 810 R8
4 APC-R700 R4
5 Alexa Fluor 647 R2
6 BUV395 UV2
7 BUV496 UV7
8 BUV563 UV9
9 BUV615 UV10
10 BUV661 UV11
# ℹ 19 more rowsWe then need to relocate our previously written code containing the instrument detector (to avoid re-typing the vector of detector names in correct sequence by hand again). Alternatively, if you have a raw .fcs file, you could load it into a GatingSet, run colnames(), and copy and paste the output into the code-chunk, adding the “” and commas as needed to save on typing time.
DetectorOrder <- c("UV1", "UV2", "UV3", "UV4", "UV5", "UV6", "UV7", "UV8",
"UV9", "UV10", "UV11", "UV12", "UV13", "UV14", "UV15", "UV16",
"V1", "V2", "V3", "V4", "V5", "V6", "V7", "V8",
"V9", "V10", "V11","V12", "V13", "V14", "V15", "V16",
"B1", "B2", "B3", "B4", "B5", "B6", "B7", "B8",
"B9", "B10", "B11", "B12", "B13", "B14",
"YG1", "YG2", "YG3", "YG4", "YG5", "YG6", "YG7", "YG8", "YG9", "YG10",
"R1", "R2", "R3", "R4", "R5", "R6", "R7", "R8")Now that we have the vector with the correct sequence, we can factor() the Detector column peak detectors data.frame, setting levels equal to our vector. Once this is done, we can use rearrange() to place the rows in the correct order on the basis of our designated levels. Once this is done, we can pull() the Fluorophore column to retrieve a vector of Fluorophore names in their now corrected sequence across lasers and wavelengths (thanks to the rearrangement we carried out on the Detector column).
PeakDetectors$Detector <- factor(PeakDetectors$Detector, levels=DetectorOrder)
FluorophoreOrder <- PeakDetectors |> arrange(Detector) |>
unique() |> pull(Fluorophore)
FluorophoreOrder [1] "BUV395" "BUV496" "BUV563" "BUV615"
[5] "BUV661" "BUV737" "BUV805" "BV421"
[9] "Pacific Blue" "BV480" "BV510" "BV605"
[13] "BV650" "BV711" "BV750" "BV786"
[17] "FITC" "Spark Blue 550" "PerCP-Cy5.5" "PE"
[21] "PE-Dazzle 594" "PE-Cy5" "PE-Vio 770" "APC"
[25] "Alexa Fluor 647" "APC-R700" "Zombie NIR" "APC-Fire 750"
[29] "APC-Fire 810" With this vector of fluorophores in the correct sequence now retrieved, we just need to rearrange our individual data.frame objects to also match. We can factor() their individual Fluorophore columns, providing the Fluorophore vector as their levels argument, before following up with the dplyr arrange() function
BeadData$Fluorophore <- factor(BeadData$Fluorophore, levels=FluorophoreOrder)
BeadData <- BeadData |> arrange(Fluorophore)
head(BeadData,3) name
1 DTR_2023_ILT_00-Reference Group-DR_CD62L BUV395 (Beads).1237114.fcs
2 DTR_2023_ILT_00-Reference Group-DR_CD8 BUV496 (Beads).1237117.fcs
3 DTR_2023_ILT_00-Reference Group-DR_CD69 BUV563 (Beads).1237115.fcs
Fluorophore Antigen Type Detector Negative UV1-A UV2-A UV3-A UV4-A
1 BUV395 CD62L Beads UV2 BeadUnstained 0.258 1.000 0.494 0.360
2 BUV496 CD8 Beads UV7 BeadUnstained 0.009 0.039 0.032 0.027
3 BUV563 CD69 Beads UV9 BeadUnstained 0.007 0.029 0.023 0.020
UV5-A UV6-A UV7-A UV8-A UV9-A UV10-A UV11-A UV12-A UV13-A UV14-A UV15-A
1 0.325 0.315 0.177 0.056 0.022 0.004 0.001 0.001 0.000 0.001 0.000
2 0.036 0.187 1.000 0.534 0.256 0.065 0.017 0.007 0.004 0.004 0.002
3 0.019 0.020 0.017 0.327 1.000 0.264 0.062 0.025 0.015 0.011 0.007
UV16-A V1-A V2-A V3-A V4-A V5-A V6-A V7-A V8-A V9-A V10-A V11-A
1 0.000 0.002 0.003 0.003 0.002 0.002 0.002 0.002 0.002 0.001 0.001 0.001
2 0.001 0.001 0.002 0.012 0.066 0.216 0.164 0.170 0.069 0.036 0.029 0.008
3 0.004 0.000 0.001 0.001 0.001 0.001 0.002 0.027 0.073 0.035 0.033 0.008
V12-A V13-A V14-A V15-A V16-A B1-A B2-A B3-A B4-A B5-A B6-A B7-A B8-A
1 0.000 0.000 0.000 0.000 0 0.001 0.001 0.001 0.001 0.000 0.000 0.000 0.000
2 0.003 0.002 0.001 0.001 0 0.068 0.043 0.034 0.010 0.006 0.003 0.001 0.001
3 0.003 0.002 0.001 0.001 0 0.000 0.006 0.109 0.218 0.118 0.076 0.025 0.015
B9-A B10-A B11-A B12-A B13-A B14-A YG1-A YG2-A YG3-A YG4-A YG5-A YG6-A YG7-A
1 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
2 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
3 0.012 0.006 0.003 0.002 0.001 0.001 0.488 0.246 0.146 0.044 0.026 0.019 0.017
YG8-A YG9-A YG10-A R1-A R2-A R3-A R4-A R5-A R6-A R7-A R8-A
1 0.000 0.000 0.000 0 0 0 0 0 0 0 0
2 0.000 0.000 0.000 0 0 0 0 0 0 0 0
3 0.006 0.004 0.002 0 0 0 0 0 0 0 0CellsData$Fluorophore <- factor(CellsData$Fluorophore, levels=FluorophoreOrder)
CellsData <- CellsData |> arrange(Fluorophore)
head(CellsData,3) name
1 DTR_2023_ILT_01-Reference Group-DR_CD62L BUV395 (Cells).1235696.fcs
2 DTR_2023_ILT_01-Reference Group-DR_CD8 BUV496 (Cells).1235699.fcs
3 DTR_2023_ILT_01-Reference Group-DR_CD69 BUV563 (Cells).1235697.fcs
Fluorophore Antigen Type Detector Negative UV1-A UV2-A UV3-A
1 BUV395 CD62L Cells UV2 PBMC_Unstained 0.259 1.000 0.487
2 BUV496 CD8 Cells UV7 CBMC_Unstained 0.009 0.039 0.031
3 BUV563 CD69 Cells UV9 PBMC_Activated_Unstained 0.007 0.030 0.024
UV4-A UV5-A UV6-A UV7-A UV8-A UV9-A UV10-A UV11-A UV12-A UV13-A UV14-A UV15-A
1 0.357 0.321 0.316 0.198 0.073 0.041 0.015 0.008 0.006 0.005 0.006 0.007
2 0.027 0.036 0.186 1.000 0.531 0.253 0.064 0.018 0.008 0.005 0.005 0.003
3 0.021 0.021 0.023 0.020 0.342 1.000 0.260 0.061 0.026 0.015 0.012 0.007
UV16-A V1-A V2-A V3-A V4-A V5-A V6-A V7-A V8-A V9-A V10-A V11-A
1 0.005 0.010 0.028 0.048 0.049 0.068 0.066 0.093 0.078 0.060 0.071 0.041
2 0.002 0.001 0.004 0.017 0.072 0.223 0.170 0.180 0.077 0.043 0.038 0.014
3 0.004 0.000 -0.002 -0.003 -0.004 -0.004 -0.001 0.025 0.069 0.033 0.031 0.008
V12-A V13-A V14-A V15-A V16-A B1-A B2-A B3-A B4-A B5-A B6-A B7-A B8-A
1 0.025 0.023 0.021 0.019 0.010 0.018 0.029 0.042 0.036 0.032 0.031 0.023 0.018
2 0.007 0.005 0.004 0.003 0.001 0.070 0.047 0.039 0.014 0.010 0.007 0.004 0.003
3 0.003 0.002 0.001 0.001 0.000 0.000 0.008 0.107 0.201 0.108 0.070 0.024 0.015
B9-A B10-A B11-A B12-A B13-A B14-A YG1-A YG2-A YG3-A YG4-A YG5-A YG6-A YG7-A
1 0.016 0.012 0.009 0.009 0.007 0.008 0.015 0.023 0.027 0.024 0.015 0.016 0.016
2 0.003 0.002 0.001 0.001 0.001 0.001 0.002 0.004 0.006 0.006 0.004 0.004 0.004
3 0.012 0.006 0.003 0.003 0.002 0.002 0.436 0.212 0.126 0.039 0.022 0.017 0.014
YG8-A YG9-A YG10-A R1-A R2-A R3-A R4-A R5-A R6-A R7-A R8-A
1 0.008 0.004 0.002 0.003 0.004 0.002 0.000 0.005 0.003 0.004 0.002
2 0.002 0.001 0.001 0.002 0.003 0.002 0.002 0.002 0.001 0.001 0.001
3 0.004 0.003 0.001 -0.001 -0.001 -0.001 -0.001 0.000 0.000 0.000 0.000ReferenceData$Fluorophore <- factor(ReferenceData$Fluorophore, levels=FluorophoreOrder)
ReferenceData <- ReferenceData |> arrange(Fluorophore)
head(ReferenceData,3)# A tibble: 3 × 65
Fluorophore UV1 UV2 UV3 UV4 UV5 UV6 UV7 UV8 UV9
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 BUV395 0.284 1 0.522 0.368 0.339 0.318 0.181 0.0509 0.0158
2 BUV496 0.01 0.0372 0.0307 0.0268 0.0360 0.193 1 0.527 0.244
3 BUV563 0.01 0.0379 0.0307 0.0260 0.0259 0.0267 0.0238 0.367 1
# ℹ 55 more variables: UV10 <dbl>, UV11 <dbl>, UV12 <dbl>, UV13 <dbl>,
# UV14 <dbl>, UV15 <dbl>, UV16 <dbl>, V1 <dbl>, V2 <dbl>, V3 <dbl>, V4 <dbl>,
# V5 <dbl>, V6 <dbl>, V7 <dbl>, V8 <dbl>, V9 <dbl>, V10 <dbl>, V11 <dbl>,
# V12 <dbl>, V13 <dbl>, V14 <dbl>, V15 <dbl>, V16 <dbl>, B1 <dbl>, B2 <dbl>,
# B3 <dbl>, B4 <dbl>, B5 <dbl>, B6 <dbl>, B7 <dbl>, B8 <dbl>, B9 <dbl>,
# B10 <dbl>, B11 <dbl>, B12 <dbl>, B13 <dbl>, B14 <dbl>, YG1 <dbl>,
# YG2 <dbl>, YG3 <dbl>, YG4 <dbl>, YG5 <dbl>, YG6 <dbl>, YG7 <dbl>, …And with that, most of our data tidying steps are complete. Lets start comparing signatures!
Before going to far, lets refresh our memory for this Cytek Aurora 5-laser panel of how these fluorophores are distributed across the various lasers and wavelengths.
TableData <- CellsData |> select(Fluorophore, Antigen) |> mutate(Day=1)
rownames(TableData) <- NULL
Table <- FluorophoreMatrix(TableData, NumberDetectors = 64)Table| D | WV | Fluorophore | Antigen | D | Fluorophore | Antigen | D | Fluorophore | Antigen | D | Fluorophore | Antigen | D | Fluorophore | Antigen |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| UV1 | 373 | ||||||||||||||
| UV2 | 388 | BUV395 | CD62L | ||||||||||||
| UV3 | 428 | V1 | BV421 | CD127 | |||||||||||
| UV4 | 443 | V2 | |||||||||||||
| UV5 | 458 | V3 | Pacific Blue | Dump_CD14 | |||||||||||
| UV6 | 473 | V4 | |||||||||||||
| 508 | V5 | BV480 | CD161 | B1 | |||||||||||
| UV7 | 514 | BUV496 | CD8 | ||||||||||||
| 525 | V6 | B2 | FITC | CD3 | |||||||||||
| UV8 | 542 | V7 | BV510 | CD45RA | B3 | Spark Blue 550 | CD3 | ||||||||
| 577 | YG1 | PE | NKG2D | ||||||||||||
| 581 | V8 | B4 | |||||||||||||
| UV9 | 582 | BUV563 | CD69 | ||||||||||||
| 598 | V9 | B5 | YG2 | ||||||||||||
| UV10 | 613 | BUV615 | CCR4 | ||||||||||||
| 615 | V10 | BV605 | CD56 | B6 | YG3 | PE-Dazzle 594 | TNFa | ||||||||
| UV11 | 662 | BUV661 | VD2 | V11 | BV650 | CCR7 | B7 | YG4 | R1 | APC | CD16 | ||||
| 679 | B8 | YG5 | PE-Cy5 | CD25 | R2 | Alexa Fluor 647 | CD3 | ||||||||
| UV12 | 695 | V12 | B9 | PerCP-Cy5.5 | CD26 | YG6 | R3 | ||||||||
| UV13 | 718 | V13 | BV711 | CD7 | B10 | YG7 | R4 | APC-R700 | CD107a | ||||||
| 738 | B11 | R5 | |||||||||||||
| UV14 | 750 | BUV737 | CXCR3 | V14 | BV750 | IFNg | YG8 | ||||||||
| 760 | B12 | R6 | Zombie NIR | Viability | |||||||||||
| UV15 | 781 | V15 | BV786 | CCR6 | B13 | YG9 | PE-Vio 770 | PD1 | R7 | APC-Fire 750 | CD27 | ||||
| UV16 | 812 | BUV805 | CD4 | V16 | B14 | YG10 | R8 | APC-Fire 810 | CD38 |
One of the main reasons I like visualizing a panel in this way, is that it is helpful when determining how various combinations of fluorophores will behave together when placed in a panel.
Speaking generally, fluorophores typically have three common types of interactions
While all the above occur partly due to the nature of the flurophores themselves, we have to remember that the extent that they cause issues in our panel can be partially mediated depending on the co-expression pattern of your cell populations for the respective antigens. If you have fluorophores that interact, but they are found on two separate cell populations, then this issue will be minimized vs. if both markers are expressed on the same cell.
Let’s go ahead and use VisualizeSignatures() to see our retrieved signatures.
For the bead single-color controls, we get the following signatures
BeadData2 <- BeadData |>
select(-c("name", "Antigen", "Type", "Detector", "Negative"))
Plots <- VisualizeSignatures(data=BeadData2,
columnname="Fluorophore", Normalize = FALSE)
plotly::ggplotly(Plots)For the cell single-color controls, we get the following signatures
CellsData2 <- CellsData |>
select(-c("name", "Antigen", "Type", "Detector", "Negative"))
Plots <- VisualizeSignatures(data=CellsData2,
columnname="Fluorophore", Normalize = FALSE)
plotly::ggplotly(Plots)And for the reference signatures, we see the following
Plots <- VisualizeSignatures(data=ReferenceData,
columnname="Fluorophore", Normalize = FALSE)
plotly::ggplotly(Plots)The simplest way of comparing our retrieved signatures would be to visualize those we want to compare in the same plot. To do this for our datasets, let us carry out some of the required data tidying first, to avoid accidental confusion (or overwriting) when there are 3 identically named “APC” fluorophores in the same data.frame.
Within the ‘Fluorophore’ column of each data.frame, we can combine the existing name with an additional identifying character string using the paste0() function, before assigning back to the original (in the process overwriting it with the modified version)
BeadData2$Fluorophore <- paste0(BeadData2$Fluorophore, "_Beads")
CellsData2$Fluorophore <- paste0(CellsData2$Fluorophore, "_Cells")
ReferenceData$Fluorophore <- paste0(ReferenceData$Fluorophore, "_Reference")Separately, the ‘Reference’ data.frame has different column names for the detector columns compared to the Cell and Bead data.frames.
setdiff(colnames(CellsData), colnames(ReferenceData)) [1] "name" "Antigen" "Type" "Detector" "Negative" "UV1-A"
[7] "UV2-A" "UV3-A" "UV4-A" "UV5-A" "UV6-A" "UV7-A"
[13] "UV8-A" "UV9-A" "UV10-A" "UV11-A" "UV12-A" "UV13-A"
[19] "UV14-A" "UV15-A" "UV16-A" "V1-A" "V2-A" "V3-A"
[25] "V4-A" "V5-A" "V6-A" "V7-A" "V8-A" "V9-A"
[31] "V10-A" "V11-A" "V12-A" "V13-A" "V14-A" "V15-A"
[37] "V16-A" "B1-A" "B2-A" "B3-A" "B4-A" "B5-A"
[43] "B6-A" "B7-A" "B8-A" "B9-A" "B10-A" "B11-A"
[49] "B12-A" "B13-A" "B14-A" "YG1-A" "YG2-A" "YG3-A"
[55] "YG4-A" "YG5-A" "YG6-A" "YG7-A" "YG8-A" "YG9-A"
[61] "YG10-A" "R1-A" "R2-A" "R3-A" "R4-A" "R5-A"
[67] "R6-A" "R7-A" "R8-A" setdiff(colnames(ReferenceData), colnames(CellsData)) [1] "UV1" "UV2" "UV3" "UV4" "UV5" "UV6" "UV7" "UV8" "UV9" "UV10"
[11] "UV11" "UV12" "UV13" "UV14" "UV15" "UV16" "V1" "V2" "V3" "V4"
[21] "V5" "V6" "V7" "V8" "V9" "V10" "V11" "V12" "V13" "V14"
[31] "V15" "V16" "B1" "B2" "B3" "B4" "B5" "B6" "B7" "B8"
[41] "B9" "B10" "B11" "B12" "B13" "B14" "YG1" "YG2" "YG3" "YG4"
[51] "YG5" "YG6" "YG7" "YG8" "YG9" "YG10" "R1" "R2" "R3" "R4"
[61] "R5" "R6" "R7" "R8" So that they can match, we will need to add the “-A” to the colnames() of ReferenceData. From there, we can combine the rows of all three data.frames via either the rbind() or bind_rows() function, ending up with a large data.frame object.
colnames(ReferenceData)[2:ncol(ReferenceData)] <- paste0(colnames(ReferenceData)[2:ncol(ReferenceData)], "-A")
Dataset <- rbind(BeadData2, CellsData2, ReferenceData)Now that we have assembled Dataset, we can compare all three sets of signatures using VisualizeSignatures() and ggplotly() to open an interactive view.
Plots <- VisualizeSignatures(data=Dataset,
columnname="Fluorophore", Normalize = FALSE)
plotly::ggplotly(Plots)One of the nice things about ggplotly() is it allows you to export an image of what you are seeing directly, by clicking on the camera icon, and selecting the desired storage location.
As we start off, we can see that for BUV395, the reference and bead controls are fairly similar, although the cell control appears to have some autofluorescence residuals left over
By contrast, for PerCP-Cy5.5, all three are discordant in their own unique ways.
Whether the differences we see between our beads and the reference controls PerCP-Cy5.5 is due to acquisition on different instruments (different lasers, gains, etc.), or differences at the fluorophore manufacturer/batch level is something to be explored on another day.
Utilmately, this approach works really well when it comes to seeing the differences, and we are able to describe how they are different, so it can be both quantitative and qualitative. Hoeverer, it is clearly not a high-throughput approach. Ultimately, many of the metrics we use when comparing signatures are attempts at putting a number to these differences, allowing for less painful implementation in a typical computational workflow.
If you use a Cytek Aurora, you may be familiar with this screen in SpectroFlo, which displays a “similarity matrix”.
When being trained on the instrument, you may have also been taught that the values displayed are meant to describe how similar two fluorophores are to each other.
Behind the curtain, these values are just the cosine similarity values between the two vectors. If both fluorophore vectors are identical to each other, they will have a value of 1 (‘cos(theta)=1’) and be considered parallel. By contrast, if they are not similar, they will have a value of 0 (‘cos(theta)=0’), being perpendicular or orthagonal. It is fareasier to describe this by simply showing it on a plot.
cos_vals <- c(1, 0)
theta <- acos(cos_vals)
df <- data.frame(
label = paste0("cos = ", cos_vals),
angle_deg = round(theta * 180 / pi, 1),
x = sin(theta),
y = cos(theta)
)
ggplot(df) +
geom_segment(aes(x = 0, y = 0, xend = x, yend = y, color = label),
arrow = arrow(length = unit(0.25, "cm")),
linewidth = 1) +
coord_fixed(xlim = c(-0.1, 1), ylim = c(-0.1, 1.1)) +
geom_hline(yintercept = 0, color = "grey80") +
geom_vline(xintercept = 0, color = "grey80") +
labs(x = NULL, y = NULL, color = "Vector",
title = NULL) +
theme_bw()
So unlike the common misperception among new users, the values are not percentages. In actuality, it is fascinating seeing the variation, considering I was often told that 0.98 is the cutoff of “practically identical” when contrasting two fluorophores or autofluorescence signatures.
cos_vals <- c(1, 0.999, 0.99, 0.98, 0.97, 0.96, 0.95, 0.90, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1)
theta <- acos(cos_vals)
df <- data.frame(
label = paste0("cos = ", cos_vals),
angle_deg = round(theta * 180 / pi, 1),
x = sin(theta),
y = cos(theta)
)
ggplot(df) +
geom_segment(aes(x = 0, y = 0, xend = x, yend = y, color = label),
arrow = arrow(length = unit(0.25, "cm")),
linewidth = 1) +
coord_fixed(xlim = c(-0.1, 1), ylim = c(-0.1, 1.1)) +
geom_hline(yintercept = 0, color = "grey80") +
geom_vline(xintercept = 0, color = "grey80") +
labs(x = NULL, y = NULL, color = "Vector",
title = NULL) +
theme_bw()
With this bit of context about what these cosine similarity values are, let’s figure out how to calculate them for our own signatures. One thing we will need to facilitate this is to install the lsa package from the CRAN repository (via install.packages()), and then attach it to our local environment (via library())
install.packages("lsa")library(lsa)Loading required package: SnowballCWith this done, rather than jumping in to the deep end of the pool, lets start with the simpler of the examples. Go ahead and filter() for the ‘BUV395’ fluorophores from our Dataset object.
BUV395 <- Dataset |> filter(str_detect(Fluorophore, "BUV395"))
BUV395 Fluorophore UV1-A UV2-A UV3-A UV4-A UV5-A UV6-A
1 BUV395_Beads 0.2580000 1 0.4940000 0.3600000 0.3250000 0.3150000
2 BUV395_Cells 0.2590000 1 0.4870000 0.3570000 0.3210000 0.3160000
3 BUV395_Reference 0.2838991 1 0.5220986 0.3675975 0.3386984 0.3181995
UV7-A UV8-A UV9-A UV10-A UV11-A UV12-A UV13-A UV14-A
1 0.1770000 0.05600000 0.02200000 0.00400000 0.0010000000 0.001 0.000 0.001
2 0.1980000 0.07300000 0.04100000 0.01500000 0.0080000000 0.006 0.005 0.006
3 0.1811009 0.05090023 0.01580275 0.01060206 0.0002981651 0.000 0.000 0.000
UV15-A UV16-A V1-A V2-A V3-A V4-A V5-A V6-A V7-A
1 0.000 0.000 0.002000000 0.003000000 0.003000000 0.002 0.002 0.002 0.002
2 0.007 0.005 0.010000000 0.028000000 0.048000000 0.049 0.068 0.066 0.093
3 0.000 0.000 0.003102064 0.003400229 0.001198394 0.000 0.000 0.000 0.000
V8-A V9-A V10-A V11-A V12-A V13-A V14-A V15-A V16-A B1-A B2-A B3-A B4-A
1 0.002 0.001 0.001 0.001 0.000 0.000 0.000 0.000 0.00 0.001 0.001 0.001 0.001
2 0.078 0.060 0.071 0.041 0.025 0.023 0.021 0.019 0.01 0.018 0.029 0.042 0.036
3 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.00 0.000 0.000 0.000 0.000
B5-A B6-A B7-A B8-A B9-A B10-A B11-A B12-A B13-A B14-A YG1-A YG2-A YG3-A
1 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
2 0.032 0.031 0.023 0.018 0.016 0.012 0.009 0.009 0.007 0.008 0.015 0.023 0.027
3 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
YG4-A YG5-A YG6-A YG7-A YG8-A YG9-A YG10-A R1-A R2-A R3-A R4-A R5-A R6-A
1 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0 0.000 0.000
2 0.024 0.015 0.016 0.016 0.008 0.004 0.002 0.003 0.004 0.002 0 0.005 0.003
3 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0 0.000 0.000
R7-A R8-A
1 0.000 0.000
2 0.004 0.002
3 0.000 0.000For lsa, we will need to provide our input as just the columns containing the numeric values (not metadata, etc.). Consequently, lets separate the character from the numeric columns for now.
Names <- BUV395 |> select(!where(is.numeric))
Numbers <- BUV395 |> select(where(is.numeric))We will next need to transpose (using base R’s t() function) our rows into a columns (similar to a pivot_longer() operation)
NumbersTransposed <- t(Numbers)
head(NumbersTransposed, 4) [,1] [,2] [,3]
UV1-A 0.258 0.259 0.2838991
UV2-A 1.000 1.000 1.0000000
UV3-A 0.494 0.487 0.5220986
UV4-A 0.360 0.357 0.3675975We can then remove the rownames() (by setting them to NULL), and swap out the colnames() for the respective identity values.
TheNames <- Names |> pull(Fluorophore)
rownames(NumbersTransposed) <- NULL
colnames(NumbersTransposed) <- TheNames
head(NumbersTransposed, 4) BUV395_Beads BUV395_Cells BUV395_Reference
[1,] 0.258 0.259 0.2838991
[2,] 1.000 1.000 1.0000000
[3,] 0.494 0.487 0.5220986
[4,] 0.360 0.357 0.3675975At this point, we are now ready, lets go ahead and run cosine() on our matrix.
CosineValues <- cosine(NumbersTransposed)
CosineValues BUV395_Beads BUV395_Cells BUV395_Reference
BUV395_Beads 1.0000000 0.9850615 0.9996056
BUV395_Cells 0.9850615 1.0000000 0.9838872
BUV395_Reference 0.9996056 0.9838872 1.0000000And tada, we have created our first cosine matrix!
But let’s scale upfrom just three signatures to the level of an entire panel, let’s do this by looking at the signatures we extracted from the bead single-color unmixing controls.
JustBeads <- Dataset |> filter(str_detect(Fluorophore, "_Bead"))And rather than have to rewrite/copy-paste the code from above every time for every new matrix we want to check, we can assemble a small function to help us out.
#' Takes a data.frame, returns the Cosine Matrix
#'
#' @param data The data.frame object with a single character column
#'
CosineCalculation <- function(data){
Names <- data |> select(!where(is.numeric))
Numbers <- data |> select(where(is.numeric))
NumbersTransposed <- t(Numbers)
TheNames <- Names |> pull(Fluorophore)
rownames(NumbersTransposed) <- NULL
colnames(NumbersTransposed) <- TheNames
CosineValues <- cosine(NumbersTransposed)
CosineValues <- round(CosineValues, 3)
return(CosineValues)
}After running the code-block (and seeing the created function object appear in the Session window of our right secondary side-bar), we can run the function and check the output.
BeadCosines <- CosineCalculation(data=JustBeads)
head(BeadCosines, 3) BUV395_Beads BUV496_Beads BUV563_Beads BUV615_Beads BUV661_Beads
BUV395_Beads 1.000 0.220 0.065 0.035 0.021
BUV496_Beads 0.220 1.000 0.317 0.079 0.021
BUV563_Beads 0.065 0.317 1.000 0.398 0.059
BUV737_Beads BUV805_Beads BV421_Beads Pacific Blue_Beads
BUV395_Beads 0.028 0.099 0.072 0.020
BUV496_Beads 0.016 0.045 0.091 0.129
BUV563_Beads 0.018 0.022 0.013 0.012
BV480_Beads BV510_Beads BV605_Beads BV650_Beads BV711_Beads
BUV395_Beads 0.047 0.055 0.005 0.004 0.003
BUV496_Beads 0.424 0.536 0.062 0.023 0.010
BUV563_Beads 0.084 0.236 0.219 0.042 0.010
BV750_Beads BV786_Beads FITC_Beads Spark Blue 550_Beads
BUV395_Beads 0.002 0.002 0.005 0.005
BUV496_Beads 0.008 0.007 0.109 0.117
BUV563_Beads 0.005 0.004 0.092 0.205
PerCP-Cy5.5_Beads PE_Beads PE-Dazzle 594_Beads PE-Cy5_Beads
BUV395_Beads 0.000 0.001 0.000 0.000
BUV496_Beads 0.005 0.024 0.012 0.002
BUV563_Beads 0.029 0.532 0.278 0.059
PE-Vio 770_Beads APC_Beads Alexa Fluor 647_Beads APC-R700_Beads
BUV395_Beads 0.000 0.000 0.000 0.000
BUV496_Beads 0.001 0.003 0.000 0.001
BUV563_Beads 0.015 0.025 0.008 0.009
APC-Fire 750_Beads APC-Fire 810_Beads
BUV395_Beads 0.000 0.000
BUV496_Beads 0.001 0.003
BUV563_Beads 0.002 0.006Likewise, we can run the Cells signature matrix
JustCells <- Dataset |> filter(str_detect(Fluorophore, "_Cell"))
CellCosines <- CosineCalculation(data=JustCells)
head(CellCosines, 3) BUV395_Cells BUV496_Cells BUV563_Cells BUV615_Cells BUV661_Cells
BUV395_Cells 1.000 0.277 0.115 0.150 0.060
BUV496_Cells 0.277 1.000 0.330 0.197 0.050
BUV563_Cells 0.115 0.330 1.000 0.426 0.071
BUV737_Cells BUV805_Cells BV421_Cells Pacific Blue_Cells
BUV395_Cells 0.085 0.135 0.194 0.118
BUV496_Cells 0.067 0.076 0.228 0.173
BUV563_Cells 0.070 0.042 0.111 0.028
BV480_Cells BV510_Cells BV605_Cells BV650_Cells BV711_Cells
BUV395_Cells 0.182 0.199 0.134 0.084 0.051
BUV496_Cells 0.439 0.555 0.129 0.054 0.028
BUV563_Cells 0.103 0.241 0.231 0.061 0.018
BV750_Cells BV786_Cells FITC_Cells Spark Blue 550_Cells
BUV395_Cells 0.037 0.046 0.059 0.109
BUV496_Cells 0.017 0.032 0.120 0.154
BUV563_Cells 0.008 0.020 0.093 0.211
PerCP-Cy5.5_Cells PE_Cells PE-Dazzle 594_Cells PE-Cy5_Cells
BUV395_Cells 0.116 0.094 0.067 0.072
BUV496_Cells 0.135 0.101 0.027 0.051
BUV563_Cells 0.112 0.512 0.257 0.082
PE-Vio 770_Cells APC_Cells Alexa Fluor 647_Cells APC-R700_Cells
BUV395_Cells 0.023 0.029 0.018 0.021
BUV496_Cells 0.017 0.018 0.017 0.018
BUV563_Cells 0.050 0.026 0.014 0.015
Zombie NIR_Cells APC-Fire 750_Cells APC-Fire 810_Cells
BUV395_Cells 0.026 0.036 0.035
BUV496_Cells 0.020 0.041 0.037
BUV563_Cells 0.015 0.034 0.033And finally, we can run the reference matrix
JustReferences <- Dataset |> filter(str_detect(Fluorophore, "_Reference"))
ReferencesCosines <- CosineCalculation(data=JustReferences)
head(ReferencesCosines, 3) BUV395_Reference BUV496_Reference BUV563_Reference
BUV395_Reference 1.000 0.220 0.077
BUV496_Reference 0.220 1.000 0.336
BUV563_Reference 0.077 0.336 1.000
BUV615_Reference BUV661_Reference BUV737_Reference
BUV395_Reference 0.036 0.021 0.028
BUV496_Reference 0.077 0.020 0.014
BUV563_Reference 0.380 0.058 0.017
BUV805_Reference BV421_Reference Pacific Blue_Reference
BUV395_Reference 0.109 0.074 0.018
BUV496_Reference 0.045 0.085 0.114
BUV563_Reference 0.021 0.015 0.009
BV480_Reference BV510_Reference BV605_Reference
BUV395_Reference 0.048 0.054 0.003
BUV496_Reference 0.420 0.525 0.054
BUV563_Reference 0.097 0.257 0.205
BV650_Reference BV711_Reference BV750_Reference
BUV395_Reference 0.003 0.002 0.001
BUV496_Reference 0.015 0.004 0.003
BUV563_Reference 0.042 0.010 0.006
BV786_Reference FITC_Reference Spark Blue 550_Reference
BUV395_Reference 0.003 0.005 0.004
BUV496_Reference 0.007 0.103 0.107
BUV563_Reference 0.005 0.089 0.195
PerCP-Cy5.5_Reference PE_Reference PE-Dazzle 594_Reference
BUV395_Reference 0.001 0.001 0.001
BUV496_Reference 0.002 0.027 0.010
BUV563_Reference 0.028 0.486 0.241
PE-Cy5_Reference PE-Vio 770_Reference APC_Reference
BUV395_Reference 0.000 0.000 0.000
BUV496_Reference 0.001 0.000 0.001
BUV563_Reference 0.057 0.017 0.019
Alexa Fluor 647_Reference APC-R700_Reference
BUV395_Reference 0.000 0.002
BUV496_Reference 0.000 0.003
BUV563_Reference 0.005 0.009
Zombie NIR_Reference APC-Fire 750_Reference
BUV395_Reference 0.000 0.000
BUV496_Reference 0.006 0.000
BUV563_Reference 0.009 0.002
APC-Fire 810_Reference
BUV395_Reference 0.000
BUV496_Reference 0.000
BUV563_Reference 0.003Comparing between them, we can see how the little bits of residual autofluorescence in the cell controls increase the similarity with the autofluorescence similar fluorophores (BUV496, BV510, Spark Blue 550, etc.)
Now that we can extract the Cosine Similarity values, lets see about plotting them with ggplot2. One thing when generating matrix style plots with geom_tile(), you will frequently encounter the need to tidy the dataset first using tidyr’s pivot_longer() function. The example tidying code can be seen below.
Melted <- ReferencesCosines |> data.frame(check.names=FALSE) |>
tibble::rownames_to_column("Fluorophore1") |>
tidyr::pivot_longer(-Fluorophore1, names_to = "Fluorophore2", values_to = "Value")head(Melted, 5)# A tibble: 5 × 3
Fluorophore1 Fluorophore2 Value
<chr> <chr> <dbl>
1 BUV395_Reference BUV395_Reference 1
2 BUV395_Reference BUV496_Reference 0.22
3 BUV395_Reference BUV563_Reference 0.077
4 BUV395_Reference BUV615_Reference 0.036
5 BUV395_Reference BUV661_Reference 0.021To avoid accidental reordering of the columns on the basis of alphabetical name, lets extract the Fluorophore vector order via pull() and unique(), before calling factor() on the two long columns (“Fluorophore1”, “Fluorophore2”), setting the Fluorophore vector as the respective levels.
CorrectFluorOrder <- Melted |> pull(Fluorophore1) |> unique()
Melted$Fluorophore1 <- factor(Melted$Fluorophore1, levels=CorrectFluorOrder)
Melted$Fluorophore2 <- factor(Melted$Fluorophore2, levels=CorrectFluorOrder)With this tidying and factoring set, we can then proceed to run the ggplot2() code. This mainly consist of the geom_tile() layer, as well as various additional layers defining our final appearance of the plot. Feel free to tinker and adjust the various pallettes to your own preference!
Plot <- ggplot(Melted, aes(Fluorophore1, Fluorophore2, fill = Value)) +
geom_tile(color = "white") +
scale_fill_gradient(low = "lightblue", high = "orange", limit = c(0,1),
space = "Lab", name="Cosine") +
theme_bw() +
geom_text(aes(label = round(Value, 2)), color = "black", size = 2) +
coord_fixed(ratio = 1) +
theme(axis.title.x = element_blank(), axis.title.y = element_blank(),
panel.grid.major = element_blank(), panel.border = element_blank(),
panel.background = element_blank(), axis.ticks = element_blank(),
legend.direction = "vertical",
axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, size = 6),
axis.text.y = element_text(size = 6),
legend.key.size = unit(0.4, "cm"))Plot
Additionally, one of the constants throughout today will be when displaying matrices, we are kind of constrained by the plots viewers within Positron/Rstudio on the display settings. It is often easier just to save out as a .png or .jpg and avoid the issue entirely.
SaveAsThis <- file.path(OutputLocation, "CosineSimilarityMatrix.png")
ggsave(SaveAsThis, units ="in", width = 10, height = 10)
Alternatively, rather than plotting according to the fluorophore order, we can rearrange based on how similar the fluorophores are to each other. This relies on this smaller helper function for the classification.
ReorderedCosine <- function(CosineMatrix){
Day <- as.dist((1-CosineMatrix)/2)
Night <- hclust(Day)
Twilight <- CosineMatrix[Night$order, Night$order]
}To be extra chaotic (and demonstrate the helper functions full potential), lets use the full Dataset containing Beads, Cells and Reference Signatures. We can then pipe through CosineCalculation() and ReorderedCosine() before repeating the factoring
DatasetCosines <- CosineCalculation(data=Dataset)
Melted <- DatasetCosines |> ReorderedCosine() |>
data.frame(check.names=FALSE) |>
tibble::rownames_to_column("Fluorophore1") |>
tidyr::pivot_longer(-Fluorophore1,
names_to = "Fluorophore2", values_to = "Value")Granted, all the fluorophores in a single small plot are likely to not be interpretable, so lets throw it into plotly ggplotly() function and make use of both the zoom and pan options to navigate
Plot <- ggplot(Melted, aes(Fluorophore1, Fluorophore2, fill = Value)) +
geom_tile(color = "white") +
scale_fill_gradient(low = "lightblue", high = "orange", limit = c(0,1),
space = "Lab", name="Cosine") +
theme_bw() +
geom_text(aes(label = round(Value, 2)), color = "black", size = 2) +
coord_fixed(ratio = 1) +
theme(axis.title.x = element_blank(), axis.title.y = element_blank(),
panel.grid.major = element_blank(), panel.border = element_blank(),
panel.background = element_blank(), axis.ticks = element_blank(),
legend.direction = "vertical",
axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, size = 6),
axis.text.y = element_text(size = 6),
legend.key.size = unit(0.4, "cm"))plotly::ggplotly(Plot)
All in all, we have a better idea of what these cosine similarity matrices are depicting, namely, comparing two fluorophores to each other, and then reordering to form a wider picture.
One of the challenges working with larger spectral flow cytometry panels (>30 fluorophores), is that with each additional fluorophore added into the unmixing matrix, additional uncertainty is present at the time of the unmixing calculation. This is in part due to the collinearity that is present between the individual reference signatures.
While spillover-spreading matrices (which we will cover in the bonus material alongside staining index reduction matrices) are focused on spread for events positive for a fluorophore on the other fluorophore channels, until recently there wasn’t a good method to measure spread on the negative population (which also impacts the staining index values).
This was very nicely covered in the recent Cytometry Part A paper by Peter Mage, Andrew Konecny and Florian Mair. In this paper, they also described how to get at a new metric to identify unmixing-dependent hotspots (UDS) that are susceptible due to their collinearity.
Since they thankfully included the math in their paper, we can take advantage of existing library of math functions in R to implement the calculation. Likewise, we can use ggplot2 to generate a visual. One thing to remember is that since these matrices can get quite large, it might be better to save them out as an image or as a .pdf rather than try to zoom in and out of the Plots window in the right secondary side-bar to see clearly.
To get started, we are starting off using one of the cosine values matrices we generated above. To make sure any edits made don’ affect the original, lets save it as new object (‘Similarity’). Also while here, lets save the rownames() containing the fluorophore names as a vector, allowing us to append them back in later on.
Similarity <- BeadCosines
TheseFluorophores <- rownames(Similarity)At this point, we can carry out the math section (see the original paper for details on why)
PseudoInverse <- MASS::ginv(Similarity)
Absolute <- abs(PseudoInverse)
TheSQRT <- sqrt(Absolute)At this point, our data ends up in a matrix that looks like this.
head(TheSQRT,3) [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 1.03995096 0.5380363 0.02107977 0.1853547 0.1840496 0.16277020 0.33763667
[2,] 0.53803632 1.3217082 0.68128795 0.2278586 0.1605087 0.08762240 0.17076990
[3,] 0.02107977 0.6812880 1.41754771 0.9909944 0.5553859 0.05244391 0.02270526
[,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 0.4413300 0.4344207 0.3961918 0.4747274 0.2228414 0.1180457 0.06666983
[2,] 0.3757432 0.4066160 0.4598203 1.1160348 0.6603164 0.2721735 0.13520823
[3,] 0.3552901 0.5465031 0.9396733 0.8998147 0.5752141 0.1783780 0.11321141
[,15] [,16] [,17] [,18] [,19] [,20] [,21]
[1,] 0.14362987 0.18500528 0.22877211 0.2021512 0.07497643 0.03462293 0.2096037
[2,] 0.13256218 0.12289875 0.56905192 0.5709505 0.11059585 0.45727247 0.4398728
[3,] 0.08468453 0.05143187 0.08256397 0.3547783 0.27701311 1.01381700 0.6939659
[,22] [,23] [,24] [,25] [,26] [,27] [,28]
[1,] 0.07395246 0.02545620 0.1478010 0.01823616 0.0743032 0.05601450 0.16745018
[2,] 0.10393920 0.06867923 0.1062398 0.21465973 0.1471791 0.01504604 0.07983378
[3,] 0.19858641 0.09368217 0.2899106 0.35495115 0.1705844 0.07623208 0.05886265So to make interpretability easier, lets reappend the vector of Fluorophore names to the original rownames position. We can also round the values to make interpretation a bit easier in the larger plot.
row.names(TheSQRT) <- TheseFluorophores
Hotspots <- round(TheSQRT, 2)head(Hotspots, 3) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
BUV395_Beads 1.04 0.54 0.02 0.19 0.18 0.16 0.34 0.44 0.43 0.40 0.47 0.22
BUV496_Beads 0.54 1.32 0.68 0.23 0.16 0.09 0.17 0.38 0.41 0.46 1.12 0.66
BUV563_Beads 0.02 0.68 1.42 0.99 0.56 0.05 0.02 0.36 0.55 0.94 0.90 0.58
[,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23]
BUV395_Beads 0.12 0.07 0.14 0.19 0.23 0.20 0.07 0.03 0.21 0.07 0.03
BUV496_Beads 0.27 0.14 0.13 0.12 0.57 0.57 0.11 0.46 0.44 0.10 0.07
BUV563_Beads 0.18 0.11 0.08 0.05 0.08 0.35 0.28 1.01 0.69 0.20 0.09
[,24] [,25] [,26] [,27] [,28]
BUV395_Beads 0.15 0.02 0.07 0.06 0.17
BUV496_Beads 0.11 0.21 0.15 0.02 0.08
BUV563_Beads 0.29 0.35 0.17 0.08 0.06For the unmixing-dependent spread matrix, the same comparison values are mirrored across the diagonal as triangles. We can use base R notation to set the upper half of the triangle values to NA, so that they show up empty in the final plot.
final_matrix <- Hotspots
lower_tri <- final_matrix
lower_tri[upper.tri(final_matrix)] <- NA At this point, if we wanted, we could save the file out as a .csv for later use. For easier reloading back into R in the future, I’d recommend converting the matrix back to a data.frame object, copy the row.names() to the columns, and switch the rownames to a proper column using the tibble packages rownames_to_column() function.
ReturnTri <- data.frame(lower_tri)
colnames(ReturnTri) <- row.names(ReturnTri)
ReturnTri <- ReturnTri |>
tibble::rownames_to_column("Fluorophore")
# write.csv(ReturnTri, "UnmixingDependentHotspot.csv", row.names=FALSE)As we may start anticipating, when working with ggplot2 geom_tile() we often need to pivot_longer() our existing data, so lets proceed through these tidying steps
Longer <- as.data.frame(lower_tri) |>
rownames_to_column("row") |>
pivot_longer(cols = -row, names_to = "col_num",
names_prefix = "V",values_to = "value") |>
mutate(row = factor(row, levels=rev(TheseFluorophores)),
col = as.numeric(col_num),
col_label = factor(TheseFluorophores[col],
levels = TheseFluorophores)) |>
filter(!is.na(value)) |> arrange(row)Which results in a “longer” data.frame that would look as follows.
head(Longer, 3)# A tibble: 3 × 5
row col_num value col col_label
<fct> <chr> <dbl> <dbl> <fct>
1 APC-Fire 810_Beads 1 0.17 1 BUV395_Beads
2 APC-Fire 810_Beads 2 0.08 2 BUV496_Beads
3 APC-Fire 810_Beads 3 0.06 3 BUV563_BeadsNow that we have a “longer” data.frame, we can pass through to ggplot2 to create the geom_tile() style plot.
Plot <- ggplot(Longer, aes(x = col_label, y = row, fill = value)) +
geom_tile(color = "white", linewidth = 0.3) +
geom_text(aes(label = round(value, 2)),
color = ifelse(Longer$value > max(Longer$value)/2, "white", "black"),
size = 3, fontface = "bold") +
scale_fill_gradient(low = "white",high = "red2",
na.value = "white",name="Hotspot") +
scale_x_discrete(position = "bottom") + labs(x = NULL, y = NULL) +
coord_fixed() + theme_minimal(base_size = 12) +
theme(panel.grid = element_blank(),legend.position = "right",
axis.text.x = element_text(angle = 45,hjust = 1,vjust = 1,),
axis.text.y = element_text(hjust = 1))When we try to run it, the dimensions may not look quite ideal.
Plot
Which is why often, I will save directly to a .png or .jpg via ggsave(), and just change the dimensions as needed depending on the total number of fluorophores.
SaveAsThis <- file.path(OutputLocation, "UnmixingDependentSpreadMatrix.png")
ggsave(SaveAsThis, units ="in", width = 10, height = 10)
For interpretation, values on the diagonal correspond to the individual fluorophore itself, with the greater values being those more impacted by collinearity. The values off the diagonal (ex. APC vs Alexa Fluor 647) designate those that are interacting and contributing towards that neighborhood hotspot.
We will deep-dive the nuances more in the bonus material and next week when evaluating our unmixing results.
Wrapping up the “simpler to derive metrics”, we have one additional value that the Cytek Aurora users see just hanging around below the cosine complexity matrix.
This “Complexity Index” basically refers to the condition number of the reference matrix. The rough interpretation basically being, the higher the number, the more likely any little bit of variance is going to cause additional spread and additional issues.
The way I have explained the concept to those starting off is basically think of it in terms of the most complicated interaction of two fluorophores in your panel, divided by the least complicated interaction of two fluorophores in your panel.
Fortunately, once we have our signature matrix, we can get at this directly via the kappa() function, remembering to set ‘exact’ argument to TRUE.
ReferenceMatrix <- ReferenceData |> select(where(is.numeric)) |> as.matrix()
#nrow(CellsMatrix)
kappa(ReferenceMatrix, exact=TRUE)[1] 10.88719Lets check the equivalent for beads
BeadsMatrix <- BeadData2 |> select(where(is.numeric)) |> as.matrix()
#nrow(CellsMatrix)
kappa(BeadsMatrix, exact=TRUE)[1] 11.17886Lets check the equivalent for cells
CellsMatrix <- CellsData2 |> select(where(is.numeric)) |> as.matrix()
#nrow(CellsMatrix)
kappa(CellsMatrix, exact=TRUE)[1] 14.80896As we can see, our bead signatures are comparable to the references, although the cell signatures with the occasional fluorophore with additional bits of autofluorescence residuals adds a bit more complication to the overall.
The challenge in interpretation, namely, because its the most complicated divided by the least complicated interaction leaves a lot of middle ground, which is why panels with substantially more fluorophores that are well spaced can end up having lower values than smaller panels with heavily overlapped markers.
And with that, in the interest of weekly consistency, we will close for now.
In this section, we were able to take a look at several of the common metrics used to assess fluorescent signatures and their interactions for spectral flow cytometry. We have seen how these can be calculated in R, and explored some of their underlying components.
As I have prefaced in the past, I am not a mathematician/statistician, nor do I pretend to be one. If you are so inclined, I gladly encourage you to dive deeper, and come back and update this walk-through so that everyone can benefit from a more-thorough explanation of the underlying math that is currently abstracted away within the R functions we rely on for the calculations.
Additionally, two other common metrics (Spillover Spreading Matrices and Staining Index Reduction) are a little more involved (as they involve unmixed samples), so these will be covered separately as part of a bonus walkthrough.
Next up, having derived the fluorescent signatures, quality checked them and made decisions about what looks reasonable, in the next primary walk-through, we will take these and actually go unmix the full-stained samples for the SDY3080 dataset. In the process, we will crack open the black box behind many of our workflows, and figure out some ways to efficiently profile our unmixed .fcs files to determine whether something is off.
We assembled “Dataset”, consisting of the three signature matrices. Use dplyr and stringr to filter for any “APC” fluorophore. Pass this smaller dataset through the various steps from today (visualize, cosine, UDS, kappa), saving images or taking screenshots as appropiate. Provide a summary blurb about the interactions you might encounter if you designed this “all-Red” panel.
Randomly subset the rows of Reference matrix into 4 smaller matrices, consisting of the same number of fluorophores (just with different fluorophores). Pass each through kappa() and see how the rank value changes. Summarize/Elaborate on your findings as a short blurb, does the value track with fluorophore overlap present?
In our unmixing-dependent spread (UDS) hotspot, we used only our panel fluorophores as inputs. But what if you wanted to evaluate spread for this panel if 1 or 2 autofluorescences signatures were present? Retrieve the unstained signature(s) from ‘Week 12’, and modify/re-run the cosine and unmixing-dependent spread matrices. After saving the plots as images or taking screenshots, briefly summarize what fluorophores would have been most impacted had you originally also accounted for autofluorescences.
