Should be RT-qPCR values standardized before PCA analysis?

Question

I come from this question of gene exression analysis: Should PCA be standardized for gene expression?

My experiment is based in 151 x 51 (individuals/samples x genes), in which, patients are subjected to 3 possible groups (~ balanced groups), so we have 2 measurements of gene expression: one in the baseline, and one after the intervention. The gene expression is actually summarised in an scalar (number) called fold-change. Fold-change is calculated through the Ct (cycle threshold), which is the number of the cycle in which you are able to detect a certain amount of the target gen. This value or cycle is normalised using Ct values of 3 reference genes which are supposed to remain stable, and each patient is measured using oneself baseline value, in the same run, to remove run-to-run variation. As you can see I am working with RT-qPCR and relative quantification values.

As you can see, "normalization" or "standardization" is performed, not in the strict meaning of mathematical foundation, but to normalize quantitation of targets in the experiments. I was wondering if, taking into account that normalization is done, when I try to perform PCA to see how these genes behave, if they form cluster, if the intervention received has influence of the expression, patterns...

Reading about, I found quite divergent opinions about scaling or standardizing the results. The conflict, as far as I can see in my own experiment, lies in all my variables are genes, measured with the same high-throughput analyzer, different batches, but normalized with reference genes, and their own baseline values.

On one hand I found quite homogeneous working with scaled values, but on the other hand, I am working with values that seem to be working on the same scale, and have been through several filters to normalize results. Not sure if executing more mathematical operations could modify in the wrong way. What is your opinion about it? (all is about using the covariance matrix or correlation matrix based PCA, as far as I understand)

Furthermore, the dataset contain missing values, so I am working with the missMDA package to impute them. In case I scaled the data, what should I do first: scale or impute the missing values?

The steps to calculate fold-change:

Update 1: About the technology used, is a OpenArray plate (analyzed with QuantStudio™ 12K Flex Real-Tim) which contains subarrays with customised target genes.It is a 3,072-reaction high-throughput real-time PCR. The reaction are individually performed (as far as I know, getting Ct for each sample and assay). About the 2nd point, it is 2x time series study (baseline and post-intervention)- But if you see the equations to calculate the fold-change, I am going to obtain a single value of gene expression combining the 2 Ct or measurements for each patient (as I told: baseline and post-intervention).

About the "genes aim" you can see that my interest is purely human genes, after an intervention (inflammation-related genes).

About the pippeting issue. As you know is never discarded, but the internal controls or reference genes (as I am adressing them) goal is to normalize results and remove potential quantification errors. However, as you said the RNA concentration (A260) and purity were calculated spectrophotometrically from peripheral blood samples origin

As far as I know, the terms correlation and covariance related to PCA, are meant to discrimanate the standardization process (in correlation matrix) and not standardized matrix (scale function in R). Not sure if I understood this is what you imply?

My guess is nowadays everything is RNAseq and microarrays, whose mechanisms, and therefore, statistical preprocesses are different from RT-qPCR. But there are papers of PCA working wiht data coming from RT-qPCR. Shouldn't it work though?

Update2

Covariance and correlation: "Covariance analysis of PCA: PCA is dimensionality reduction based on maximising multi-dimensional variance so that would be unusual, but correlation matrix of PCA is understandable. Normally covariance analysis is performed directly on standardised data." Let me insist on this topic a little bit more.I am not sure if I understood correctly what you stated, but the information I've reviewed seems to say otherwise. As far as I was able to find, https://aedin.github.io/PCAworkshop/articles/b_PCA.html : "PCA was computed as a singular value decomposition (SVD) of a column centered, scaled matrix. This was PCA of the correlation matrix. If the matrix is centered but not scaled, it is PCA of the covariance matrix." I am not sure if I understood correctly what you stated, but the information I've reviewed seems to say otherwise

"The correlation matrix is the standardized version of the covariance matrix. Analysing the correlation matrix is a useful default method because it takes the standardized form of the matrix; therefore, if variables have been measured using different scales this will not affect the analysis. Often you will want to analyse variables that use different measurement scales. Analysing the correlation matrix ensures that differences in measurement scales are accounted for. In addition, even variables measured using the same scale can have very different variances and this too creates problems for principal component analysis. Using the correlation matrix eliminates this problem also. There are statistical reasons for preferring to analyse the covariance matrix (the reason being that correlation coefficients are insensitive to variations in the dispersion of data whereas covariance is and so produces better-defined factor structures (Tinsley & Tinsley, 1987)) and generally the results will differ from analysis on the correlation matrix. However, the covariance matrix should be analysed only when your variables are commensurable." Conclusion: shouldn't it be the covariance matrix the one without standardized values, because variables share same scale?

Missing values: this is an issue, at least with several functions in R (prcomp(), princomp())which don't work with NAs in the matrix.

About this statement: "The key question, which is standardisation and I would say it is essential for RT-qPCR derived data at a minimum to at least investigate consistency of signal against a control. That minimum is what the second equation describes (why it works is complicated)."

I don't really understand what you mean standardize at a minimum. The endogenous control genes we are using are individualized for every sample. Let say we are investigating the fold-change in the individual i and the variable k. We are using endogenous controls in the i individual to "standardize" values. The

Answer 1

This is as good an answer I can give based the information. The issue here is that the experiment is derived via RT-qPCR and thats not the same as RNAseq type approach that the normal basis for questions on this site. Generally, we sort of assume RNAseq for gene expression where all the gene expression within a sample are automatically standardised against each other. Thus PCA can be performed directly on that data type. I will summarise this at the end because it's the key point I'm trying to make.

In your experiment there are 51 separate reactions per patient (or is 51 the total time series e.g. for one gene??) and it's not clear whether this constitutes a time series (but there must be some time-series aspect to the study: I assume the 'baseline' is at least one point in a time series). That is a lot of opportunity for variation within one patient sample, even if it is the same time point. Maybe the tech doing the RT-qPCR makes a pipetting error for 1 in 10 samples - the results could collapse. The other issue is whether the RT-qPCR is internally standardised, i.e. are you assessing two amplicons per PCR and the same control amplicon is used for each different gene under investigation? This traps the 1 in 10 hypothetical error rate because that is the first point of standardisation. The question is whether a further standardisation is required beyond that, i.e. per patient - thats tricky but basically its difficult to avoid because of variation between different samples from the same patient (again you must be doing some sort of time series). So for example the RNA extraction could have different amounts of blood etc ... However, it needs some thought because here you are standardisation a standardised result against a backdrop of an unstated hypothesis (comments).

Summary The single point I want to make is that large scale RT-qPCR is quite different from the standard RNA-seq approach on this site from a statistical perspective, i.e. it appear to be 51 reactions per sample whereas normally the approach is 1 reaction per sample which generates all the expression data across thousands of genes, thus PCA is a fantastic analysis for this data-type. This difference creates a layer of statistical complexity and I what I would do is carefully consider my hypothesis and formulate an analytical plan in accordance with the aims.

The question is whether a simple standardisation could be performed and then shift straight into PCA? Yes. Is that sufficient? I dunno. I think you might be avoiding time-series analysis by standardisation and then subtraction of the baseline which is then placed into PCA. At a minimum I would perform two PCA analyses one for baseline and one for post-treatment.

Generally, a complex standardisation (like I'm inferring) is not cool if the next stage of the analysis is PCA.

Final thoughts Although PCA is considered the primary analysis here. I think it takes serious thinking given all the available information to formulate an analytical plan from my experience in this type of data, but as a rough guide (given I dunno the experimental aims) a complex standardisation and then I personally would consider covariance thereafter for a time series rather than PCA as general approach here. However, two PCAs before and after using a simple standardisation and compare presence/absence, that could work and very clear results would be accepted.

---

Note Just to note that "correlation matrix based PCA" is basically PCA, it's just trying to sort out the groups. However "covariance" is a different analysis to PCA and statistically and the aims are very different. Also just to mention, but I assume you know this, "correlation" is different from "correlation matrix based PCA".

Missing values Missing values for PCA isn't an issue - you just miss them out. Missing values in a time series can have lots of different solutions, the key is justification.

---

Update

Covariance analysis of PCA: PCA is dimensionality reduction based on maximising multi-dimensional variance so that would be unusual, but correlation matrix of PCA is understandable. Normally covariance analysis is performed directly on standardised data.

Overall, I am now clearer one part of the calculation: this is essentially an ad hoc statistic to assess fold-change of the raw Ct via a single /minimal PCAs following one or more standardisations.

What you are measuring is the optimised multi-dimensional variance of the fold-change of raw Ct after at least one or else numerous standardisations. So this is not the actual fold-change of Ct, but variance of the fold-change derived via Eigen values.

If you are clear about that and what this represents and the reader is clear about this then it would be up to the reviewers.

The key question, which is standardisation and I would say it is essential for RT-qPCR derived data at a minimum to at least investigate consistency of signal against a control. That minimum is what the second equation describes (why it works is complicated).