14 - Unmixing in R
2026-08-05
For the YouTube livestream recordings, see here
Background
.
Welcome to the fourteenth week of the Cytometry in R course. In this walk-through, we will gather the various coding knowledge building blocks we have been tinkering with since Week 09, and attempt to unmix our raw full-stained .fcs files using the spectral signature matrices that we derived from our unmixing controls (both single-color and unstained). We will also make some necessary formatting changes to the unmixed .fcs files, to avoid visualizing our data as tiny distant dots if they are later opened using commercial software.
.
More importantly, we hope to crack open what to many in the commmunity is seen as a black box. By understanding the underlying mechanics operating behind the curtain, we can gain a greater appreciation of the elements involved in the unmixing process (and hopefully assist with contextualizing/troubleshooting any unmixing errors we might encounter in the future).
.
The goal of this primary course walkthrough is focused on providing broad detailed look at the unmixing process. Additional details on unmixed .fcs file re-formatting are provided in a separate bonus walk-through for those that are interested. In the upcoming week’s primary walkthrough, we will cover how to systematically evaluate our unmixed .fcs files for unmixing errors, and make decisions on what signature variants to use for unmixing reattempts, etc.
.
Additionally, there are two separate community walk-throughs (for AutoSpectral and TRU-OLS respectively) focusing on these alternative unmixing approaches already available.
.
And with that, lets get started!
Walk Through
.
For unmixing, we work with raw full-stained .fcs files, so we will need to retrieve these from the SDY3080 ImmPort repository dataset that we have been using these last couple weeks. Since these original files can be massive in size, these were downsampled to allow for their sharing via GitHub without running into file size limit. These downsampled .fcs files can be found in this week’s “data” folder.
.
In addition to our full-stained samples, we will need the signature matrices we derived from our unmixing controls (both single-colors and unstained) from the previous walkthroughs. An important thing to note, we have not fully validated these signatures, so any unmixed .fcs files from today will likely have some unmixing errors. Since we plan to follow up in the next session on ways to evaluate how the different signatures included impacted our unmixing results, using these intermediate matrices should be fine for now.
.
As always, you are also welcome to follow along using your own .fcs files. Main thing to watch for is that choice of keywords may be different from those shown if your .fcs files originated from a different manufacturers spectral flow cytometer. Similarly, the detector column name endings (ending in “-A” for this dataset) may also differ (showing as “-H”, “-W”, none at all, etc.).
.
Consequently, small changes may need to be carried out to the existing code to account for these differences. If you get stuck, or want to share your solutions for different manufacturers instruments so that others can benefit, please post on the Discussions page, and we can help you navigate for your own respective instruments configuration.
.
And with that, let’s begin this journey into the blackhole that is spectral unmixing.
Set Up
.
As always, we can start by attaching the various R packages that we will need throughout today to our local environment. Since we have already derived the signature matrices and will not be creating new gates, the number of required packages is fewer compared to other weeks.
.
With this done, we can designate the file.path() to the location of both our storage and output folders.
.
With the storage location of our files established, we can then identify the .csv and .fcs files included in this week’s dataset using list.files().
[1] "data/DTR_2023_ILT_01-INF052-Ctrl_Antibody.1235515.fcs"
[2] "data/DTR_2023_ILT_01-INF052-PMA_Antibody.1235651.fcs"
[3] "data/DTR_2023_ILT_01-ND050_v1-Ctrl_Antibody.1235673.fcs"
[4] "data/DTR_2023_ILT_01-ND050_v1-PMA_Antibody.1235676.fcs" .
Lets go ahead and load in the individual .csv files containing the signature matrices for both Beads and Cells using read.csv(), saving each to its own object/variable. In the process, we can remove some of the unecessary metadata columns leftover from processing, since in addition to the detector columns, we will only need the Fluorophore and Antigen information.
[1] "Fluorophore" "Antigen" "UV1-A" "UV2-A" "UV3-A"
[6] "UV4-A" "UV5-A" "UV6-A" "UV7-A" "UV8-A"
[11] "UV9-A" "UV10-A" "UV11-A" "UV12-A" "UV13-A"
[16] "UV14-A" "UV15-A" "UV16-A" "V1-A" "V2-A"
[21] "V3-A" "V4-A" "V5-A" "V6-A" "V7-A"
[26] "V8-A" "V9-A" "V10-A" "V11-A" "V12-A"
[31] "V13-A" "V14-A" "V15-A" "V16-A" "B1-A"
[36] "B2-A" "B3-A" "B4-A" "B5-A" "B6-A"
[41] "B7-A" "B8-A" "B9-A" "B10-A" "B11-A"
[46] "B12-A" "B13-A" "B14-A" "YG1-A" "YG2-A"
[51] "YG3-A" "YG4-A" "YG5-A" "YG6-A" "YG7-A"
[56] "YG8-A" "YG9-A" "YG10-A" "R1-A" "R2-A"
[61] "R3-A" "R4-A" "R5-A" "R6-A" "R7-A"
[66] "R8-A" Fluorophore Antigen UV1-A UV2-A UV3-A UV4-A UV5-A UV6-A UV7-A UV8-A UV9-A
1 BUV395 CD62L 0.258 1.000 0.494 0.360 0.325 0.315 0.177 0.056 0.022
2 BUV496 CD8 0.009 0.039 0.032 0.027 0.036 0.187 1.000 0.534 0.256
3 BUV563 CD69 0.007 0.029 0.023 0.020 0.019 0.020 0.017 0.327 1.000
UV10-A UV11-A UV12-A UV13-A UV14-A UV15-A UV16-A V1-A V2-A V3-A V4-A
1 0.004 0.001 0.001 0.000 0.001 0.000 0.000 0.002 0.003 0.003 0.002
2 0.065 0.017 0.007 0.004 0.004 0.002 0.001 0.001 0.002 0.012 0.066
3 0.264 0.062 0.025 0.015 0.011 0.007 0.004 0.000 0.001 0.001 0.001
V5-A V6-A V7-A V8-A V9-A V10-A V11-A V12-A V13-A V14-A V15-A V16-A B1-A
1 0.002 0.002 0.002 0.002 0.001 0.001 0.001 0.000 0.000 0.000 0.000 0 0.001
2 0.216 0.164 0.170 0.069 0.036 0.029 0.008 0.003 0.002 0.001 0.001 0 0.068
3 0.001 0.002 0.027 0.073 0.035 0.033 0.008 0.003 0.002 0.001 0.001 0 0.000
B2-A B3-A B4-A B5-A B6-A B7-A B8-A B9-A B10-A B11-A B12-A B13-A B14-A
1 0.001 0.001 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
2 0.043 0.034 0.010 0.006 0.003 0.001 0.001 0.001 0.000 0.000 0.000 0.000 0.000
3 0.006 0.109 0.218 0.118 0.076 0.025 0.015 0.012 0.006 0.003 0.002 0.001 0.001
YG1-A YG2-A YG3-A YG4-A YG5-A YG6-A YG7-A YG8-A YG9-A YG10-A R1-A R2-A R3-A
1 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0 0 0
2 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0 0 0
3 0.488 0.246 0.146 0.044 0.026 0.019 0.017 0.006 0.004 0.002 0 0 0
R4-A R5-A R6-A R7-A R8-A
1 0 0 0 0 0
2 0 0 0 0 0
3 0 0 0 0 0 [1] "Fluorophore" "Antigen" "UV1-A" "UV2-A" "UV3-A"
[6] "UV4-A" "UV5-A" "UV6-A" "UV7-A" "UV8-A"
[11] "UV9-A" "UV10-A" "UV11-A" "UV12-A" "UV13-A"
[16] "UV14-A" "UV15-A" "UV16-A" "V1-A" "V2-A"
[21] "V3-A" "V4-A" "V5-A" "V6-A" "V7-A"
[26] "V8-A" "V9-A" "V10-A" "V11-A" "V12-A"
[31] "V13-A" "V14-A" "V15-A" "V16-A" "B1-A"
[36] "B2-A" "B3-A" "B4-A" "B5-A" "B6-A"
[41] "B7-A" "B8-A" "B9-A" "B10-A" "B11-A"
[46] "B12-A" "B13-A" "B14-A" "YG1-A" "YG2-A"
[51] "YG3-A" "YG4-A" "YG5-A" "YG6-A" "YG7-A"
[56] "YG8-A" "YG9-A" "YG10-A" "R1-A" "R2-A"
[61] "R3-A" "R4-A" "R5-A" "R6-A" "R7-A"
[66] "R8-A" Fluorophore Antigen UV1-A UV2-A UV3-A UV4-A UV5-A UV6-A UV7-A UV8-A UV9-A
1 BUV395 CD62L 0.259 1.000 0.487 0.357 0.321 0.316 0.198 0.073 0.041
2 BUV496 CD8 0.009 0.039 0.031 0.027 0.036 0.186 1.000 0.531 0.253
3 BUV563 CD69 0.007 0.030 0.024 0.021 0.021 0.023 0.020 0.342 1.000
UV10-A UV11-A UV12-A UV13-A UV14-A UV15-A UV16-A V1-A V2-A V3-A V4-A
1 0.015 0.008 0.006 0.005 0.006 0.007 0.005 0.010 0.028 0.048 0.049
2 0.064 0.018 0.008 0.005 0.005 0.003 0.002 0.001 0.004 0.017 0.072
3 0.260 0.061 0.026 0.015 0.012 0.007 0.004 0.000 -0.002 -0.003 -0.004
V5-A V6-A V7-A V8-A V9-A V10-A V11-A V12-A V13-A V14-A V15-A V16-A
1 0.068 0.066 0.093 0.078 0.060 0.071 0.041 0.025 0.023 0.021 0.019 0.010
2 0.223 0.170 0.180 0.077 0.043 0.038 0.014 0.007 0.005 0.004 0.003 0.001
3 -0.004 -0.001 0.025 0.069 0.033 0.031 0.008 0.003 0.002 0.001 0.001 0.000
B1-A B2-A B3-A B4-A B5-A B6-A B7-A B8-A B9-A B10-A B11-A B12-A B13-A
1 0.018 0.029 0.042 0.036 0.032 0.031 0.023 0.018 0.016 0.012 0.009 0.009 0.007
2 0.070 0.047 0.039 0.014 0.010 0.007 0.004 0.003 0.003 0.002 0.001 0.001 0.001
3 0.000 0.008 0.107 0.201 0.108 0.070 0.024 0.015 0.012 0.006 0.003 0.003 0.002
B14-A YG1-A YG2-A YG3-A YG4-A YG5-A YG6-A YG7-A YG8-A YG9-A YG10-A R1-A
1 0.008 0.015 0.023 0.027 0.024 0.015 0.016 0.016 0.008 0.004 0.002 0.003
2 0.001 0.002 0.004 0.006 0.006 0.004 0.004 0.004 0.002 0.001 0.001 0.002
3 0.002 0.436 0.212 0.126 0.039 0.022 0.017 0.014 0.004 0.003 0.001 -0.001
R2-A R3-A R4-A R5-A R6-A R7-A R8-A
1 0.004 0.002 0.000 0.005 0.003 0.004 0.002
2 0.003 0.002 0.002 0.002 0.001 0.001 0.001
3 -0.001 -0.001 -0.001 0.000 0.000 0.000 0.000.
And with that, we have the main elements we will need for our unmixing.
Planning
.
Before diving straight in, it is often worthwhile to strategize what we are trying to accomplish first.
.
Over the last couple weeks, we have derived out signature matrices for our single-color controls (from both beads and cells). Similarly, we have an additional matrix that contains unstained signatures. These matrices in addition to the detector columns also contain metadata, in the form of the Fluorophore and Antigen columns.
.
Separately, we have file.paths to the location of our individual raw .fcs files. Once these are loaded into R, we would have access to the exprs() matrices where the raw measurement values are stored. After these are extracted, we need to provide these and the equivalent signature matrix columns to the function that carries out unmixing calculation and returns the respective fluorophore contributions.
.
We consequently end up with a new matrix containing the calculated abundances of each fluorophore for a given cell. These in turn need to be swapped in for the exprs() slot within the “flowFrame” object. At this point, any corresponding changes to the parameters() and keyword() slots would need to be carried out, converting over the formatting from one typical of a raw .fcs file to that of an unmixed .fcs file. Some of these steps may overlap with the Week 10 bonus material, so we may be able to reuse some of that code instead of writing from scratch.
.
Lastly, we would then need to save our new .fcs file to a designated location. At which point, we can visualize it (to make sure no formatting mistakes were made in the coding process). If none are encountered, we would then be set so that next time we can start off by evaluating the unmixed .fcs file for unmixing errors that might affect the downstream analysis process.
Loading FCS files
.
To get started, lets create a basic function skeleton. We can designate the function name as OldFashionedUnmix(). As always when building functions, individual arguments go inside the “()”, and code being run goes within the “{}”.
.
In terms of what should be our OldFashionedUnmix() function’s first argument, at the end of the day, we would want to simply provide to our function a vector containing the file.path()’s to where our .fcs files, and then stand back and let the function handle the rest of the process.
.
So we will need to set it up to iterate (using either purrr’s map() or walk() functions, or base R equivalents). Consequently, lets simply use “x” as our first argument. We can add it to our function by placing the new argument name between the “()”.
.
Additionally, to help our future-selves out when we forget what each argument does next month (or who are we kidding, maybe tomorrow), lets also write out the arguments documentation in the roxygen2 skeleton, by providing the name after “@param”, followed by a short description of purpose.
.
We can select “Run Cell” on our code-block to run the function, which results in it being created in our local environment (appearing in Positron listed under Session tab in the right secondary sidebar).
.
With our OldFashionedUnmix() function now created, we can evaluate the output of the code inside by running the code block below using the purrr packages walk() function to iterate through our vector of .fcs file paths.
[1] "data/DTR_2023_ILT_01-INF052-Ctrl_Antibody.1235515.fcs"
[1] "data/DTR_2023_ILT_01-INF052-PMA_Antibody.1235651.fcs"
[1] "data/DTR_2023_ILT_01-ND050_v1-Ctrl_Antibody.1235673.fcs"
[1] "data/DTR_2023_ILT_01-ND050_v1-PMA_Antibody.1235676.fcs".
As we proceed and build out the rest of the function, remember to re-run the code block to refresh OldFashionedUnmix() in the local environment, so that you can see the effect of the changes that are made.
.
Next up, we will need to load the designated .fcs file into our local environment from where it is being stored (as being designated by the iterated in file path stored in “x”).
.
As we encountered in Week 05, there are two ways we can approach this, loading the file in as either a ‘flowFrame’ or a ‘cytoframe’ object. The main difference between these two versions are whether the object is loaded immediately into rapid access memory (RAM) or not (via a pointer to the location).
.
In the case of unmixing, we need access to the internals, are not planning to pre-gate, and are unmixing individual objects at a time. Because of this, lets use the flowCore packages read.FCS() function to load our raw .fcs file in directly as a flowFrame().
.
Since our function now depends on an external package, lets make a note of this by adding the “@importFrom” line to our roxygen2 skeleton, listing the parent package and function. Since we are not connecting to a Description/Namespace file at this point of the course, lets also add “flowCore::” before the functions name inside the code, so that this dependency is also noted in different context. From here, refresh the function and test the output.
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#'
#' @importFrom flowCore read.FCS
#'
OldFashionedUnmix <- function(x){
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
return(TheRawFCS)
}[[1]]
flowFrame object 'DTR_2023_ILT_01-INF052-Ctrl_Antibody.1235515.fcs'
with 10000 cells and 74 observables:
name desc range minRange maxRange
$P1 Time NA 1428432 0 1428431
$P2 UV1-A NA 4194304 -111 4194304
$P3 UV2-A NA 4194304 -111 4194304
$P4 UV3-A NA 4194304 -111 4194304
$P5 UV4-A NA 4194304 -111 4194304
... ... ... ... ... ...
$P70 R4-A NA 4194304 -111 4194304
$P71 R5-A NA 4194304 -111 4194304
$P72 R6-A NA 4194304 -111 4194304
$P73 R7-A NA 4194304 -111 4194304
$P74 R8-A NA 4194304 -111 4194304
719 keywords are stored in the 'description' slot.
Alright, we are getting back our flowFrame object as expected. We can now start gathering the actual components for the unmixing.
Retrieving exprs values
.
Next order of business is retrieving the exprs() slot content (containing the raw MFI values). This will contain our detector columns, as well as additional columns for “Time”, “SSC” (off the violet laser), “FSC”, and “SSC-B” (off the blue laser).
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#'
#' @importFrom flowCore read.FCS exprs
#'
OldFashionedUnmix <- function(x){
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
return(TheRawMatrix)
} Time UV1-A UV2-A UV3-A UV4-A UV5-A UV6-A UV7-A
[1,] 887309 964.4953 7116.870 6668.760 5442.837 6575.271 29927.244 148718.08
[2,] 1393476 573.2829 2249.769 1380.251 1108.857 1924.750 4140.159 14098.38
[3,] 1258794 3870.6978 13829.658 9413.271 6581.815 6863.400 8449.966 10353.08
UV8-A UV9-A UV10-A UV11-A UV12-A UV13-A UV14-A
[1,] 83858.398 49607.305 13840.966 7501.835 8714.148 18463.369 16111.409
[2,] 22587.686 42985.773 11908.258 5140.280 2746.446 1408.366 2365.944
[3,] 7921.755 7253.571 3959.576 5106.959 10668.052 27839.826 24112.010
UV15-A UV16-A SSC-W SSC-H SSC-A V1-A V2-A
[1,] 15055.137 11021.854 683436.8 843643 960961.0 12381.396 14378.654
[2,] 3713.767 7759.097 658488.6 1080396 1185714.0 1923.075 5794.869
[3,] 31397.920 45514.152 756646.2 714195 900654.8 20689.764 22286.967
V3-A V4-A V5-A V6-A V7-A V8-A V9-A V10-A
[1,] 16559.46 19534.35 43587.71 38181.98 51864.67 35652.65 21050.21 22895.95
[2,] 10455.71 12758.49 24218.83 23529.07 33043.38 24329.87 15209.02 16929.69
[3,] 22783.75 13194.50 14704.46 16057.25 29045.98 18846.78 11669.49 16529.01
V11-A V12-A V13-A V14-A V15-A V16-A FSC-W
[1,] 25843.68 42143.965 107268.664 70643.219 62196.21 32139.05 696000.9
[2,] 12336.57 5332.318 5721.375 6142.263 10714.00 14052.50 665755.6
[3,] 25978.15 64206.992 170656.406 109576.125 93684.66 49434.56 723983.9
FSC-H FSC-A SSC-B-W SSC-B-H SSC-B-A B1-A B2-A B3-A
[1,] 1473050 1708740 676320.5 608196 685559.1 13091.453 25392.572 54251.03
[2,] 1411870 1566601 660120.8 782257 860640.1 1477.710 2176.785 6405.36
[3,] 1276900 1540758 715180.6 554999 661540.9 3955.835 19526.842 48903.08
B4-A B5-A B6-A B7-A B8-A B9-A B10-A
[1,] 54079.10 30304.557 17815.393 8657.025 6662.045 10566.27 9565.594
[2,] 9868.17 6189.625 4117.425 2333.955 1899.690 1821.17 1807.325
[3,] 20192.06 14296.358 8729.176 6668.154 6366.620 11939.72 12764.375
B11-A B12-A B13-A B14-A YG1-A YG2-A YG3-A YG4-A
[1,] 7321.2754 6748.429 7185.166 7885.149 58635.73 26546.73 15374.80 11338.110
[2,] 977.0795 1137.500 1444.950 3106.935 18065.18 8613.78 6474.16 6698.371
[3,] 9258.7959 7538.699 7746.180 9104.485 2733.22 1766.38 3066.42 8749.650
YG5-A YG6-A YG7-A YG8-A YG9-A YG10-A R1-A R2-A
[1,] 8914.430 10226.93 22174.04 10758.09 19897.71 21207.76 6629.140 8471.260
[2,] 5722.640 4309.27 7127.54 3617.74 6346.34 16151.24 8631.559 9068.079
[3,] 8039.431 10995.11 28742.22 12658.80 18004.07 19135.83 9863.699 14070.698
R3-A R4-A R5-A R6-A R7-A R8-A
[1,] 23756.39 35678.38 26992.84 27509.79 62521.95 56157.36
[2,] 10793.09 11519.06 9824.64 7929.60 22465.80 53715.19
[3,] 36907.98 55936.11 42062.24 31035.13 60446.53 58800.29.
We get back from exprs() a “matrix” style object, we can convert this to a “data.frame” object to allow for easier selection of columns using the dplyr package.
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#'
#' @importFrom flowCore read.FCS exprs
#'
OldFashionedUnmix <- function(x){
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
return(TheRawDataFrame)
} Time UV1-A UV2-A UV3-A UV4-A UV5-A UV6-A UV7-A
1 887309 964.4953 7116.870 6668.760 5442.837 6575.271 29927.244 148718.08
2 1393476 573.2829 2249.769 1380.251 1108.857 1924.750 4140.159 14098.38
3 1258794 3870.6978 13829.658 9413.271 6581.815 6863.400 8449.966 10353.08
UV8-A UV9-A UV10-A UV11-A UV12-A UV13-A UV14-A
1 83858.398 49607.305 13840.966 7501.835 8714.148 18463.369 16111.409
2 22587.686 42985.773 11908.258 5140.280 2746.446 1408.366 2365.944
3 7921.755 7253.571 3959.576 5106.959 10668.052 27839.826 24112.010
UV15-A UV16-A SSC-W SSC-H SSC-A V1-A V2-A V3-A
1 15055.137 11021.854 683436.8 843643 960961.0 12381.396 14378.654 16559.46
2 3713.767 7759.097 658488.6 1080396 1185714.0 1923.075 5794.869 10455.71
3 31397.920 45514.152 756646.2 714195 900654.8 20689.764 22286.967 22783.75
V4-A V5-A V6-A V7-A V8-A V9-A V10-A V11-A
1 19534.35 43587.71 38181.98 51864.67 35652.65 21050.21 22895.95 25843.68
2 12758.49 24218.83 23529.07 33043.38 24329.87 15209.02 16929.69 12336.57
3 13194.50 14704.46 16057.25 29045.98 18846.78 11669.49 16529.01 25978.15
V12-A V13-A V14-A V15-A V16-A FSC-W FSC-H FSC-A
1 42143.965 107268.664 70643.219 62196.21 32139.05 696000.9 1473050 1708740
2 5332.318 5721.375 6142.263 10714.00 14052.50 665755.6 1411870 1566601
3 64206.992 170656.406 109576.125 93684.66 49434.56 723983.9 1276900 1540758
SSC-B-W SSC-B-H SSC-B-A B1-A B2-A B3-A B4-A B5-A
1 676320.5 608196 685559.1 13091.453 25392.572 54251.03 54079.10 30304.557
2 660120.8 782257 860640.1 1477.710 2176.785 6405.36 9868.17 6189.625
3 715180.6 554999 661540.9 3955.835 19526.842 48903.08 20192.06 14296.358
B6-A B7-A B8-A B9-A B10-A B11-A B12-A B13-A
1 17815.393 8657.025 6662.045 10566.27 9565.594 7321.2754 6748.429 7185.166
2 4117.425 2333.955 1899.690 1821.17 1807.325 977.0795 1137.500 1444.950
3 8729.176 6668.154 6366.620 11939.72 12764.375 9258.7959 7538.699 7746.180
B14-A YG1-A YG2-A YG3-A YG4-A YG5-A YG6-A YG7-A
1 7885.149 58635.73 26546.73 15374.80 11338.110 8914.430 10226.93 22174.04
2 3106.935 18065.18 8613.78 6474.16 6698.371 5722.640 4309.27 7127.54
3 9104.485 2733.22 1766.38 3066.42 8749.650 8039.431 10995.11 28742.22
YG8-A YG9-A YG10-A R1-A R2-A R3-A R4-A R5-A
1 10758.09 19897.71 21207.76 6629.140 8471.260 23756.39 35678.38 26992.84
2 3617.74 6346.34 16151.24 8631.559 9068.079 10793.09 11519.06 9824.64
3 12658.80 18004.07 19135.83 9863.699 14070.698 36907.98 55936.11 42062.24
R6-A R7-A R8-A
1 27509.79 62521.95 56157.36
2 7929.60 22465.80 53715.19
3 31035.13 60446.53 58800.29.
With this done, we will need to isolate the detector columns (since it will be their values we will need for unmixing) from the other columns (“Time”, “SSC”, “FSC” and “SSC-B”). We will still want to retain these columns, as following unmixing we will need to add them back to our unmixed .fcs file so we can gate based off scatter.
.
Rather than directly specifying which columns, we can use the dplyr packages select() and matches() functions in combination to grab any column whose column name matches a particular character string (or strings when separated by “|”). So for the case above, we could within matches() place “FSC|SSC|Time” to enable these columns to be grabbed.
.
However, for certain instruments or cases, there could be other columns we might want to retain. Rather than needing to modify the function directly, we can set up a second argument (“retainThese”), and provide a default value (“FSC|SSC|Time”). That way we retain the ability to easily modify the pattern to handle the edge cases, without needing to specify the values each time when the default suffices.
.
Lets go ahead and add the argument and it’s default value between the “()”, and update the roxygen2 documentation for this new argument (and what it does). We also need to list the two dplyr functions in the roxygen, and add the “dplyr::” before the respective functions in the code to avoid dependency issues.
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#'
#' @importFrom flowCore read.FCS exprs
#' @importFrom dplyr select matches
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time"){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
TheColNames <- colnames(StashedColumns)
return(TheColNames)
}.
Alright, we have a generalized way to designate which columns need to kept outside of the unmixing (but later placed back). By adding an “!” to that line of codes select(), we can also gather the non-stashed columns (which will mainly be our detector columns of interest).
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#'
#' @importFrom flowCore read.FCS exprs
#' @importFrom dplyr select matches
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time"){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
WorkingColumns <- TheRawDataFrame |>
dplyr::select(!dplyr::matches(retainThese))
TheColNames <- colnames(WorkingColumns)
return(TheColNames)
}[[1]]
[1] "UV1-A" "UV2-A" "UV3-A" "UV4-A" "UV5-A" "UV6-A" "UV7-A" "UV8-A"
[9] "UV9-A" "UV10-A" "UV11-A" "UV12-A" "UV13-A" "UV14-A" "UV15-A" "UV16-A"
[17] "V1-A" "V2-A" "V3-A" "V4-A" "V5-A" "V6-A" "V7-A" "V8-A"
[25] "V9-A" "V10-A" "V11-A" "V12-A" "V13-A" "V14-A" "V15-A" "V16-A"
[33] "B1-A" "B2-A" "B3-A" "B4-A" "B5-A" "B6-A" "B7-A" "B8-A"
[41] "B9-A" "B10-A" "B11-A" "B12-A" "B13-A" "B14-A" "YG1-A" "YG2-A"
[49] "YG3-A" "YG4-A" "YG5-A" "YG6-A" "YG7-A" "YG8-A" "YG9-A" "YG10-A"
[57] "R1-A" "R2-A" "R3-A" "R4-A" "R5-A" "R6-A" "R7-A" "R8-A" .
For the Cytek Aurora .fcs files in this example, this last line of code was sufficient to isolate the detector columns. However, this was in part to the instrument settings not having also acquired the “-H” or “-W” values for the various detectors. Since for some manufacturers instruments, these are also acquired by default, we maay need an additional step to designate whether the detector columns being selected for end in “-A”, “-H”, or “-W”.
.
To do this, we could use another argument where we either select the variant we want, or designate the variants we want to exclude. For now, lets go with the exclude approach, externalizing this in the form of a third argument (“detectorExclude”) with a default value of “-H|-W”.
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#' @param detectorExclude Default is "-H|-W", intended to remove
#' additional detector columns other than -A, adjust as needed for
#' your own instruments configuration
#'
#' @importFrom flowCore read.FCS exprs
#' @importFrom dplyr select matches
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time",
detectorExclude="-H|-W"){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
WorkingColumns <- TheRawDataFrame |>
dplyr::select(!dplyr::matches(retainThese))
WorkingColumns <- WorkingColumns |>
dplyr::select(!dplyr::matches(detectorExclude))
TheColNames <- colnames(WorkingColumns)
return(TheColNames)
}[[1]]
[1] "UV1-A" "UV2-A" "UV3-A" "UV4-A" "UV5-A" "UV6-A" "UV7-A" "UV8-A"
[9] "UV9-A" "UV10-A" "UV11-A" "UV12-A" "UV13-A" "UV14-A" "UV15-A" "UV16-A"
[17] "V1-A" "V2-A" "V3-A" "V4-A" "V5-A" "V6-A" "V7-A" "V8-A"
[25] "V9-A" "V10-A" "V11-A" "V12-A" "V13-A" "V14-A" "V15-A" "V16-A"
[33] "B1-A" "B2-A" "B3-A" "B4-A" "B5-A" "B6-A" "B7-A" "B8-A"
[41] "B9-A" "B10-A" "B11-A" "B12-A" "B13-A" "B14-A" "YG1-A" "YG2-A"
[49] "YG3-A" "YG4-A" "YG5-A" "YG6-A" "YG7-A" "YG8-A" "YG9-A" "YG10-A"
[57] "R1-A" "R2-A" "R3-A" "R4-A" "R5-A" "R6-A" "R7-A" "R8-A" .
Having pre-empted a potential issue, and still retaining the detector columns of interest, lets swap the return() content to get back the “data.frame” rather than just the colnames() output, so that we can see the underlying exprs() values.
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#' @param detectorExclude Default is "-H|-W", intended to remove
#' additional detector columns other than -A, adjust as needed for
#' your own instruments configuration
#'
#' @importFrom flowCore read.FCS exprs
#' @importFrom dplyr select matches
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time",
detectorExclude="-H|-W"){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
WorkingColumns <- TheRawDataFrame |>
dplyr::select(!dplyr::matches(retainThese))
WorkingColumns <- WorkingColumns |>
dplyr::select(!dplyr::matches(detectorExclude))
# TheColNames <- colnames(WorkingColumns)
return(WorkingColumns)
} UV1-A UV2-A UV3-A UV4-A UV5-A UV6-A UV7-A UV8-A
1 964.4953 7116.870 6668.760 5442.837 6575.271 29927.244 148718.08 83858.398
2 573.2829 2249.769 1380.251 1108.857 1924.750 4140.159 14098.38 22587.686
3 3870.6978 13829.658 9413.271 6581.815 6863.400 8449.966 10353.08 7921.755
UV9-A UV10-A UV11-A UV12-A UV13-A UV14-A UV15-A
1 49607.305 13840.966 7501.835 8714.148 18463.369 16111.409 15055.137
2 42985.773 11908.258 5140.280 2746.446 1408.366 2365.944 3713.767
3 7253.571 3959.576 5106.959 10668.052 27839.826 24112.010 31397.920
UV16-A V1-A V2-A V3-A V4-A V5-A V6-A V7-A
1 11021.854 12381.396 14378.654 16559.46 19534.35 43587.71 38181.98 51864.67
2 7759.097 1923.075 5794.869 10455.71 12758.49 24218.83 23529.07 33043.38
3 45514.152 20689.764 22286.967 22783.75 13194.50 14704.46 16057.25 29045.98
V8-A V9-A V10-A V11-A V12-A V13-A V14-A V15-A
1 35652.65 21050.21 22895.95 25843.68 42143.965 107268.664 70643.219 62196.21
2 24329.87 15209.02 16929.69 12336.57 5332.318 5721.375 6142.263 10714.00
3 18846.78 11669.49 16529.01 25978.15 64206.992 170656.406 109576.125 93684.66
V16-A B1-A B2-A B3-A B4-A B5-A B6-A B7-A
1 32139.05 13091.453 25392.572 54251.03 54079.10 30304.557 17815.393 8657.025
2 14052.50 1477.710 2176.785 6405.36 9868.17 6189.625 4117.425 2333.955
3 49434.56 3955.835 19526.842 48903.08 20192.06 14296.358 8729.176 6668.154
B8-A B9-A B10-A B11-A B12-A B13-A B14-A YG1-A
1 6662.045 10566.27 9565.594 7321.2754 6748.429 7185.166 7885.149 58635.73
2 1899.690 1821.17 1807.325 977.0795 1137.500 1444.950 3106.935 18065.18
3 6366.620 11939.72 12764.375 9258.7959 7538.699 7746.180 9104.485 2733.22
YG2-A YG3-A YG4-A YG5-A YG6-A YG7-A YG8-A YG9-A
1 26546.73 15374.80 11338.110 8914.430 10226.93 22174.04 10758.09 19897.71
2 8613.78 6474.16 6698.371 5722.640 4309.27 7127.54 3617.74 6346.34
3 1766.38 3066.42 8749.650 8039.431 10995.11 28742.22 12658.80 18004.07
YG10-A R1-A R2-A R3-A R4-A R5-A R6-A R7-A
1 21207.76 6629.140 8471.260 23756.39 35678.38 26992.84 27509.79 62521.95
2 16151.24 8631.559 9068.079 10793.09 11519.06 9824.64 7929.60 22465.80
3 19135.83 9863.699 14070.698 36907.98 55936.11 42062.24 31035.13 60446.53
R8-A
1 56157.36
2 53715.19
3 58800.29Loading signature matrix
.
From our raw .fcs file, we have now retrieved exprs() values we will use for unmixing. Within this “data.frame”, each row contains the detector measurements recorded for an individual cell. During unmixing, each row is unmixed vs. the reference signature matrix, which allows us to estimate the presence/abscence of each fluorophore, as well as its relative abundance. When unmixing is complete, we have a “matrix” object with the same number of cells, but instead of the detector columns, new columns for each fluorophore in the provided signature matrix. This “matrix” object is ultimately the one added back to the exprs() slot.
.
Beyond providing the signature values that distinguish each fluorophore from each other, the signature matrix can also provides the source for the fluorophore and antigen names. Lets go ahead and load our signature matrix into our function via new argument (“SignatureData”).
.
We can also retain the flexibility of either providing a “data.frame” or a file.path to a “.csv” file through the use of a conditional statement (“if”, “else”, “else if”) depending on the type of object provided to “SignatureData”
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#' @param detectorExclude Default is "-H|-W", intended to remove
#' additional detector columns other than -A, adjust as needed for
#' your own instruments configuration
#' @param SignatureData A data.frame containing a Fluorophore,
#' Antigen and Detector columns.
#'
#' @importFrom flowCore read.FCS exprs
#' @importFrom dplyr select matches
#' @importFrom utils read.csv
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time",
detectorExclude="-H|-W", SignatureData){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
WorkingColumns <- TheRawDataFrame |>
dplyr::select(!dplyr::matches(retainThese))
WorkingColumns <- WorkingColumns |>
dplyr::select(!dplyr::matches(detectorExclude))
# TheColNames <- colnames(WorkingColumns)
# Load the signature matrix
if(is.data.frame(SignatureData)){
Signatures <- SignatureData
} else {
Signatures <-read.csv(SignatureData, check.names=FALSE)
}
return(Signatures)
} Fluorophore Antigen UV1-A UV2-A UV3-A UV4-A UV5-A UV6-A UV7-A UV8-A UV9-A
1 BUV395 CD62L 0.258 1.000 0.494 0.360 0.325 0.315 0.177 0.056 0.022
2 BUV496 CD8 0.009 0.039 0.032 0.027 0.036 0.187 1.000 0.534 0.256
3 BUV563 CD69 0.007 0.029 0.023 0.020 0.019 0.020 0.017 0.327 1.000
UV10-A UV11-A UV12-A UV13-A UV14-A UV15-A UV16-A V1-A V2-A V3-A V4-A
1 0.004 0.001 0.001 0.000 0.001 0.000 0.000 0.002 0.003 0.003 0.002
2 0.065 0.017 0.007 0.004 0.004 0.002 0.001 0.001 0.002 0.012 0.066
3 0.264 0.062 0.025 0.015 0.011 0.007 0.004 0.000 0.001 0.001 0.001
V5-A V6-A V7-A V8-A V9-A V10-A V11-A V12-A V13-A V14-A V15-A V16-A B1-A
1 0.002 0.002 0.002 0.002 0.001 0.001 0.001 0.000 0.000 0.000 0.000 0 0.001
2 0.216 0.164 0.170 0.069 0.036 0.029 0.008 0.003 0.002 0.001 0.001 0 0.068
3 0.001 0.002 0.027 0.073 0.035 0.033 0.008 0.003 0.002 0.001 0.001 0 0.000
B2-A B3-A B4-A B5-A B6-A B7-A B8-A B9-A B10-A B11-A B12-A B13-A B14-A
1 0.001 0.001 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
2 0.043 0.034 0.010 0.006 0.003 0.001 0.001 0.001 0.000 0.000 0.000 0.000 0.000
3 0.006 0.109 0.218 0.118 0.076 0.025 0.015 0.012 0.006 0.003 0.002 0.001 0.001
YG1-A YG2-A YG3-A YG4-A YG5-A YG6-A YG7-A YG8-A YG9-A YG10-A R1-A R2-A R3-A
1 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0 0 0
2 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0 0 0
3 0.488 0.246 0.146 0.044 0.026 0.019 0.017 0.006 0.004 0.002 0 0 0
R4-A R5-A R6-A R7-A R8-A
1 0 0 0 0 0
2 0 0 0 0 0
3 0 0 0 0 0Separating Metadata
.
With this accomplished, we can separate our character columns (containing our metadata for Fluorophore and Antigen) from the numeric columns (i.e. the detector columns used during unmixing).
.
We can then similarly separate out the character columns (being used for labeling) from the numeric columns (being used in the unmixing process). This can be accomplished through use of select(), where() and is.numeric() as we have previously encountered.
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#' @param detectorExclude Default is "-H|-W", intended to remove
#' additional detector columns other than -A, adjust as needed for
#' your own instruments configuration
#' @param SignatureData A data.frame containing a Fluorophore,
#' Antigen and Detector columns.
#'
#' @importFrom flowCore read.FCS exprs
#' @importFrom dplyr select matches where
#' @importFrom utils read.csv
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time",
detectorExclude="-H|-W", SignatureData){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
WorkingColumns <- TheRawDataFrame |>
dplyr::select(!dplyr::matches(retainThese))
WorkingColumns <- WorkingColumns |>
dplyr::select(!dplyr::matches(detectorExclude))
# TheColNames <- colnames(WorkingColumns)
# Load the signature matrix
if(is.data.frame(SignatureData)){
Signatures <- SignatureData
} else {
Signatures <-read.csv(SignatureData, check.names=FALSE)
}
Metadata <- Signatures |>
dplyr::select(!dplyr::where(is.numeric))
Numerics <- Signatures |>
dplyr::select(dplyr::where(is.numeric))
return(colnames(Numerics))
} [1] "UV1-A" "UV2-A" "UV3-A" "UV4-A" "UV5-A" "UV6-A" "UV7-A" "UV8-A"
[9] "UV9-A" "UV10-A" "UV11-A" "UV12-A" "UV13-A" "UV14-A" "UV15-A" "UV16-A"
[17] "V1-A" "V2-A" "V3-A" "V4-A" "V5-A" "V6-A" "V7-A" "V8-A"
[25] "V9-A" "V10-A" "V11-A" "V12-A" "V13-A" "V14-A" "V15-A" "V16-A"
[33] "B1-A" "B2-A" "B3-A" "B4-A" "B5-A" "B6-A" "B7-A" "B8-A"
[41] "B9-A" "B10-A" "B11-A" "B12-A" "B13-A" "B14-A" "YG1-A" "YG2-A"
[49] "YG3-A" "YG4-A" "YG5-A" "YG6-A" "YG7-A" "YG8-A" "YG9-A" "YG10-A"
[57] "R1-A" "R2-A" "R3-A" "R4-A" "R5-A" "R6-A" "R7-A" "R8-A" Normalizing Check
.
Additionally, our signature retrieval functions can return either the raw values or normalized/scaled values that range from 0 to 1. For unmixing, we need the later. We can add an internal check through the use of another conditional statement, where if it is trigerred by the presence of a value greater than 1, it will proceed to normalize the values in the matrix before proceeding.
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#' @param detectorExclude Default is "-H|-W", intended to remove
#' additional detector columns other than -A, adjust as needed for
#' your own instruments configuration
#' @param SignatureData A data.frame containing a Fluorophore,
#' Antigen and Detector columns.
#'
#' @importFrom flowCore read.FCS exprs
#' @importFrom dplyr select matches where
#' @importFrom utils read.csv
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time",
detectorExclude="-H|-W", SignatureData){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
WorkingColumns <- TheRawDataFrame |>
dplyr::select(!dplyr::matches(retainThese))
WorkingColumns <- WorkingColumns |>
dplyr::select(!dplyr::matches(detectorExclude))
# TheColNames <- colnames(WorkingColumns)
# Load the signature matrix
if(is.data.frame(SignatureData)){
Signatures <- SignatureData
} else {
Signatures <-read.csv(SignatureData, check.names=FALSE)
}
Metadata <- Signatures |>
dplyr::select(!dplyr::where(is.numeric))
Numerics <- Signatures |>
dplyr::select(dplyr::where(is.numeric))
if (any(Numerics > 1)) {
message("Signature values greater than 1 detected, normalizing")
n <- Numerics
# n[n < 0] <- 0
A <- do.call(pmax, n)
Normalized <- n/A
Numerics <- Normalized
}
return(Numerics)
} UV1-A UV2-A UV3-A UV4-A UV5-A UV6-A UV7-A UV8-A UV9-A UV10-A UV11-A UV12-A
1 0.258 1.000 0.494 0.360 0.325 0.315 0.177 0.056 0.022 0.004 0.001 0.001
2 0.009 0.039 0.032 0.027 0.036 0.187 1.000 0.534 0.256 0.065 0.017 0.007
3 0.007 0.029 0.023 0.020 0.019 0.020 0.017 0.327 1.000 0.264 0.062 0.025
UV13-A UV14-A UV15-A UV16-A V1-A V2-A V3-A V4-A V5-A V6-A V7-A V8-A
1 0.000 0.001 0.000 0.000 0.002 0.003 0.003 0.002 0.002 0.002 0.002 0.002
2 0.004 0.004 0.002 0.001 0.001 0.002 0.012 0.066 0.216 0.164 0.170 0.069
3 0.015 0.011 0.007 0.004 0.000 0.001 0.001 0.001 0.001 0.002 0.027 0.073
V9-A V10-A V11-A V12-A V13-A V14-A V15-A V16-A B1-A B2-A B3-A B4-A B5-A
1 0.001 0.001 0.001 0.000 0.000 0.000 0.000 0 0.001 0.001 0.001 0.001 0.000
2 0.036 0.029 0.008 0.003 0.002 0.001 0.001 0 0.068 0.043 0.034 0.010 0.006
3 0.035 0.033 0.008 0.003 0.002 0.001 0.001 0 0.000 0.006 0.109 0.218 0.118
B6-A B7-A B8-A B9-A B10-A B11-A B12-A B13-A B14-A YG1-A YG2-A YG3-A YG4-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.003 0.001 0.001 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
3 0.076 0.025 0.015 0.012 0.006 0.003 0.002 0.001 0.001 0.488 0.246 0.146 0.044
YG5-A YG6-A YG7-A 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.000 0.000 0.000 0 0 0 0 0 0 0 0
2 0.000 0.000 0.000 0.000 0.000 0.000 0 0 0 0 0 0 0 0
3 0.026 0.019 0.017 0.006 0.004 0.002 0 0 0 0 0 0 0 0Dimensions Check
.
Before diving into unmixing, we should also make sure the number of detector columns present in our retrieved exprs() object matches the number of numeric columns retrieved from our signature matrix, to avoid having the code fail out. The simplest implementation would be with a conditional statement, using all() and colnames() as shown here to set up the check.
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#' @param detectorExclude Default is "-H|-W", intended to remove
#' additional detector columns other than -A, adjust as needed for
#' your own instruments configuration
#' @param SignatureData A data.frame containing a Fluorophore,
#' Antigen and Detector columns.
#'
#' @importFrom flowCore read.FCS exprs
#' @importFrom dplyr select matches where
#' @importFrom utils read.csv
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time",
detectorExclude="-H|-W", SignatureData){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
WorkingColumns <- TheRawDataFrame |>
dplyr::select(!dplyr::matches(retainThese))
WorkingColumns <- WorkingColumns |>
dplyr::select(!dplyr::matches(detectorExclude))
# TheColNames <- colnames(WorkingColumns)
# Load the signature matrix
if(is.data.frame(SignatureData)){
Signatures <- SignatureData
} else {
Signatures <-read.csv(SignatureData, check.names=FALSE)
}
Metadata <- Signatures |>
dplyr::select(!dplyr::where(is.numeric))
Numerics <- Signatures |>
dplyr::select(dplyr::where(is.numeric))
if (any(Numerics > 1)) {
message("Signature values greater than 1 detected, normalizing")
n <- Numerics
# n[n < 0] <- 0
A <- do.call(pmax, n)
Normalized <- n/A
Numerics <- Normalized
}
# Make sure dimensions match each other
if (!all(colnames(Numerics) == colnames(WorkingColumns))){
stop("colnames of SignatureData due not match the internal colnames of exprs")
}
return(Numerics)
}
.
With this check now implemented and our starting components collected, we are ready to proceed with the actual unmixing.
Ordinary Least Squares
.
When unmixing, we need to go from the total signal measurements from the individual cells, and referencing the signature reference matrix, calculate out the abundance for each fluorophore that best fits the signal seen on the cell.
.
This process is often mediated via various least squares methods, the more frequently encountered for spectral flow cytometry data being ordinary least squares. In R, this is implemented in the lsfit() function from base R stats package.
.
To carry out its task, both the signature values and the sample values be in the “long” (more rows than columns) style format. Since both our data.frames are in the “wide” (more columns than rows) style format currently, we will need to swing/pivot them.
.
To avoid issues, we can back up the colnames() corresponding to each detectors as a precaution, and then transpose (t()) each to swing to a long format matrix.
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#' @param detectorExclude Default is "-H|-W", intended to remove
#' additional detector columns other than -A, adjust as needed for
#' your own instruments configuration
#' @param SignatureData A data.frame containing a Fluorophore,
#' Antigen and Detector columns.
#'
#' @importFrom flowCore read.FCS exprs
#' @importFrom dplyr select matches where
#' @importFrom utils read.csv
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time",
detectorExclude="-H|-W", SignatureData){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
WorkingColumns <- TheRawDataFrame |>
dplyr::select(!dplyr::matches(retainThese))
WorkingColumns <- WorkingColumns |>
dplyr::select(!dplyr::matches(detectorExclude))
# TheColNames <- colnames(WorkingColumns)
# Load the signature matrix
if(is.data.frame(SignatureData)){
Signatures <- SignatureData
} else {
Signatures <-read.csv(SignatureData, check.names=FALSE)
}
Metadata <- Signatures |>
dplyr::select(!dplyr::where(is.numeric))
Numerics <- Signatures |>
dplyr::select(dplyr::where(is.numeric))
if (any(Numerics > 1)) {
message("Signature values greater than 1 detected, normalizing")
n <- Numerics
# n[n < 0] <- 0
A <- do.call(pmax, n)
Normalized <- n/A
Numerics <- Normalized
}
# Make sure dimensions match each other
if (!all(colnames(Numerics) == colnames(WorkingColumns))){
stop("colnames of SignatureData due not match the internal colnames of exprs")
}
DetectorNameBackups <- colnames(Numerics)
TransposedSignatureValues <- t(Numerics)
TransposedSampleValues <- t(WorkingColumns)
return(TransposedSampleValues)
} [,1] [,2] [,3] [,4] [,5] [,6] [,7]
UV1-A 964.4953 573.2829 3870.698 554.3912 1329.527 -328.5889 468.6372
UV2-A 7116.8696 2249.7693 13829.658 2200.5334 3142.493 388.6842 2027.7600
UV3-A 6668.7598 1380.2511 9413.271 2342.8127 3085.299 679.0435 2014.5955
[,8] [,9] [,10]
UV1-A -18.14732 1866.813 1273.151
UV2-A 370.98251 5948.810 2584.233
UV3-A 750.88953 5262.106 2588.696.
With both our signature and sample values now properly oriented, lets run lsfit() to carry out the unmixing.
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#' @param detectorExclude Default is "-H|-W", intended to remove
#' additional detector columns other than -A, adjust as needed for
#' your own instruments configuration
#' @param SignatureData A data.frame containing a Fluorophore,
#' Antigen and Detector columns.
#'
#' @importFrom flowCore read.FCS exprs
#' @importFrom dplyr select matches where
#' @importFrom utils read.csv
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time",
detectorExclude="-H|-W", SignatureData){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
WorkingColumns <- TheRawDataFrame |>
dplyr::select(!dplyr::matches(retainThese))
WorkingColumns <- WorkingColumns |>
dplyr::select(!dplyr::matches(detectorExclude))
# TheColNames <- colnames(WorkingColumns)
# Load the signature matrix
if(is.data.frame(SignatureData)){
Signatures <- SignatureData
} else {
Signatures <-read.csv(SignatureData, check.names=FALSE)
}
Metadata <- Signatures |>
dplyr::select(!dplyr::where(is.numeric))
Numerics <- Signatures |>
dplyr::select(dplyr::where(is.numeric))
if (any(Numerics > 1)) {
message("Signature values greater than 1 detected, normalizing")
n <- Numerics
# n[n < 0] <- 0
A <- do.call(pmax, n)
Normalized <- n/A
Numerics <- Normalized
}
# Make sure dimensions match each other
if (!all(colnames(Numerics) == colnames(WorkingColumns))){
stop("colnames of SignatureData due not match the internal colnames of exprs")
}
DetectorNameBackups <- colnames(Numerics)
TransposedSignatureValues <- t(Numerics)
TransposedSampleValues <- t(WorkingColumns)
# Unmixing
LeastSquares <- lsfit(x = TransposedSignatureValues,
y = TransposedSampleValues, intercept = FALSE)
return(LeastSquares)
}.
The returned object is a named list, so using names() we can see what the returned elements consist of.
[1] "coefficients" "residuals" "intercept" "qr" Coefficients
.
From the unmixing returns we get back from lsfit(), the first entry in the list is the “coefficients”. If we run dim() to check the dimensions, we get back
.
So rows correspond to the calculated fluorophores, while the columns are the individual cells. Not the usual presentation for our data in an .fcs file, but we can transpose it back to regular shape in a bit. Also, this was a 29-fluorophore panel, so it looks like we missed a fluorophore (likely Zombie NIR as it wasn’t present in the bead reference controls).
.
We can also quickly see a few rows via head() to just confirm everything looks normal, and that we didn’t just get back a matrix with entirely NULL values.
Y1 Y2 Y3 Y4 Y5 Y6 Y7
X1 1651.841 1104.071 12205.677 825.1315 3597.724 777.57145 2235.356
X2 143294.054 1320.456 1234.847 25280.4914 1575.461 -19.94797 1717.416
X3 5202.604 34150.543 1123.961 15173.0698 4062.194 2308.03111 4552.457
Y8 Y9 Y10 Y11 Y12 Y13 Y14
X1 747.78943 4003.4381 3214.621 863.7965 1558.463 1262.8752 417.8141
X2 92.65133 610.1143 2072.842 24839.5463 1492.749 536.1644 1392.5397
X3 1703.87298 3447.0353 5358.980 7929.3683 12000.004 10380.7488 16603.0314
Y15 Y16 Y17 Y18 Y19 Y20
X1 875.4954 1280.8764 672.3464 840.6662 4629.2804 685.4916
X2 3000.2793 320.8317 17470.7487 1173.5742 275.5238 67540.2830
X3 8820.0658 9523.3735 14357.8813 19086.3260 4882.0914 3030.7180Residuals
.
Next up, we have the “residuals”. We can start by checking with dim()
.
In this case, we have 64 rows, and 10000 columns, so these must track with the original detector columns that were provided. We will circle back and visualize these residuals shortly to get a better understanding of how they work.
Intercept
.
The next named entry is “intercept”, which in the actual code line we set to FALSE. This is similarly case here, where it is storing that information for reference.
qr
.
The final list entry is “qr”, which is also a list. When we check the names of this nested list, we get back the following values, which are outside our range to follow up with for the time being.
Residuals
.
So overall, its the “coefficient” values that we will need to grab, since these correspond to our unmixed fluorophore columns that need to go into the .fcs file. However, there is some additional useful information in the residuals that are worth taking a look at, since more complicated unmixing methods take advantage of these same values in their own calculations.
.
First off, we need to a helper function ResidualPlots() to plot the underlying data. This can take the lsfit() returned list and work from there, so we can set the only argument as “LeastSquaresFit”.
#' Takes lsfit output, retrieves and plots the residuals
#'
#' @param LeastSquaresList The returned list from lsfit
#'
ResidualPlots <- function(LeastSquaresList){
library(ggplot2)
residuals <- LeastSquaresList$residual
ResidualsDF <- data.frame(
obs = seq_len(nrow(residuals)), # Vector length residual rows ("Detectors")
rms = sqrt(rowMeans(residuals^2)) # root-mean-square of the fit residuals
)
Plot <- ggplot(ResidualsDF, aes(x = obs, y = rms)) +
geom_segment(aes(xend = obs, yend = 0)) +
labs(title = "RMS residual per observation",
x = "Observation index", y = "RMS residual") +
theme_minimal()
Plot
}.
At which point, we can add ResidualPlots() as the next line in OldFashionedUnmix()
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#' @param detectorExclude Default is "-H|-W", intended to remove
#' additional detector columns other than -A, adjust as needed for
#' your own instruments configuration
#' @param SignatureData A data.frame containing a Fluorophore,
#' Antigen and Detector columns.
#'
#' @importFrom flowCore read.FCS exprs
#' @importFrom dplyr select matches where
#' @importFrom utils read.csv
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time",
detectorExclude="-H|-W", SignatureData){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
WorkingColumns <- TheRawDataFrame |>
dplyr::select(!dplyr::matches(retainThese))
WorkingColumns <- WorkingColumns |>
dplyr::select(!dplyr::matches(detectorExclude))
# TheColNames <- colnames(WorkingColumns)
# Load the signature matrix
if(is.data.frame(SignatureData)){
Signatures <- SignatureData
} else {
Signatures <-read.csv(SignatureData, check.names=FALSE)
}
Metadata <- Signatures |>
dplyr::select(!dplyr::where(is.numeric))
Numerics <- Signatures |>
dplyr::select(dplyr::where(is.numeric))
if (any(Numerics > 1)) {
message("Signature values greater than 1 detected, normalizing")
n <- Numerics
# n[n < 0] <- 0
A <- do.call(pmax, n)
Normalized <- n/A
Numerics <- Normalized
}
# Make sure dimensions match each other
if (!all(colnames(Numerics) == colnames(WorkingColumns))){
stop("colnames of SignatureData due not match the internal colnames of exprs")
}
DetectorNameBackups <- colnames(Numerics)
TransposedSignatureValues <- t(Numerics)
TransposedSampleValues <- t(WorkingColumns)
# Unmixing
LeastSquares <- lsfit(x = TransposedSignatureValues,
y = TransposedSampleValues, intercept = FALSE)
Plot <- ResidualPlots(LeastSquaresList = LeastSquares)
return(Plot)
}.
And lets go ahead and return the plots for all the fcs files for comparison
[[1]]
[[2]]
[[3]]
[[4]]
.
Alright, so in this case (where we forgot to include our Zombie NIR viability dye and unstained signatures in the reference matrix), the place where we see the most leftover residuals (signal that didn’t fit into the calculation) occurs in roughly the same locations where those fluorophores are found.
.
So, lets add them back in and see what happens. We can first retrieve the Unstained and Zombie Signatures to have them on hand.
AlternateCellSignatures <- read.csv(csv_files[3], check.names=FALSE)
AlternateCellSignatures <- AlternateCellSignatures |>
select(-c(name, Type, Detector, Negative))
Unstained <- AlternateCellSignatures |>
dplyr::filter(stringr::str_detect(Fluorophore, "PBMC_Unstained"))
Zombie <- AlternateCellSignatures |>
dplyr::filter(stringr::str_detect(Fluorophore, "Zombie")).
And lets first add back the Zombie signature to the beads signature matrix, and compare vs. the previous residual plot when it wasn’t included
.
So before
.
And after
.
Alright, we no longer have the peak around what would have been the R8 detector. What if we also add back an unstained cell signature?
.
Before with Zombie
.
After with both Unstained and Zombie
.
Some shifts across the board, but the big dramatic shift in residuals was when we accounted for the previously missing Zombie NIR.
.
Just to confirm, what if we removed something important?
.
As you can start to see, even though not the main output we need from lsfit(), the residual output can be quite useful in determining whether all the signal was properly accounted for in the signature matrix or not (which is why other more complicated unmixing methods take advantage of it).
.
Before we turn our focus back to the coefficients, (i.e., our unmixed outputs that were of our main reason for being here), lets add one more argument that would allow us to get back the ResidualPlots() if we wanted to in the future. We can call this argument “returnType”, setting the default to “fcs”, but also specifying a “residual” condition that return these plots if requested.
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#' @param detectorExclude Default is "-H|-W", intended to remove
#' additional detector columns other than -A, adjust as needed for
#' your own instruments configuration
#' @param SignatureData A data.frame containing a Fluorophore,
#' Antigen and Detector columns.
#' @param returnType Default "fcs", for residual plots use "residual"
#'
#' @importFrom flowCore read.FCS exprs
#' @importFrom dplyr select matches where
#' @importFrom utils read.csv
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time",
detectorExclude="-H|-W", SignatureData, returnType="fcs"){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
WorkingColumns <- TheRawDataFrame |>
dplyr::select(!dplyr::matches(retainThese))
WorkingColumns <- WorkingColumns |>
dplyr::select(!dplyr::matches(detectorExclude))
# TheColNames <- colnames(WorkingColumns)
# Load the signature matrix
if(is.data.frame(SignatureData)){
Signatures <- SignatureData
} else {
Signatures <-read.csv(SignatureData, check.names=FALSE)
}
Metadata <- Signatures |>
dplyr::select(!dplyr::where(is.numeric))
Numerics <- Signatures |>
dplyr::select(dplyr::where(is.numeric))
if (any(Numerics > 1)) {
message("Signature values greater than 1 detected, normalizing")
n <- Numerics
# n[n < 0] <- 0
A <- do.call(pmax, n)
Normalized <- n/A
Numerics <- Normalized
}
# Make sure dimensions match each other
if (!all(colnames(Numerics) == colnames(WorkingColumns))){
stop("colnames of SignatureData due not match the internal colnames of exprs")
}
DetectorNameBackups <- colnames(Numerics)
TransposedSignatureValues <- t(Numerics)
TransposedSampleValues <- t(WorkingColumns)
# Unmixing
LeastSquares <- lsfit(x = TransposedSignatureValues,
y = TransposedSampleValues, intercept = FALSE)
# Residual Plot Fork
if (returnType == "residuals"){
Plot <- ResidualPlots(LeastSquaresList = LeastSquares)
return(Plot)
}
return(LeastSquares)
}
Y1 Y2 Y3 Y4 Y5 Y6 Y7
X1 1651.841 1104.071 12205.677 825.1315 3597.724 777.57145 2235.356
X2 143294.054 1320.456 1234.847 25280.4914 1575.461 -19.94797 1717.416
X3 5202.604 34150.543 1123.961 15173.0698 4062.194 2308.03111 4552.457
Y8 Y9 Y10 Y11 Y12 Y13 Y14
X1 747.78943 4003.4381 3214.621 863.7965 1558.463 1262.8752 417.8141
X2 92.65133 610.1143 2072.842 24839.5463 1492.749 536.1644 1392.5397
X3 1703.87298 3447.0353 5358.980 7929.3683 12000.004 10380.7488 16603.0314
Y15 Y16 Y17 Y18 Y19 Y20
X1 875.4954 1280.8764 672.3464 840.6662 4629.2804 685.4916
X2 3000.2793 320.8317 17470.7487 1173.5742 275.5238 67540.2830
X3 8820.0658 9523.3735 14357.8813 19086.3260 4882.0914 3030.7180Reassembly
.
Alright, back to our main reason for being here, the coefficients, i.e. our unmixed data, that now that we have derived we need to change back to a wider format more typical of an .fcs file. We can reverse down this route we earlier traveled by once again transposing (t())
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#' @param detectorExclude Default is "-H|-W", intended to remove
#' additional detector columns other than -A, adjust as needed for
#' your own instruments configuration
#' @param SignatureData A data.frame containing a Fluorophore,
#' Antigen and Detector columns.
#' @param returnType Default "fcs", for residual plots use "residual"
#'
#' @importFrom flowCore read.FCS exprs
#' @importFrom dplyr select matches where
#' @importFrom utils read.csv
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time",
detectorExclude="-H|-W", SignatureData, returnType="fcs"){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
WorkingColumns <- TheRawDataFrame |>
dplyr::select(!dplyr::matches(retainThese))
WorkingColumns <- WorkingColumns |>
dplyr::select(!dplyr::matches(detectorExclude))
# TheColNames <- colnames(WorkingColumns)
# Load the signature matrix
if(is.data.frame(SignatureData)){
Signatures <- SignatureData
} else {
Signatures <-read.csv(SignatureData, check.names=FALSE)
}
Metadata <- Signatures |>
dplyr::select(!dplyr::where(is.numeric))
Numerics <- Signatures |>
dplyr::select(dplyr::where(is.numeric))
if (any(Numerics > 1)) {
message("Signature values greater than 1 detected, normalizing")
n <- Numerics
# n[n < 0] <- 0
A <- do.call(pmax, n)
Normalized <- n/A
Numerics <- Normalized
}
# Make sure dimensions match each other
if (!all(colnames(Numerics) == colnames(WorkingColumns))){
stop("colnames of SignatureData due not match the internal colnames of exprs")
}
DetectorNameBackups <- colnames(Numerics)
TransposedSignatureValues <- t(Numerics)
TransposedSampleValues <- t(WorkingColumns)
# Unmixing
LeastSquares <- lsfit(x = TransposedSignatureValues,
y = TransposedSampleValues, intercept = FALSE)
# Residual Plot Fork
if (returnType == "residuals"){
Plot <- ResidualPlots(LeastSquaresList = LeastSquares)
return(Plot)
}
# Reassembly
TransposedLeastSquares <- t(LeastSquares$coefficients)
return(TransposedLeastSquares)
}.
And before we go further, lets make sure we keep the updated signature matrix with both Zombie NIR and Unstained (ideally reorganized in correct sequence, and unstained renamed to AF).
[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" "APC-Fire 750" "APC-Fire 810"
[29] "Zombie NIR" "PBMC_Unstained" RelocateElements <- function(data, from, after) {
index <- seq_len(nrow(data))
index <- index[index != from]
insert_at <- which(index == after)
index <- append(index, from, after = insert_at)
data[index, ]
}
UpdatedBeadReferences <- RelocateElements(data=UpdatedCellReferences, from=29, after=27)
UpdatedBeadReferences <- UpdatedBeadReferences |>
mutate(Fluorophore=case_when(
Fluorophore == "PBMC_Unstained" ~ "AF",
TRUE ~ Fluorophore))
tail(UpdatedBeadReferences, 5) Fluorophore Antigen UV1-A UV2-A UV3-A UV4-A UV5-A UV6-A UV7-A UV8-A UV9-A
26 APC-R700 CD107a 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
27 APC-Fire 750 CD27 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
29 Zombie NIR Viability 0.001 0.003 0.002 0.003 0.004 0.005 0.009 0.007 0.007
28 APC-Fire 810 CD38 0.000 0.000 0.000 0.000 0.000 0.000 0.001 0.001 0.001
30 AF 0.046 0.087 0.075 0.095 0.146 0.251 0.537 0.405 0.372
UV10-A UV11-A UV12-A UV13-A UV14-A UV15-A UV16-A V1-A V2-A V3-A V4-A
26 0.000 0.002 0.016 0.042 0.034 0.020 0.012 0.000 0.000 0.000 0.000
27 0.000 0.001 0.000 0.001 0.011 0.057 0.053 0.000 0.000 0.001 0.001
29 0.004 0.003 0.002 0.009 0.066 0.075 0.034 0.002 0.006 0.010 0.011
28 0.000 0.004 0.002 0.002 0.003 0.021 0.093 0.000 0.001 0.001 0.001
30 0.151 0.092 0.056 0.044 0.054 0.041 0.035 0.076 0.258 0.434 0.505
V5-A V6-A V7-A V8-A V9-A V10-A V11-A V12-A V13-A V14-A V15-A V16-A
26 0.000 0.000 0.001 0.000 0.000 0.000 0.006 0.041 0.112 0.064 0.042 0.017
27 0.001 0.001 0.001 0.001 0.001 0.001 0.003 0.001 0.002 0.024 0.142 0.087
29 0.017 0.016 0.024 0.020 0.014 0.016 0.010 0.008 0.041 0.223 0.280 0.081
28 0.002 0.002 0.002 0.002 0.001 0.002 0.012 0.006 0.006 0.006 0.063 0.202
30 0.768 0.730 1.000 0.746 0.527 0.609 0.350 0.204 0.195 0.180 0.164 0.097
B1-A B2-A B3-A B4-A B5-A B6-A B7-A B8-A B9-A B10-A B11-A B12-A
26 0.000 0.000 0.000 0.000 0.000 0.000 0.001 0.003 0.018 0.025 0.015 0.009
27 0.000 0.000 0.000 0.000 0.000 0.000 0.001 0.000 0.000 0.000 0.001 0.005
29 0.009 0.011 0.016 0.011 0.010 0.008 0.007 0.005 0.009 0.030 0.114 0.276
28 0.001 0.001 0.001 0.001 0.001 0.001 0.002 0.001 0.001 0.001 0.001 0.001
30 0.248 0.302 0.418 0.284 0.248 0.201 0.148 0.110 0.118 0.075 0.058 0.054
B13-A B14-A YG1-A YG2-A YG3-A YG4-A YG5-A YG6-A YG7-A YG8-A YG9-A YG10-A
26 0.006 0.006 0.000 0.000 0.000 0.013 0.036 0.158 0.500 0.168 0.108 0.055
27 0.014 0.017 0.000 0.000 0.001 0.007 0.004 0.003 0.005 0.043 0.232 0.185
29 0.231 0.153 0.006 0.006 0.005 0.006 0.006 0.007 0.041 0.112 0.126 0.048
28 0.006 0.035 0.001 0.001 0.003 0.032 0.021 0.016 0.021 0.012 0.074 0.298
30 0.045 0.058 0.128 0.128 0.158 0.145 0.112 0.101 0.121 0.056 0.050 0.031
R1-A R2-A R3-A R4-A R5-A R6-A R7-A R8-A
26 0.028 0.121 0.640 1.000 0.700 0.374 0.349 0.155
27 0.015 0.013 0.008 0.010 0.043 0.291 1.000 0.604
29 0.009 0.012 0.030 0.151 0.595 1.000 0.970 0.316
28 0.063 0.061 0.052 0.043 0.037 0.051 0.372 1.000
30 0.039 0.049 0.051 0.054 0.031 0.028 0.028 0.017 1 2 3 4 5 6
Y1 698.9385 142103.07416 4693.3116 40.70749 504.9595 -105.37370
Y2 431.6759 494.90676 33796.1594 -837.73170 609.7961 -63.70081
Y3 11175.8241 -56.28136 572.2067 1169.13066 -512.1527 -566.20791
7 8 9 10 11 12 13
Y1 -288.1342 3795.2743 -394.3967 1440.6184 4549.111 947.5550 14926.849
Y2 1813.3706 440.2116 4273.7052 247.4025 26076.313 207.4692 4482.672
Y3 30539.7372 7388.5224 -623.5300 577.2415 5688.807 -678.7335 16255.307
14 15 16 17 18 19 20
Y1 97542.8745 1733.8254 -211.4355 -604.68497 45958.6701 2098.7171 55249.1340
Y2 -270.8378 1423.2353 3333.3268 68.41591 693.3397 424.2276 717.9945
Y3 162022.4787 -770.5102 2684.0677 -97.15930 45425.0009 2690.2426 552.9612
21 22 23 24 25 26 27
Y1 -515.4923 386.9981 4409.8038 1867.690 -1011.577 53.23036 34909.681
Y2 766.2239 194.9767 762.6181 2758.993 1524.624 6348.41146 -2686.721
Y3 -388.2691 357.2628 2599.5774 3812.097 2123.882 -2193.30755 18450.067
29 28 30
Y1 707.1137 27332.27 6806.192
Y2 1887.3333 52433.72 4700.989
Y3 392.5792 29704.08 7383.022.
So we are now in the wide format, which is expected by exprs(). However, we don’t have any colnames(), so its time to overwrite this a vector of the Fluorophore names. For Cytek Aurora instruments, the fluorophore names typically have a designated letter (“-A”) appended at the end, giving the characteristic display seen in many unmixed .fcs files (FITC-A, etc.).
.
We can try to generalize this process based on our kept column names, to avoid yet another argument to specify.
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#' @param detectorExclude Default is "-H|-W", intended to remove
#' additional detector columns other than -A, adjust as needed for
#' your own instruments configuration
#' @param SignatureData A data.frame containing a Fluorophore,
#' Antigen and Detector columns.
#' @param returnType Default "fcs", for residual plots use "residual"
#'
#' @importFrom flowCore read.FCS exprs
#' @importFrom dplyr select matches where pull
#' @importFrom utils read.csv
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time",
detectorExclude="-H|-W", SignatureData, returnType="fcs"){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
WorkingColumns <- TheRawDataFrame |>
dplyr::select(!dplyr::matches(retainThese))
WorkingColumns <- WorkingColumns |>
dplyr::select(!dplyr::matches(detectorExclude))
# TheColNames <- colnames(WorkingColumns)
# Load the signature matrix
if(is.data.frame(SignatureData)){
Signatures <- SignatureData
} else {
Signatures <-read.csv(SignatureData, check.names=FALSE)
}
Metadata <- Signatures |>
dplyr::select(!dplyr::where(is.numeric))
Numerics <- Signatures |>
dplyr::select(dplyr::where(is.numeric))
if (any(Numerics > 1)) {
message("Signature values greater than 1 detected, normalizing")
n <- Numerics
# n[n < 0] <- 0
A <- do.call(pmax, n)
Normalized <- n/A
Numerics <- Normalized
}
# Make sure dimensions match each other
if (!all(colnames(Numerics) == colnames(WorkingColumns))){
stop("colnames of SignatureData due not match the internal colnames of exprs")
}
DetectorNameBackups <- colnames(Numerics)
TransposedSignatureValues <- t(Numerics)
TransposedSampleValues <- t(WorkingColumns)
# Unmixing
LeastSquares <- lsfit(x = TransposedSignatureValues,
y = TransposedSampleValues, intercept = FALSE)
# Residual Plot Fork
if (returnType == "residuals"){
Plot <- ResidualPlots(LeastSquaresList = LeastSquares)
return(Plot)
}
# Reassembly
TransposedLeastSquares <- t(LeastSquares$coefficients)
FluorophoreNames <- Metadata |> dplyr::pull(Fluorophore)
TheDetectorColNames <- colnames(WorkingColumns)
AppendThisLetter <- sub("^[^-]*", "", TheDetectorColNames) |> unique()
FluorophoreNames <- paste0(FluorophoreNames, AppendThisLetter)
colnames(TransposedLeastSquares) <- FluorophoreNames
return(TransposedLeastSquares)
} BUV395-A BUV496-A BUV563-A BUV615-A BUV661-A BUV737-A
Y1 698.9385 142103.07416 4693.3116 40.70749 504.9595 -105.37370
Y2 431.6759 494.90676 33796.1594 -837.73170 609.7961 -63.70081
Y3 11175.8241 -56.28136 572.2067 1169.13066 -512.1527 -566.20791
BUV805-A BV421-A Pacific Blue-A BV480-A BV510-A BV605-A BV650-A
Y1 -288.1342 3795.2743 -394.3967 1440.6184 4549.111 947.5550 14926.849
Y2 1813.3706 440.2116 4273.7052 247.4025 26076.313 207.4692 4482.672
Y3 30539.7372 7388.5224 -623.5300 577.2415 5688.807 -678.7335 16255.307
BV711-A BV750-A BV786-A FITC-A Spark Blue 550-A PerCP-Cy5.5-A
Y1 97542.8745 1733.8254 -211.4355 -604.68497 45958.6701 2098.7171
Y2 -270.8378 1423.2353 3333.3268 68.41591 693.3397 424.2276
Y3 162022.4787 -770.5102 2684.0677 -97.15930 45425.0009 2690.2426
PE-A PE-Dazzle 594-A PE-Cy5-A PE-Vio 770-A APC-A Alexa Fluor 647-A
Y1 55249.1340 -515.4923 386.9981 4409.8038 1867.690 -1011.577
Y2 717.9945 766.2239 194.9767 762.6181 2758.993 1524.624
Y3 552.9612 -388.2691 357.2628 2599.5774 3812.097 2123.882
APC-R700-A APC-Fire 750-A Zombie NIR-A APC-Fire 810-A AF-A
Y1 53.23036 34909.681 707.1137 27332.27 6806.192
Y2 6348.41146 -2686.721 1887.3333 52433.72 4700.989
Y3 -2193.30755 18450.067 392.5792 29704.08 7383.022.
With our unmixed fluorophore columns now set, we can put back the stashed columns for “Time”, “SSC”, “FSC” and “SSC-B”.
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#' @param detectorExclude Default is "-H|-W", intended to remove
#' additional detector columns other than -A, adjust as needed for
#' your own instruments configuration
#' @param SignatureData A data.frame containing a Fluorophore,
#' Antigen and Detector columns.
#' @param returnType Default "fcs", for residual plots use "residual"
#'
#' @importFrom flowCore read.FCS exprs
#' @importFrom dplyr select matches where pull bind_cols
#' @importFrom utils read.csv
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time",
detectorExclude="-H|-W", SignatureData, returnType="fcs"){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
WorkingColumns <- TheRawDataFrame |>
dplyr::select(!dplyr::matches(retainThese))
WorkingColumns <- WorkingColumns |>
dplyr::select(!dplyr::matches(detectorExclude))
# TheColNames <- colnames(WorkingColumns)
# Load the signature matrix
if(is.data.frame(SignatureData)){
Signatures <- SignatureData
} else {
Signatures <-read.csv(SignatureData, check.names=FALSE)
}
Metadata <- Signatures |>
dplyr::select(!dplyr::where(is.numeric))
Numerics <- Signatures |>
dplyr::select(dplyr::where(is.numeric))
if (any(Numerics > 1)) {
message("Signature values greater than 1 detected, normalizing")
n <- Numerics
# n[n < 0] <- 0
A <- do.call(pmax, n)
Normalized <- n/A
Numerics <- Normalized
}
# Make sure dimensions match each other
if (!all(colnames(Numerics) == colnames(WorkingColumns))){
stop("colnames of SignatureData due not match the internal colnames of exprs")
}
DetectorNameBackups <- colnames(Numerics)
TransposedSignatureValues <- t(Numerics)
TransposedSampleValues <- t(WorkingColumns)
# Unmixing
LeastSquares <- lsfit(x = TransposedSignatureValues,
y = TransposedSampleValues, intercept = FALSE)
# Residual Plot Fork
if (returnType == "residuals"){
Plot <- ResidualPlots(LeastSquaresList = LeastSquares)
return(Plot)
}
# Reassembly
TransposedLeastSquares <- t(LeastSquares$coefficients)
FluorophoreNames <- Metadata |> dplyr::pull(Fluorophore)
TheDetectorColNames <- colnames(WorkingColumns)
AppendThisLetter <- sub("^[^-]*", "", TheDetectorColNames) |> unique()
FluorophoreNames <- paste0(FluorophoreNames, AppendThisLetter)
colnames(TransposedLeastSquares) <- FluorophoreNames
# Bind the stashed Time, SSC and FSC columns to the new fluorophore columns
UnmixedData <- bind_cols(StashedColumns, TransposedLeastSquares)
return(UnmixedData)
} Time SSC-W SSC-H SSC-A FSC-W FSC-H FSC-A SSC-B-W SSC-B-H
1 887309 683436.8 843643 960961.0 696000.9 1473050 1708740 676320.5 608196
2 1393476 658488.6 1080396 1185714.0 665755.6 1411870 1566601 660120.8 782257
3 1258794 756646.2 714195 900654.8 723983.9 1276900 1540758 715180.6 554999
SSC-B-A BUV395-A BUV496-A BUV563-A BUV615-A BUV661-A BUV737-A
1 685559.1 698.9385 142103.07416 4693.3116 40.70749 504.9595 -105.37370
2 860640.1 431.6759 494.90676 33796.1594 -837.73170 609.7961 -63.70081
3 661540.9 11175.8241 -56.28136 572.2067 1169.13066 -512.1527 -566.20791
BUV805-A BV421-A Pacific Blue-A BV480-A BV510-A BV605-A BV650-A
1 -288.1342 3795.2743 -394.3967 1440.6184 4549.111 947.5550 14926.849
2 1813.3706 440.2116 4273.7052 247.4025 26076.313 207.4692 4482.672
3 30539.7372 7388.5224 -623.5300 577.2415 5688.807 -678.7335 16255.307
BV711-A BV750-A BV786-A FITC-A Spark Blue 550-A PerCP-Cy5.5-A
1 97542.8745 1733.8254 -211.4355 -604.68497 45958.6701 2098.7171
2 -270.8378 1423.2353 3333.3268 68.41591 693.3397 424.2276
3 162022.4787 -770.5102 2684.0677 -97.15930 45425.0009 2690.2426
PE-A PE-Dazzle 594-A PE-Cy5-A PE-Vio 770-A APC-A Alexa Fluor 647-A
1 55249.1340 -515.4923 386.9981 4409.8038 1867.690 -1011.577
2 717.9945 766.2239 194.9767 762.6181 2758.993 1524.624
3 552.9612 -388.2691 357.2628 2599.5774 3812.097 2123.882
APC-R700-A APC-Fire 750-A Zombie NIR-A APC-Fire 810-A AF-A
1 53.23036 34909.681 707.1137 27332.27 6806.192
2 6348.41146 -2686.721 1887.3333 52433.72 4700.989
3 -2193.30755 18450.067 392.5792 29704.08 7383.022.
And with that, congratulations, we now have our unmixed data in a format that can be returned to the exprs() slot.
UnmixInternal Extravaganza
.
At this point we have assembled a “data.frame” containing the unmixed data for each cell in our sample. However, if you remember the Week 3 and Week 10 walk-throughs, you may remember how interconnected the elements between the different flowFrame slots (exprs(), parameters(), keyword()) slots actually are.
.
In this case, instead of creating a new .fcs file entirely from scratch, we need to retain certain aspects of the original raw .fcs files metadata. However, we need to also remove the metadata associated with the detector columns (which are no longer present), and add new metadata for our newly derived unmixed fluorophore columns.
.
Since there are multiple moving pieces, this is the job for a nested internal function. For anyone who nerds out about .fcs file internals, you can find the additional details here in this bonus walkthrough. For the majority who would rather not, the final function can be found in the code-chunk below (as well as in a “.R” file in this week’s R folder for future use).
Code
#' Internal function for OldFashinedUnmix, handles fixing the
#' formatting going from raw to unmixed .fcs file
#'
#' @param ff The original raw flowFrame (used for the initial metadata)
#' @param data The updated data.frame following unmixing
#' @param panel Metadata from the signature matrix, containing
#' Fluorophore and Antigen
#'
#' @importFrom flowCore parameters keyword
#' @importFrom dplyr filter mutate row_number relocate pull
#' select
#' @importFrom tibble rownames_to_column column_to_rownames
#
UnmixInternal <- function(ff, data, panel){
# Identify the retained columns
TheOriginalColumns <- colnames(ff)
TheNewColumns <- colnames(data)
RetainedColumns <- intersect(TheOriginalColumns, TheNewColumns)
# Retrieve original parameters data with "$P" rownames
TheOriginalParameters <- flowCore::parameters(ff)@data
# Identify "$P" rownames to eliminate or keep
KeepThese <- TheOriginalParameters |>
dplyr::filter(name %in% RetainedColumns)
GetRidThese <- TheOriginalParameters |>
dplyr::filter(!name %in% RetainedColumns) |> rownames()
# Identify existing keyword names
TheOriginalDescription <- flowCore::keyword(ff)
OriginalKeywordNames <- names(TheOriginalDescription)
# Identification of keywords containing "$P"
Escaped <- gsub("\\$", "\\\\$", GetRidThese)
RegexFormatted <- paste0("^", Escaped, "($|[^0-9])")
RegexCombinatorial <- paste(RegexFormatted, collapse = "|")
IdentifiedElimination <- OriginalKeywordNames[
grepl(RegexCombinatorial, OriginalKeywordNames)]
IdentifiedRetention <- OriginalKeywordNames[
!grepl(RegexCombinatorial, OriginalKeywordNames)]
# Identification of keywords containing "flowCore_$P"
SecondRegexPattern <- paste0("(?<!\\d)", Escaped, "(?!\\d)")
SecondRegexCombinatorial <- paste(SecondRegexPattern, collapse = "|")
SecondElimination <- IdentifiedRetention[
grepl(SecondRegexCombinatorial, IdentifiedRetention, perl = TRUE)]
SecondRetention <- IdentifiedRetention[
!grepl(SecondRegexCombinatorial, IdentifiedRetention, perl = TRUE)]
# Identification of keywords containing "PDisplay"
NoDollars <- sub("^\\$", "", GetRidThese)
NoDollarsRegex <- paste0("(?<!\\d)\\$?", NoDollars, "(?!\\d)")
NoDollarsCombined <- paste(NoDollarsRegex, collapse = "|")
ThirdElimination <- SecondRetention[
grepl(NoDollarsCombined, SecondRetention, perl = TRUE)]
ThirdRetention <- SecondRetention[
!grepl(NoDollarsCombined, SecondRetention, perl = TRUE)]
# Subset original description for keyword names we want to retain
IntermediateDescription <- TheOriginalDescription[ThirdRetention]
# Renumbering retained SSC, FSC, SSC-B columns in parameters
IntermediateParameters <- KeepThese |>
tibble::rownames_to_column("OriginalRowNumber") |>
dplyr::mutate(NewRowNumber=paste0("$P", dplyr::row_number())) |>
relocate(NewRowNumber, .before=1)
# Renumbering the retained keywords with new "$P" row numbers
NewRowNumbers <- IntermediateParameters |>
dplyr::pull(NewRowNumber)
OriginalRowNumbers <- IntermediateParameters |>
dplyr::pull(OriginalRowNumber)
ForLoopDescription <- IntermediateDescription
# For-loop to renumber the retained keywords
for (i in seq_along(NewRowNumbers)) {
NewX <- NewRowNumbers[i]
NewXNum <- sub("^\\$P", "", NewX)
OldX <- OriginalRowNumbers[i]
OldXNum <- sub("^\\$P", "", OldX)
# Rename matching keywords with new row numbers
InternalEscaped <- gsub("\\$", "\\\\$", OldX)
InternalRegexFormatted <- paste0(
"(?<!\\d)", InternalEscaped, "(?!\\d)")
InternalRegexCombinatorial <- paste(
InternalRegexFormatted, collapse = "|")
InternalIdentifiedRename <- names(ForLoopDescription)[
grepl(InternalRegexCombinatorial,
names(ForLoopDescription), perl = TRUE)]
## Actual renaming in action
RenamePattern <- paste0("^(flowCore_)?(\\$)?P", OldXNum, "(.*)$")
InternalRenamed <- sub(RenamePattern,
paste0("\\1\\2P", NewXNum, "\\3"), InternalIdentifiedRename)
names(ForLoopDescription)[match(InternalIdentifiedRename,
names(ForLoopDescription))] <- InternalRenamed
# Matching for rename the "PDisplay" keywords
InternalNoDollars <- sub("^\\$", "", OldX) # "P18"
InternalNoDollarsRegex <- paste0(
"(?<!\\d)\\$?", InternalNoDollars, "(?!\\d)")
ThirdRename <- names(ForLoopDescription)[
grepl(InternalNoDollarsRegex, names(ForLoopDescription),
perl = TRUE)]
# Final renaming in action for "PDisplay"
DisplayRenamePattern <- paste0("^(\\$)?P", OldXNum, "(.*)$")
ThirdRenamed <- sub(DisplayRenamePattern,
paste0("\\1P", NewXNum, "\\2"), ThirdRename)
names(ForLoopDescription)[match(ThirdRename,
names(ForLoopDescription))] <- ThirdRenamed
}
# Add the new parameter rows for the Unmixed Fluorophores
IntermediateParameters <- IntermediateParameters |>
dplyr::select(-OriginalRowNumber) |>
tibble::column_to_rownames("NewRowNumber")
UpdatedParameters <- UnmixedParameterUpdate(OldParameters=IntermediateParameters, NewExprs=data)
TheAntigens <- panel |> pull(Antigen)
UpdatedParameters$desc <- TheAntigens
pd <- rbind(IntermediateParameters, UpdatedParameters)
# Generate new keywords for the new fluorophore rows in parameters
new_pid <- rownames(UpdatedParameters)
new_kw <- ForLoopDescription
for (i in new_pid){
NoDollarCode <- sub("^\\$", "", i) # For "PDisplay keyword"
new_kw[paste0(i,"B")] <- new_kw["$P1B"] # Bytes
new_kw[paste0(i,"E")] <- "0,0"
new_kw[paste0(i,"N")] <- pd[[i,1]] # Fluorophore Name
new_kw[paste0(i,"R")] <- pd[[i,5]] # Range Default
new_kw[paste0(i,"S")] <- pd[[i,2]] # Antigen Name
new_kw[paste0(i,"TYPE")] <- "Unmixed_Fluorescence"
new_kw[paste0(i,"V")] <- "0" # Voltage/Gain, unmixed default 0
new_kw[paste0("flowCore_", i,"Rmax")] <- pd[[i,5]] # maxRange
new_kw[paste0("flowCore_", i,"Rmin")] <- pd[[i,4]] # minRange
new_kw[paste0(NoDollarCode,"DISPLAY")] <- "LOG"
}
# Order Keywords by default sequence
new_kw <- new_kw[order(names(new_kw))]
# Overwrite old parameters "data" with the new parameters
OriginalParameterSlot <- flowCore::parameters(ff)
OriginalParameterSlot@data <- pd
# Convert unmixed data.frame to unmixed matrix
UnmixedMatrix <- as.matrix(data)
# Last keyword overrides
new_kw$`CREATOR` <- "CytometryInR version 1.0.0"
TheColumnNames <- colnames(data)
TheSpilloverNames <- TheColumnNames[!grepl("Time|FSC|SSC", TheColumnNames)]
MatrixSize <- length(TheSpilloverNames)
NewMatrix <- matrix(0, nrow = MatrixSize, ncol = MatrixSize, byrow = TRUE)
diag(NewMatrix) <- 1
colnames(NewMatrix) <- TheSpilloverNames
new_kw$`$SPILLOVER` <- NewMatrix
new_fcs <- new("flowFrame", exprs=UnmixedMatrix, parameters=OriginalParameterSlot,
description=new_kw)
return(new_fcs)
}
#' Internal for UnmixInternal, creates the new parameter
#' data rows needed to properly integrate new fluorophore columns
#' in exprs matrix
#'
#' @param OldParameters The parameter data.frame with
#' modified row numbers
#' @param NewExprs The unmixed data that will eventually
#' be placed back into exprs
#'
#' @importFrom Biobase pData
#' @importFrom dplyr pull select
#' @importFrom tidyselect all_of
#'
UnmixedParameterUpdate <- function(OldParameters, NewExprs){
# Remove the Retaind
OldNames <- OldParameters |> dplyr::pull(name) |> unname()
Overlapped <- intersect(OldNames, colnames(NewExprs))
NewExprs <- NewExprs |> select(-all_of(Overlapped))
# Create new rows for the unmixed fluorophore columns
NewColumnLength <- ncol(NewExprs)
NewColumnNames <- colnames(NewExprs)
NewParameter <- max(as.integer(gsub("\\$P", "", rownames(OldParameters)))) + 1
NewParameter <- seq(NewParameter, length.out = NewColumnLength)
NewParameter <- paste0("$P", NewParameter)
# Provide Hard Coded Numbers for Range, minRange and maxRange
SSCRange <- OldParameters[2,3] # Hard-Coded based on Cytek Aurora SpectroFlo unmixed value
MinRange <- -111.0001 # Hard-Coded based on Cytek Aurora SpectroFlo unmixed value
MaxRange <- 4192505.7500 # Hard-Coded based on Cytek Aurora SpectroFlo unmixed value
UpdatedParameters <- do.call(rbind, lapply(NewColumnNames, function(i){
vec <- NewExprs[,i]
rg <- range(vec)
data.frame(name = i,
desc = NA,
range = SSCRange,
minRange = MinRange,
maxRange = MaxRange)
}))
rownames(UpdatedParameters) <- NewParameter
return(UpdatedParameters)
}.
And with that, add the line for UnmixInternal() into OldFashionedUnmix(), and see if it returns a properly formatted flowFrame.
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#' @param detectorExclude Default is "-H|-W", intended to remove
#' additional detector columns other than -A, adjust as needed for
#' your own instruments configuration
#' @param SignatureData A data.frame containing a Fluorophore,
#' Antigen and Detector columns.
#' @param returnType Default "fcs", for residual plots use "residual"
#'
#' @importFrom flowCore read.FCS exprs
#' @importFrom dplyr select matches where pull bind_cols
#' @importFrom utils read.csv
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time",
detectorExclude="-H|-W", SignatureData, returnType="fcs"){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
WorkingColumns <- TheRawDataFrame |>
dplyr::select(!dplyr::matches(retainThese))
WorkingColumns <- WorkingColumns |>
dplyr::select(!dplyr::matches(detectorExclude))
# TheColNames <- colnames(WorkingColumns)
# Load the signature matrix
if(is.data.frame(SignatureData)){
Signatures <- SignatureData
} else {
Signatures <-read.csv(SignatureData, check.names=FALSE)
}
Metadata <- Signatures |>
dplyr::select(!dplyr::where(is.numeric))
Numerics <- Signatures |>
dplyr::select(dplyr::where(is.numeric))
if (any(Numerics > 1)) {
message("Signature values greater than 1 detected, normalizing")
n <- Numerics
# n[n < 0] <- 0
A <- do.call(pmax, n)
Normalized <- n/A
Numerics <- Normalized
}
# Make sure dimensions match each other
if (!all(colnames(Numerics) == colnames(WorkingColumns))){
stop("colnames of SignatureData due not match the internal colnames of exprs")
}
DetectorNameBackups <- colnames(Numerics)
TransposedSignatureValues <- t(Numerics)
TransposedSampleValues <- t(WorkingColumns)
# Unmixing
LeastSquares <- lsfit(x = TransposedSignatureValues,
y = TransposedSampleValues, intercept = FALSE)
# Residual Plot Fork
if (returnType == "residuals"){
Plot <- ResidualPlots(LeastSquaresList = LeastSquares)
return(Plot)
}
# Reassembly
TransposedLeastSquares <- t(LeastSquares$coefficients)
FluorophoreNames <- Metadata |> dplyr::pull(Fluorophore)
TheDetectorColNames <- colnames(WorkingColumns)
AppendThisLetter <- sub("^[^-]*", "", TheDetectorColNames) |> unique()
FluorophoreNames <- paste0(FluorophoreNames, AppendThisLetter)
colnames(TransposedLeastSquares) <- FluorophoreNames
# Bind the stashed Time, SSC and FSC columns to the new fluorophore columns
UnmixedData <- bind_cols(StashedColumns, TransposedLeastSquares)
# Return properly formatted flowFrame
new_fcs <- UnmixInternal(ff=TheRawFCS, data=UnmixedData, panel=Metadata)
return(new_fcs)
}[[1]]
flowFrame object 'DTR_2023_ILT_01-INF052-Ctrl_Antibody.1235515.fcs'
with 10000 cells and 40 observables:
name desc range minRange maxRange
$P1 Time NA 1428432 0 1428431
$P2 SSC-W NA 4194304 0 4194303
$P3 SSC-H NA 4194304 0 4194303
$P4 SSC-A NA 4194304 0 4194303
$P5 FSC-W NA 4194304 0 4194303
... ... ... ... ... ...
$P36 APC-R700-A CD107a 4194304 -111 4192506
$P37 APC-Fire 750-A CD27 4194304 -111 4192506
$P38 Zombie NIR-A Viability 4194304 -111 4192506
$P39 APC-Fire 810-A CD38 4194304 -111 4192506
$P40 AF-A 4194304 -111 4192506
443 keywords are stored in the 'description' slot.
And with that, we are in the final stretch!
Final Push
.
With UnmixInternal() now taking our “Unmixed Fluorophore” “data.frame”, and returning a properly formatted unmixed “flowFrame”, all we need to do is any last minute formatting changes to our file name keywords, as well as designate where we want the .fcs file to be saved to on our computers.
Renaming
.
For renaming, lets give the option to rename on the basis of up to two existing keywords (handled via a new “sample.name” argument), as well as an “addon” argument that can be used to append “_Unmixed” before “.fcs” in the name.
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#' @param detectorExclude Default is "-H|-W", intended to remove
#' additional detector columns other than -A, adjust as needed for
#' your own instruments configuration
#' @param SignatureData A data.frame containing a Fluorophore,
#' Antigen and Detector columns.
#' @param returnType Default "fcs", for residual plots use "residual"
#' @param sample.name Keywords to pull for the new name, Cytek Aurora
#' default keywords are c("GROUPNAME", "TUBENAME").
#' @param addon Default "Unmixed", added right before the ".fcs"
#'
#' @importFrom flowCore read.FCS exprs keyword
#' @importFrom dplyr select matches where pull bind_cols
#' @importFrom utils read.csv
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time",
detectorExclude="-H|-W", SignatureData, returnType="fcs",
sample.name=c("GROUPNAME", "TUBENAME"), addon="Unmixed"){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
WorkingColumns <- TheRawDataFrame |>
dplyr::select(!dplyr::matches(retainThese))
WorkingColumns <- WorkingColumns |>
dplyr::select(!dplyr::matches(detectorExclude))
# TheColNames <- colnames(WorkingColumns)
# Load the signature matrix
if(is.data.frame(SignatureData)){
Signatures <- SignatureData
} else {
Signatures <-read.csv(SignatureData, check.names=FALSE)
}
Metadata <- Signatures |>
dplyr::select(!dplyr::where(is.numeric))
Numerics <- Signatures |>
dplyr::select(dplyr::where(is.numeric))
if (any(Numerics > 1)) {
message("Signature values greater than 1 detected, normalizing")
n <- Numerics
# n[n < 0] <- 0
A <- do.call(pmax, n)
Normalized <- n/A
Numerics <- Normalized
}
# Make sure dimensions match each other
if (!all(colnames(Numerics) == colnames(WorkingColumns))){
stop("colnames of SignatureData due not match the internal colnames of exprs")
}
DetectorNameBackups <- colnames(Numerics)
TransposedSignatureValues <- t(Numerics)
TransposedSampleValues <- t(WorkingColumns)
# Unmixing
LeastSquares <- lsfit(x = TransposedSignatureValues,
y = TransposedSampleValues, intercept = FALSE)
# Residual Plot Fork
if (returnType == "residuals"){
Plot <- ResidualPlots(LeastSquaresList = LeastSquares)
return(Plot)
}
# Reassembly
TransposedLeastSquares <- t(LeastSquares$coefficients)
FluorophoreNames <- Metadata |> dplyr::pull(Fluorophore)
TheDetectorColNames <- colnames(WorkingColumns)
AppendThisLetter <- sub("^[^-]*", "", TheDetectorColNames) |> unique()
FluorophoreNames <- paste0(FluorophoreNames, AppendThisLetter)
colnames(TransposedLeastSquares) <- FluorophoreNames
# Bind the stashed Time, SSC and FSC columns to the new fluorophore columns
UnmixedData <- dplyr::bind_cols(StashedColumns, TransposedLeastSquares)
# Return properly formatted flowFrame
new_fcs <- UnmixInternal(ff=TheRawFCS, data=UnmixedData, panel=Metadata)
# Provide a new name to GUID and FIL
if (length(sample.name) == 2){
first <- sample.name[[1]]
second <- sample.name[[2]]
first <- keyword(new_fcs, first)
second <- keyword(new_fcs, second)
name <- paste(first, second, sep="_")
} else {name <- keyword(new_fcs, sample.name)}
if (!is.null(addon)){name <- paste0(name, "_", addon)}
AssembledName <- paste0(name, ".fcs")
new_fcs@description$GUID <- AssembledName
new_fcs@description$`$FIL` <- AssembledName
return(new_fcs)
}[[1]]
flowFrame object 'INF052_Ctrl_Antibody_Unmixed.fcs'
with 10000 cells and 40 observables:
name desc range minRange maxRange
$P1 Time NA 1428432 0 1428431
$P2 SSC-W NA 4194304 0 4194303
$P3 SSC-H NA 4194304 0 4194303
$P4 SSC-A NA 4194304 0 4194303
$P5 FSC-W NA 4194304 0 4194303
... ... ... ... ... ...
$P36 APC-R700-A CD107a 4194304 -111 4192506
$P37 APC-Fire 750-A CD27 4194304 -111 4192506
$P38 Zombie NIR-A Viability 4194304 -111 4192506
$P39 APC-Fire 810-A CD38 4194304 -111 4192506
$P40 AF-A 4194304 -111 4192506
443 keywords are stored in the 'description' slotSaving
.
And last, but definitely not least, we need to designate where to store the file. We can copy in some of the previous code we used as part of Week 10 to simplify this process.
#' This function unmixes our raw full-stained .fcs files using the
#' signature matrix we provide.
#'
#' @param x A file.path to a raw full-stained .fcs file we want to
#' unmix.
#' @param retainThese Default "FSC|SSC|Time", used to separate out
#' columns not used for unmixing, but that should be retained for
#' the final unmixed .fcs files.
#' @param detectorExclude Default is "-H|-W", intended to remove
#' additional detector columns other than -A, adjust as needed for
#' your own instruments configuration
#' @param SignatureData A data.frame containing a Fluorophore,
#' Antigen and Detector columns.
#' @param returnType Default "fcs", for residual plots use "residual"
#' @param sample.name Keywords to pull for the new name, Cytek Aurora
#' default keywords are c("GROUPNAME", "TUBENAME").
#' @param addon Default "Unmixed", added right before the ".fcs"
#' @param outpath File.path to where we want to store the .fcs files,
#' Default is NULL, which results in being saved to working directory
#'
#' @importFrom flowCore read.FCS exprs write.FCS keyword
#' @importFrom dplyr select matches where pull bind_cols
#' @importFrom utils read.csv
#'
OldFashionedUnmix <- function(x, retainThese="FSC|SSC|Time",
detectorExclude="-H|-W", SignatureData, returnType="fcs",
sample.name=c("GROUPNAME", "TUBENAME"), addon="Unmixed",
outpath=NULL){
# Retrieve the underlying MFI values from exprs slot
TheRawFCS <- flowCore::read.FCS(filename=x,
transformation=FALSE, truncate_max_range = FALSE)
TheRawMatrix <- flowCore::exprs(TheRawFCS)
TheRawDataFrame <- data.frame(TheRawMatrix, check.names=FALSE)
# Identify the detector columns
StashedColumns <- TheRawDataFrame |>
dplyr::select(dplyr::matches(retainThese))
WorkingColumns <- TheRawDataFrame |>
dplyr::select(!dplyr::matches(retainThese))
WorkingColumns <- WorkingColumns |>
dplyr::select(!dplyr::matches(detectorExclude))
# TheColNames <- colnames(WorkingColumns)
# Load the signature matrix
if(is.data.frame(SignatureData)){
Signatures <- SignatureData
} else {
Signatures <-read.csv(SignatureData, check.names=FALSE)
}
Metadata <- Signatures |>
dplyr::select(!dplyr::where(is.numeric))
Numerics <- Signatures |>
dplyr::select(dplyr::where(is.numeric))
if (any(Numerics > 1)) {
message("Signature values greater than 1 detected, normalizing")
n <- Numerics
# n[n < 0] <- 0
A <- do.call(pmax, n)
Normalized <- n/A
Numerics <- Normalized
}
# Make sure dimensions match each other
if (!all(colnames(Numerics) == colnames(WorkingColumns))){
stop("colnames of SignatureData due not match the internal colnames of exprs")
}
DetectorNameBackups <- colnames(Numerics)
TransposedSignatureValues <- t(Numerics)
TransposedSampleValues <- t(WorkingColumns)
# Unmixing
LeastSquares <- lsfit(x = TransposedSignatureValues,
y = TransposedSampleValues, intercept = FALSE)
# Residual Plot Fork
if (returnType == "residuals"){
Plot <- ResidualPlots(LeastSquaresList = LeastSquares)
return(Plot)
}
# Reassembly
TransposedLeastSquares <- t(LeastSquares$coefficients)
FluorophoreNames <- Metadata |> dplyr::pull(Fluorophore)
TheDetectorColNames <- colnames(WorkingColumns)
AppendThisLetter <- sub("^[^-]*", "", TheDetectorColNames) |> unique()
FluorophoreNames <- paste0(FluorophoreNames, AppendThisLetter)
colnames(TransposedLeastSquares) <- FluorophoreNames
# Bind the stashed Time, SSC and FSC columns to the new fluorophore columns
UnmixedData <- dplyr::bind_cols(StashedColumns, TransposedLeastSquares)
# Return properly formatted flowFrame
new_fcs <- UnmixInternal(ff=TheRawFCS, data=UnmixedData, panel=Metadata)
# Provide a new name to GUID and FIL
if (length(sample.name) == 2){
first <- sample.name[[1]]
second <- sample.name[[2]]
first <- flowCore::keyword(new_fcs, first)
second <- flowCore::keyword(new_fcs, second)
name <- paste(first, second, sep="_")
} else {name <- flowCore::keyword(new_fcs, sample.name)}
if (!is.null(addon)){name <- paste0(name, "_", addon)}
AssembledName <- paste0(name, ".fcs")
new_fcs@description$GUID <- AssembledName
new_fcs@description$`$FIL` <- AssembledName
if (is.null(outpath)) {outpath <- getwd()}
fileSpot <- file.path(outpath, AssembledName)
if (returnType == "fcs") {
flowCore::write.FCS(new_fcs, filename = fileSpot, delimiter="#")
} else {return(new_fcs)}
}.
And if we check our designated output location…
.
Wooh! We have the .fcs file.
Evaluation
.
Alright, moment of truth. Let’s open it in either Floreada.io or if you have a license your propietary software of choice, and see if it actually behaves like a unmixed .fcs file.
.
Alright, no immediate red flags! It appears that we have been successful in our quest to unmix the raw .fcs files we retrieved from ImmPort entirely within R from start to finish.
.
As mentioned earlier, we used a signature matrix for unmixing that we haven’t fully optimized whether the signatures included were the best candidates, so we still need to evaluate how our unmixing looks, and based on any encountered unmixing errors update our decision making about what signatures we provide to the unmixing matrix.
.
We will get to this in the next primary walkthrough, but for now, the main tool in our unmixing arsenal is now set up to handle whichever optimized signature matrix we ultimately decide to provide it with.
.
Likewise, beyond using OLS, there are other methods via which we can also unmix. Some like non-negative least squares involve simply swapping out lsfit() for another R function. Others are more involved, and have different inputs, some of which was covered in the Autospectral as well as TRU-OLS community walkthroughs.
Take Away
.
In this walk-through, we have culminated some of the elements we have been building out over the last several sessions, and actually unmixed out our .fcs files. This is not something everyone can say that they have done, so congratulations on having made it this far.
.
As you saw from the process, the main elements involved are actually processing the data needed to unmix, and after unmix, formatting to prepare the unmixed .fcs file with the new columns. The actual unmixing function at the heart is rather swapable. Likewise, it is possible to swap this step out entirely for faster implementations in lower-level programming languages like Fortran, C++ and Rust which is what several R packages do behind the scenes.
.
Next time in the primary course material, we will look at how to evaluate our unmixing outputs, with the following week set to look at various algorithms that can help us screen our unmixed .fcs files for instrumental issues arising from fluidic and laser instability. We will also explore some additional comparison metrics in the form of spreading matrices and staining index in the bonus content over the next couple weeks.
.
Until then, cheers! We are now half way through the original conception of our course.
Additional Resources
Mathematics of spectral unmixing │Peter Mage │ Babraham Institute Spectral Symposium 2022
CytoBytes: Spectral cytometry have you feeling all mixed up? Let’s get unmixed!
Take-home Problems
Problem 1
Modify OldFashionedUnmix()’s returnType argument to instead of returning a “.fcs” file or a “residual” plot, to also returned the Unmixed data.frame, as well as the flowFrame (since both can be useful!). Add any conditional statements or arguments as you deem necessary.
Problem 2
We briefly looked at the residual plots, including a case where we left out BUV805 CD4 and saw corresponding increase on those detectors in the residual plot. Rerun this example with a signature matrix missing one of the fluorophores, selecting returnType = “residual” to get back the residual plot. At this point, rerun with returnType = “fcs”, and generate a badly unmixed .fcs file. Then repeat this process with a signature matrix that contains all the signature references. Visualizing the two .fcs files, which fluorophores ended up having messed up unmixing due to the exclusion, and which ones were mostly unaffected?
Problem 3
One of the first packages that I encountered that implemented unmixing in R was Christopher Hall’s flowUnmix package. Since, several additional R packages have also implemented ways to unmix using various methods. While how they process their unmixing controls to get their signature matrices often vary wildly, at the heart the basic unmixing steps are similar, with choice of unmixing method often just needing to swtich out lsfit() for another function that implements the other method. Since you may be interested in implementing one of these methods on your own in the future, lets go on a treasure hunt.
For these 2 packages (flowUnmix, and Autospectral), dive into their GitHub repositories and open their R folder. Check the .R files for any with an “unmixing-something-or-order” style name and examine it. Try reading your way through the code until you encounter the function that actually carries out the unmixing. Look up that functions original package, and check the help file. Take a screenshot, and report back. Feel free to repeat for any unmixing methods of interest (OLS, NNLS, WLS, etc.)
