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:
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",
macs2.path = "/home/stuartt/miniconda3/envs/signac/bin/macs2"
)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 | 565196 | 565405 | 210 | * | CD14+_Monocytes,CD8_Naive,B_cell_progenitor |
| chr1 | 569272 | 569498 | 227 | * | CD4_Naive,CD4_Memory,Double_negative_T_cell,CD8_Naive,B_cell_progenitor,CD14+_Monocytes,CD8_effector,pre-B_cell |
| chr1 | 713484 | 714543 | 1060 | * | CD14+_Monocytes,CD8_effector,CD4_Naive,pre-B_cell,Double_negative_T_cell,CD8_Naive,CD4_Memory,B_cell_progenitor,NK_dim,Dendritic_cell,CD16+_Monocytes,pDC,NK_bright |
| chr1 | 752307 | 752820 | 514 | * | CD14+_Monocytes,CD16+_Monocytes,NK_dim,CD4_Memory,NK_bright |
| chr1 | 762087 | 763172 | 1086 | * | CD14+_Monocytes,CD4_Memory,CD4_Naive,pre-B_cell,CD8_effector,CD8_Naive,CD16+_Monocytes,NK_dim,Double_negative_T_cell,B_cell_progenitor,NK_bright |
| chr1 | 779510 | 780112 | 603 | * | CD4_Naive,CD4_Memory,CD8_effector |
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"
)
Session Info
## R version 4.0.1 (2020-06-06)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
##
## 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.0.0 Seurat_4.0.1.9005 Signac_1.2.0
##
## loaded via a namespace (and not attached):
## [1] fastmatch_1.1-0 systemfonts_1.0.1 plyr_1.8.6
## [4] igraph_1.2.6 lazyeval_0.2.2 splines_4.0.1
## [7] BiocParallel_1.24.1 listenv_0.8.0 SnowballC_0.7.0
## [10] scattermore_0.7 GenomeInfoDb_1.26.5 ggplot2_3.3.3
## [13] digest_0.6.27 htmltools_0.5.1.1 fansi_0.4.2
## [16] magrittr_2.0.1 memoise_2.0.0 tensor_1.5
## [19] cluster_2.1.2 ROCR_1.0-11 globals_0.14.0
## [22] Biostrings_2.58.0 matrixStats_0.58.0 docopt_0.7.1
## [25] pkgdown_1.6.1 spatstat.sparse_2.0-0 colorspace_2.0-0
## [28] ggrepel_0.9.1 textshaping_0.3.3 xfun_0.22
## [31] dplyr_1.0.5 sparsesvd_0.2 crayon_1.4.1
## [34] RCurl_1.98-1.3 jsonlite_1.7.2 spatstat.data_2.1-0
## [37] survival_3.2-11 zoo_1.8-9 glue_1.4.2
## [40] polyclip_1.10-0 gtable_0.3.0 zlibbioc_1.36.0
## [43] XVector_0.30.0 leiden_0.3.7 future.apply_1.7.0
## [46] BiocGenerics_0.36.0 abind_1.4-5 scales_1.1.1
## [49] DBI_1.1.1 miniUI_0.1.1.1 Rcpp_1.0.6
## [52] viridisLite_0.4.0 xtable_1.8-4 reticulate_1.19
## [55] spatstat.core_2.1-2 stats4_4.0.1 htmlwidgets_1.5.3
## [58] httr_1.4.2 RColorBrewer_1.1-2 ellipsis_0.3.1
## [61] ica_1.0-2 pkgconfig_2.0.3 farver_2.1.0
## [64] ggseqlogo_0.1 sass_0.3.1 uwot_0.1.10
## [67] deldir_0.2-10 utf8_1.2.1 labeling_0.4.2
## [70] tidyselect_1.1.0 rlang_0.4.10 reshape2_1.4.4
## [73] later_1.2.0 munsell_0.5.0 tools_4.0.1
## [76] cachem_1.0.4 generics_0.1.0 ggridges_0.5.3
## [79] evaluate_0.14 stringr_1.4.0 fastmap_1.1.0
## [82] yaml_2.2.1 ragg_1.1.2 goftest_1.2-2
## [85] knitr_1.33 fs_1.5.0 fitdistrplus_1.1-3
## [88] purrr_0.3.4 RANN_2.6.1 pbapply_1.4-3
## [91] future_1.21.0 nlme_3.1-152 mime_0.10
## [94] slam_0.1-48 RcppRoll_0.3.0 compiler_4.0.1
## [97] plotly_4.9.3 png_0.1-7 spatstat.utils_2.1-0
## [100] tibble_3.1.1 tweenr_1.0.2 bslib_0.2.4
## [103] stringi_1.5.3 highr_0.9 desc_1.3.0
## [106] lattice_0.20-41 Matrix_1.3-2 vctrs_0.3.7
## [109] pillar_1.6.0 lifecycle_1.0.0 spatstat.geom_2.1-0
## [112] lmtest_0.9-38 jquerylib_0.1.4 RcppAnnoy_0.0.18
## [115] data.table_1.14.0 cowplot_1.1.1 bitops_1.0-7
## [118] irlba_2.3.3 httpuv_1.6.0 patchwork_1.1.1
## [121] GenomicRanges_1.42.0 R6_2.5.0 promises_1.2.0.1
## [124] lsa_0.73.2 KernSmooth_2.23-18 gridExtra_2.3
## [127] IRanges_2.24.1 parallelly_1.24.0 codetools_0.2-18
## [130] MASS_7.3-53.1 assertthat_0.2.1 rprojroot_2.0.2
## [133] qlcMatrix_0.9.7 sctransform_0.3.2 Rsamtools_2.6.0
## [136] S4Vectors_0.28.1 GenomeInfoDbData_1.2.4 mgcv_1.8-33
## [139] parallel_4.0.1 grid_4.0.1 rpart_4.1-15
## [142] tidyr_1.1.3 rmarkdown_2.7 Rtsne_0.15
## [145] ggforce_0.3.3 shiny_1.6.0