RunPCA step by step

[hyeonjin-cho]

Hi, I'm running PCA step-by-step with Seurat code published in GitHub. However, with smaller dataset (1000 cells and 1290 genes) it worked fine, which means my step-by-step PCA matched the one automatically produced with Seurat. However, when I increased the number of genes (1000 cells and 11K genes), it didn't match. Does anyone have any clue which part of the PCA step is different when you use larger dataset?

Thank you!

[timoast]

I don't understand what you mean exactly. Could you provide the code you're running, and explain the expected result?

[hyeonjin-cho]

So I have a count matrix with 1000 cells and 1290 genes, and I'm running

simData.seurat <- RunPCA(simData.seurat, approx=F)

where simData.seurat is my Seurat object and I turned off the approx parameter to use prcomp() instead of irlba().

Then, I used prcomp() manually to compare the results.

Specifically, I used:

object <- as.matrix(simData.seurat[["RNA"]]@scale.data)
npcs <- min(50, nrow(object)-1)
pca.results <- prcomp(t(object), rank.=npcs)

And got the same result as RunPCA (except for the +/- signs).
To be more specific, I get:

> simData.seurat[["pca"]]@cell.embeddings[1:10,1:2]
                PC_1       PC_2
cell_1_0   2.0204074 -2.5447096
cell_2_0   8.2292489  2.0277615
cell_3_0   0.7066860 -2.2447125
cell_4_0   4.0846358 -1.3404078
cell_5_0  -0.1353387 -2.6935367
cell_6_0   2.3407943  0.2888769
cell_7_0   1.3386422  1.3447453
cell_8_0   0.7712575 -2.7390502
cell_9_0   1.0018514 -2.3214214
cell_10_0 -6.6129032  3.0455656
> pca.results$x[1:10,1:2]
                 PC1        PC2
cell_1_0   2.0204074  2.5447096
cell_2_0   8.2292489 -2.0277615
cell_3_0   0.7066860  2.2447125
cell_4_0   4.0846358  1.3404078
cell_5_0  -0.1353387  2.6935367
cell_6_0   2.3407943 -0.2888769
cell_7_0   1.3386422 -1.3447453
cell_8_0   0.7712575  2.7390502
cell_9_0   1.0018514  2.3214214
cell_10_0 -6.6129032 -3.0455656

However, when I used 11K genes instead of 1290, and run both codes above (RunPCA and prcomp), I get different results
Again, simData.seurat is my Seurat object, but this time with 1000 cells and 11K genes

> simData.seurat[["pca"]]@cell.embeddings[1:10,1:2]
                PC_1       PC_2
cell_1_0    2.612392 -5.8370546
cell_2_0    6.897207 -0.2805351
cell_3_0   -1.030449  1.7576380
cell_4_0    7.044210  0.6586270
cell_5_0    1.661446  3.2977973
cell_6_0    3.746600  3.9372971
cell_7_0   -5.007205 -4.4025976
cell_8_0    4.505317  6.0685604
cell_9_0   -1.343331 -3.1441043
cell_10_0 -11.287677 -5.0507961
> pca.results$x[1:10,1:2]
                  PC1        PC2
cell_1_0    0.1847801  -7.106062
cell_2_0   14.2456962   6.314401
cell_3_0   -3.7191144  -2.394562
cell_4_0    1.7788215 -13.924483
cell_5_0   -3.1858047  -3.604660
cell_6_0   -0.3557700  -6.090833
cell_7_0   -1.3147928   3.908382
cell_8_0    0.1687961  -7.194880
cell_9_0   -3.4474973  -4.444379
cell_10_0 -10.1837476  10.842258

So I'm just wondering if any steps in Seurat is different when using larger dataset.

Thank you!

[timoast]

RunPCA by default uses the variable features, so the two are likely using different feature sets as input here

[hyeonjin-cho]

Yes it now matches! Thank you Tim