forked from ropensci/tidyqpcr
-
Notifications
You must be signed in to change notification settings - Fork 0
/
amp_melt_curve_functions.R
188 lines (184 loc) · 6.4 KB
/
amp_melt_curve_functions.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
#' Remove baseline from amplification curves (BETA)
#'
#' Remove baseline from qPCR amplification curves
#' by subtracting median of initial cycles.
#'
#' BETA function version because:
#'
#' - assumes Roche Lightcycler format,
#' we should ideally replace "program_no == 2" by something more generic?
#'
#' - the rule-of thumb "baseline is median of initial 10 cycles"
#' has not been tested robustly
#'
#' @param plateamp data frame with plate amplification data, including variables
#' well, cycle, fluor_raw (raw fluorescence value), and program_no. Assume
#' program 2 for amplification curves from Roche Lightcycler format data.
#' @param maxcycle maximum cycle value to use for baseline, before
#' amplification.
#'
#' @return platemap with additional columns per well:
#' \tabular{ll}{
#' fluor_base \tab baseline /background value \cr
#' fluor_signal \tab normalized fluorescence signal,
#' i.e. fluor_raw - fluor_base
#' }
#'
#' @export
#' @importFrom rlang .data
#'
#' @examples
#' # create simple dataset of raw fluorescence
#' # with two samples over 15 cycles
#' raw_fluor_tibble <- tibble(sample_id = rep(c("S1", "S2"), each = 15),
#' target_id = "T1",
#' well_row = "A",
#' well_col = rep(c(1, 2), each = 15),
#' well = rep(c("A1", "A2"), each = 15),
#' cycle = rep(1:15,2),
#' fluor_raw = c(1:15, 6:20),
#' program_no = 2)
#'
#' # remove base fluorescence from dataset
#' raw_fluor_tibble |>
#' debaseline()
#'
debaseline <- function(plateamp, maxcycle = 10) {
assertthat::assert_that(
assertthat::has_name(
plateamp,
c("well", "program_no", "cycle", "fluor_raw"))
)
assertthat::assert_that(
!assertthat::has_name(plateamp,"fluor_base"),
msg = "plateamp contains a variable called fluor_base, which breaks debaseline"
)
baseline <-
plateamp |>
dplyr::group_by(.data$well) |>
dplyr::filter(.data$program_no == 2,
.data$cycle <= maxcycle) |>
dplyr::summarize(fluor_base = stats::median(.data$fluor_raw))
plateamp |>
dplyr::left_join(baseline, by = "well") |>
dplyr::mutate(fluor_signal = .data$fluor_raw - .data$fluor_base)
}
#' Calculate dy/dx vector from vectors y and x
#'
#' Used in tidyqpcr to calculate dR/dT for a melt curve of fluorescence signal R
#' vs temperature T.
#'
#' @param x input variable, numeric vector, assumed to be temperature
#' @param y output variable, numeric vector of same length as x, assumed
#' to be fluorescence signal.
#' @param method to use for smoothing:
#'
#' "spline" default, uses smoothing spline stats::smooth.spline.
#'
#' "diff" base::diff for lagged difference
#'
#' @param ... other arguments to pass to smoothing method.
#'
#' @return estimated first derivative of y with respect to x, numeric vector of
#' same length as y.
#'
#' @family melt_curve_functions
#'
#' @export
#'
#' @examples
#' # create simple curve
#' x = 1:5
#' y = x^2
#'
#' # calculate gradient of curve
#' #----- use case 1 : using splines
#' calculate_dydx(x, y)
#'
#' # optional arguments are passed to smooth.splines function
#' calculate_dydx(x, y, spar = 0.5)
#'
#' #----- use case 2 : using difference between adjacent points
#' calculate_dydx(x, y, method = "diff")
#'
calculate_dydx <- function(x, y, method = "spline", ...) {
assertthat::assert_that(is.numeric(x))
assertthat::assert_that(is.numeric(y))
assertthat::assert_that(length(x) == length(y))
if (method == "diff") {
return(-c(diff(y) / diff(x), NA))
} else if (method == "spline") {
fit <- stats::smooth.spline(x = x, y = y, ...)
return(-1 * stats::predict(object = fit, x = x, deriv = 1)$y)
}
}
#' Calculate dR/dT of melt curves for of every well in a plate.
#'
#' dR/dT, the derivative of the melt curve (of fluorescence signal R vs
#' temperature T), has a maximum at the melting temperature Tm. A single peak in
#' this suggests a single-length PCR product is present in the well.
#'
#' Note that this function does not group by plate, only by well.
#' The function will give strange results if you pass it data from
#' more than one plate. Avoid this by analysing one plate at a time.
#'
#' @param platemelt data frame describing melt curves, including variables
#' well, temperature, fluor_raw (raw fluorescence value).
#' @param method to use for smoothing:
#'
#' "spline" default, uses smoothing spline stats::smooth.spline.
#'
#' "diff" base::diff for lagged difference
#'
#' @param ... other arguments to pass to smoothing method.
#'
#' @return platemelt with additional column dRdT.
#'
#' @family melt_curve_functions
#'
#' @export
#'
#' @examples
#' # create simple curve
#' # create simple dataset of raw fluorescence with two samples
#' temp_tibble <- tibble(sample_id = rep(c("S1", "S2"), each = 10),
#' target_id = "T1",
#' well_row = "A",
#' well_col = rep(c(1, 2), each = 10),
#' well = rep(c("A1", "A2"), each = 10),
#' temperature = rep(56:65,2),
#' fluor_raw = c(1:10, 6:15))
#'
#' # calculate drdt of all melt curves
#' #----- use case 1 : using splines
#' temp_tibble |>
#' calculate_drdt_plate()
#'
#' # optional arguments are passed to smooth.splines function
#' temp_tibble |>
#' calculate_drdt_plate(spar = 0.5)
#'
#' #----- use case 2 : using difference between adjacent points
#' temp_tibble |>
#' calculate_drdt_plate(method = "diff")
#'
calculate_drdt_plate <- function(platemelt, method = "spline", ...) {
assertthat::assert_that(
assertthat::has_name(
platemelt,
c("well", "temperature", "fluor_raw"))
)
if (assertthat::has_name(platemelt,"plate") ) {
warning("platemelt has a plate column, but calculate_drdt_plate works only for single plates.")
}
platemelt |>
dplyr::arrange(.data$well, .data$temperature) |>
dplyr::group_by(.data$well) |>
dplyr::mutate(dRdT =
calculate_dydx(x = .data$temperature,
y = .data$fluor_raw,
method = method,
...)
) |>
dplyr::ungroup()
}