In this tutorial, we demonstrate how to call peaks on a single-cell ATAC-seq dataset using MACS2.

To use the peak calling functionality in Signac you will first need to install MACS2. This can be done using pip or conda, or by building the package from source.

In this demonstration we use scATAC-seq data for human PBMCs. See our vignette for the code used to generate this object, and links to the raw data. First, load the required packages and the pre-computed Seurat object:

library(Signac)
library(Seurat)

pbmc <- readRDS("../vignette_data/pbmc.rds")
DimPlot(pbmc)

Peak calling can be performed using the CallPeaks() function, and can either be done separately for different groups of cells, or performed using data from all the cells. To call peaks on each annotated cell type, we can use the group.by argument:

peaks <- CallPeaks(
  object = pbmc,
  group.by = "predicted.id"
)

The results are returned as a GRanges object, with an additional metadata column listing the cell types that each peak was identified in:

seqnames start end width strand peak_called_in
chr1 10146 10433 288 * CD4_Memory
chr1 565204 565404 201 * CD8_Naive
chr1 569273 569499 227 * CD4_Naive,CD4_Memory,pre-B_cell,Double_negative_T_cell,CD8_Naive,B_cell_progenitor,CD14+_Monocytes,CD8_effector
chr1 713485 714543 1059 * CD14+_Monocytes,CD8_effector,CD4_Naive,CD8_Naive,CD4_Memory,pre-B_cell,Double_negative_T_cell,B_cell_progenitor,NK_dim,CD16+_Monocytes,Dendritic_cell,pDC
chr1 752308 752820 513 * CD14+_Monocytes,Double_negative_T_cell,CD16+_Monocytes,NK_dim,CD4_Memory
chr1 762088 763182 1095 * CD14+_Monocytes,CD4_Naive,pre-B_cell,CD8_effector,CD4_Memory,CD16+_Monocytes,CD8_Naive,NK_dim,Double_negative_T_cell,Dendritic_cell,B_cell_progenitor

To quantify counts in each peak, you can use the FeatureMatrix() function.

We can visualize the cell-type-specific MACS2 peak calls alongside the 10x Cellranger peak calls (currently being used in the pbmc object) with the CoveragePlot() function. Here the Cellranger peaks are shown in grey and the MACS2 peaks in red:

CoveragePlot(
  object = pbmc,
  region = "CD8A",
  ranges = peaks,
  ranges.title = "MACS2"
)
## Warning: Removed 1 rows containing missing values (`geom_segment()`).
Session Info
## R version 4.2.2 (2022-10-31)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux 8.6 (Ootpa)
## 
## Matrix products: default
## BLAS/LAPACK: /home/users/astar/gis/stuartt/mambaforge/envs/renv/lib/libopenblasp-r0.3.21.so
## 
## 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       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] SeuratObject_4.1.3 Seurat_4.3.0       Signac_1.10.0     
## 
## loaded via a namespace (and not attached):
##   [1] fastmatch_1.1-3        systemfonts_1.0.4      plyr_1.8.8            
##   [4] igraph_1.4.3           lazyeval_0.2.2         sp_1.6-0              
##   [7] splines_4.2.2          BiocParallel_1.32.6    listenv_0.9.0         
##  [10] scattermore_0.8        GenomeInfoDb_1.34.9    ggplot2_3.4.2         
##  [13] digest_0.6.31          htmltools_0.5.5        fansi_1.0.4           
##  [16] magrittr_2.0.3         memoise_2.0.1          tensor_1.5            
##  [19] cluster_2.1.4          ROCR_1.0-11            globals_0.16.2        
##  [22] Biostrings_2.66.0      matrixStats_0.63.0     pkgdown_2.0.7         
##  [25] spatstat.sparse_3.0-1  colorspace_2.1-0       ggrepel_0.9.3         
##  [28] textshaping_0.3.6      xfun_0.39              dplyr_1.1.2           
##  [31] crayon_1.5.2           RCurl_1.98-1.12        jsonlite_1.8.4        
##  [34] progressr_0.13.0       spatstat.data_3.0-1    survival_3.5-5        
##  [37] zoo_1.8-11             glue_1.6.2             polyclip_1.10-4       
##  [40] gtable_0.3.3           zlibbioc_1.44.0        XVector_0.38.0        
##  [43] leiden_0.4.3           future.apply_1.10.0    BiocGenerics_0.44.0   
##  [46] abind_1.4-5            scales_1.2.1           spatstat.random_3.1-4 
##  [49] miniUI_0.1.1.1         Rcpp_1.0.10            viridisLite_0.4.2     
##  [52] xtable_1.8-4           reticulate_1.28        stats4_4.2.2          
##  [55] htmlwidgets_1.6.2      httr_1.4.6             RColorBrewer_1.1-3    
##  [58] ellipsis_0.3.2         ica_1.0-3              farver_2.1.1          
##  [61] pkgconfig_2.0.3        sass_0.4.6             uwot_0.1.14           
##  [64] deldir_1.0-9           utf8_1.2.3             labeling_0.4.2        
##  [67] tidyselect_1.2.0       rlang_1.1.1            reshape2_1.4.4        
##  [70] later_1.3.0            munsell_0.5.0          tools_4.2.2           
##  [73] cachem_1.0.8           cli_3.6.1              generics_0.1.3        
##  [76] ggridges_0.5.4         evaluate_0.21          stringr_1.5.0         
##  [79] fastmap_1.1.1          yaml_2.3.7             ragg_1.2.4            
##  [82] goftest_1.2-3          knitr_1.43             fs_1.6.2              
##  [85] fitdistrplus_1.1-8     purrr_1.0.1            RANN_2.6.1            
##  [88] pbapply_1.7-0          future_1.32.0          nlme_3.1-162          
##  [91] mime_0.12              RcppRoll_0.3.0         compiler_4.2.2        
##  [94] plotly_4.10.1          png_0.1-8              spatstat.utils_3.0-2  
##  [97] tibble_3.2.1           bslib_0.4.2            stringi_1.7.12        
## [100] highr_0.10             desc_1.4.2             lattice_0.21-8        
## [103] Matrix_1.5-4           vctrs_0.6.2            pillar_1.9.0          
## [106] lifecycle_1.0.3        spatstat.geom_3.1-0    lmtest_0.9-40         
## [109] jquerylib_0.1.4        RcppAnnoy_0.0.20       data.table_1.14.8     
## [112] cowplot_1.1.1          bitops_1.0-7           irlba_2.3.5.1         
## [115] httpuv_1.6.9           patchwork_1.1.2        GenomicRanges_1.50.2  
## [118] R6_2.5.1               promises_1.2.0.1       KernSmooth_2.23-20    
## [121] gridExtra_2.3          IRanges_2.32.0         parallelly_1.35.0     
## [124] codetools_0.2-19       MASS_7.3-58.3          rprojroot_2.0.3       
## [127] withr_2.5.0            sctransform_0.3.5      Rsamtools_2.14.0      
## [130] S4Vectors_0.36.2       GenomeInfoDbData_1.2.9 parallel_4.2.2        
## [133] grid_4.2.2             tidyr_1.3.0            rmarkdown_2.21        
## [136] Rtsne_0.16             spatstat.explore_3.1-0 shiny_1.7.4