In LucasFVoges/DoubleqpcR: DMAS-qPCR data analysis

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Library DoubleqpcR

library(DoubleqpcR)
library(ggplot2)
library(investr)
calculateEverything = FALSE
# Set this to true if all things should be calculated (time consuming)

InputCq <- system.file("extdata", "rps4_8fach_cq.csv", package = "DoubleqpcR")
InputWells <- system.file("extdata", "rps4_8fach.csv", package = "DoubleqpcR")
bitterData <- system.file("extdata", "DMASqPCR_bitter_marzipan.csv", package = "DoubleqpcR")
sweetData <- system.file("extdata", "DMASqPCR_sweet_marzipan.csv", package = "DoubleqpcR")

Data Handling

To import the data the following can be used. The input.cq, input.raw, input.raw.melt, input.raw.melt.usedOnly are necessary and will be used globally! But they can be created without these functions, of course.

The example data used here is a dilution series and therefore the samples are numeric concentrations of the target genotype. Some methods for the regression will expect, that the sample column are indeed numbers! Sample 100 means a 100% target genome concentration (This can be confusing, because it is not the actual concentration).

read.cqTable(csv = InputCq) # mistake in name on purpose to demonstrate
read.fluorescenceTable(csv = InputWells)

Cq value determination

To add more Cq values from fluorescence curves any method can be used. DoubleqpcR has two implemented methods and a wrapper for the qpcR package. Both use the first and second derivative maxima to determine the Cq value. (FD = First Derivative = TP = Tourning Point) and (SD = Second Derivative = exponential slope)

if(calculateEverything){
# Calculate new Cq values (needs fluorescence table!)
calc.Cq(method = "TP", cq.new = "TP")
calc.Cq(method = "SD", cq.new = "SD")

calc.Cq.qpcR(method = "TP", cq.new = "TP.qpcR")
calc.Cq.qpcR(method = "SD", cq.new = "SD.qpcR")
}

Get the Cq values

This step is optional and more for internal function use. But if the Cq values (with outlier removal) for a sample are needed in a data format, one can use the following.

myCqValues <- table.Cq(sample = 8.8, target = "bitter", CqType = c("Regression"), outliers = TRUE, alpha = 0.05, format = "data")

Plot the fluorescence curves {.tabset .tabset-fade}

To change the visuals you can overwrite the ggplot theme settings and details. This is commented in the plot for all curves.

p <- plot.curve(sample = "all")
p
# p + theme_tq() + scale_color_tq() + theme(legend.position = "none")

An Example:

We choose Dixon outlier test here, as we have only 4 measurements per individual subsample (or no outlier removal). One can later choose outlier detection again when combining the subsamples.

# Here is the command in simpler form:
plot.curve(sample = "81.4.a", target = "bitter")  # sample can be also "81.4". then no percentage is shown.
table.Cq(sample = "81.4.a", target = "bitter", CqType = c("Regression", "Autocalculated"), outliers = TRUE, outliers.method = "Dixon", alpha = 0.05, format = "kable")

# plot all samples with Headings, ideal for tabsets.
if(calculateEverything){
for (sampleX in unique(input.cq$sample)) {
   cat('\n\n#### Sample `', sampleX, '`\n\n')
   plot(plot.curve(sample = sampleX, target = "bitter"))
   cat('\n\n')
   show(table.Cq(sample = sampleX, target = "bitter", CqType = c("Regression","Autocalculated","TP.qpcR","SD.qpcR","SD","TP"), outliers = TRUE, alpha = 0.05, format = "kable"))
}
}

Combine data for further analysis

At this point a new list of dataframes (data.cq) will be created via make.data.cq() which looses the connection to the fluorescence table. The steps above can be repeated and added to data.cq via the option add = TRUE. This procedure allows to combine multiple experiments. Cq value types that are wanted to use needs to be selected.

If all values are already present in the input.cq like it is the case in this example we only need to make the data.cq list.

make.Cq.data(add = FALSE, target = "bitter", CqType = c("Regression", "Autocalculated"), outliers = TRUE, outliers.method = "Dixon", alpha = 0.05, silent = TRUE)

In this example we have seperated two experiments with an indicator "a" and "b". It is now time to combine the two individual subsamples to one.

combineSubsamples(delimiter = ".", outliers = TRUE, outliers.method = "Grubbs")

Summary of Cq values

A summary of the cq values can be given via a data.cq.sum data frame. Samples which are names numeric can be filtered.

Cq.data.mean(CqType = "Regression", onlyNumeric = TRUE)
data.Cq.sum

Plot Cq values

Here all samples that are in the data.cq list will be plotted for one Cq type.

plot.Cq(CqType = "Regression", onlyNumeric = TRUE)

Here all samples that are in the data.cq list will be plotted for one Cq type against the log of the concentration. Therefore the sample name has to be a numeric value which corresponds to the target concentration.

plot <- plot.Cq.efficiency(CqType = "Regression")

plot

Delta Cq Values

To get the delta Cq values the function delta.Cq.data() can be used. It has different methods to calculate the differences...

Here we use all combinations of Cq values from one sample.

delta.Cq.data(CqType = "Regression", method = "c", onlyNumeric = TRUE)
head(delta.Cq)

Plot delta Cq values

A typical box plot representation of delta Cq values. This is of cause only for demonstration purposes as these plots can be made in various ways from the data with standard visualisation methods. (Obviously this makes no sense with "delta.Cq.data(method="mean")".

plot <- plot.delta.Cq.box(points = TRUE) 
plot + theme_bw()

Regression analysis (polynomial)

At this point all samples need to be numeric and represent a percentage!

All samples should be between 0 % and 100 %. Samples that are not numeric and between 0 to 100 will be skipped.

An overview can be plotted for different Cq types from data.Cq.

plot.delta.Cq.overview(CqType = c("Regression", "Autocalculated"))

Model generation

A model (linear, poly3, ...) can be used to describe the values from DMAS-qPCR. The R package investR is used for prediction and plotting. And the package caret for cross validation.

regression.delta.Cq(CqType = "Regression", fit = "poly3", cv = TRUE, cvComplete = TRUE)

Model analysis

The model can be analysed by standard methods:

# Model Summary:
summary(model.delta.Cq)

# Plot the residuals:
plot(summary(model.delta.Cq)$residuals, main = "Residuals for fit", ylab = "residuals")

Model prediction

Values can be predicted with investr package. The mean.response should be set to TRUE. But one could argue that the intrinsic bias from the combinatoric data for model creation is not fit for a mean response and a individual response should be considered (mean.response = FALSE). This leads to enormous uncertainty around delta CQ of 0.

predictionValue <- 0

res <- invest(model.delta.Cq, y0 = predictionValue, interval = "inversion", level = 0.95, mean.response = FALSE)

plotFit(model.delta.Cq, interval = "both",
         level = 0.95, shade = TRUE,
         col.pred = "whitesmoke",
         col.conf = "skyblue",
         extend.range = TRUE,
         xlab = "proportion of target genotype",
         ylab = bquote(Delta*Cq))
abline(h = predictionValue, v = c(res$lower, res$estimate, res$upper), lty = 2, lwd = 0.5, col = "salmon")

For the prediction y = r predictionValue the model predicts a response of r res$estimate with upper r res$upper and lower r res$lower interval.

Regression analysis (Linearisation)

At this point all samples need to be numeric and represent a percentage!

All samples should be between 0 % and 100 %. Samples that are not numeric and between 0 to 100 will be skipped.

An overview can be plotted for different Cq types from data.Cq.

plot.delta.Cq.overview(CqType = c("Regression", "Autocalculated"), linSqrtTrans = TRUE, xlab = "shifted square root transformation of proportion")

Model generation

The shifted square root transformation to linearise the data can be used with the parameter linSqrtTrans = TRUE

regression.delta.Cq(CqType = "Regression", linSqrtTrans = TRUE, fit = "linear", cv = TRUE, cvComplete = TRUE)

Model analysis

The model can be analysed by standard methods:

# Model Summary:
summary(model.delta.Cq)

Model prediction

Values can be predicted with investr package. The mean.response should be set to TRUE. But one could argue that the intrinsic bias from the combinatoric data for model creation is not fit for a mean response and a individual response should be considered (mean.response = FALSE). This leads to enormous uncertainty around delta CQ of 0.

predictionValue <- 0

res <- invest(model.delta.Cq, y0 = predictionValue, interval = "inversion", level = 0.95, mean.response = FALSE)

plotFit(model.delta.Cq, interval = "both",
         level = 0.95, shade = TRUE,
         col.pred = "whitesmoke",
         col.conf = "skyblue",
         extend.range = TRUE,
         xlab = "proportion of target genotype",
         ylab = bquote(Delta*Cq))
abline(h = predictionValue, v = c(res$lower, res$estimate, res$upper), lty = 2, lwd = 0.5, col = "salmon")

For the prediction y = r predictionValue the model predicts a response of r reTransform(res$estimate) with upper r reTransform(res$upper) and lower r reTransform(res$lower) interval.

Example

We look at an example from the study where we examined bitter and sweet almond samples. Here, two datasets from sweet M0 marzipan and from self-made bitter marzipan are used as an example. We use the regression analysis from before. The data in this part has no concentration and we do not look at the fluorescence data.

Import

Let use first import the bitter marzipan data. The .csv table is manually saved from the .xlxs excel supporting data file. We can see that there are two errors. First the column names are not as described by the package which are automatically renamed and we can ignore this error. Secondly a warning about the genotype is presented. This is from the NTC measurement and we can delete this row.

Next we can add this data to the data.cq dataframe. Actually we make a new one, because we do not want to add the data to the regression data we did before.

read.cqTable(csv = "DMASqPCR_bitter_marzipan.csv")

input.cq <- input.cq[1:(nrow(input.cq)-1),]

make.Cq.data(add = FALSE, target = "bitter", CqType = "Regression")

Additionally we can think about removing the "debittered 6" experiment as only values for one genotype are present. But here we just leave them be.

Now we add the values from the sweet Marzipan experiment. Here we have a problem with the sample names as they are called "Marzipan1" and so on. we have to manually change it to "Marzipan 1" to be coherent with the other naming...

read.cqTable(csv = "DMASqPCR_sweet_marzipan.csv")

input.cq$sample <- sub("([0-9])"," \\1",input.cq$sample)

make.Cq.data(add = TRUE, target = "bitter", CqType = "Regression")

Now we can combine our subsamples and have two entries in the data.Cq!

combineSubsamples(delimiter = " ", useMeans = TRUE)
data.Cq

Note that in this case we have searched for outliers two times. When we make or add the data.Cq in one subsample and when combining all subsamples. If the means are used, then the means will be compared for outliers! So in this case it makes sense. If we work with useMeans = FALSE than the outlier test in the make.Cq.data() function should be omitted!

Plot delta Cq values with prediction

Again the typical box plot representation of delta Cq values but now we include the predictions from the previously generated regression model.

First we have to generate the delta Cq values (this will overwrite the delta.Cq dataframe from regression analysis) As the model in the environment we have to tell the prediction to reTransform the values. We also plot the P value from a t.test for the two boxplots.

delta.Cq.data(method = "combinatorial")
plot.delta.Cq.box(useModel = TRUE, reTransform = TRUE, predictionRange = seq(-6,2,2), pVal = TRUE)

This should be an example for possible Boxplot analysis. These functions do not have the goal to be complete or to cover all possibilities.

Happy R-coding have a nice day!

LucasFVoges/DoubleqpcR documentation built on Feb. 19, 2024, 7:21 p.m.