PCA with DESeq2

Installation

Install using BiocManager::install().

if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("DESeq2")

Package version.

packageVersion("DESeq2")

[1] '1.46.0'

Count table

https://zenodo.org/records/13970886.

my_url <- 'https://zenodo.org/records/13970886/files/rsem.merged.gene_counts.tsv?download=1'
my_file <- 'rsem.merged.gene_counts.tsv'

if(file.exists(my_file) == FALSE){
  download.file(url = my_url, destfile = my_file)
}

gene_counts <- read_tsv("rsem.merged.gene_counts.tsv", show_col_types = FALSE)
head(gene_counts)

# A tibble: 6 × 10
  gene_id   `transcript_id(s)` ERR160122 ERR160123 ERR160124 ERR164473 ERR164550
  <chr>     <chr>                  <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
1 ENSG0000… ENST00000373020,E…        2         6         5      374       1637 
2 ENSG0000… ENST00000373031,E…       19        40        28        0          1 
3 ENSG0000… ENST00000371582,E…      268.      274.      429.     489        637 
4 ENSG0000… ENST00000367770,E…      360.      449.      566.     363.       606.
5 ENSG0000… ENST00000286031,E…      156.      185.      265.      85.4      312.
6 ENSG0000… ENST00000374003,E…       24        23        40     1181        423 
# ℹ 3 more variables: ERR164551 <dbl>, ERR164552 <dbl>, ERR164554 <dbl>

Metadata.

tibble::tribble(
  ~sample, ~run_id, ~group,
  "C2_norm", "ERR160122", "normal",
  "C3_norm", "ERR160123", "normal",
  "C5_norm", "ERR160124", "normal",
  "C1_norm", "ERR164473", "normal",
  "C1_cancer", "ERR164550", "cancer",
  "C2_cancer", "ERR164551", "cancer",
  "C3_cancer", "ERR164552", "cancer",
  "C5_cancer", "ERR164554", "cancer"
) -> my_metadata

my_metadata$group <- factor(my_metadata$group, levels = c('normal', 'cancer'))

Matrix.

gene_counts |>
  dplyr::select(starts_with("ERR")) |>
  mutate(across(everything(), as.integer)) |>
  as.matrix() -> gene_counts_mat

row.names(gene_counts_mat) <- gene_counts$gene_id

idx <- match(colnames(gene_counts_mat), my_metadata$run_id)
colnames(gene_counts_mat) <- my_metadata$sample[idx]

tail(gene_counts_mat)

                C2_norm C3_norm C5_norm C1_norm C1_cancer C2_cancer C3_cancer
ENSG00000293594       0       0       0       0         0         0         0
ENSG00000293595       3       5       3       0         0         0         0
ENSG00000293596       0       0       0       0         0         0         0
ENSG00000293597       1       2      11       1         2         3         1
ENSG00000293599       2       0       1       0         1         2         0
ENSG00000293600      45      59      85     561       789      1099       701
                C5_cancer
ENSG00000293594         0
ENSG00000293595         0
ENSG00000293596         0
ENSG00000293597         2
ENSG00000293599         0
ENSG00000293600       845

PCA

Create DESeqDataSet object.

lung_cancer <- DESeqDataSetFromMatrix(
  countData = gene_counts_mat,
  colData   = my_metadata,
  design    = ~ group
)
lung_cancer

class: DESeqDataSet 
dim: 63140 8 
metadata(1): version
assays(1): counts
rownames(63140): ENSG00000000003 ENSG00000000005 ... ENSG00000293599
  ENSG00000293600
rowData names(0):
colnames(8): C2_norm C3_norm ... C3_cancer C5_cancer
colData names(3): sample run_id group

Quickly estimate dispersion trend and apply a variance stabilizing transformation.

pas_vst <- vst(lung_cancer)
pas_vst

class: DESeqTransform 
dim: 63140 8 
metadata(1): version
assays(1): ''
rownames(63140): ENSG00000000003 ENSG00000000005 ... ENSG00000293599
  ENSG00000293600
rowData names(4): baseMean baseVar allZero dispFit
colnames(8): C2_norm C3_norm ... C3_cancer C5_cancer
colData names(4): sample run_id group sizeFactor

Plot PCA.

plotPCA(pas_vst, intgroup = "group") +
  theme_minimal()

using ntop=500 top features by variance

Past versions of plot_PCA-1.png

PCA data.

pca_data <- plotPCA(pas_vst, intgroup = "group", returnData = TRUE)

using ntop=500 top features by variance

pca_data

                 PC1         PC2  group group.1      name
C2_norm   -117.37808   7.0149036 normal  normal   C2_norm
C3_norm   -116.03913 -13.2708426 normal  normal   C3_norm
C5_norm   -115.34425   6.2593420 normal  normal   C5_norm
C1_norm     69.82117   0.8856676 normal  normal   C1_norm
C1_cancer   67.53316  -2.1774785 cancer  cancer C1_cancer
C2_cancer   70.13468   0.7640921 cancer  cancer C2_cancer
C3_cancer   70.50438   0.1610337 cancer  cancer C3_cancer
C5_cancer   70.76807   0.3632821 cancer  cancer C5_cancer

Source code.

plotPCA_copied = function(object, intgroup="condition", ntop=500, returnData=FALSE, pcsToUse=1:2){
  message(paste0("using ntop=",ntop," top features by variance"))
  
  # calculate the variance for each gene
  rv <- rowVars(assay(object))

  # select the ntop genes by variance
  select <- order(rv, decreasing=TRUE)[seq_len(min(ntop, length(rv)))]

  # perform a PCA on the data in assay(x) for the selected genes
  pca <- prcomp(t(assay(object)[select,]))

  # the contribution to the total variance for each component
  percentVar <- pca$sdev^2 / sum( pca$sdev^2 )

  if (!all(intgroup %in% names(colData(object)))) {
    stop("the argument 'intgroup' should specify columns of colData(dds)")
  }

  # add the intgroup factors together to create a new grouping factor
  group <- if (length(intgroup) > 1) {
    intgroup.df <- as.data.frame(colData(object)[, intgroup, drop=FALSE])
    factor(apply( intgroup.df, 1, paste, collapse=":"))
  } else {
    colData(object)[[intgroup]]
  }

  # assembly the data for the plot
  pcs <- paste0("PC", pcsToUse)
  d <- data.frame(V1=pca$x[,pcsToUse[1]],
                  V2=pca$x[,pcsToUse[2]],
                  group=group, name=colnames(object), colData(object))
  colnames(d)[1:2] <- pcs
  
  if (returnData) {
    attr(d, "percentVar") <- percentVar[pcsToUse]
    return(d)
  }

  ggplot(data=d, aes_string(x=pcs[1], y=pcs[2], color="group")) +
    geom_point(size=3) + 
    xlab(paste0(pcs[1],": ",round(percentVar[pcsToUse[1]] * 100),"% variance")) +
      ylab(paste0(pcs[2],": ",round(percentVar[pcsToUse[2]] * 100),"% variance")) +
        coord_fixed()
}

Plot using the copied function.

plotPCA_copied(pas_vst, intgroup = "group") +
  theme_minimal()

using ntop=500 top features by variance

Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
ℹ Please use tidy evaluation idioms with `aes()`.
ℹ See also `vignette("ggplot2-in-packages")` for more information.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.

Past versions of plot_PCA_copied-1.png

Perform PCA as per plotPCA().

# calculate the variance for each gene
rv <- rowVars(assay(pas_vst))
head(rv)

ENSG00000000003 ENSG00000000005 ENSG00000000419 ENSG00000000457 ENSG00000000460 
      6.3038555       3.4528327       0.1236461       0.5966854       1.0711842 
ENSG00000000938 
      3.6559005 

# select the ntop genes by variance
ntop <- 500
topgenes <- order(rv, decreasing=TRUE)[seq_len(min(ntop, length(rv)))]
head(topgenes)

[1] 12679 12446 12359 20116 53333  8411

# perform a PCA on the data in assay(x) for the selected genes
pca <- prcomp(t(assay(pas_vst)[topgenes,]))

Loadings are the coefficients that define how strongly each genes contributes to a principal component (PC).

loadings <- pca$rotation
loadings[1:10, 1:2]

                        PC1          PC2
ENSG00000169876 -0.06639949  0.008294131
ENSG00000168878  0.06304686  0.005076416
ENSG00000168484  0.05002038  0.117877848
ENSG00000205277 -0.05947887  0.023446665
ENSG00000275969 -0.05923209  0.013935605
ENSG00000143520 -0.05881607  0.007613250
ENSG00000087086  0.05824667  0.016638253
ENSG00000204849 -0.05843500  0.001321622
ENSG00000276040 -0.05730797  0.013343573
ENSG00000185775 -0.05689974 -0.001037869

Squared loadings represent how much variance of each PC is attributable to each variable.

contribution <- loadings^2
contribution[1:10, 1:2]

                        PC1          PC2
ENSG00000169876 0.004408893 6.879261e-05
ENSG00000168878 0.003974907 2.577000e-05
ENSG00000168484 0.002502039 1.389519e-02
ENSG00000205277 0.003537736 5.497461e-04
ENSG00000275969 0.003508441 1.942011e-04
ENSG00000143520 0.003459330 5.796158e-05
ENSG00000087086 0.003392674 2.768315e-04
ENSG00000204849 0.003414650 1.746686e-06
ENSG00000276040 0.003284204 1.780509e-04
ENSG00000185775 0.003237580 1.077171e-06

Gene contributing most to PC1.

sort(contribution[, 1], decreasing = TRUE) |> head()

ENSG00000169876 ENSG00000168878 ENSG00000205277 ENSG00000275969 ENSG00000143520 
    0.004408893     0.003974907     0.003537736     0.003508441     0.003459330 
ENSG00000204849 
    0.003414650 

Note that the order is similar, because plotPCA() already sorted the genes by their variance!

Session information

sessionInfo()

R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] DESeq2_1.46.0               SummarizedExperiment_1.36.0
 [3] Biobase_2.66.0              MatrixGenerics_1.18.1      
 [5] matrixStats_1.5.0           GenomicRanges_1.58.0       
 [7] GenomeInfoDb_1.42.3         IRanges_2.40.1             
 [9] S4Vectors_0.44.0            BiocGenerics_0.52.0        
[11] lubridate_1.9.3             forcats_1.0.0              
[13] stringr_1.5.1               dplyr_1.1.4                
[15] purrr_1.0.2                 readr_2.1.5                
[17] tidyr_1.3.1                 tibble_3.2.1               
[19] ggplot2_3.5.1               tidyverse_2.0.0            
[21] workflowr_1.7.1            

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.1        farver_2.1.2            fastmap_1.2.0          
 [4] promises_1.3.2          digest_0.6.37           timechange_0.3.0       
 [7] lifecycle_1.0.4         processx_3.8.4          magrittr_2.0.3         
[10] compiler_4.4.1          rlang_1.1.4             sass_0.4.9             
[13] tools_4.4.1             utf8_1.2.4              yaml_2.3.10            
[16] knitr_1.48              S4Arrays_1.6.0          labeling_0.4.3         
[19] bit_4.5.0               DelayedArray_0.32.0     abind_1.4-8            
[22] BiocParallel_1.40.0     withr_3.0.2             grid_4.4.1             
[25] git2r_0.35.0            colorspace_2.1-1        scales_1.3.0           
[28] cli_3.6.3               rmarkdown_2.28          crayon_1.5.3           
[31] generics_0.1.3          rstudioapi_0.17.1       httr_1.4.7             
[34] tzdb_0.4.0              cachem_1.1.0            zlibbioc_1.52.0        
[37] parallel_4.4.1          XVector_0.46.0          vctrs_0.6.5            
[40] Matrix_1.7-0            jsonlite_1.8.9          callr_3.7.6            
[43] hms_1.1.3               bit64_4.5.2             locfit_1.5-9.12        
[46] jquerylib_0.1.4         glue_1.8.0              codetools_0.2-20       
[49] ps_1.8.1                stringi_1.8.4           gtable_0.3.6           
[52] later_1.3.2             UCSC.utils_1.2.0        munsell_0.5.1          
[55] pillar_1.10.1           htmltools_0.5.8.1       GenomeInfoDbData_1.2.13
[58] R6_2.5.1                rprojroot_2.0.4         vroom_1.6.5            
[61] evaluate_1.0.1          lattice_0.22-6          highr_0.11             
[64] httpuv_1.6.15           bslib_0.8.0             Rcpp_1.0.13            
[67] SparseArray_1.6.2       whisker_0.4.1           xfun_0.48              
[70] fs_1.6.4                getPass_0.2-4           pkgconfig_2.0.3