# Help with performing PCA analysis on data from 3 different datasets of normalized counts data for RNA-seq experiment

So, I have limited knowledge of R but I need to do a PCA analysis of 3 different datasets of gene expression as a result of combined growth or mono-culture growth.

The 3 different datasets I performed DESeq2 analysis on are as follows:

Dataset 1: This was sequenced by Ion Torrent and it was single-end reads. Samples that were compared were biological quadruplicates of yeast species A grown as pure culture versus species A grown in direct contact with yeast species B. This was all done for 1 time point @ 24h.

Dataset 2: This was sequenced by Illumina MiSeq and paired-end reads were generated, but this dataset had two parts to it.

Samples that were compared were biological quadruplicates of yeast species A grown as pure culture versus species A grown in direct contact with yeast species B. This dataset contains 2 time points @ 2h and 24h.

For the other part of this, samples that were compared were biological quadruplicates of yeast species A grown as pure culture versus species A grown via indirect (metabolic exchange possible) contact with yeast species B. This dataset contains 2 time points @ 2h and 24h.

Dataset 3: This was sequenced by Illumina MiSeq and paired-end reads were generated, but this dataset had two parts to it.

Here, two biological repeats of samples of yeast species A grown as pure culture versus two biological repeats of species A grown in direct contact with yeast species B, under AEROBIC conditions were performed. This dataset contains 1 time point @ 48h.

Then one biological repeat of sample of yeast species A grown as pure culture versus two biological repeats of species A grown in direct contact with yeast species B, under ANAEROBIC conditions were performed. This dataset contains 1 time point @ 48h.

I wrote the results out of the normalized counts per sample combined with the DESeq2 results as follows:

res <- results(dds) table(res$padj<0.05) res <- res[order(res$padj), ] resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)), by="row.names", sort=FALSE) names(resdata)[1] <- "Gene" head(resdata) write.csv(resdata, file="diffexpr-results.csv")

So the issue I am having is, because these are different datasets and conditions differ across them there are different numbers of genes expressed in each dataset that are not present across all datasets. So I need to combine all the datasets to only compare the genes found in all the datasets mentioned and then perform a PCA analysis on those genes only.

An example of the setup of the excel files I have, without the other info that you get with DESeq2 analysis e.g. columns "baseMean", "log2FoldChange", "lfcSE", "stat", "pvalue" and "padj", is as follows:

Gene_ID Sample1_PC Sample2_PC Sample3_PC Sample4_PC Sample5_Mixed Sample6_Mixed Sample7_Mixed
YL09980 1240,016 759,683 985,322 1025,855 5002,654 4200,963 4002,855
YLW9980 589,456 452,236 630,228 523,821 2325,556 2566,987 1800,633

So these tables are just the normalized counts of the genes per sample.

I was told that the data in the table has to be transposed basically, so that the genes are the variables and the samples are the rows. I am just not sure how to do this in R.

I was also told there is a good chance that the noise resulting from the use of different platforms would likely cause these datasets to separate from each other based on that, making this analysis meaningless or not allowing for comparison of gene expression that is conserved across the datasets as a result of mixing the two species during growth.

However, I just need to do this analysis and find that out as soon as possible, as this is a bottleneck to my progress at this point. Please any assistance would be greatly appreciated. I am currently teaching myself more about "tidyverse" in an effort to try and do this, but I am on a deadline with this.

Edit: So I did the PCA plot using the following function:

data_total.pca <- prcomp(data_total[,c(2:32)], center = TRUE, scale. = TRUE)

and this was the output I got:

Importance of components:
PC1    PC2    PC3    PC4     PC5     PC6     PC7     PC8     PC9    PC10    PC11    PC12    PC13    PC14
Standard deviation     4.2180 2.5400 1.9188 1.2866 0.65393 0.47799 0.40826 0.38879 0.35585 0.34915 0.27036 0.18549 0.14335 0.11830
Proportion of Variance 0.5739 0.2081 0.1188 0.0534 0.01379 0.00737 0.00538 0.00488 0.00408 0.00393 0.00236 0.00111 0.00066 0.00045
Cumulative Proportion  0.5739 0.7820 0.9008 0.9542 0.96800 0.97537 0.98074 0.98562 0.98971 0.99364 0.99600 0.99711 0.99777 0.99822
PC15    PC16    PC17    PC18    PC19    PC20    PC21    PC22    PC23    PC24    PC25    PC26    PC27    PC28
Standard deviation     0.10153 0.09239 0.08407 0.07765 0.07536 0.05978 0.04883 0.04556 0.04406 0.03960 0.03764 0.03378 0.03167 0.02897
Proportion of Variance 0.00033 0.00028 0.00023 0.00019 0.00018 0.00012 0.00008 0.00007 0.00006 0.00005 0.00005 0.00004 0.00003 0.00003
Cumulative Proportion  0.99855 0.99883 0.99906 0.99925 0.99943 0.99955 0.99963 0.99969 0.99976 0.99981 0.99985 0.99989 0.99992 0.99995
PC29    PC30    PC31
Standard deviation     0.02695 0.02126 0.02085
Proportion of Variance 0.00002 0.00001 0.00001
Cumulative Proportion  0.99997 0.99999 1.00000


So according to the tutorial I am following a general rule of thumb is that the first 3 PCs should explain at least 70% of the variance. The issue is my PCA plot looks weird all of the points are squished up at the x = 0; y = 0 axes of my PCA plot. So I think I need to transpose the data so that the genes are the columns (variables) and the samples are the rows (objects). So I tried this:

data_total_t <- transpose(data_total) rownames(data_total_t) <- colnames(data_total) colnames(data_total_t) <- rownames(data_total)

And this worked, but now I need to exclude row 1 from the function I used to plot the PCA as the PCA cannot have character (gene name) info in it.

Any help on how I can do this please?

## 1 Answer

I think you are making this harder than you have to. Why can't you follow the code used in the vignette for making PCA plots? It doesn't even involve getting all the way to the DE genes.

• Hi there @swbarnes2, I found a tutorial online to follow but now I am uncertain of the plot I am getting in terms of if it is correct or not. All of my data is now squished into the x = 0; y = 0 point of the PCA axes. Is this because I should rather be plotting the genes as columns and the sample names should make up the rows? Oct 19 '21 at 8:31
• IF you are already using DESeq, why not just do what they do in the vignette? Oct 19 '21 at 23:38
• I will give it a look thanks @swbarnes2, I have also generated PCA plots already now for the data but I am not sure if they are very meaningful or if the data can actually be compared in this way. Oct 20 '21 at 11:13