Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-13T04:44:28.170Z
Answers.Text: Im using prcomp function would that work as you mentioned princomp ?
Answers.UserDisplayName: kcm
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-13T10:47:06.230Z
Answers.Text: `Z%*%pZ$rotation` should be the same as `pZ$x` (that’s in fact how `pZ$x` is calculated by `prcomp`).
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-13T10:50:15.167Z
Answers.Text: 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.
Answers.UserDisplayName: llrs
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-13T11:41:44.650Z
Answers.Text: @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
Answers.UserDisplayName: kcm
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-13T11:53:46.963Z
Answers.Text: @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
Answers.UserDisplayName: llrs
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-13T13:15:40.477Z
Answers.Text: okay i will try to find it out
Answers.UserDisplayName: kcm
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-13T14:39:23.897Z
Answers.Text: @Llopis i tried that i get error like this Error in Z %*% pZ$rotation : 
    requires numeric/complex matrix/vector arguments
Answers.UserDisplayName: kcm
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-14T07:49:39.713Z
Answers.Text: @krushnachChandra See the new edit. Could you explain what code are you using and what is Z?
Answers.UserDisplayName: llrs
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-12-12T16:50:18.493Z
Answers.LastActivityDate: 2017-12-14T07:48:48.733Z
Answers.LastEditDate: 2017-12-14T07:48:48.733Z
Answers.UserDisplayName: llrs
Body: <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>
  
Tags: <r><rna-seq>
Title: Extracting features or gene from PCA after calculating PCA for downstream analysis
ThreadText: # 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](https://www.nature.com/articles/nbt.2967#s1):
  
  > 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](http://stat.ethz.ch/R-manual/R-devel/library/stats/html/prcomp.html) 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**: @username_2 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?  
  
  ---
  
  