Is it ok to interprete PCA plot this way?

I have several samples (C2, C4, C5) and want to check if they are at a certain stage. I included some known samples (D0 - D77) which were generated at different stages by another lab. In the PCA plot, my samples cluster together on the left and the known samples are dispersed on the right. I think the major difference among all the samples is different experimental protocols (PC1) and the second is different stages (PC2). So my samples are at the same stage. Is that right? And can we say my samples are at a stage between D12 and D19 (when projected to known samples, my samples are located between D12 and D19)? I have no strong mathematical background. Hope someone with math background can give some explanation. Thanks!

[ UPDATE1 ] I did this analysis using prcomp in R. The input is a 2D numeric matrix with 17436 rows and 42 columns. Each row represents a gene and each column represents a sample. The number is the gene expression level for a gene in a sample. The gene expression level is normalized using DESeq2 and thus the numbers are comparable across genes (rows) and samples (columns). For the 42 columns, 18 are from my experiments and the rest from published datasets. Besides different protocols, it is possible there are other differences. Generally, the table is a combination of two sources of data. In my data , C2, C4 and C5 are three cell lines which we processed in parallel in experiments. In published datasets, they sampled the cells at different time points (Day 0 to Day 77).

[ UPDATE2 ] R Code for PCA and plotting

# normCounts: normalized count from DESeq2, with 17436 rows and 42 columns
normCounts0 <- normCounts[ rowSums(normCounts) > 0, ]
tab <- t(normCounts0)
pca <- prcomp(tab, scale = TRUE)
tmp.x <- as.data.frame(pca$x) tmp.x$sample<-c(rep("C2_Con",3),rep("C2_KD",3),rep("C5_Con",3),rep("C5_KD",3),rep("C4_Con",3),rep("C4_KD",3), rep("D0",4),rep("D7",4),rep("D12",2),rep("D19",4),rep("D26",2),rep("D33",2),rep("D49",2),rep("D63",2),rep("D77",2))

require("ggplot2")
p <- ggplot(tmp.x, aes(x=PC1, y=PC2, color=sample))
p + geom_point() + scale_color_discrete(breaks=c("C2_Con","C2_KD","C5_Con","C5_KD","C4_Con","C4_KD", "D0","D7", "D12","D19","D26","D33","D49","D63","D77"))


• It appears you are mis-applying PCA with regard to this data. It is not as simple as making the association that PC2-->stage. When you look at the data for your known stages, they do not align with the y (PC2) axis. And your samples appear to be extreme outliers with respect to at least one of your input variables. What is the dimensionality of your inputs? You might consider generating some scatter plots to figure out which input feature is blowing up your variance in the first PC. – bogatron Jun 17 '15 at 20:39
• Thanks, bogatron. I updated the question with more details about the data. The differences between my 18 samples and 24 published samples are expected. Some points overlap on the PCA plot. On the left, there are 18 points. Sorry, this might mislead you. – Dejian Jun 18 '15 at 15:07
• If I understand your update to the question, you have a data set with 17436 variables (gene expressions) and 42 samples. So I'm curious how you're even doing PCA because your covariance matrix is definitely singular. I wonder if you have your data matrix transposed (inadvertently treating it as 42 variables with 17436 samples). Someone could help you better if you post the relevant portion of your code. – bogatron Jun 18 '15 at 20:16
• I transposed the matrix. The code for PCA is pasted in the update. Thanks. – Dejian Jun 19 '15 at 2:22

Given the information provided, I do not think it's enough to say that your samples are between stages D12 and D19. I admit knowing nothing about your domain, but with what you're showing it is very conceivable that other factors are causing the spread of your samples from the others on PC1 and it's too much of a stretch to say that your samples "project" into the others in between D12 and D19. Why does it have to be a projection orthogonal to the axis described by PC1? Remeber that PC1 is a linear combination of all your inputs. That is, it describes an arbitrary line through your input space that is probably not aligned with any one thing in particular. For example, all three of the following projections could be envisioned for your problem:

Another red flag for me is that you say PC1 is related to differences in experimental protocols. You haven't posted a scree plot or listed the eigenvalues for the different PCs but obviously PC1 soaks up the most variance so it's a pretty big deal. If you're saying that PC1 is correlated with experimental protocol, then that would seem good because hopefully the impact of experimental protocol is minimized in other PCs. However, that's never completely the case, really.

Again, not knowing your domain it's hard to tell what the impact of experimental protocol is on the distribution of your data in the PC plot, but I'd be VERY wary of comparing these data in this way if your data was produced using a different experimental design and, as you noted, "there may be other differences."

My hunch is that you would need to continue minimizing differences in the protocols between your experiments and the published experiments before you could really use this type of analysis. Others might disagree but if you posted the sources of your data, the eigenvectors, and possibly a subset of the data, experts in your field may be able to chime in on this methodology.

The short answer is: you could be right in your interpretation but at this point you don't have enough evidence to suggest your samples are "in between" D12 and D19 in terms of the stage. And I worry about differences in methodology completely swamping your ability to interpret things this way.

Don't just mix the data of two different sources.

It looks as if the two sources may simply be using different scales. But that makes the whole PCA analysis meaningless.

All your plot may be visualizing is that the other samples use a different way of measuring things. It seems all the variance is in the other samples. In particular this means that you cannot draw any conclusions about how your samples relate to the reference samples.

All what I'd say is that you are most probably right!

Your plot is compatible to the explanation you have for your data so most probably that is the case. I assume you know what features to use (see example) but still better to post some info about input data -e.g. the dimensionality of the data- to be sure everything is right there. The only reason which my ruin the explanation is presence of a confounder.

• Example

I'm not an expert of your field but just imagine there is a difference between materials you and the other lab used to get the data. Then I'll ask you how would you know that the difference on PC1 is not according to the difference of materials but difference of protocoles? you can confidently say that PC1 is capturing the difference of protocols if you are sure there is nothing else affecting it.

Hope it helped!

• Thanks, kasramsh. I update my questions, adding more info about the input dataset. You are right. It is possible there are other factors leading to the difference on PC1. Generally, that difference is due to two sources of data, not limiting to protocol itself. – Dejian Jun 18 '15 at 14:57
• What I am interested in is the difference on PC2. Based on the public datasets (points on the right), the difference on PC2 is due to time. What upsets me is whether I can project my samples (points on the left) to public data (points on the right) to conclude my samples are on stages between D12 and D19. I tried to find the answer by learning PCA algorithm but I'm still not sure. Thanks! – Dejian Jun 18 '15 at 14:57
• Regarding your update I'd say be careful about interpretation! you do not know the procedure of getting those data and it would be problematic for interpretation. The point is that those datasets are from different stages but also from different sources! so we can not make decision here. The other point is about using prcomp function. I've never used it but usually in machine learning tools the setting is like samples in rows and features in columns so I recommend to have a quick look at it again. – Kasra Manshaei Jun 18 '15 at 18:34
• I transposed the matrix, thus rows are for samples and columns for genes. The code I use is in the 2nd update. Thanks! – Dejian Jun 19 '15 at 2:28