I think there are two different definitions of "working":
1) Checking if there is a bug in the code - The code doesn't work as described (in terms of carrying out the individual steps that it is supposed to, without crashing or having unexpected errors like "NA" values). This can be identified and fixed, and I think that is the focus of the responses above (in other words, if you know a package as a function, check that you are using the correct commands for the function that you want to use.
2) Assessment of the normalization/adjustment (concerns about over-fitting, etc.) - I think this is harder to define in a way that applies equally well to all projects. However, I think checking expression before and after the extra step may be important.
This can be a little tricky in situations where you include co-variates that affect clustering. In other words, if the concern is that the normalization makes pre-defined groups artificially cluster together more than can be robustly observed in other studies (which I would consider "over-fitting" and/or "over-correction"), then the overall clustering in a PCA plot may be difficult to assess.
If I had to give some possibly useful recommendations (that may or may not help for a particular project), this is what I can think of:
a) Visually check clustering of replicates in the heatmap. I think this matters more for the p-value calculation (checking that good concordance among replicates was indeed defined), but I admittedly use multi-variate p-value models (tested with different methods for each project) more often than a different normalization (or, just for the heatmap, I might do something like center by batch/tissue/cell prior to setting the overall mean to 0 and standard deviation of 1).
b) Visually inspect some individual candidate genes (possibly positive controls from your experiment). For example, if you are correcting an effect, make sure that the groups are randomized and you can observe the same trend if you separate groups by the variable you are trying to correct. For example, if you have 3 batches and 2 groups, I would check i) is one sample from each group in each batch and ii) assuming your candidate is up-regulated for each treatment, can you confirm that expression is increased for each treatment relative to each control for each separate batch? The adjustment may make this easier to see, but I think you still need the original expression (possibly in a supplemental figure) do confirm you haven't done something weird (like change the relative trend in expression within batches).
If there are other ideas of how to address #2, I would love to hear them!
Please edit your original post and use these instructions to add images inline: How to add images to a Biostars post
Thanks for showing how to do it!
Strange. Which package is the "PCAplot" function from? Why is there a circle? and how much of the variation does the PCA 1 and 2 explain? I assume the "MDSplot" function is edgeR's own right?
How did you use the "removeBatchEffect"? What input did you give?
Hi again! Thanks for answering. I constructed the PCAplot from a short pipeline I found online in a blog. . Maybe someone can suggest a better way to create the PCA plot, I would be very happy if someone could share it with me.
In the before removing batch effect PCAplot the variation explain with PCA1 and PCA2 are:
In the after removing batch effect PCAplot the variation explain with PCA1 and PCA2 is:
On the other way, yes MDSplot function is from edgeR. And now I am going to write the pipeline used for edgeR normalization and plots generation to see if I am correct in the followed steps, and on the same way answerd your question about "removeBatchEffect" use:
I am going to try to add the images in another way following the guides "genomax" mentions.
If you need more information about my followed steps or something more tell me. Any appreciation, recomendation or help will be welcomed! Thansk,
Iraia
For the PCA bi-plots, you should be plotting the
xcomponent, not therotationcomponent. So, it should beTMM.pca$x, notTMM.pca$rotation. DESeq2's PCA function and my own PCAtools package plot thexcomponent.To read more:
This is not true. In RNA-seq and expression studies, control and treatment do not necessarily have to segregate on a bi-plot.
Perhaps you should explain more about the laboratory preparatory work that you performed: How was RNA isolation and extraction performed?; which kits were used?; how was sequencing performed and to which depth?
I find that many people worry too much about batch effects. The key to tackling batch is in your experimental design, which you need to describe in detail here.
If the poster didn't transpose his count matrix before doing prcomp, then
$rotationmight be what he wants.Not sure about that. It looks like
prcomp()was used. The variable names in prcomp are confusing.rotationis just the component loadings, whereasxis the original data (possibly centered and scaled) multiplied by the loadings, which is what people usually plot. I have very rarely seen people make bi-plots of component loadings.Let's do an example. Create random data-matrix with genes as rows and samples as columns:
Perform PCA:
Transpose the input data (
t(a)) to see what happens:So, iraia.munoa, if your input data has genes as rows and samples as columns, you need to run prcomp like this:
Then plot the
xvariable ofpcaThe important thing is that you use the x-component and that you use log transformed data - and I would also suggest that you use scale = TRUE in prcomp.
so you are agree with Kevin Blighe in doing the transposition? I would try it right know and tell you all about the result. Thanks everybody!
Hi again! I think I have manage to follow your recomendations. I write them right here to see if I am correct (after TMM normalization, log transformation and removing batch effect):
Defining prcomp with scale=TRUE, and transposing to use the x values:
Plotting the results:
Conclusions:
Is this correct? I mean, I am working with this data for my phd project, would it be correct to use it in this way? or someone recomends me to follow other steps? And something more important, after this steps and with that PCAplots would you follow the proccess and do the DE, believing in the results? I would expect that after removing batch effect, all control samples would be grouped in a corner and all treatment samples in the other corner, I mean with a good separation, and I only see this clasification if I draw a line from the upper left corner to the bottom right corner, and I don't know if it is enough to explain correctly sample distribution. Maybe I have an incorrect understanding of how to understand the PCA plot significance in analyzing RNA-Seq data, so any explanaition would be welcomed!
Thanks in advance
The PCA bi-plot looks as expected now, i.e., using the
xvalue. The batch correction seems to have clearly identified a batch effect and adjusted for it, too. You are still not showing the percent explained variation on the axes on the bi-plots, though (?).Oh! sorry, with my coments I did not want you to say if everything I have done is correct or not. I only want recomendation of which steps to follow or detection if something I have shown here is correct or not. I mean, we don't have bioinformaticians in our research group, so I am trying to understand and learn everything by myself (as I think a lot of researchers do with sequencing analysis) and trying to do the job in the better possible way. So I only wanted to know if after the batch effect adjust it is a good idea to follow the DE analysis, I mean if after looking at the second PCA plot, experts in the field would say "Well, now the sample clusters are better, and it is worthy to follow with the DE analysis", or on the other hand "Well, the adjust is correct, but maybe there are another steps worthy to follow before DE analysis".
Even so, I will try your idea of doing the same DE analysis, with data before adjustment and after adjustment and compare both DEG lists (and of course the comparisons of all the plots you mentioned above).
Answering your questions, the variation percent is:
68% and 77% of cumulative variance is a small proportion?
And finally, my goal of using removeBatchEffect() and doing before and after MDS and PCAplots was to see if after adjusting the batch effect, the clusterization of samples was better. But as I have read, I can't use the data after removing the batch effect for DEG analysis right? I need to do the edgeR protocol (or DESeq) but taken into account the "batch" in the matrix design right?
Thanks again
I see what you mean. For your differential expression analysis, just include
batchin the design formula and use the un-adjusted normalised counts. The batch effect looks consistent in your study, so, the regression model should be easy to fit by EdgeR's internal engine.Please take time to read responses here: https://support.bioconductor.org/p/87718/
Thanks again for all your answers, you have helped me a lot! If someone more wants to add some other points of view they are also welcome to writ in this post.
What about the variance values? Are they enough to focud on PC1 vs PC2 or do I need to change to PC3 or something like that?
Regarding the explained variation, nothing seems out of the ordinary to me.
Just regarding these lines:
This should be:
..or:
Take a look here: https://support.bioconductor.org/p/116614/#116621
I didn't know about "prior.count". I have read what Gordon says in the post you linked, and I have used it in my samples, and the MDSplot looks the same. I don't understand really well why I should use it. It is to avoid that 0 counts become Inf when doing the log transformation? And how do I know which number to chose there? why 2? Thanks,
To follow with RamRS answer, I didn't know the need to transpose the count matrix, and I simply followed a pipeline posted in my first reply. I though that it was correct to use it as a global value for each variables, as they explain the rotations value as coefficients of the linear combinations of the continuous variables.