BMR 617: Statistical Techniques for the Biomedical Sciences

Analyzing Real-Time Quantitiative PCR Data

Real-time quantitative PCR (RT-qPCR) is used to quantify the amount of DNA present in a sample. It is almost always used to quantify RNA, by using reverse-transcriptase to create cDNA from the RNA, and then to quanitfy the cDNA.

RT-qPCR works by replicating the DNA present in each machine cycle. The assay starts with a quantity of DNA that is too small to detect, and repeatedly amplifies it until it reaches some threshold, which is comfortably above the minimum detection level. The number of cycles to reach the threshold is recorded. The basic idea here is that the more DNA was present at the beginning, the fewer cycles it will take to reach the threshold.

Suppose the experiment starts with some number of DNA molecules, which we'll call \(A\). On each cycle of the experiment, each DNA molecule is replicated. So the number of molecules after one cycle is \(2\times A\). After two cycles it is \(2\times 2\times A\), after three it is \(2\times 2\times 2\times A\), and in general after \(t\) cycles it is \(2^tA\).

If the threshold is \(T\) and the number of cycles to reach the threshold is \(C_T\), then we must have \[T=2^{C_T}A\] and \[\frac{A}{T}=2^{-C_T}\]

One general problem with RT-qPCR is that it's difficult or impossible to control the total amount of RNA extracted from a sample. To account for this, from each sample the cDNA from a gene of interest is amplified, along with cDNA from a control gene, which is assumed to have constant expression across all samples and experimental conditions. Such control genes are called "housekeeping genes". Common choices for these genes are GAPDH and β-actin.

If \(g\) is the amount of cDNA from the gene of interest, \(h\) is the amount of cDNA from the housekeeping gene, \(C_{T,g}\) the threshold cycle for the gene of interest, and \(C_{T,h}\) the threshold cycle from the housekeeping gene, we have \[\frac{g}{T}=2^{-C_{T,g}}\] and \[\frac{h}{T}=2^{-C_{T,h}}\] Dividing one of these by the other gives \[\frac{g}{h}=\frac{2^{-C_{T,g}}}{2^{-C_{T,h}}}\] or \[\frac{g}{h}=2^{-(C_{T,g}-C_{T,h})}\] If we denote the difference in \(C_T\) values between the gene of interest and the housekeeping gene by \[\Delta C_T=C_{T,g}-C_{T,h}\] then \[\frac{g}{h}=2^{-\Delta C_T}\] In other words, \(2^{-\Delta C_T}\) represents the amount of cDNA from the gene of interest, relative to the amount of cDNA from the housekeeping gene.

While this is rarely noted, we can also write this in log form: \[-\Delta C_T=\log_2\left(\frac{g}{h}\right)\] So \(-\Delta C_T\) is the log (base 2) of the amount of cDNA from the gene of interest relative to the housekeeping gene.

Let's grab some sample data to work with. These data are kindly provided by Dr. Chris Risher at Marshall University.

These data comprise a simple experiment to compare the expression of a gene called α2δ-1 in the brains of male and female rats (i.e. we are comparing expression between male rat brains and female rat brains). 


library(tidyverse)
pcr_data <- read_csv("https://denvirlab.marshall.edu/BMR617-2024/data/M-F_a2d1_qPCR_quant.csv")
The data are presented as CT values, for the target gene and the housekeeping gene, which in this case is GAPDH. There are nine male rat brains and nine female rat brains. 

pcr_data

As described above, it's difficult or impossible to control the total amount of mRNA extracted from each brain. To account for this, we use a housekeeping gene (GAPDH), for which we assume the mRNA expression is constant across all experimental samples. The quantity \(\Delta C_T\) (the difference in \(C_T\) values between the gene of interest, α2δ-1 and the housekeeping gene) is the negative log (base 2) of the mRNA expression of α2δ-1 relative to the expression of GAPDH.

We can start by adding a column to the data table for \(\Delta C_T\): 


pcr_data <- pcr_data %>% mutate(DeltaCT = a2d1_CT - GAPDH_CT)
View pcr_data (by clicking on it in the "Environment" tab). Pick a couple of rows at random and verify that the DeltaCT value really is the \(C_T\) value for the target gene, minus the \(C_T\) value for the housekeeping gene.

Note that we have a \(\Delta C_T\) value for each sample, which was computed by subtracting a different value from the target \(C_T\) for each sample (the housekeeping \(C_T\) for that sample). We did this because we are correcting for the fact that there is a different total amount of mRNA in each sample.

Remember from above that \(2^{-\Delta C_T}\) is the amount of cDNA (representing the amount of mRNA) in the target gene, relative to the housekeeping gene. Let's just look at what these values look like (but we won't keep them): 


pcr_data %>% mutate(ExpRelGAPDH = 2^-DeltaCT)
This gives 

# A tibble: 18 × 6
   SampleID Sex   a2d1_CT GAPDH_CT DeltaCT ExpRelGAPDH
                        
 1 Sample1  M        34.7     25.9    8.81    0.00223 
 2 Sample2  M        35.6     25.8    9.79    0.00113 
 3 Sample3  M        35.6     25.7    9.88    0.00106 
 4 Sample4  M        34.6     25.7    8.93    0.00206 
 5 Sample5  M        34.6     25.8    8.85    0.00217 
 6 Sample6  M        34.4     25.7    8.75    0.00232 
 7 Sample7  M        36.4     26.9    9.50    0.00138 
 8 Sample8  M        36.8     26.9    9.85    0.00109 
 9 Sample9  M        36.9     26.8   10.1     0.000912
10 Sample10 F        36.9     26.4   10.5     0.000674
11 Sample11 F        36.1     26.4    9.78    0.00114 
12 Sample12 F        36.3     26.4    9.91    0.00104 
13 Sample13 F        36.0     24.8   11.3     0.000408
14 Sample14 F        36.1     24.8   11.3     0.000395
15 Sample15 F        35.3     24.7   10.7     0.000618
16 Sample16 F        36.3     25.9   10.4     0.000720
17 Sample17 F        34.9     26.0    8.93    0.00205 
18 Sample18 F        35.8     25.9    9.89    0.00105
We could actually work with these values, but they're not very intuitive. The units don't make much sense. The only interpretation here is that, for example, sample 1 expresses 0.00223 times as much α2δ-1 as it does GAPDH.

To remedy this, we will change units so that the mean \(\Delta C_T\) in one of our groups is zero. Typically, we use the control group for this; in this case we can arbitrarily choose one group, and we'll choose the female rat brains.

Here all we are doing is effectively changing units. To do this in a valid way, we must subtract the exact same value from all the \(\Delta C_T\)s. The value we will subtract is the mean \(\Delta C_T\) from one of the groups. This value is called a callibrator in the PCR literature. We first filter to get the female samples, and then pull the \(\Delta C_T\) values from them. Then we'll take the mean: 


f_delta_CTs <- pcr_data %>% filter(Sex == "F") %>% pull(DeltaCT)
callibrator <- mean(f_delta_CTs)
Verify that the callibrator is just a single value (not a collection of values): 

callibrator
Now we subtract the callibrator from each \(\Delta C_T\). Since this is a difference between the \(\Delta C_T\) for the sample, and an average \(\Delta C_T\) it's the difference between two \(\Delta C_T\)s, and it's known as the \(\Delta\Delta C_T\): 

pcr_data <- pcr_data %>% mutate(DeltaDeltaCT = DeltaCT - callibrator)
We saw before that \(2^{-\Delta C_T}\) was the expression of the target gene in each sample, relative to the housekeeping gene. We called this quantity \(\frac{g}{h}\). Since the callibrator is the average \(\Delta C_T\) for one of our groups, \(2^{-\text{callibrator}}\) is essentially the average expression of the target gene relative to the housekeeping gene in the female rat brains. (Technically, it is the geometric mean of these values.) So now if we take \(2^{-\Delta\Delta C_{T_i}}\) for each of our samples, we have \[ \begin{aligned} 2^{-\Delta\Delta C_{T_i}}&=2^{-\left(\Delta C_{T_i} - \text{callibrator}\right)} \\ &=\frac{2^{-\Delta C_{T_i}}}{2^{-\text{callibrator}}} \\ &=\frac{\frac{g_i}{h}}{\text{Average}_\text{Female}\left(\frac{g}{h}\right)} \\ &=\frac{g_i}{\text{Average}_\text{Female}g} \end{aligned} \] In other words, \(2^{-\Delta\Delta C_T}\) is the expression of the target gene in the sample, relative to the average (the geometric mean) of the expression of the target gene in the Female group. Note here that we assumed the expression of the housekeeping gene was constant.

Let's add these values to our table: 


pcr_data <- pcr_data %>% mutate(RelativeExpression = 2^-DeltaDeltaCT)
We can summarize our values. Note the mean \(\Delta\Delta C_T\) value should be zero. The mean relative expression will not be 1, because we effectively used the geometric mean of the expression for the female group as the normalizer, but it should be reasonably close: 

pcr_data %>% group_by(Sex) %>% 
  summarise(DeltaDeltaCT = mean(DeltaDeltaCT), RelativeExpression = mean(RelativeExpression))
which gives 

# A tibble: 2 × 3
  Sex   DeltaDeltaCT RelativeExpression
                        
1 F         5.92e-16               1.14
2 M        -9.17e- 1               2.01

The preceding work is published in a famous paper by Livak and Schmittgen. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method . Methods. 2001 Dec;25(4):402-8. doi: 10.1006/meth.2001.1262. PMID: 11846609.

The following are important to remember when using this method:

The first Δ is the difference between the \(C_T\) for the target gene and the housekeeping gene. It is calculated on a per-sample basis; that is, the value subtracted will be different for each sample. This is because we are adjusting for the fact that a different amount of total RNA may have been extracted for each sample.
The second Δ is the difference between the \(\Delta C_T\) for each sample and the average (mean) \(\Delta C_T\) for the control group. In this subtraction, the same value is subtracted from every \(\Delta C_T\). This is because we are simply changing units to give a more meaningful scale.

Hypothesis testing for PCR data

In the example data, we want to test the hypothesis that the expression of α2δ-1 is different in male rat brains to female rat brains. There are two schools of thought on how to do this:

Livak and Schmittgen, in the paper cited above, say:

The endpoint of real-time PCR analysis is the threshold cycle or \(C_T\). The \(C_T\) is determined from a log-linear plot of the PCR signal versus the cycle number. Thus, \(C_T\) is an exponential and not a linear term. For this reason, any statistical presentation using the raw \(C_T\) values should be avoided.

On the other hand, in their paper Statistical Analysis of real-time PCR data, BMC Bioinformatics, 2006, Yuan et al. say

The primary assumption with this approach is that the additive effect of concentration, gene, and replicate can be adjusted by subtracting \(C_t\) number of target gene from that of reference gene, which will provide \(\Delta C_t\) as shown in Table 2. The \(\Delta C_t\) for treatment and control can therefore be subject to simple t-test, which will yield the estimation of \(\Delta\Delta C_t\).
Livak and Schmittgen here are arguing that because the \(C_T\) values are used as exponents in the equation for the relative quantity of DNA in the experiment, that variation in the \(C_T\) values will be considerably smaller than variation in the quantity of DNA. While this is true, the argument is unconvincing; the quantities themselves will also be considerably larger, and so will any differences between them.

Yuan et al. are arguing that the \(\Delta C_T\) values are subject to a number of random effects which are additive, and therefore are likely to be normally distributed. This is a slightly more convincing argument, but it lacks some justification.

Really, the question should boil down to which values are sampled from a normal distribution. The way to really test this would be to perform real-time PCR on a large number of biologocial samples and then examine the distribution of the \(C_T\) values and the \(2^{-\Delta\Delta C_T}\) values. To my knowledge, no-one has done this.

Let's try both approaches here. 


	t.test(DeltaCT ~ Sex, data=pcr_data)
This gives a p-value of 0.01037. The mean \(\Delta C_T\) for the female group is 10.300471 and for the male group is 9.383845. Note the difference in means is 9.383845-10.300471=0.916626, which is the \(\Delta\Delta C_T\) value for the male group, as expected. 


	t.test(RelativeExpression ~ Sex, data=pcr_data)
This gives a reasonably close p-value of 0.01673.

Either of these approaches is acceptable, but be sure to provide enough information showing what you did. E.g. for the first approach, you could say "differences between expression of α2δ-1 between male and female rat brains were tested using a T-test on the \(\Delta C_T\) values, following Yuan et al. (2006), giving p=0.01037. The average fold change was calculated as 1.88 fold, using the \(\Delta\Delta C_T\) method of Livak and Schmittgen (2001).

For the second approach, you might say "relative expression was calculated using the \(\Delta\Delta C_T\) method of Livak and Schmittgen (2001) and compared between males and females using a T-test. The resulting p-value was 0.01673, with a fold change of 1.88 fold.

Finally, we can create plots using the relative expression. For example: 


ggplot(pcr_data, aes(x=Sex, y=RelativeExpression)) + 
  geom_boxplot(outlier.shape=NA) + 
  geom_point(position=position_jitter(width=0.1))
or 

pcr_data %>% group_by(Sex) %>% 
  summarise(REsd=sd(RelativeExpression), RelativeExpression = mean(RelativeExpression), n=n(), sem=REsd/sqrt(n))  %>%
  ggplot(aes(x=Sex, y=RelativeExpression)) + 
    geom_col(fill='#00b140') + 
    geom_errorbar(aes(ymin=RelativeExpression-sem, ymax=RelativeExpression+sem), width=0.2)