<p>The <a href="http://stat.ethz.ch/R-manual/R-devel/library/stats/html/prcomp.html" rel="nofollow noreferrer">help page</a> explains that there is a component which:</p>

<blockquote>
  <p>rotation<br>
  the matrix of variable loadings (i.e., a matrix whose columns contain the eigenvectors). The function princomp returns this in the element loadings.</p>
</blockquote>

<p>So using using the rotation one can use it to calculate the importance of each feature in the input:</p>

<pre><code>C &lt;- chol(S &lt;- toeplitz(.9 ^ (0:31))) # Cov.matrix and its root
set.seed(17)
X &lt;- matrix(rnorm(32000), 1000, 32)
Z &lt;- X %*% C  ## ==&gt;  cov(Z) ~=  C'C = S
pZ &lt;- prcomp(Z)
genesImportance &lt;- Z %*% pZ$rotation
genesImportance2 &lt;- pZ$x
all.equal(genesImportance, genesImportance2)
## [1] "Mean relative difference: 0.03680457"
</code></pre>

<p>This new matrix has the loading for each gene for each component. So you can now select them according to how much do they participate in each component or whatever you want. Or you can use it to predict the position of samples you didn't set in the initial data. </p>

CC BY-SA 3.0
2017-12-13T04:44:28.170Z
Im using prcomp function would that work as you mentioned princomp ?
kcm
CC BY-SA 3.0
2017-12-13T10:47:06.230Z
`Z%*%pZ$rotation` should be the same as `pZ$x` (that’s in fact how `pZ$x` is calculated by `prcomp`).
Konrad Rudolph
CC BY-SA 3.0
2017-12-13T10:50:15.167Z
Totally true, I didn't went further in the help page, but with `rotation` one can predict the position of new samples not included in the original data. I'll edit it.
llrs
CC BY-SA 3.0
2017-12-13T11:41:44.650Z
@KonradRudolph I'm still not getting where is this coming from  _italic_ pZ$x   im little ignorant as I see examples and do it so i have to do a back referencing  to in R methods and tutorials , is it in my data frame as well ?with context to my data frame are you suggesting this  _italic_"pca$x"  and what is that  Z  actually
kcm
CC BY-SA 3.0
2017-12-13T11:53:46.963Z
@krushnachChandra Z is the initial data for the PCA. pZ is the output of `prcomp`. See the edit, I might not explain it clearly enough
llrs
CC BY-SA 3.0
2017-12-13T13:15:40.477Z
okay i will try to find it out
kcm
CC BY-SA 3.0
2017-12-13T14:39:23.897Z
@Llopis i tried that i get error like this Error in Z %*% pZ$rotation : 
  requires numeric/complex matrix/vector arguments
kcm
CC BY-SA 3.0
2017-12-14T07:49:39.713Z
@krushnachChandra See the new edit. Could you explain what code are you using and what is Z?
llrs
CC BY-SA 3.0
2017-12-12T16:50:18.493Z
2017-12-14T07:48:48.733Z
2017-12-14T07:48:48.733Z
llrs
<p>The quote below is from <a href="https://www.nature.com/articles/nbt.2967#s1" rel="nofollow noreferrer">this paper</a>:</p>

<blockquote>
  <p>We performed principal component analysis (PCA) of low-coverage sequencing data to identify genes explaining variation across cells. PCA separated the cells into groups corresponding to the source populations (Fig. 2c and Supplementary Figs. 3–5), and the genes distinguishing each group reflected the biological properties of the cell types (Supplementary Fig. 5 and Supplementary Table 3). PCA of low- and high-coverage sequencing data revealed a remarkably similar graphical distribution of analyzed cells, and the majority (78%) of the top 500 genes determined by PCA were shared between PCA performed on low- and high-coverage data (Supplementary Figs. 4 and 6 and Supplementary Table 4). </p>
</blockquote>

<p>How are they performing PCA and then finding out the components which contribute to the strongest PC and then taking out the genes for downstream analysis? Let's say there are genes involved in PC1 and PC2 which define a certain cell type and that distinguish it. I think for a biologist, like me, I just want to get the genes or gene list that are getting involved in determining lineage or cell type. How can I do this?</p>

<p>As a to test what they might be doing i did this </p>

<pre><code>test&lt;- read.csv('NON_CODING.csv',header = T,row.names = 1)

# preform PCA
pca = prcomp(t(test), center=TRUE, scale=TRUE)
</code></pre>

<p>then after calculating PCA i do see this in my console when I type pca</p>

<pre><code>&gt; pca
Standard deviations (1, .., p=16):
 [1] 3.826851e+01 2.080405e+01 1.568739e+01 1.349256e+01 1.119348e+01 9.980547e+00 8.365034e+00
 [8] 8.098841e+00 7.519507e+00 6.184505e+00 5.880260e+00 5.139609e+00 4.851091e+00 4.335606e+00
[15] 3.870918e+00 3.173450e-14

Rotation (n x k) = (2899 x 16):
                                       PC1           PC2           PC3           PC4
5S_rRNA                      -1.090574e-02 -2.412665e-03  1.637689e-02 -3.603865e-02
AB019441.29                  -1.928250e-02  1.821083e-03 -9.713724e-03 -1.978737e-03
ABBA01017803.1               -1.823266e-02 -3.727144e-03 -9.790131e-04 -3.937062e-02
ABC14-1080714F14.1            2.438024e-02  3.816019e-03  1.657784e-04 -7.141480e-03
ABC7-481722F1.1               2.467403e-02  2.873432e-03  3.781153e-03  3.790824e-03
AC000036.4                   -1.066777e-02  2.838305e-02  2.642673e-02 -7.998608e-03
AC000089.3                   -1.249026e-02 -7.870864e-04 -2.569602e-02 -3.165073e-03
AC000120.7                   -1.313696e-02  5.923245e-03  1.633255e-02  4.077602e-02
AC000123.4                    2.415996e-02  1.029732e-02  8.930792e-03 -1.118660e-02
AC000403.4                    8.692835e-03 -3.060170e-02  3.746839e-02 -6.968971e-03
AC002064.4                   -1.799590e-02  2.147160e-02 -2.303742e-02  7.647937e-03
</code></pre>

<p>Now how do i find which are the list or set of genes that gives me PC1 or PC2 so on </p>

<p>Any suggestion or help how to get the genes that make most difference between two major component that would be very helpful </p>

<r><rna-seq>
Extracting features or gene from PCA after calculating PCA for downstream analysis
# Extracting features or gene from PCA after calculating PCA for downstream analysis

(Asked by: username_0 on 2017-12-12)

The quote below is from this paper:

We performed principal component analysis (PCA) of low-coverage sequencing data to identify genes explaining variation across cells. PCA separated the cells into groups corresponding to the source populations (Fig. 2c and Supplementary Figs. 3–5), and the genes distinguishing each group reflected the biological properties of the cell types (Supplementary Fig. 5 and Supplementary Table 3). PCA of low- and high-coverage sequencing data revealed a remarkably similar graphical distribution of analyzed cells, and the majority (78%) of the top 500 genes determined by PCA were shared between PCA performed on low- and high-coverage data (Supplementary Figs. 4 and 6 and Supplementary Table 4). 

How are they performing PCA and then finding out the components which contribute to the strongest PC and then taking out the genes for downstream analysis? Let's say there are genes involved in PC1 and PC2 which define a certain cell type and that distinguish it. I think for a biologist, like me, I just want to get the genes or gene list that are getting involved in determining lineage or cell type. How can I do this?
As a to test what they might be doing i did this 
test<- read.csv('NON_CODING.csv',header = T,row.names = 1)

# preform PCA
pca = prcomp(t(test), center=TRUE, scale=TRUE)

then after calculating PCA i do see this in my console when I type pca
> pca
Standard deviations (1, .., p=16):
 [1] 3.826851e+01 2.080405e+01 1.568739e+01 1.349256e+01 1.119348e+01 9.980547e+00 8.365034e+00
 [8] 8.098841e+00 7.519507e+00 6.184505e+00 5.880260e+00 5.139609e+00 4.851091e+00 4.335606e+00
[15] 3.870918e+00 3.173450e-14

Rotation (n x k) = (2899 x 16):
                                       PC1           PC2           PC3           PC4
5S_rRNA                      -1.090574e-02 -2.412665e-03  1.637689e-02 -3.603865e-02
AB019441.29                  -1.928250e-02  1.821083e-03 -9.713724e-03 -1.978737e-03
ABBA01017803.1               -1.823266e-02 -3.727144e-03 -9.790131e-04 -3.937062e-02
ABC14-1080714F14.1            2.438024e-02  3.816019e-03  1.657784e-04 -7.141480e-03
ABC7-481722F1.1               2.467403e-02  2.873432e-03  3.781153e-03  3.790824e-03
AC000036.4                   -1.066777e-02  2.838305e-02  2.642673e-02 -7.998608e-03
AC000089.3                   -1.249026e-02 -7.870864e-04 -2.569602e-02 -3.165073e-03
AC000120.7                   -1.313696e-02  5.923245e-03  1.633255e-02  4.077602e-02
AC000123.4                    2.415996e-02  1.029732e-02  8.930792e-03 -1.118660e-02
AC000403.4                    8.692835e-03 -3.060170e-02  3.746839e-02 -6.968971e-03
AC002064.4                   -1.799590e-02  2.147160e-02 -2.303742e-02  7.647937e-03

Now how do i find which are the list or set of genes that gives me PC1 or PC2 so on 
Any suggestion or help how to get the genes that make most difference between two major component that would be very helpful

**Tags:** r, rna-seq

---

## 1 answers

### Answer by username_1 on 2017-12-12. Score: 7 | Accepted answer

The help page explains that there is a component which:

rotation
  the matrix of variable loadings (i.e., a matrix whose columns contain the eigenvectors). The function princomp returns this in the element loadings.

So using using the rotation one can use it to calculate the importance of each feature in the input:
C <- chol(S <- toeplitz(.9 ^ (0:31))) # Cov.matrix and its root
set.seed(17)
X <- matrix(rnorm(32000), 1000, 32)
Z <- X %*% C  ## ==>  cov(Z) ~=  C'C = S
pZ <- prcomp(Z)
genesImportance <- Z %*% pZ$rotation
genesImportance2 <- pZ$x
all.equal(genesImportance, genesImportance2)
## [1] "Mean relative difference: 0.03680457"

This new matrix has the loading for each gene for each component. So you can now select them according to how much do they participate in each component or whatever you want. Or you can use it to predict the position of samples you didn't set in the initial data.

> **Comments:**
> **username_0**: Im using prcomp function would that work as you mentioned princomp ?  
> **username_2**: `Z%*%pZ$rotation` should be the same as `pZ$x` (that’s in fact how `pZ$x` is calculated by `prcomp`).  
> **username_1**: Totally true, I didn't went further in the help page, but with `rotation` one can predict the position of new samples not included in the original data. I'll edit it.  
> **username_0**: @KonradRudolph I'm still not getting where is this coming from  _italic_ pZ$x   im little ignorant as I see examples and do it so i have to do a back referencing  to in R methods and tutorials , is it in my data frame as well ?with context to my data frame are you suggesting this  _italic_"pca$x"  and what is that  Z  actually  
> **username_1**: @krushnachChandra Z is the initial data for the PCA. pZ is the output of `prcomp`. See the edit, I might not explain it clearly enough  
> **username_0**: okay i will try to find it out  
> **username_0**: @Llopis i tried that i get error like this Error in Z %*% pZ$rotation :    requires numeric/complex matrix/vector arguments  
> **username_1**: @krushnachChandra See the new edit. Could you explain what code are you using and what is Z?  

---

