<p>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.</p>

CC BY-SA 4.0
2021-10-19T08:31:27.763Z
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?
Justin1609
CC BY-SA 4.0
2021-10-19T23:38:25.300Z
IF you are already using DESeq, why not just do what they do in the vignette?
swbarnes2
CC BY-SA 4.0
2021-10-20T11:13:41.283Z
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.
Justin1609
CC BY-SA 4.0
2021-10-18T18:00:01.200Z
2021-10-18T18:00:01.200Z
swbarnes2
<blockquote>
<p>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.</p>
</blockquote>
<p>You can merge all your datasets for common genes in all datasets and perform <code>PCA</code>.</p>
<blockquote>
<p>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.</p>
</blockquote>
<p>Your choice of transposing data will depend on whether you want to show variability between genes or samples. However you can get results for both individuals (genes) as well as variables (samples), using the <code>prcomp()</code>, <code>factoMineR</code> and <code>factoexta</code> packages in 'R'. Usually I like to use <code>$rotation</code> values ​​to plot the variability between samples using <code>ggplot2</code> package. You can use this <a href="http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/118-principal-component-analysis-in-r-prcomp-vs-princomp/#:%7E:text=The%20function%20princomp()%20uses,preferred%20compared%20to%20princomp()." rel="nofollow noreferrer">PCA Tutorial</a> to learn more.</p>
<blockquote>
<p>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.</p>
</blockquote>
<p>This is true as the datasets are from different platforms, there may be some noise or batch effects. However, you can plot the 'PCA' and decide for yourself whether it is good or not.</p>

CC BY-SA 4.0
2022-08-18T07:08:26.897Z
2022-08-18T07:08:26.897Z
Nitesh Shriwash
<p>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.</p>
<p>The 3 different datasets I performed DESeq2 analysis on are as follows:</p>
<p>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.</p>
<p>Dataset 2: This was sequenced by Illumina MiSeq and paired-end reads were generated, but this dataset had two parts to it.</p>
<p>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.</p>
<p>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.</p>
<p>Dataset 3: This was sequenced by Illumina MiSeq and paired-end reads were generated, but this dataset had two parts to it.</p>
<p>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.</p>
<p>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.</p>
<p>I wrote the results out of the normalized counts per sample combined with the DESeq2 results as follows:</p>
<pre><code>res &lt;- results(dds)
table(res$padj&lt;0.05)
res &lt;- res[order(res$padj), ]
resdata &lt;- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)), by=&quot;row.names&quot;, sort=FALSE)
names(resdata)[1] &lt;- &quot;Gene&quot;
head(resdata)
write.csv(resdata, file=&quot;diffexpr-results.csv&quot;)
</code></pre>
<p>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.</p>
<p>An example of the setup of the excel files I have, without the other info that you get with DESeq2 analysis e.g. columns <code>baseMean</code>, <code>log2FoldChange</code>, <code>lfcSE</code>, <code>stat</code>, <code>pvalue</code> and <code>padj</code>, is as follows:</p>
<div class="s-table-container">
<table class="s-table">
<thead>
<tr>
<th>Gene_ID</th>
<th>Sample1_PC</th>
<th>Sample2_PC</th>
<th>Sample3_PC</th>
<th>Sample4_PC</th>
<th>Sample5_Mixed</th>
<th>Sample6_Mixed</th>
<th>Sample7_Mixed</th>
</tr>
</thead>
<tbody>
<tr>
<td>YL09980</td>
<td>1240,016</td>
<td>759,683</td>
<td>985,322</td>
<td>1025,855</td>
<td>5002,654</td>
<td>4200,963</td>
<td>4002,855</td>
</tr>
<tr>
<td>YLW9980</td>
<td>589,456</td>
<td>452,236</td>
<td>630,228</td>
<td>523,821</td>
<td>2325,556</td>
<td>2566,987</td>
<td>1800,633</td>
</tr>
</tbody>
</table>
</div>
<p>So these tables are just the normalized counts of the genes per sample.</p>
<p>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.</p>
<p>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.</p>
<p>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 &quot;tidyverse&quot; in an effort to try and do this, but I am on a deadline with this.</p>
<p>Edit: So I did the PCA plot using the following function:</p>
<pre><code>data_total.pca &lt;- prcomp(data_total[,c(2:32)], center = TRUE, scale. = TRUE)
</code></pre>
<p>and this was the output I got:</p>
<pre><code>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
</code></pre>
<p>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:</p>
<pre><code>data_total_t &lt;- transpose(data_total)
rownames(data_total_t) &lt;- colnames(data_total)
colnames(data_total_t) &lt;- rownames(data_total)
</code></pre>
<p>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.</p>
<p>Any help on how I can do this please?</p>

<r><rna-seq><deseq2><pca><meta-analysis>
Help with performing PCA analysis on data from 3 different datasets of normalized counts data for RNA-seq experiment
# Help with performing PCA analysis on data from 3 different datasets of normalized counts data for RNA-seq experiment

(Asked by: username_0 on 2021-10-18)

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?

**Tags:** r, rna-seq, deseq2, pca, meta-analysis

---

## 2 answers

### Answer by username_1 on 2021-10-18. Score: 1

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.

> **Comments:**
> **username_0**: 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?  
> **username_1**: IF you are already using DESeq, why not just do what they do in the vignette?  
> **username_0**: 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.  

---

### Answer by username_2 on 2022-08-18. Score: 1

> 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.

You can merge all your datasets for common genes in all datasets and perform `PCA`.

> 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.

Your choice of transposing data will depend on whether you want to show variability between genes or samples. However you can get results for both individuals (genes) as well as variables (samples), using the `prcomp()`, `factoMineR` and `factoexta` packages in 'R'. Usually I like to use `$rotation` values ​​to plot the variability between samples using `ggplot2` package. You can use this [PCA Tutorial](http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/118-principal-component-analysis-in-r-prcomp-vs-princomp/#:%7E:text=The%20function%20princomp()%20uses,preferred%20compared%20to%20princomp().) to learn more.

> 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.

This is true as the datasets are from different platforms, there may be some noise or batch effects. However, you can plot the 'PCA' and decide for yourself whether it is good or not.


---

