In januarharianto/respR: Import, Process, Analyse, and Calculate Rates from Respirometry Data
library(knitr) # load knitr to enable options
library(respR) # load respR
opts_chunk$set(collapse = TRUE,
comment = "#>",
cache = FALSE,
tidy = FALSE,
highlight = TRUE,
fig.width = 10,
fig.height = 5,
fig.align = "center",
R.options = list(
scipen = 999,
digits = 3))
This vignette shows a intermittent-flow analysis from a lengthy experiment, with over 100 replicates and a background rate that changes over the course of the experiment. This type of experiment is one many fish physiologists in particular will be familiar with.
The example data (zeb_intermittent.rd) was kindly provided by Dr. Davide Thambithurai, and is from an experiment on a zebrafish. It is comprised of one replicate of 14 minutes duration (12 minutes recording, 2 minutes flushing), followed by 105 replicates of 11 minutes duration (9 minutes recording, 2 minutes flushing). In addition there are pre- and post-experiment background recordings. See help(zeb_intermittent.rd) for row locations of all stages.
See vignette("intermittent_short") for an example with a smaller dataset, which covers in more detail the basics of analysing intermittent-flow data and is a better starting point if you are new to respR or analysing these types of experiment. Even if not, we would recommend reading both before starting your analysis.
Experiment overview
From this experiment dataset we want to extract several metrics:
-
Pre- and post-experiment background rates. Background rates were recorded at the start and end of the experiment (initial and end sections with shallow slopes). This was a lengthy experiment in high temperatures, and so the microbial background rate increases substantially over the course of the experiment. We will assume this increase is constant (i.e. linear), and each replicate rate will be adjusted with a background rate calculated according its time within the experiment. That is, later replicates will be adjusted by a greater amount.
-
Maximum metabolic rate (MMR). Before this experiment the fish had been exercised to exhaustion, then immediately placed in the respirometer. The first replicate during which it is recovering from the accrued oxygen debt will represent MMR.
-
Routine metabolic rate (RMR). This will be calculated using the later replicates once the animal has recovered from the exercise. This is generally defined as the routine use of oxygen to fuel basal metabolism and minor spontaneous movements for maintaining station or posture.
-
Standard metabolic rate (SMR). Sometimes also known as basal or minimum metabolic rate. This is also calculated using the later replicates once the animal has recovered. This is generally defined as the lowest oxygen use rates observed, which are representative of the requirements of basal metabolism only.
Which of these last two metrics is valid or appropriate to report depends on the organism and the circumstances and aims of the experiments. Their definitions are widely debated, and there is a huge literature on the differences between these physiological metrics. See here for some papers discussing these. Here, we will show how to extract what may be regarded as both RMR and SMR from this dataset.
Analysis methods {#analysis}
Extracting the background rates and MMR are relatively straightforward and will be covered in their own sections.
After this, we will show two different approaches to extract RMR and SMR using respR functions. The first uses the dedicated calc_rate.int() function to run calc_rate on each replicate to extract a rate from the same data region of each, before filtering these to get the RMR. The second uses the dedicated auto_rate.int() function to run auto_rate on each replicate to extract the lowest rate from each of a specified width, before filtering these to get the SMR.
Please note: You should not simply copy and repeat the following code on your own data. They are examples to show how the respR functions are adaptable and flexible and how rates can be extracted, filtered and selected to arrive at a final reportable rate in a variety of ways. They are not recommendations that you should use these specific approaches for particular outputs. There is a massive literature around respirometry, the different physiological metrics, and how to quantify them. Before analysis of your own data you should familiarise yourself with it, with similar experiments in the literature, and the current state of best practices in the field. Many factors such as appropriate regions or widths to extract rates from, method used to extract rates, and final selection criteria differ between experiments for very good reasons, in that every study is different. It is important to select appropriate analysis factors objectively and report them in your methods. See especially Killen et al. (2021) and Prinzing et al. (2021).
Note we use the native pipe operator (|>) introduced in R v4.1 for much of the following. These can be substituted for dpylr pipes (%>%) if you have not yet updated, or alternatively save output objects and enter them into the next function.
Inspecting the data {#inspecting}
inspect() is the general respR function for visualising and inspecting data prior to analysis, and saving it to an object that can be used in subsequent functions. This is a very large dataset (nearly 80,000 rows), so if running the code in this vignette note that operations may take some time.
Before progressing, we'll inspect a portion of the data at the start without saving the result to show the structure.
zeb_intermittent.rd |>
subset_data(from = 1, to = 18000, by = "row") |>
inspect()
The top plot shows the initial background recording up to row 5000, the first replicate of 14 minutes duration to determine MMR, followed by replicates of 11 minutes duration. We can see in these oxygen decreasing, separated by regular flushes where it increases back to ambient levels. The bottom plot, a rolling rate across a window of 10% of the data, is not particularly useful as this window will include multiple replicates.
We can similarly have a quick look at the end of the dataset.
zeb_intermittent.rd |>
subset_data(from = 65000, by = "row") |>
inspect(width = 0.01)
We can see that by the end of the experiment the replicates are much more regular, suggesting rates have stabilised, and we can also see the end background recording. This time we passed a different width to inspect, so the rolling rate plot is across a 1% window. This makes it slightly more useful, and we can see within each replicate rates seem to be consistent at around -0.0025, while the background rate is much lower and closer to zero (dashed line).
Now that we are fairly satisfied the experiment looks as expected we can inspect the entire dataset and this time save the result to an object we will use for the rest of the analysis.
zeb <- inspect(zeb_intermittent.rd)
Here because the dataset is so large the plots are of limited use, but we can at least see that much more oxygen was used in the early replicates, as would be expected because the specimen was exercised prior to the experiment, after which they stabilise and are more regular.
Background rates {#background}
The initial and end sections of the data will be used to calculate background rates and how they change over the course of the experiment. Timepoints of these will typically have been recorded as part of the experiment, so can be used to subset the relevant sections (see ?zeb_intermittent.rd for row locations of all stages).
We subset these regions from the inspect object we just saved above, and save them as separate calc_rate.bg objects.
bg_pre <- zeb |>
subset_data(from = 1, to = 4999, by = "row") |>
calc_rate.bg()
bg_post <- zeb |>
subset_data(from = 75140, to = 79251, by = "row") |>
calc_rate.bg()
The two different background rates (only the post-experiment plot is shown above) have been saved to these two objects. We can see the actual values by printing them:
bg_pre
bg_post
We can see the background rate increases by around 70% over the course of the experiment. We will use these later to apply a dynamic background adjustment to specimen rates based on their time during the experiment.
MMR {#mmr}
Here we calculate maximum metabolic rate (MMR) from the first replicate. This replicate is longer than the others; 14 mins (840 rows) vs. 11 mins (660 rows) for all others.
Subset {#mmrsubset}
Here, we subset the first replicate which starts at row 5000 and is 14 minutes long (840 rows) to a separate inspect object, excluding the two minutes of flush (120 rows) at the end, and plot it to see how the rate changes over the replicate.
# subset rep 1
zeb_rep_1 <- subset_data(zeb,
from = 5000,
to = 5000 + 840 - 120,
by = "row") |>
plot(width = 0.2)
In this replicate, the rolling rate plot shows that calculated at this particular width, the rate declines rapidly from around -0.012 to around -0.008, then declines more slowly over the course of the replicate (note we pass a higher width than the default 0.1 to smooth the rolling rate a little more). This is as we would expect from an animal recovering from exercise. However the very initial stages of replicates are often unstable and excluded from rate analysis (usually known as the "wait" phase). Therefore we will probably want to extract our MMR rate after this initial stage but before it declines too much.
Extract rates
We will show three different ways of extracting what could be defined as an MMR. First, we will use calc_rate() to extract a rate from a single specific region. Secondly, we will use the auto_rate "highest" method to automatically extract the highest rate of a specific width. Lastly, we will use the auto_rate "linear" method with the default inputs to identify linear regions of the data and choose the highest from amongst these. In this last example we will also show how to adjust and convert rates, then filter them using select_rate() to arrive at a final reportable result.
1. calc_rate
This function (see vignette("calc_rate")) simply allows you to extract a rate (or multiple rates) from a specified row, time, or oxygen region. In the later example for extracting rates from the other replicates using calc_rate.int() we use the same three minute wait phase and five minute measure phase in each replicate. We will use the same here to allow for comparison with those rates. This is fairly common practice in respirometry studies and perfectly acceptable, in that it is objective and consistent.
zeb_mmr <- calc_rate(zeb_rep_1,
from = 180, # three minutes 'wait' after start
to = 480, # five minutes later = 'measure' phase
by = "row") |>
summary()
We can see the rate as the final column in summary, and this would be perfectly acceptable to report as the MMR as long as the analysis criteria were also included. See option 3 for how to adjust and convert the rate.
2. auto_rate method = "highest"
auto_rate also allows you to also extract rates of a specific time or row width, but in contrast to calc_rate it calculates every rate of this width across the dataset and then orders them in various ways. One method is to order by highest absolute value to get the highest uptake rate.
Here we run this method and choose a width of two minutes (120 rows), that is the function will calculate every rate of this width and order them from highest to lowest. We use the same three minute wait phase to exclude the initial part of the replicate, and for that we use subset_data() to remove the first 180 rows, then pipe the output to auto_rate (we could also have done this when we first subset and inspected the data above).
zeb_mmr <- zeb_rep_1 |>
subset_data(from = 180, by = "row") |>
auto_rate(method = "highest",
width = 120,
by = "row") |>
summary()
We have already seen that rates are high at the start and then decline across the replicate, so it's not a surprise that the highest rates are from early in the replicate. We can see our MMR under these criteria would be the first row of this summary table. See next section for how these multiple rates can be adjusted, converted, and filtered to get a final output rate.
3. auto_rate method = "linear" {#opt3}
Here we simply run auto_rate and let the default inputs be applied to find the linear regions of the data. Be sure to read help("auto_rate") and vignette("auto_rate") to understand what this method does, appropriate inputs and how to modify the defaults. We again subset the data to remove the initial 3 minutes.
zeb_mmr <- zeb_rep_1 |>
subset_data(from = 180, by = "row") |>
auto_rate() |>
summary()
The function has found seven linear regions. All of these rates are valid results depending on the question being asked, in that they have been found to be linear regions, and so of consistently sustained rates. In this case, since we are interested in the maximum rate, the fifth ranked result has the highest rate value, and comes from early in the replicate, so is probably the result we are most interested in. We can plot it using pos.
plot(zeb_mmr, pos = 5)
Many aspects will feed into rate selection criteria for a specific experiment or study, and it depends on these whether this in an acceptable, reportable result. Maybe at less than 2 minutes this is not a long enough duration for us to be comfortable that this was a sustained MMR. In which case perhaps one of the other results is better reported. Or maybe we could report an average of some of the rates, such as the top three. The important thing is to apply these selection criteria consistently and report them alongside results.
Changing the width {#width}
These are all valid options, as long as they are reported. However, we can also change the analysis parameters. The width input in auto_rate determines the width of the rolling regressions used to find linear regions. See auto_rate() and other vignettes for full details, but briefly in the case of the "linear" method this is the starting point for the analysis. Resulting regressions will vary in width, but in general the minimum width will be increased by increasing the width. The default value is 0.2, representing 20% of the data length. Here, with 10 minutes of data this would represent 2 minutes, which is quite a brief time window. So instead we will increase this to 0.5 (it can also be specified in the units of the "by" input).
See Prinzing et al. (2021) for an excellent discussion of appropriate widths in rolling regressions to determine MMR.
zeb_mmr <- zeb_rep_1 |>
subset_data(from = 180, by = "row") |>
auto_rate(width = 0.5) |>
summary()
Now we have six results varying between 4 and 5 minutes duration, and they are much more consistent in rate value. In this case the fifth result has the highest rate, is early in the replicate, and so is likely a good estimation of MMR. Note how the rolling rate plot is smoother with this higher width.
plot(zeb_mmr, pos = 5)
With non-linear data like this, such as from an animal recovering from exercise, increasing the width over which a rate is determined will necessarily decrease the rate values. The researcher needs to keep this in mind, and find a balance between a width that provides a good estimation of the specimen's physiological state, but also is not affected by other factors such as data noise or anomalies. Most importantly, these decisions and criteria should be reported.
High width values can lead to overfitting and mis-estimation of rates, and low widths underfitting and rate values which are too sensitive to data noise. See here for discussion of this in the context of auto_rate, and how to choose an appropriate width. We recommend reading Prinzing et al. (2021) to understand the nuances of determining MMR and regression widths. See also Killen et al. (2021).
Adjust
The adjust_rate() function allows rates to be background adjusted in a number of ways. See vignette("adjust_rate") for examples. Here we will use it to apply a dynamic linear correction using the start and end background rates we determined earlier. Essentially this calculates the appropriate background adjustment for a rate using the time at which it was determined, assuming the background rate increases linearly from the pre- to post-experiment values.
For MMR, we will adjust the results from option 3 above, but this will work with the other options also. We have our saved auto_rate object, so it's relatively easy to adjust. We just need to enter the two background rate objects we saved earlier and specify the method.
zeb_mmr_adj <- adjust_rate(zeb_mmr,
by = bg_pre, # first background rate
by2 = bg_post, # second background rate
method = "linear")
summary(zeb_mmr_adj)
There are several things to note here. We can see the ultimate adjustment value is small in comparison to the rate of the specimen, as it should be in most experiments. Also, the adjustment values differ slightly; they are lower for rates that come from the start of the replicate. This is what we would expect because we are assuming the background rate is increasing over time. adjust_rate uses the midpoint time of each rate regression to calculate the appropriate adjustment value, so earlier rates will have a lower background adjustment. Lastly, the adjustments are very close in value to the pre-experiment background rate of -0.0000742 we established earlier. Again, this is what we would expect because this is the first replicate and conducted just after the first background recording. If we are assuming there is a linear increase in background rate over time this would mean they should be very close in value.
See later section for how adjustment values differ in other replicates.
Convert
Now that the rates have been adjusted we can convert them to units. This will convert every rate in the object.
The convert_rate function can output rates as absolute, that is of the whole animal or chamber, mass-specific if a mass is entered, or area-specific if an area is entered. The output.unit should be correctly formatted for whichever of these are chosen, that is in the order Oxygen/Time, Oxygen/Time/Mass, or Oxygen/Time/Area.
Note the volume input is volume of water in the respirometer, not the volume of the respirometer. That is, it represents the effective volume. A specimen will displace some of the water, therefore this is the volume of the respirometer minus the volume of the specimen. There are several approaches to calculate the effective volume or specimen volume: geometrically; through displacement in a separate vessel; or calculated from the mass and density. For example, fish are often assumed to have an equal density as water (~1000 kg/m^3), so their mass is measured, converted to volume and subtracted from the respirometer volume. Volume could also be determined directly by pouring out and measuring the water at the end of the experiment, or by weighing the respirometer. See the respfun respirometry utilities package for functions to calculate the effective volume and convert water mass to volume.
Here we use the adjust_rate object saved in the previous step, and we'll convert to a mass-specific rate.
zeb_mmr_adj_conv <- convert_rate(zeb_mmr_adj,
oxy.unit = "mg/L", # oxygen units of the original raw data
time.unit = "secs", # time units of the original raw data
output.unit = "mg/h/g", # desired output unit
volume = 0.12, # effective volume of the respirometer in L
mass = 0.0009) # mass of the specimen in kg
summary(zeb_mmr_adj_conv)
The $summary table contains all rate regression parameters and data locations, adjustments (if applied), units, and more. In fact it contains all relevant information about the analysis and is a great way of saving the results, especially using the export = TRUE input to save it as a data frame. One handy tip: you may want to enter S, t, and P even if they are not required for conversions to the output rate unit because they are saved in the $summary table, and this can help in keeping track of results across different experiments. Note the row and endrow columns will refer to rows of the most recent subset of the data, however the time and endtime will always refer back to the time values of original raw data.
Select
The last analysis step is to select our final reportable MMR from the multiple rates we determined using the linear method. We saw earlier that after changing some analysis parameters one of the MMR rates, the 5th in the summary table, was higher than the others.
summary(zeb_mmr_adj_conv)
In this case, we could fairly confidently report this as our MMR of -3.05 mg/h/g. However, the select_rate() function allows you to filter rates by applying selection criteria. This example will select and extract the highest n rates to a new object, in this case the highest single rate. This is quite a simplistic example, but see later and vignette("select_rate") for more advanced examples including how to apply multiple criteria.
MMR <- select_rate(zeb_mmr_adj_conv,
method = "highest",
n = 1)
If we need to save any of the other regression data for this rate, summary has additional inputs that are useful: pos allows specific rows of summary tables to be extracted, and export allows it to be saved as a separate data frame. Here, we'll use this on the final object to export the coefficients and other data for this rate.
mmr_results <- summary(MMR, export = TRUE)
mmr_results
Complete MMR analysis
This is the complete analysis we have just conducted to get MMR.
mmr_results <-
zeb_rep_1 |> # using the inspected replicate 1 data...
subset_data(from = 180,
by = "row") |> # subset to apply a 'wait' period
auto_rate(width = 0.5) |> # use auto_rate to get most linear regions
adjust_rate(by = bg_pre,
by2 = bg_post,
method = "linear") |> # adjust
convert_rate(oxy.unit = "mg/L",
time.unit = "secs",
output.unit = "mg/h/g",
volume = 0.12,
mass = 0.0009) |> # convert
select_rate(method = "highest",
n = 1) |> # select highest rate
summary(export = TRUE) # final reportable result
mmr_results
This allows you to document all analysis parameters and rate selection criteria so that they can be consistently applied and reported alongside the results. See Reporting section for how analyses might be reported concisely in text.
RMR
calc_rate.int {#crint}
To extract a rate from every replicate and calculate a routine metabolic rate (RMR) from amongst the results we will use calc_rate.int(). To repeat the note from above, this is not a recommendation that this is how you should extract an RMR from your own data; it is simply an example of one approach and how the respR functions are flexible and adaptable.
calc_rate.int was introduced in respR v2.1 and allows you to specify a replicate structure and run calc_rate on each replicate to extract a rate from the same time or row region in each. See vignette("calc_rate.int") for details of the functionality.
The difference between this and the approach below using auto_rate.int() is that you specify a region over which to calculate a rate, and this will be the same in each replicate. This takes away some of the objectivity in determining rates that auto_rate allows for, but as long as selection criteria are reported and applied consistently is a perfectly valid way of processing intermittent-flow data.
We've shown above how to get background rates and MMR, so we will use calc_rate.int to extract a rate from the same region of each replicate and filter these rates to report an RMR. In these data, after the first replicate of 14 minutes duration for MMR there is a regular structure of 105 replicates of 11 minutes duration including the flush.
Subset data
In the examples here we had to specify the start and end location of every replicate in calc_rate.int. However, if replicates cycle at a regular interval we can simply enter this as starts. This tells the function that starting from row 1 replicates cycle at this regular interval, and it can be entered as a row interval or a duration in time, as specified using the by input. We only need to subset the data so that the first replicate starts at row 1 before passing it to calc_rate.int.
We start by subsetting the data so that the 105 regularly-spaced replicates start at row 1, that is we exclude the starting background recording, the first replicate for MMR, and the end background recording. We then inspect it.
zeb_insp <- zeb_intermittent.rd |>
subset_data(from = 5840,
to = 75139,
by = "row") |>
inspect()
Extract rates
Now we run calc_rate.int and specify a 660 row interval between replicates using starts. We'll calculate a rate across the same five minute time period in each replicate, from three minutes to eight minutes. To do this we specify a 180 row wait phase and 300 row measure phase.
zeb_cr.int <- calc_rate.int(zeb_insp,
starts = 660,
wait = 180,
measure = 300,
by = "row")
By default the function plots the first 20 replicates (pos can be used to choose others up to a maximum of 20). Note how the lower time axis has the actual raw data time values, but the top row axis refers to the rows of each replicate subset. In these plots the rate data is the yellow points, the red shaded region is the wait phase, and the green shaded region the measure phase.
We can use summary() to view the results.
summary(zeb_cr.int)
As we saw from our inspection of the data earlier, rates are higher at the start of the dataset because the specimen was exercised before being placed in the respirometer. We will filter these out shortly. For now we can see we have a rate from every replicate. Our entered time range has been applied within each replicate, and the replicate each rate originates from can be seen in the $rep column.
Adjust {#crintadj}
calc_rate.int objects also work with adjust_rate, including with dynamic adjustment methods. Here, we'll adjust using the pre- and post-experiment background recordings assuming the background rate increases linearly from the initial to ending background rates.
zeb_cr.int_adj <- adjust_rate(zeb_cr.int,
by = bg_pre,
by2 = bg_post,
method = "linear")
If the adjustments were correctly applied we would expect the adjustment value of the early replicates to be close to the value of the initial background rate (-0.000074), and those of the later ones close to that of the ending background rate (-0.000120). Let's look at the summary to check, and use pos to pick out one from the start, middle and end.
summary(zeb_cr.int_adj, pos = c(1, 50, 105))
As expected the $adjustment column values increase over the experiment between these values, so we can proceed to conversion.
Convert
Now we will convert the rates, then plot the values to decide how to select a final RMR rate.
zeb_cr.int_conv <- convert_rate(zeb_cr.int_adj,
oxy.unit = "mg/L", # oxygen units of the original raw data
time.unit = "secs", # time units of the original raw data
output.unit = "mg/h/g", # desired output unit
volume = 0.12, # effective volume of the respirometer in L
mass = 0.0009) # mass of the specimen in kg
The converted rates can be seen in the summary table.
summary(zeb_cr.int_conv)
We can see the converted rates as the last column. The $summary table contains the all rate regression parameters and data locations, adjustments (if applied), units, and more. It can be exported as a separate data frame by using export = TRUE in summary(). This is a great way of exporting all the relevant data for your final results.
Select {#crintselect}
It will depend on the aims and circumstances of an experiment what rate results to extract, summarise, and report. For instance we might report RMR as a mean of a representative selection of rates after the specimens uptake rate has stabilised, or a mean of a number or percentile of lowest rates found. The important thing is to apply these selection criteria consistently and report them alongside results.
As of respR v2.1 convert_rate has new plotting functionality to help visualise and explore rate results. For intermittent-flow results the most useful of these is type = "rate".
plot(zeb_cr.int_conv, type = "rate")
This shows how the calculated rate values vary across the dataset. We can see from this plot they are high at the start, stabilise after around timepoint 20000, then are a little unstable again if not slightly elevated after timepoint 60000. Between these however they appear to be very stable.
The select_rate() function allows convert_rate results to be filtered or reordered in various ways. See help("select_rate") and vignette("select_rate") for more details. For this analysis, lets only use rates from replicates between timepoints 20000 and 60000. Then, simply to show an example of how rates can be further selected, let's only select those with r-squared of 0.97 and above, select the lowest 20, and finally take the mean of them.
RMR <- zeb_cr.int_conv |>
select_rate(method = "time", n = c(20000, 60000)) |>
select_rate(method = "rsq", n = c(0.97, 1)) |>
select_rate(method = "lowest", n = 20) |>
summary() |>
mean()
That's it; this is our final routine metabolic rate (at least, as we have chosen to define it).
SMR
auto_rate.int {#arint}
To extract a rate from every replicate and determine a standard metabolic rate (SMR) from amongst the results we will use auto_rate.int(). To repeat the note from above, this is not a recommendation that this is how you should extract an SMR from your own data; it is simply an example of one approach and how the respR functions are flexible and adaptable.
Introduced in respR v2.1, auto_rate.int() is a dedicated function for intermittent-flow respirometry data, in which you specify a replicate structure and run auto_rate() on each replicate. This function extracts rates in a number of ways including the most linear or most consistently maintained rates, plus highest and lowest rates of a specified width. See vignette("auto_rate") for more details, and vignette("auto_rate.int") for an overview of this function.
The difference between this and the approach above using calc_rate.int() is that the specified measure phase is not the same as the region over which a rate is calculated. Instead a measure phase is specified and the auto_rate analysis is done within this region using the method and width inputs. This can result in rates of different widths coming from within this region depending on the inputs.
Subset data
We again subset the data so that the 105 regularly-spaced replicates start at row 1.
zeb_insp <- zeb_intermittent.rd |>
subset_data(from = 5840,
to = 75139,
by = "row") |>
inspect()
Extract rates {#arintextract}
To get SMR we necessarily need to define it, in particular the duration over which rates are calculated within each replicate. These should be chosen appropriately and with reference to best practices and previous studies. Here, let's say we want to extract the lowest rate across 2 minutes within each replicate, and we will define our SMR for this specimen as the mean of the lowest 10th percentile of all 105 extracted rates.
We run auto_rate.int on the inspected data and use method = "lowest" and width = 120. This means the function calculates every rate of this width within each replicate, orders them according to the input method, and by default (n = 1) returns the top ranked result, in this case the lowest, for each replicate.
We again specify a 660 row interval between replicates using starts. Just to be different from the above example, this time we will specify a 2 minute (120 row) wait phase, and six minute (360 row) measure phase.
zeb_ar.int <- auto_rate.int(zeb_insp,
starts = 660,
wait = 120,
measure = 360,
method = "lowest",
width = 120,
by = "row") |>
summary()
Note how the plots show the rates coming from different regions within each replicate measure phase. The summary contains all the results, that is the lowest 2-minute rate from each replicate. Before selecting from amongst these to get our final SMR we first need to adjust and convert them. This is exactly the same process we followed in the example above.
Adjust
zeb_ar.int_adj <- adjust_rate(zeb_ar.int,
by = bg_pre,
by2 = bg_post,
method = "linear")
Convert
zeb_ar.int_conv <- convert_rate(zeb_ar.int_adj,
oxy.unit = "mg/L", # oxygen units of the original raw data
time.unit = "secs", # time units of the original raw data
output.unit = "mg/h/g", # desired output unit
volume = 0.12, # effective volume of the respirometer in L
mass = 0.0009) # mass of the specimen in kg
summary(zeb_ar.int_conv)
Quality control
We can plot the convert_rate object for a quick look at how rates vary across the data.
plot(zeb_ar.int_conv, type = "rate")
These rates seem to be a bit lower than those we extracted for RMR above, as we would expect for an SMR.
This is an opportunity to do some quality control before we select our final rates to determine SMR. Looking at the plot above there is one result that seems to be uncharacteristically low compared to the others with a value just over 0.5. Looking at the summary table we can see this is from replicate 48. We can plot a selection of rates from this region using pos and the y-axis range will adapt to provide a more detailed view of the range of rates.
plot(zeb_ar.int_conv, type = "rate", pos = 30:80)
Replicate 48 does seem to be anomalously low. If we go back to the auto_rate.int results object we can use pos to plot this replicate for a closer look. We'll do this in two different ways using type. Firstly the rate in the context of the whole replicate, secondly the auto_rate output plot for this replicate which has more diagnostic plots, but only shows the measure phase.
plot(zeb_ar.int, type = "rep", pos = 48)
plot(zeb_ar.int, type = "ar", pos = 48)
Looking at these plots this result does seem to be anomalous, and in this replicate the rate extracted is not really representative. It also has a lower r-squared which otherwise varies from 0.72 to 0.98. There seems to be an anomaly in the oxygen trace, maybe a sensor glitch or slight temperature change. Issues like this causing glitches in the trace are very common.
Therefore, it would be perfectly acceptable and justifiable to manually exclude this result (there are also 1 or 2 others we might want to investigate). To apply such selection criteria objectively and systematically would require examining multiple experiments. This is why definitions of metrics like SMR and analysis parameters should usually be decided upon after thoroughly exploring all your results. After this you should be able to arrive at a set of objective and reasonable selection criteria.
Select
To get our final SMR we just need to use select_rate again with our chosen criteria, and pipe the result to mean. We decided above our definition of SMR would be the mean of the lowest 10th percentile of the rates extracted from each replicate. We also want to exclude the result from replicate 48.
Excluding this replicate could be achieved in a number of ways in select_rate. We could use method = "rsq" to select only results with r-squared above a certain value which excludes this one. We could use the time_omit method to exclude this rate based on its timepoint (found in the summary table), or use the rep method simply to select all replicates except this one. There is a specific rep_omit method however to remove particular replicates.
SMR <- zeb_ar.int_conv |>
select_rate(method = "rep_omit", n = 48) |>
select_rate(method = "lowest_percentile", n = 0.1) |>
summary() |>
mean()
We can see the summary of the final 11 selected rates, and then the mean of these. This is our final SMR.
Summary {#ratessumm}
Now we have our final results for each metric.
MMR: -3.05 mg/h/g
RMR: -0.94 mg/h/g
SMR: -0.70 mg/h/g
Typically these rates would now be reported as positive values.
Reporting {#report}
The above is a lengthy description of an analysis, so this vignette may not give a good impression of just how succinct a respR analysis may be. This code block shows the entire analysis for SMR.
Complete analysis
# Inspect data ------------------------------------------------------------
zeb <- inspect(zeb_intermittent.rd)
# Background --------------------------------------------------------------
bg_pre <- subset_data(zeb, from = 0, to = 4999, by = "time") |>
calc_rate.bg()
bg_post <- subset_data(zeb, from = 75140, to = 79251, by = "time") |>
calc_rate.bg()
# SMR Analysis ------------------------------------------------------------
# Subset ------------------------------------------------------------------
zeb |>
subset_data(from = 5840,
to = 75139,
by = "row") |>
# Extract rates -----------------------------------------------------------
auto_rate.int(starts = 660,
wait = 120,
measure = 360,
method = "lowest",
width = 120,
by = "row") |>
# Adjust ------------------------------------------------------------------
adjust_rate(by = bg_pre,
by2 = bg_post,
method = "linear") |>
# Convert -----------------------------------------------------------------
convert_rate(oxy.unit = "mg/L",
time.unit = "secs",
output.unit = "mg/h/g",
volume = 0.12,
mass = 0.0009) |>
# Select -----------------------------------------------------------------
select_rate(method = "rep_omit", n = 48) |>
select_rate(method = "lowest_percentile", n = 0.1) |>
summary() |>
mean()
How to report analysis methods
The ultimate aim of respR was to enable clear reporting and reproducible analyses. This is an example of how we might report the above analysis in a manuscript.
"Data was analysed using the R package respR (Harianto et al. 2019). The auto_rate.int function was used to extract the lowest rate across a two minute window from each replicate. Rates were adjusted for background using pre- and post-experiment blank recordings with no specimen, assuming a linear increase in background rate over time. Obviously anomalous rates were excluded, and SMR was defined as the mean of the lowest 10th percentile of remaining adjusted rates."
Alternative approaches
Previous versions of this vignette demonstrated use of different methods, such as looping and subsetting of data, to extract individual replicates for processing through calc_rate or auto_rate. These approaches still work and are archived here in case they are of use in some cases, for example if you want to extract rates from the same oxygen range in each replicate. However for most analyses use of calc_rate.int or auto_rate.int is the recommended approach.
januarharianto/respR documentation built on April 20, 2024, 4:34 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(knitr) # load knitr to enable options library(respR) # load respR opts_chunk$set(collapse = TRUE, comment = "#>", cache = FALSE, tidy = FALSE, highlight = TRUE, fig.width = 10, fig.height = 5, fig.align = "center", R.options = list( scipen = 999, digits = 3))
This vignette shows a intermittent-flow analysis from a lengthy experiment, with over 100 replicates and a background rate that changes over the course of the experiment. This type of experiment is one many fish physiologists in particular will be familiar with.
The example data (zeb_intermittent.rd) was kindly provided by Dr. Davide Thambithurai, and is from an experiment on a zebrafish. It is comprised of one replicate of 14 minutes duration (12 minutes recording, 2 minutes flushing), followed by 105 replicates of 11 minutes duration (9 minutes recording, 2 minutes flushing). In addition there are pre- and post-experiment background recordings. See help(zeb_intermittent.rd) for row locations of all stages.
See vignette("intermittent_short") for an example with a smaller dataset, which covers in more detail the basics of analysing intermittent-flow data and is a better starting point if you are new to respR or analysing these types of experiment. Even if not, we would recommend reading both before starting your analysis.
Experiment overview
From this experiment dataset we want to extract several metrics:
-
Pre- and post-experiment background rates. Background rates were recorded at the start and end of the experiment (initial and end sections with shallow slopes). This was a lengthy experiment in high temperatures, and so the microbial background rate increases substantially over the course of the experiment. We will assume this increase is constant (i.e. linear), and each replicate rate will be adjusted with a background rate calculated according its time within the experiment. That is, later replicates will be adjusted by a greater amount.
-
Maximum metabolic rate (MMR). Before this experiment the fish had been exercised to exhaustion, then immediately placed in the respirometer. The first replicate during which it is recovering from the accrued oxygen debt will represent MMR.
-
Routine metabolic rate (RMR). This will be calculated using the later replicates once the animal has recovered from the exercise. This is generally defined as the routine use of oxygen to fuel basal metabolism and minor spontaneous movements for maintaining station or posture.
-
Standard metabolic rate (SMR). Sometimes also known as basal or minimum metabolic rate. This is also calculated using the later replicates once the animal has recovered. This is generally defined as the lowest oxygen use rates observed, which are representative of the requirements of basal metabolism only.
Which of these last two metrics is valid or appropriate to report depends on the organism and the circumstances and aims of the experiments. Their definitions are widely debated, and there is a huge literature on the differences between these physiological metrics. See here for some papers discussing these. Here, we will show how to extract what may be regarded as both RMR and SMR from this dataset.
Analysis methods {#analysis}
Extracting the background rates and MMR are relatively straightforward and will be covered in their own sections.
After this, we will show two different approaches to extract RMR and SMR using respR functions. The first uses the dedicated calc_rate.int() function to run calc_rate on each replicate to extract a rate from the same data region of each, before filtering these to get the RMR. The second uses the dedicated auto_rate.int() function to run auto_rate on each replicate to extract the lowest rate from each of a specified width, before filtering these to get the SMR.
Please note: You should not simply copy and repeat the following code on your own data. They are examples to show how the respR functions are adaptable and flexible and how rates can be extracted, filtered and selected to arrive at a final reportable rate in a variety of ways. They are not recommendations that you should use these specific approaches for particular outputs. There is a massive literature around respirometry, the different physiological metrics, and how to quantify them. Before analysis of your own data you should familiarise yourself with it, with similar experiments in the literature, and the current state of best practices in the field. Many factors such as appropriate regions or widths to extract rates from, method used to extract rates, and final selection criteria differ between experiments for very good reasons, in that every study is different. It is important to select appropriate analysis factors objectively and report them in your methods. See especially Killen et al. (2021) and Prinzing et al. (2021).
Note we use the native pipe operator (|>) introduced in R v4.1 for much of the following. These can be substituted for dpylr pipes (%>%) if you have not yet updated, or alternatively save output objects and enter them into the next function.
Inspecting the data {#inspecting}
inspect() is the general respR function for visualising and inspecting data prior to analysis, and saving it to an object that can be used in subsequent functions. This is a very large dataset (nearly 80,000 rows), so if running the code in this vignette note that operations may take some time.
Before progressing, we'll inspect a portion of the data at the start without saving the result to show the structure.
zeb_intermittent.rd |> subset_data(from = 1, to = 18000, by = "row") |> inspect()
The top plot shows the initial background recording up to row 5000, the first replicate of 14 minutes duration to determine MMR, followed by replicates of 11 minutes duration. We can see in these oxygen decreasing, separated by regular flushes where it increases back to ambient levels. The bottom plot, a rolling rate across a window of 10% of the data, is not particularly useful as this window will include multiple replicates.
We can similarly have a quick look at the end of the dataset.
zeb_intermittent.rd |> subset_data(from = 65000, by = "row") |> inspect(width = 0.01)
We can see that by the end of the experiment the replicates are much more regular, suggesting rates have stabilised, and we can also see the end background recording. This time we passed a different width to inspect, so the rolling rate plot is across a 1% window. This makes it slightly more useful, and we can see within each replicate rates seem to be consistent at around -0.0025, while the background rate is much lower and closer to zero (dashed line).
Now that we are fairly satisfied the experiment looks as expected we can inspect the entire dataset and this time save the result to an object we will use for the rest of the analysis.
zeb <- inspect(zeb_intermittent.rd)
Here because the dataset is so large the plots are of limited use, but we can at least see that much more oxygen was used in the early replicates, as would be expected because the specimen was exercised prior to the experiment, after which they stabilise and are more regular.
Background rates {#background}
The initial and end sections of the data will be used to calculate background rates and how they change over the course of the experiment. Timepoints of these will typically have been recorded as part of the experiment, so can be used to subset the relevant sections (see ?zeb_intermittent.rd for row locations of all stages).
We subset these regions from the inspect object we just saved above, and save them as separate calc_rate.bg objects.
bg_pre <- zeb |> subset_data(from = 1, to = 4999, by = "row") |> calc_rate.bg() bg_post <- zeb |> subset_data(from = 75140, to = 79251, by = "row") |> calc_rate.bg()
The two different background rates (only the post-experiment plot is shown above) have been saved to these two objects. We can see the actual values by printing them:
bg_pre
bg_post
We can see the background rate increases by around 70% over the course of the experiment. We will use these later to apply a dynamic background adjustment to specimen rates based on their time during the experiment.
MMR {#mmr}
Here we calculate maximum metabolic rate (MMR) from the first replicate. This replicate is longer than the others; 14 mins (840 rows) vs. 11 mins (660 rows) for all others.
Subset {#mmrsubset}
Here, we subset the first replicate which starts at row 5000 and is 14 minutes long (840 rows) to a separate inspect object, excluding the two minutes of flush (120 rows) at the end, and plot it to see how the rate changes over the replicate.
# subset rep 1 zeb_rep_1 <- subset_data(zeb, from = 5000, to = 5000 + 840 - 120, by = "row") |> plot(width = 0.2)
In this replicate, the rolling rate plot shows that calculated at this particular width, the rate declines rapidly from around -0.012 to around -0.008, then declines more slowly over the course of the replicate (note we pass a higher width than the default 0.1 to smooth the rolling rate a little more). This is as we would expect from an animal recovering from exercise. However the very initial stages of replicates are often unstable and excluded from rate analysis (usually known as the "wait" phase). Therefore we will probably want to extract our MMR rate after this initial stage but before it declines too much.
Extract rates
We will show three different ways of extracting what could be defined as an MMR. First, we will use calc_rate() to extract a rate from a single specific region. Secondly, we will use the auto_rate "highest" method to automatically extract the highest rate of a specific width. Lastly, we will use the auto_rate "linear" method with the default inputs to identify linear regions of the data and choose the highest from amongst these. In this last example we will also show how to adjust and convert rates, then filter them using select_rate() to arrive at a final reportable result.
1. calc_rate
This function (see vignette("calc_rate")) simply allows you to extract a rate (or multiple rates) from a specified row, time, or oxygen region. In the later example for extracting rates from the other replicates using calc_rate.int() we use the same three minute wait phase and five minute measure phase in each replicate. We will use the same here to allow for comparison with those rates. This is fairly common practice in respirometry studies and perfectly acceptable, in that it is objective and consistent.
zeb_mmr <- calc_rate(zeb_rep_1, from = 180, # three minutes 'wait' after start to = 480, # five minutes later = 'measure' phase by = "row") |> summary()
We can see the rate as the final column in summary, and this would be perfectly acceptable to report as the MMR as long as the analysis criteria were also included. See option 3 for how to adjust and convert the rate.
2. auto_rate method = "highest"
auto_rate also allows you to also extract rates of a specific time or row width, but in contrast to calc_rate it calculates every rate of this width across the dataset and then orders them in various ways. One method is to order by highest absolute value to get the highest uptake rate.
Here we run this method and choose a width of two minutes (120 rows), that is the function will calculate every rate of this width and order them from highest to lowest. We use the same three minute wait phase to exclude the initial part of the replicate, and for that we use subset_data() to remove the first 180 rows, then pipe the output to auto_rate (we could also have done this when we first subset and inspected the data above).
zeb_mmr <- zeb_rep_1 |> subset_data(from = 180, by = "row") |> auto_rate(method = "highest", width = 120, by = "row") |> summary()
We have already seen that rates are high at the start and then decline across the replicate, so it's not a surprise that the highest rates are from early in the replicate. We can see our MMR under these criteria would be the first row of this summary table. See next section for how these multiple rates can be adjusted, converted, and filtered to get a final output rate.
3. auto_rate method = "linear" {#opt3}
Here we simply run auto_rate and let the default inputs be applied to find the linear regions of the data. Be sure to read help("auto_rate") and vignette("auto_rate") to understand what this method does, appropriate inputs and how to modify the defaults. We again subset the data to remove the initial 3 minutes.
zeb_mmr <- zeb_rep_1 |> subset_data(from = 180, by = "row") |> auto_rate() |> summary()
The function has found seven linear regions. All of these rates are valid results depending on the question being asked, in that they have been found to be linear regions, and so of consistently sustained rates. In this case, since we are interested in the maximum rate, the fifth ranked result has the highest rate value, and comes from early in the replicate, so is probably the result we are most interested in. We can plot it using pos.
plot(zeb_mmr, pos = 5)
Many aspects will feed into rate selection criteria for a specific experiment or study, and it depends on these whether this in an acceptable, reportable result. Maybe at less than 2 minutes this is not a long enough duration for us to be comfortable that this was a sustained MMR. In which case perhaps one of the other results is better reported. Or maybe we could report an average of some of the rates, such as the top three. The important thing is to apply these selection criteria consistently and report them alongside results.
Changing the width {#width}
These are all valid options, as long as they are reported. However, we can also change the analysis parameters. The width input in auto_rate determines the width of the rolling regressions used to find linear regions. See auto_rate() and other vignettes for full details, but briefly in the case of the "linear" method this is the starting point for the analysis. Resulting regressions will vary in width, but in general the minimum width will be increased by increasing the width. The default value is 0.2, representing 20% of the data length. Here, with 10 minutes of data this would represent 2 minutes, which is quite a brief time window. So instead we will increase this to 0.5 (it can also be specified in the units of the "by" input).
See Prinzing et al. (2021) for an excellent discussion of appropriate widths in rolling regressions to determine MMR.
zeb_mmr <- zeb_rep_1 |> subset_data(from = 180, by = "row") |> auto_rate(width = 0.5) |> summary()
Now we have six results varying between 4 and 5 minutes duration, and they are much more consistent in rate value. In this case the fifth result has the highest rate, is early in the replicate, and so is likely a good estimation of MMR. Note how the rolling rate plot is smoother with this higher width.
plot(zeb_mmr, pos = 5)
With non-linear data like this, such as from an animal recovering from exercise, increasing the width over which a rate is determined will necessarily decrease the rate values. The researcher needs to keep this in mind, and find a balance between a width that provides a good estimation of the specimen's physiological state, but also is not affected by other factors such as data noise or anomalies. Most importantly, these decisions and criteria should be reported.
High width values can lead to overfitting and mis-estimation of rates, and low widths underfitting and rate values which are too sensitive to data noise. See here for discussion of this in the context of auto_rate, and how to choose an appropriate width. We recommend reading Prinzing et al. (2021) to understand the nuances of determining MMR and regression widths. See also Killen et al. (2021).
Adjust
The adjust_rate() function allows rates to be background adjusted in a number of ways. See vignette("adjust_rate") for examples. Here we will use it to apply a dynamic linear correction using the start and end background rates we determined earlier. Essentially this calculates the appropriate background adjustment for a rate using the time at which it was determined, assuming the background rate increases linearly from the pre- to post-experiment values.
For MMR, we will adjust the results from option 3 above, but this will work with the other options also. We have our saved auto_rate object, so it's relatively easy to adjust. We just need to enter the two background rate objects we saved earlier and specify the method.
zeb_mmr_adj <- adjust_rate(zeb_mmr, by = bg_pre, # first background rate by2 = bg_post, # second background rate method = "linear")
summary(zeb_mmr_adj)
There are several things to note here. We can see the ultimate adjustment value is small in comparison to the rate of the specimen, as it should be in most experiments. Also, the adjustment values differ slightly; they are lower for rates that come from the start of the replicate. This is what we would expect because we are assuming the background rate is increasing over time. adjust_rate uses the midpoint time of each rate regression to calculate the appropriate adjustment value, so earlier rates will have a lower background adjustment. Lastly, the adjustments are very close in value to the pre-experiment background rate of -0.0000742 we established earlier. Again, this is what we would expect because this is the first replicate and conducted just after the first background recording. If we are assuming there is a linear increase in background rate over time this would mean they should be very close in value.
See later section for how adjustment values differ in other replicates.
Convert
Now that the rates have been adjusted we can convert them to units. This will convert every rate in the object.
The convert_rate function can output rates as absolute, that is of the whole animal or chamber, mass-specific if a mass is entered, or area-specific if an area is entered. The output.unit should be correctly formatted for whichever of these are chosen, that is in the order Oxygen/Time, Oxygen/Time/Mass, or Oxygen/Time/Area.
Note the volume input is volume of water in the respirometer, not the volume of the respirometer. That is, it represents the effective volume. A specimen will displace some of the water, therefore this is the volume of the respirometer minus the volume of the specimen. There are several approaches to calculate the effective volume or specimen volume: geometrically; through displacement in a separate vessel; or calculated from the mass and density. For example, fish are often assumed to have an equal density as water (~1000 kg/m^3), so their mass is measured, converted to volume and subtracted from the respirometer volume. Volume could also be determined directly by pouring out and measuring the water at the end of the experiment, or by weighing the respirometer. See the respfun respirometry utilities package for functions to calculate the effective volume and convert water mass to volume.
Here we use the adjust_rate object saved in the previous step, and we'll convert to a mass-specific rate.
zeb_mmr_adj_conv <- convert_rate(zeb_mmr_adj, oxy.unit = "mg/L", # oxygen units of the original raw data time.unit = "secs", # time units of the original raw data output.unit = "mg/h/g", # desired output unit volume = 0.12, # effective volume of the respirometer in L mass = 0.0009) # mass of the specimen in kg summary(zeb_mmr_adj_conv)
The $summary table contains all rate regression parameters and data locations, adjustments (if applied), units, and more. In fact it contains all relevant information about the analysis and is a great way of saving the results, especially using the export = TRUE input to save it as a data frame. One handy tip: you may want to enter S, t, and P even if they are not required for conversions to the output rate unit because they are saved in the $summary table, and this can help in keeping track of results across different experiments. Note the row and endrow columns will refer to rows of the most recent subset of the data, however the time and endtime will always refer back to the time values of original raw data.
Select
The last analysis step is to select our final reportable MMR from the multiple rates we determined using the linear method. We saw earlier that after changing some analysis parameters one of the MMR rates, the 5th in the summary table, was higher than the others.
summary(zeb_mmr_adj_conv)
In this case, we could fairly confidently report this as our MMR of -3.05 mg/h/g. However, the select_rate() function allows you to filter rates by applying selection criteria. This example will select and extract the highest n rates to a new object, in this case the highest single rate. This is quite a simplistic example, but see later and vignette("select_rate") for more advanced examples including how to apply multiple criteria.
MMR <- select_rate(zeb_mmr_adj_conv, method = "highest", n = 1)
If we need to save any of the other regression data for this rate, summary has additional inputs that are useful: pos allows specific rows of summary tables to be extracted, and export allows it to be saved as a separate data frame. Here, we'll use this on the final object to export the coefficients and other data for this rate.
mmr_results <- summary(MMR, export = TRUE)
mmr_results
Complete MMR analysis
This is the complete analysis we have just conducted to get MMR.
mmr_results <- zeb_rep_1 |> # using the inspected replicate 1 data... subset_data(from = 180, by = "row") |> # subset to apply a 'wait' period auto_rate(width = 0.5) |> # use auto_rate to get most linear regions adjust_rate(by = bg_pre, by2 = bg_post, method = "linear") |> # adjust convert_rate(oxy.unit = "mg/L", time.unit = "secs", output.unit = "mg/h/g", volume = 0.12, mass = 0.0009) |> # convert select_rate(method = "highest", n = 1) |> # select highest rate summary(export = TRUE) # final reportable result
mmr_results
This allows you to document all analysis parameters and rate selection criteria so that they can be consistently applied and reported alongside the results. See Reporting section for how analyses might be reported concisely in text.
RMR
calc_rate.int {#crint}
To extract a rate from every replicate and calculate a routine metabolic rate (RMR) from amongst the results we will use calc_rate.int(). To repeat the note from above, this is not a recommendation that this is how you should extract an RMR from your own data; it is simply an example of one approach and how the respR functions are flexible and adaptable.
calc_rate.int was introduced in respR v2.1 and allows you to specify a replicate structure and run calc_rate on each replicate to extract a rate from the same time or row region in each. See vignette("calc_rate.int") for details of the functionality.
The difference between this and the approach below using auto_rate.int() is that you specify a region over which to calculate a rate, and this will be the same in each replicate. This takes away some of the objectivity in determining rates that auto_rate allows for, but as long as selection criteria are reported and applied consistently is a perfectly valid way of processing intermittent-flow data.
We've shown above how to get background rates and MMR, so we will use calc_rate.int to extract a rate from the same region of each replicate and filter these rates to report an RMR. In these data, after the first replicate of 14 minutes duration for MMR there is a regular structure of 105 replicates of 11 minutes duration including the flush.
Subset data
In the examples here we had to specify the start and end location of every replicate in calc_rate.int. However, if replicates cycle at a regular interval we can simply enter this as starts. This tells the function that starting from row 1 replicates cycle at this regular interval, and it can be entered as a row interval or a duration in time, as specified using the by input. We only need to subset the data so that the first replicate starts at row 1 before passing it to calc_rate.int.
We start by subsetting the data so that the 105 regularly-spaced replicates start at row 1, that is we exclude the starting background recording, the first replicate for MMR, and the end background recording. We then inspect it.
zeb_insp <- zeb_intermittent.rd |> subset_data(from = 5840, to = 75139, by = "row") |> inspect()
Extract rates
Now we run calc_rate.int and specify a 660 row interval between replicates using starts. We'll calculate a rate across the same five minute time period in each replicate, from three minutes to eight minutes. To do this we specify a 180 row wait phase and 300 row measure phase.
zeb_cr.int <- calc_rate.int(zeb_insp, starts = 660, wait = 180, measure = 300, by = "row")
By default the function plots the first 20 replicates (pos can be used to choose others up to a maximum of 20). Note how the lower time axis has the actual raw data time values, but the top row axis refers to the rows of each replicate subset. In these plots the rate data is the yellow points, the red shaded region is the wait phase, and the green shaded region the measure phase.
We can use summary() to view the results.
summary(zeb_cr.int)
As we saw from our inspection of the data earlier, rates are higher at the start of the dataset because the specimen was exercised before being placed in the respirometer. We will filter these out shortly. For now we can see we have a rate from every replicate. Our entered time range has been applied within each replicate, and the replicate each rate originates from can be seen in the $rep column.
Adjust {#crintadj}
calc_rate.int objects also work with adjust_rate, including with dynamic adjustment methods. Here, we'll adjust using the pre- and post-experiment background recordings assuming the background rate increases linearly from the initial to ending background rates.
zeb_cr.int_adj <- adjust_rate(zeb_cr.int, by = bg_pre, by2 = bg_post, method = "linear")
If the adjustments were correctly applied we would expect the adjustment value of the early replicates to be close to the value of the initial background rate (-0.000074), and those of the later ones close to that of the ending background rate (-0.000120). Let's look at the summary to check, and use pos to pick out one from the start, middle and end.
summary(zeb_cr.int_adj, pos = c(1, 50, 105))
As expected the $adjustment column values increase over the experiment between these values, so we can proceed to conversion.
Convert
Now we will convert the rates, then plot the values to decide how to select a final RMR rate.
zeb_cr.int_conv <- convert_rate(zeb_cr.int_adj, oxy.unit = "mg/L", # oxygen units of the original raw data time.unit = "secs", # time units of the original raw data output.unit = "mg/h/g", # desired output unit volume = 0.12, # effective volume of the respirometer in L mass = 0.0009) # mass of the specimen in kg
The converted rates can be seen in the summary table.
summary(zeb_cr.int_conv)
We can see the converted rates as the last column. The $summary table contains the all rate regression parameters and data locations, adjustments (if applied), units, and more. It can be exported as a separate data frame by using export = TRUE in summary(). This is a great way of exporting all the relevant data for your final results.
Select {#crintselect}
It will depend on the aims and circumstances of an experiment what rate results to extract, summarise, and report. For instance we might report RMR as a mean of a representative selection of rates after the specimens uptake rate has stabilised, or a mean of a number or percentile of lowest rates found. The important thing is to apply these selection criteria consistently and report them alongside results.
As of respR v2.1 convert_rate has new plotting functionality to help visualise and explore rate results. For intermittent-flow results the most useful of these is type = "rate".
plot(zeb_cr.int_conv, type = "rate")
This shows how the calculated rate values vary across the dataset. We can see from this plot they are high at the start, stabilise after around timepoint 20000, then are a little unstable again if not slightly elevated after timepoint 60000. Between these however they appear to be very stable.
The select_rate() function allows convert_rate results to be filtered or reordered in various ways. See help("select_rate") and vignette("select_rate") for more details. For this analysis, lets only use rates from replicates between timepoints 20000 and 60000. Then, simply to show an example of how rates can be further selected, let's only select those with r-squared of 0.97 and above, select the lowest 20, and finally take the mean of them.
RMR <- zeb_cr.int_conv |> select_rate(method = "time", n = c(20000, 60000)) |> select_rate(method = "rsq", n = c(0.97, 1)) |> select_rate(method = "lowest", n = 20) |> summary() |> mean()
That's it; this is our final routine metabolic rate (at least, as we have chosen to define it).
SMR
auto_rate.int {#arint}
To extract a rate from every replicate and determine a standard metabolic rate (SMR) from amongst the results we will use auto_rate.int(). To repeat the note from above, this is not a recommendation that this is how you should extract an SMR from your own data; it is simply an example of one approach and how the respR functions are flexible and adaptable.
Introduced in respR v2.1, auto_rate.int() is a dedicated function for intermittent-flow respirometry data, in which you specify a replicate structure and run auto_rate() on each replicate. This function extracts rates in a number of ways including the most linear or most consistently maintained rates, plus highest and lowest rates of a specified width. See vignette("auto_rate") for more details, and vignette("auto_rate.int") for an overview of this function.
The difference between this and the approach above using calc_rate.int() is that the specified measure phase is not the same as the region over which a rate is calculated. Instead a measure phase is specified and the auto_rate analysis is done within this region using the method and width inputs. This can result in rates of different widths coming from within this region depending on the inputs.
Subset data
We again subset the data so that the 105 regularly-spaced replicates start at row 1.
zeb_insp <- zeb_intermittent.rd |> subset_data(from = 5840, to = 75139, by = "row") |> inspect()
Extract rates {#arintextract}
To get SMR we necessarily need to define it, in particular the duration over which rates are calculated within each replicate. These should be chosen appropriately and with reference to best practices and previous studies. Here, let's say we want to extract the lowest rate across 2 minutes within each replicate, and we will define our SMR for this specimen as the mean of the lowest 10th percentile of all 105 extracted rates.
We run auto_rate.int on the inspected data and use method = "lowest" and width = 120. This means the function calculates every rate of this width within each replicate, orders them according to the input method, and by default (n = 1) returns the top ranked result, in this case the lowest, for each replicate.
We again specify a 660 row interval between replicates using starts. Just to be different from the above example, this time we will specify a 2 minute (120 row) wait phase, and six minute (360 row) measure phase.
zeb_ar.int <- auto_rate.int(zeb_insp, starts = 660, wait = 120, measure = 360, method = "lowest", width = 120, by = "row") |> summary()
Note how the plots show the rates coming from different regions within each replicate measure phase. The summary contains all the results, that is the lowest 2-minute rate from each replicate. Before selecting from amongst these to get our final SMR we first need to adjust and convert them. This is exactly the same process we followed in the example above.
Adjust
zeb_ar.int_adj <- adjust_rate(zeb_ar.int, by = bg_pre, by2 = bg_post, method = "linear")
Convert
zeb_ar.int_conv <- convert_rate(zeb_ar.int_adj, oxy.unit = "mg/L", # oxygen units of the original raw data time.unit = "secs", # time units of the original raw data output.unit = "mg/h/g", # desired output unit volume = 0.12, # effective volume of the respirometer in L mass = 0.0009) # mass of the specimen in kg
summary(zeb_ar.int_conv)
Quality control
We can plot the convert_rate object for a quick look at how rates vary across the data.
plot(zeb_ar.int_conv, type = "rate")
These rates seem to be a bit lower than those we extracted for RMR above, as we would expect for an SMR.
This is an opportunity to do some quality control before we select our final rates to determine SMR. Looking at the plot above there is one result that seems to be uncharacteristically low compared to the others with a value just over 0.5. Looking at the summary table we can see this is from replicate 48. We can plot a selection of rates from this region using pos and the y-axis range will adapt to provide a more detailed view of the range of rates.
plot(zeb_ar.int_conv, type = "rate", pos = 30:80)
Replicate 48 does seem to be anomalously low. If we go back to the auto_rate.int results object we can use pos to plot this replicate for a closer look. We'll do this in two different ways using type. Firstly the rate in the context of the whole replicate, secondly the auto_rate output plot for this replicate which has more diagnostic plots, but only shows the measure phase.
plot(zeb_ar.int, type = "rep", pos = 48) plot(zeb_ar.int, type = "ar", pos = 48)
Looking at these plots this result does seem to be anomalous, and in this replicate the rate extracted is not really representative. It also has a lower r-squared which otherwise varies from 0.72 to 0.98. There seems to be an anomaly in the oxygen trace, maybe a sensor glitch or slight temperature change. Issues like this causing glitches in the trace are very common.
Therefore, it would be perfectly acceptable and justifiable to manually exclude this result (there are also 1 or 2 others we might want to investigate). To apply such selection criteria objectively and systematically would require examining multiple experiments. This is why definitions of metrics like SMR and analysis parameters should usually be decided upon after thoroughly exploring all your results. After this you should be able to arrive at a set of objective and reasonable selection criteria.
Select
To get our final SMR we just need to use select_rate again with our chosen criteria, and pipe the result to mean. We decided above our definition of SMR would be the mean of the lowest 10th percentile of the rates extracted from each replicate. We also want to exclude the result from replicate 48.
Excluding this replicate could be achieved in a number of ways in select_rate. We could use method = "rsq" to select only results with r-squared above a certain value which excludes this one. We could use the time_omit method to exclude this rate based on its timepoint (found in the summary table), or use the rep method simply to select all replicates except this one. There is a specific rep_omit method however to remove particular replicates.
SMR <- zeb_ar.int_conv |> select_rate(method = "rep_omit", n = 48) |> select_rate(method = "lowest_percentile", n = 0.1) |> summary() |> mean()
We can see the summary of the final 11 selected rates, and then the mean of these. This is our final SMR.
Summary {#ratessumm}
Now we have our final results for each metric.
MMR: -3.05 mg/h/g RMR: -0.94 mg/h/g SMR: -0.70 mg/h/g
Typically these rates would now be reported as positive values.
Reporting {#report}
The above is a lengthy description of an analysis, so this vignette may not give a good impression of just how succinct a respR analysis may be. This code block shows the entire analysis for SMR.
Complete analysis
# Inspect data ------------------------------------------------------------ zeb <- inspect(zeb_intermittent.rd) # Background -------------------------------------------------------------- bg_pre <- subset_data(zeb, from = 0, to = 4999, by = "time") |> calc_rate.bg() bg_post <- subset_data(zeb, from = 75140, to = 79251, by = "time") |> calc_rate.bg() # SMR Analysis ------------------------------------------------------------ # Subset ------------------------------------------------------------------ zeb |> subset_data(from = 5840, to = 75139, by = "row") |> # Extract rates ----------------------------------------------------------- auto_rate.int(starts = 660, wait = 120, measure = 360, method = "lowest", width = 120, by = "row") |> # Adjust ------------------------------------------------------------------ adjust_rate(by = bg_pre, by2 = bg_post, method = "linear") |> # Convert ----------------------------------------------------------------- convert_rate(oxy.unit = "mg/L", time.unit = "secs", output.unit = "mg/h/g", volume = 0.12, mass = 0.0009) |> # Select ----------------------------------------------------------------- select_rate(method = "rep_omit", n = 48) |> select_rate(method = "lowest_percentile", n = 0.1) |> summary() |> mean()
How to report analysis methods
The ultimate aim of respR was to enable clear reporting and reproducible analyses. This is an example of how we might report the above analysis in a manuscript.
"Data was analysed using the R package respR (Harianto et al. 2019). The auto_rate.int function was used to extract the lowest rate across a two minute window from each replicate. Rates were adjusted for background using pre- and post-experiment blank recordings with no specimen, assuming a linear increase in background rate over time. Obviously anomalous rates were excluded, and SMR was defined as the mean of the lowest 10th percentile of remaining adjusted rates."
Alternative approaches
Previous versions of this vignette demonstrated use of different methods, such as looping and subsetting of data, to extract individual replicates for processing through calc_rate or auto_rate. These approaches still work and are archived here in case they are of use in some cases, for example if you want to extract rates from the same oxygen range in each replicate. However for most analyses use of calc_rate.int or auto_rate.int is the recommended approach.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.