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?
