ATACCoGAPS
Rossin Erbe
Introduction
Coordinated Gene Association in Pattern Sets (CoGAPS) is a technique for latent space learning. CoGAPS is a member of the Nonnegative Matrix Factorization (NMF) class of algorithms. NMFs factorize a data matrix into two related matrices containing gene weights, the Amplitude (A) matrix, and sample weights, the Pattern (P) Matrix. Each column of A or row of P defines a feature and together this set of features defines the latent space among features and samples, respectively. In NMF, the values of the elements in the A and P matrices are constrained to be greater than or equal to zero. This constraint simultaneously reflects the non-negative nature of molecular data and enforces the additive nature of the resulting feature dimensions, generating solutions that are biologically intuitive to interpret (Seung and Lee (1999)).
This package extends CoGAPS for usage with scATACseq data, allowing for summarization to genomic peaks or motifs as features for input into CoGAPS. A count matrix - peaks or motifs by cells is used as input, yielding patterns of accessibility that distinguish biological variation between cells.
Installation
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("ATACCoGAPS")
Package Overview
To outline the ATACCoGAPS pipeline, we will use as an example data set single-cell ATAC sequencing data containing 10 cell lines, published by Schep et al, 2017. The data was downloaded from GEO accession number GSE99172 and preprocessed using dataSubsetBySparsity() from the ATACCoGAPS package to remove cells and peaks with more than 99% sparsity (more than 99% zeroes). We will use a subset of the data containing 50 random cells from each cell line and the top ~5000 most variable peaks as an example for this vignette.
data("subsetSchepData")
data("schepPeaks")
data("schepCellTypes")
We use these data to set the hyperparameters of the CoGAPS algorithm. Here we tell CoGAPS to find 7 patterns in 10000 iterations of the algorithm. We use the sparseOptimization method as our data are sparse single-cell data. We run the algorithm distributed across the genome since we have more genomic features than cells (if it was the opposite we would set the distributed pattern to “single-cell”). We then input the peak and cell type information to be returned as part of our result object.
params <- CogapsParams(nPatterns=7, nIterations=10000, seed=42, sparseOptimization=TRUE, distributed="genome-wide", geneNames=schepPeaks, sampleNames=as.character(schepCellTypes))
params
## -- Standard Parameters --
## nPatterns 7
## nIterations 10000
## seed 42
## sparseOptimization TRUE
## distributed genome-wide
##
## -- Sparsity Parameters --
## alpha 0.01
## maxGibbsMass 100
##
## -- Distributed CoGAPS Parameters --
## nSets 4
## cut 7
## minNS 2
## maxNS 6
##
## 5036 gene names provided
## first gene name: chr1-1051395-1051894
##
## 600 sample names provided
## first sample name: Fibroblasts
We now call CoGAPS via the R function. CoGAPS is a Bayesian Non-Negative Matrix Factorization algorithm (Fertig et al, 2010). It factorizes count matrices from RNA or epigenetic sequencing data and returns patterns which distinguish both features and samples, allowing for the discovery of regulatory differences between samples. In the case of scATAC-seq our features are usually peaks and our samples are indvidual cells.
schepCogapsResult <- CoGAPS(data = subsetSchepData, params = params, nThreads = 1)
Without parallelization the above takes a little while to run, so we provide the pre-computed output object
data("schepCogapsResult")
Pattern Matrix Visualization
The first quick visualization of CoGAPS results is generally plotting the Pattern Matrix (the output matrix which is patterns x cells). These plots allow us to determine which patterns differentiate which cell types.
We can either plot each pattern indvidually
#colors to plot by
col <- viridis::plasma(n=12)
cgapsPlot(cgaps_result = schepCogapsResult, sample.classifier = schepCellTypes, cols = col, ylab = "Pattern Weight")
Or all together in a heatmap
heatmapPatternMatrix(cgaps_result = schepCogapsResult, sample.classifier = schepCellTypes, cellCols = col, col = viridis::magma(n=9))
We can note which patterns differentiate which cell types (for example that pattern 9 seems to be defining the LCL cells). If any patterns are unclear, such as pattern 5, we can perform a Wilcoxon Rank Sum test to determine which cell types are most significantly associated with the pattern.
#get the pattern Matrix
patMatrix <- getSampleFactors(schepCogapsResult)
#perform a pairwise Wilcoxon test
pairwise.wilcox.test(patMatrix[,5], schepCellTypes, p.adjust.method = "BH")
##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: patMatrix[, 5] and schepCellTypes
##
## Fibroblasts LCL 1 LCL 2 Rep1 LCL 2 Rep2 ESCs
## LCL 1 1.9e-13 - - - -
## LCL 2 Rep1 7.9e-09 0.00215 - - -
## LCL 2 Rep2 4.8e-09 0.02188 0.73104 - -
## ESCs 0.38529 2.1e-11 2.4e-06 1.5e-06 -
## HL60 Leukemia 0.00182 8.6e-09 0.00105 0.00061 0.07705
## K562 Erythroleukemia 1.4e-09 9.7e-16 6.3e-15 7.7e-15 2.5e-09
## LMPP 6.3e-11 0.00718 0.67835 0.95876 4.1e-08
## Monocyte 1.2e-12 0.21600 0.09527 0.28185 1.1e-09
## AML Patient 070 1.7e-10 0.06924 0.43176 0.74080 1.3e-07
## AML Patient 353 2.2e-11 0.10374 0.20668 0.38529 5.6e-09
## Erythroblast 1.2e-12 0.43285 0.00022 0.00306 1.7e-11
## HL60 Leukemia K562 Erythroleukemia LMPP Monocyte
## LCL 1 - - - -
## LCL 2 Rep1 - - - -
## LCL 2 Rep2 - - - -
## ESCs - - - -
## HL60 Leukemia - - - -
## K562 Erythroleukemia 1.7e-12 - - -
## LMPP 6.0e-05 1.4e-15 - -
## Monocyte 6.3e-07 9.7e-16 0.22002 -
## AML Patient 070 0.00011 1.4e-15 0.76118 0.38301
## AML Patient 353 3.1e-06 1.8e-15 0.33555 0.80825
## Erythroblast 5.2e-09 1.4e-15 0.00039 0.05713
## AML Patient 070 AML Patient 353
## LCL 1 - -
## LCL 2 Rep1 - -
## LCL 2 Rep2 - -
## ESCs - -
## HL60 Leukemia - -
## K562 Erythroleukemia - -
## LMPP - -
## Monocyte - -
## AML Patient 070 - -
## AML Patient 353 0.67361 -
## Erythroblast 0.02682 0.01342
##
## P value adjustment method: BH
We see that pattern 5 is most strongly associated with the K562 cells.
If we do not have pre-established cell annotations, we can cluster cells by pattern association.
cellClass <- patternMarkerCellClassifier(schepCogapsResult)
cellClasses <- cellClass$cellClassifier
heatmapPatternMatrix(schepCogapsResult, as.factor(cellClasses), col = viridis::magma(n=9))
Finding Regulatory Differences between Cell Types
Now that we know which patterns distinguish which cell types, we can look at those same patterns in the amplitude matrix (peaks by patterns) to determine which peaks are differentially accessible between the patterns and thus which peaks are differentially accessible between the cell types.
We can use the patternMarker Statistic (Stein-O’Brien et al, 2017) to find which peaks are most differentially accessible. To show the degree of differentiation, we can plot the most pattern differentiating peaks for each pattern from the original data.
heatmapPatternMarkers(cgaps_result = schepCogapsResult, atac_data = subsetSchepData, celltypes = schepCellTypes, numregions = 5, colColors = col, col = viridis::plasma(n = 2))
Pathway Based Analysis
To make use of these pttern marker peaks, one option is to try to find genes that fall within these peaks and determine whether the accessibility of certain groups of genes suggests differential pathway activation.
data("schepGranges")
#loading TxDb of human genes
library(Homo.sapiens)
#find genes known to fall within thresholded patternMarker peaks for each pattern
genes <- genePatternMatch(cogapsResult = schepCogapsResult, generanges = schepGranges, genome = Homo.sapiens)
#download hallmark pathways using msigdbr
library(dplyr)
pathways <- msigdbr::msigdbr(species = "Homo sapiens", category =
"H") %>% dplyr::select(gs_name, gene_symbol) %>% as.data.frame()
#match these pattern Gene sets to hallmark pathways, using an adjusted p-value threshold of 0.001.
matchedPathways <- pathwayMatch(gene_list = genes, pathways = pathways, p_threshold = 0.001)
lapply(matchedPathways, function(x) {x[4]})
## [[1]]
## [[1]]$summaryTable
## NULL
##
##
## [[2]]
## [[2]]$summaryTable
## NULL
##
##
## [[3]]
## [[3]]$summaryTable
## NULL
##
##
## [[4]]
## [[4]]$summaryTable
## NULL
##
##
## [[5]]
## [[5]]$summaryTable
## NULL
##
##
## [[6]]
## [[6]]$summaryTable
## NULL
##
##
## [[7]]
## [[7]]$summaryTable
## NULL
##
##
## [[8]]
## [[8]]$summaryTable
## NULL
##
##
## [[9]]
## [[9]]$summaryTable
## NULL
Several patterns do not return Hallmark pathways at this level of significance, but those that do seem logical in the cell types those patterns differentiate.
Motif/Transcription Factor Based Analysis
The other way we can use pattern marker peak information is to match to DNA motifs and known Transcription Factor binding at those motifs.
motifResults <- simpleMotifTFMatch(cogapsResult = schepCogapsResult, generanges = schepGranges, organism = "Homo sapiens", genome = "hg19", motifsPerRegion = 1)
##
## Attaching package: 'Biostrings'
## The following object is masked from 'package:base':
##
## strsplit
##
## Attaching package: 'rtracklayer'
## The following object is masked from 'package:BiocIO':
##
## FileForFormat
We can get a summary of TF binding, generally having more confidence in those TFs that have higher numbers of motifs at which the same TF could bind.
motifResults$tfMatchSummary[[1]]
## SPI1
## 1
We can also examine the accessibility of the TF itself in a cell type of interest in order to gather information on whether the TF is expressed to bind at accessible sites. We’ll test the accessibility of EGR1 in monocytes as an example.
#get peaks overlapping with the gene
EGR1peaks <- geneAccessibility("EGR1", schepGranges, subsetSchepData, Homo.sapiens)
## 403 genes were dropped because they have exons located on both strands
## of the same reference sequence or on more than one reference sequence,
## so cannot be represented by a single genomic range.
## Use 'single.strand.genes.only=FALSE' to get all the genes in a
## GRangesList object, or use suppressMessages() to suppress this message.
## 'select()' returned 1:1 mapping between keys and columns
## 403 genes were dropped because they have exons located on both strands
## of the same reference sequence or on more than one reference sequence,
## so cannot be represented by a single genomic range.
## Use 'single.strand.genes.only=FALSE' to get all the genes in a
## GRangesList object, or use suppressMessages() to suppress this message.
## 'select()' returned 1:1 mapping between keys and columns
#make binary accessibility matrix
binaryMatrix <- (subsetSchepData > 0) + 0
#find accessibility of those peaks relative to others among monocyte cells
foldAccessibility(peaksAccessibility = EGR1peaks$EGR1, cellTypeList = schepCellTypes, cellType = "Monocyte", binaryMatrix = binaryMatrix)
## [1] 1.429779
EGR1 is 1.42 times more accessible than average in Monocytes
Session Info
## R version 4.4.0 beta (2024-04-15 r86425 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows Server 2022 x64 (build 20348)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=C
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] BSgenome.Hsapiens.UCSC.hg19_1.4.3
## [2] BSgenome_1.72.0
## [3] rtracklayer_1.64.0
## [4] BiocIO_1.14.0
## [5] Biostrings_2.72.0
## [6] XVector_0.44.0
## [7] dplyr_1.1.4
## [8] Homo.sapiens_1.3.1
## [9] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
## [10] org.Hs.eg.db_3.19.1
## [11] GO.db_3.19.1
## [12] OrganismDbi_1.46.0
## [13] GenomicFeatures_1.56.0
## [14] GenomicRanges_1.56.0
## [15] GenomeInfoDb_1.40.0
## [16] AnnotationDbi_1.66.0
## [17] IRanges_2.38.0
## [18] S4Vectors_0.42.0
## [19] Biobase_2.64.0
## [20] BiocGenerics_0.50.0
## [21] ATACCoGAPS_1.6.0
## [22] CoGAPS_3.24.0
##
## loaded via a namespace (and not attached):
## [1] later_1.3.2
## [2] bitops_1.0-7
## [3] GeneOverlap_1.40.0
## [4] filelock_1.0.3
## [5] R.oo_1.26.0
## [6] tibble_3.2.1
## [7] graph_1.82.0
## [8] XML_3.99-0.16.1
## [9] DirichletMultinomial_1.46.0
## [10] lifecycle_1.0.4
## [11] httr2_1.0.1
## [12] pwalign_1.0.0
## [13] lattice_0.22-6
## [14] magrittr_2.0.3
## [15] plotly_4.10.4
## [16] sass_0.4.9
## [17] rmarkdown_2.26
## [18] jquerylib_0.1.4
## [19] yaml_2.3.8
## [20] httpuv_1.6.15
## [21] cowplot_1.1.3
## [22] DBI_1.2.2
## [23] CNEr_1.40.0
## [24] RColorBrewer_1.1-3
## [25] abind_1.4-5
## [26] zlibbioc_1.50.0
## [27] purrr_1.0.2
## [28] R.utils_2.12.3
## [29] msigdbr_7.5.1
## [30] RCurl_1.98-1.14
## [31] pracma_2.4.4
## [32] rappdirs_0.3.3
## [33] GenomeInfoDbData_1.2.12
## [34] seqLogo_1.70.0
## [35] annotate_1.82.0
## [36] codetools_0.2-20
## [37] DelayedArray_0.30.0
## [38] DT_0.33
## [39] xml2_1.3.6
## [40] tidyselect_1.2.1
## [41] Mus.musculus_1.3.1
## [42] UCSC.utils_1.0.0
## [43] viridis_0.6.5
## [44] matrixStats_1.3.0
## [45] BiocFileCache_2.12.0
## [46] GenomicAlignments_1.40.0
## [47] jsonlite_1.8.8
## [48] motifmatchr_1.26.0
## [49] tools_4.4.0
## [50] progress_1.2.3
## [51] TFMPvalue_0.0.9
## [52] Rcpp_1.0.12
## [53] glue_1.7.0
## [54] gridExtra_2.3
## [55] SparseArray_1.4.0
## [56] xfun_0.43
## [57] MatrixGenerics_1.16.0
## [58] withr_3.0.0
## [59] BiocManager_1.30.22
## [60] fastmap_1.1.1
## [61] rhdf5filters_1.16.0
## [62] fansi_1.0.6
## [63] caTools_1.18.2
## [64] digest_0.6.35
## [65] R6_2.5.1
## [66] mime_0.12
## [67] colorspace_2.1-0
## [68] gtools_3.9.5
## [69] poweRlaw_0.80.0
## [70] chromVAR_1.26.0
## [71] biomaRt_2.60.0
## [72] RSQLite_2.3.6
## [73] R.methodsS3_1.8.2
## [74] tidyr_1.3.1
## [75] utf8_1.2.4
## [76] generics_0.1.3
## [77] data.table_1.15.4
## [78] prettyunits_1.2.0
## [79] httr_1.4.7
## [80] htmlwidgets_1.6.4
## [81] S4Arrays_1.4.0
## [82] org.Mm.eg.db_3.19.1
## [83] TFBSTools_1.42.0
## [84] pkgconfig_2.0.3
## [85] gtable_0.3.5
## [86] blob_1.2.4
## [87] SingleCellExperiment_1.26.0
## [88] htmltools_0.5.8.1
## [89] fgsea_1.30.0
## [90] RBGL_1.80.0
## [91] scales_1.3.0
## [92] tidyverse_2.0.0
## [93] png_0.1-8
## [94] knitr_1.46
## [95] tzdb_0.4.0
## [96] reshape2_1.4.4
## [97] rjson_0.2.21
## [98] JASPAR2016_1.31.0
## [99] curl_5.2.1
## [100] cachem_1.0.8
## [101] rhdf5_2.48.0
## [102] stringr_1.5.1
## [103] KernSmooth_2.23-22
## [104] parallel_4.4.0
## [105] miniUI_0.1.1.1
## [106] restfulr_0.0.15
## [107] pillar_1.9.0
## [108] grid_4.4.0
## [109] vctrs_0.6.5
## [110] gplots_3.1.3.1
## [111] promises_1.3.0
## [112] dbplyr_2.5.0
## [113] xtable_1.8-4
## [114] cluster_2.1.6
## [115] evaluate_0.23
## [116] readr_2.1.5
## [117] cli_3.6.2
## [118] compiler_4.4.0
## [119] Rsamtools_2.20.0
## [120] rlang_1.1.3
## [121] crayon_1.5.2
## [122] plyr_1.8.9
## [123] forcats_1.0.0
## [124] stringi_1.8.3
## [125] TxDb.Mmusculus.UCSC.mm10.knownGene_3.10.0
## [126] viridisLite_0.4.2
## [127] BiocParallel_1.38.0
## [128] babelgene_22.9
## [129] txdbmaker_1.0.0
## [130] munsell_0.5.1
## [131] lazyeval_0.2.2
## [132] Matrix_1.7-0
## [133] hms_1.1.3
## [134] bit64_4.0.5
## [135] ggplot2_3.5.1
## [136] Rhdf5lib_1.26.0
## [137] KEGGREST_1.44.0
## [138] shiny_1.8.1.1
## [139] SummarizedExperiment_1.34.0
## [140] highr_0.10
## [141] memoise_2.0.1
## [142] bslib_0.7.0
## [143] fastmatch_1.1-4
## [144] bit_4.0.5
## [145] BSgenome.Mmusculus.UCSC.mm10_1.4.3