/
cyto_spillover_spread_compute.R
309 lines (269 loc) · 11.7 KB
/
cyto_spillover_spread_compute.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
## CYTO_SPILLOVER_SPREAD_COMPUTE -----------------------------------------------
#' Compute Spillover Spreading Matrix
#'
#' \code{cyto_spillover_spread_compute} computes the spillover spreading matrix
#' as described by Nguyen et al. (2013) using gated and compensated compensation
#' controls.
#'
#' @param x object of class \code{\link[flowCore:flowSet-class]{flowSet}} or
#' \code{\link[flowWorkspace:GatingSet-class]{GatingSet}} containing
#' compensated, transformed and gated compensation controls.
#' @param parent name of the population to use for downstream calculations when
#' a \code{GatingSet} object is supplied, set to the last gated population by
#' default (e.g. "Single Cells").
#' @param axes_trans object of class
#' \code{\link[flowWorkspace:transformerList]{transformerList}} generated by a
#' \code{cyto_transformer} which contains the transformer definitions to be
#' used internally to transform the fluorescent channels of the samples as
#' required.
#' @param channel_match name of csv file to associate a fluorescent channel with
#' each of the compensation controls. The \code{channel_match} file should
#' contain a "name" column with the names of the compensation controls and a
#' "channel" column to associate a fluorescent channel with each compensation
#' control. Users need not generate this file by hand as it will be created
#' following the channel selection process.
#' @param compensated logical indicating whether the supplied compensation
#' controls have been compensated. Compensation controls must be compensated
#' prior to calculation of spillover spreading matrix. If \code{compensated}
#' is set to FALSE, the compensation controls will be compensated internally
#' using the supplied \code{spillover} matrix.
#' @param spillover name of the output spillover matrix csv file to be used
#' internally to compensate the compensation controls. If no spillover matrix
#' is supplied \code{spillover_spread_compute}, the spillover matrix attached
#' to compensation controls will be applied.
#' @param spillover_spread name of the csv file to which the spillover spreading
#' matrix will be saved, set to \code{"date-Spillover-Spread-Matrix.csv"} by
#' default.
#' @param axes_limits options include \code{"auto"}, \code{"data"} or
#' \code{"machine"} to use optimised, data or machine limits respectively. Set
#' to \code{"machine"} by default to use entire axes ranges. Fine control over
#' axes limits can be obtained by altering the \code{xlim} and \code{ylim}
#' arguments.
#' @param ... additional arguments passed to \code{\link{cyto_plot}}.
#'
#' @references Nguyen R, Perfetto S, Mahnke YD, Chattopadhyay P & Roederer M,
#' (2013). Quantifying spillover spreading for comparing instrument
#' performance and aiding in multicolor panel design. Cytometry A,
#' 83(3):306-15.
#'
#' @importFrom flowWorkspace pData
#' @importFrom flowCore Subset parameters exprs
#' @importFrom stats quantile na.omit
#' @importFrom grDevices graphics.off dev.new
#' @importFrom methods is
#'
#' @author Dillon Hammill, \email{Dillon.Hammill@anu.edu.au}
#'
#' @seealso \code{\link{cyto_spillover_compute}}
#' @seealso \code{\link{cyto_spillover_edit}}
#'
#' @export
cyto_spillover_spread_compute <- function(x,
parent = NULL,
axes_trans = NULL,
channel_match = NULL,
compensated = FALSE,
spillover = NULL,
spillover_spread = NULL,
axes_limits = "machine",
...){
# PREPARE DATA ----------------------------------------------------------
# COPY
cyto_copy <- cyto_copy(x)
# CHANNELS
channels <- cyto_fluor_channels(cyto_copy)
# TRANSFORMATIONS
axes_trans <- cyto_transformer_extract(cyto_copy)
# GATINGSET COMPENSATION
if(is(cyto_copy, "GatingSet")){
spill <- cyto_spillover_extract(cyto_copy)
if(length(spill) == 0){
compensated <- FALSE
}else{
compensated <- TRUE
}
}
# DATA SHOULD BE COMPENSATED & TRANSFORMED
if(compensated == FALSE){
# UNTRANSFORMED DATA
if(.all_na(axes_trans) | is.null(axes_trans)){
cyto_copy <- cyto_compensate(cyto_copy, spillover = spillover)
axes_trans <- cyto_transformer_biex(cyto_copy,
channels = channels,
type = "biex",
plot = FALSE)
cyto_copy <- suppressWarnings(cyto_transform(cyto_copy,
trans = axes_trans,
plot = FALSE))
# TRANSFORMED DATA
}else{
# INVERSE TRANSFORM
suppressWarnings(cyto_transform(cyto_extract(cyto_copy),
trans = axes_trans,
inverse = TRUE,
plot = FALSE))
# COMPENSATE
cyto_compensate(cyto_copy, spillover = spillover)
# TRANSFORM
suppressWarnings(cyto_transform(cyto_extract(cyto_copy),
trans = axes_trans,
plot = FALSE))
}
}else if(compensated == TRUE){
# TRANSFORMATIONS
if(.all_na(axes_trans) | is.null(axes_trans)){
axes_trans <- cyto_transformer_biex(cyto_copy,
channels = channels,
type = "biex",
plot = FALSE)
suppressWarnings(cyto_transform(cyto_copy,
trans = axes_trans,
plot = FALSE))
}
}
## SPILLOVER SPREAD COMPUTATION ----------------------------------------------
# NEGATIVE AND POSITIVE POPULATIONS (TRANSFORMED)
pops <- .cyto_spillover_pops(cyto_copy,
parent = parent,
channel_match = channel_match,
axes_trans = axes_trans,
axes_limits = axes_limits,
...)
neg_pops <- pops[["negative"]]
pos_pops <- pops[["positive"]]
# UPDATE DETAILS IN X
cyto_details(x) <- cyto_details(cyto_copy)
# EXPERIMENT DETAILS
pd <- cyto_details(pos_pops)
# SAMPLE NAMES
nms <- cyto_names(pos_pops)
# INVERSE TRANSFORMATIONS
neg_pops <- cyto_transform(neg_pops,
trans = axes_trans,
inverse = TRUE,
plot = FALSE)
pos_pops <- cyto_transform(pos_pops,
trans = axes_trans,
inverse = TRUE,
plot = FALSE)
# Calculate spillover spread for each compensation control
SSM <- lapply(seq_along(neg_pops), function(z){
# Extract NIL and pop
NIL <- neg_pops[[z]]
pop <- pos_pops[[z]]
# Channel
chan <- pd$channel[pd$name == cyto_names(pop)]
# Calculate 50th and 84th percentile for NIL
NIL_50 <- LAPPLY(channels, function(channel){
quantile(exprs(NIL)[,channel], 0.5)
})
names(NIL_50) <- channels
NIL_84 <- LAPPLY(channels, function(channel){
quantile(exprs(NIL)[,channel], 0.84)
})
names(NIL_84) <- channels
# Combine into a list - NIL_stats
NIL_stats <- list("50%" = NIL_50, "84%" = NIL_84)
# Calculate 50th and 84th percentile for pop
pop_50 <- LAPPLY(channels, function(channel){
quantile(exprs(pop)[,channel], 0.5)
})
names(pop_50) <- channels
pop_84 <- LAPPLY(channels, function(channel){
quantile(exprs(pop)[,channel], 0.84)
})
names(pop_84) <- channels
# Combine into a list - pop_stats
pop_stats <- list("50%" = pop_50, "84%" = pop_84)
# Calculate pop SD and VARIANCE
pop_SD <- pop_stats[[2]] - pop_stats[[1]]
pop_VAR <- pop_SD^2
# Calculate NIL SD and VARIANCE
NIL_SD <- NIL_stats[[2]] - NIL_stats[[1]]
NIL_VAR <- NIL_SD^2
# Difference in VAR of NIL and pop in other detectors
VAR_diff <- pop_VAR - NIL_VAR
# Replace negative numbers with 0 - stain has less spread than reference
if(any(VAR_diff < 0)){
VAR_diff[VAR_diff < 0] <- 0
}
# Difference in SD of NIL and pop in other detectors
SD_diff <- sqrt(VAR_diff)
# Change in fluorescence in primary detector
signal <- pop_stats[[1]][match(chan, channels)] -
NIL_stats[[1]][match(chan, channels)]
# Spillover spread values
return(SD_diff/sqrt(signal))
})
# Make matrix
SSM <- do.call("rbind", SSM)
# Add rownames (channel associated with each control)
rownames(SSM) <- pd$channel[match(cyto_names(pos_pops), pd$name)]
# Replace diagonal with NA
lapply(seq(1, nrow(SSM), 1), function(z) {
SSM[z, match(rownames(SSM)[z], colnames(SSM))] <<- NA
})
# Sort SSM by column names
ind <- na.omit(match(colnames(SSM), rownames(SSM)))
SSM <- SSM[ind, ]
# Turn off pop-up graphics device (gating)
graphics.off()
dev.new()
# Heatmap label size
ylab_width <- max(nchar(rownames(SSM))) * 0.02
xlab_width <- max(nchar(colnames(SSM))) * 0.04
# Plot heatmap
superheat::superheat(SSM,
title = "Spillover Spreading Matrix \n",
title.size = 6,
title.alignment = "center",
left.label.size = ylab_width,
left.label.col = "white",
left.label.text.alignment = "right",
bottom.label.text.angle = 90,
bottom.label.size = xlab_width,
bottom.label.col = "white",
bottom.label.text.alignment = "right",
legend.height = 0.2,
legend.vspace = 0,
row.title = "Fluorochrome \n",
row.title.size = 6,
heat.na.col = "grey",
heat.pal = c("white",
"steelblue",
"steelblue2",
"steelblue3",
"steelblue4",
"purple",
"mediumorchid",
"orchid",
"magenta",
"deeppink2",
"deeppink4",
"red"),
heat.pal.values = c(0,
1,
2,
3,
6,
10,
20,
30,
65,
100,
200,
300))
# Default spillover spread file name
if(is.null(spillover_spread)){
spillover_spread <- paste0(format(Sys.Date(), "%d%m%y"),
"-Spillover-Spread.csv")
}
# Write to csv file
if (!is(spillover_spread, "character")) {
stop("'spillover_spread' should be the name of a csv file.")
} else {
spillover_spread <- file_ext_append(spillover_spread, ".csv")
write.csv(SSM, spillover_spread)
}
return(SSM)
}