ATACseqTFEA Guide
25 April 2023
1 Introduction
ATAC-seq, an assay for Transposase-Accessible Chromatin using sequencing, is a widely used technique for chromatin accessibility analysis. Detecting differential activation of transcription factors between two different experiment conditions provides the possibility of decoding the key factors in a phenotype. Lots of tools have been developed to detect the differential activity of TFs (DATFs) for different groups of samples. Those tools can be divided into two groups. One group detects DATFs from differential accessibility analysis, such as MEME1, HOMER2, enrichr3, and ChEA4. Another group finds the DATFs by enrichment tests, such as BiFET5, diffTF6, and TFEA7. For single-cell ATAC-seq analysis, Signac and chromVar are widely used tools.
2 Motivation
All of these tools detect the DATF by only considering the open status of chromatin. None of them take the TF footprint into count. The open status provides the possibility of TF can bind to that position. The TF footprint by ATAC-seq shows the status of TF bindings.
To help researchers quickly assess the differential activity of hundreds of TFs by detecting the difference in TF footprint via enrichment score8, we have developed the ATACseqTFEA package. The ATACseqTFEA package is a robust and reliable computational tool to identify the key regulators responding to a phenotype.
schematic diagram of ATACseqTFEA
3 Quick start
Here is an example using ATACseqTFEA with a subset of ATAC-seq data.
3.1 Installation
First, install ATACseqTFEA and other packages required to run
the examples.
Please note that the example dataset used here is from zebrafish.
To run an analysis with a dataset from a different species or different assembly,
please install the corresponding Bsgenome and “TxDb”.
For example, to analyze mouse data aligned to “mm10”,
please install “BSgenome.Mmusculus.UCSC.mm10”,
and “TxDb.Mmusculus.UCSC.mm10.knownGene”.
You can also generate a TxDb object by
functions makeTxDbFromGFF from a local “gff” file,
or makeTxDbFromUCSC, makeTxDbFromBiomart, and makeTxDbFromEnsembl,
from online resources in the GenomicFeatures package.
library(BiocManager)
BiocManager::install(c("ATACseqTFEA",
"ATACseqQC",
"Rsamtools",
"BSgenome.Drerio.UCSC.danRer10",
"TxDb.Drerio.UCSC.danRer10.refGene"))3.2 Load library
library(ATACseqTFEA)
library(BSgenome.Drerio.UCSC.danRer10) ## for binding sites search
library(ATACseqQC) ## for footprint3.3 Prepare binding sites
To do TFEA, there are two inputs, the binding sites, and the change ranks.
To get the binding sites, the ATACseqTFEA package provides the
prepareBindingSites function. Users can also try to get the binding sites
list by other tools such as “fimo”9.
The prepareBindingSites function request a cluster of position weight matrix
(PWM) of TF motifs. ATACseqTFEA prepared a merged PWMatrixList for
405 motifs. The PWMatrixList is a collection of jasper2018, jolma2013 and
cisbp_1.02 from package motifDB (v 1.28.0) and merged by distance smaller than
1e-9 calculated by MotIV::motifDistances function (v 1.42.0).
The merged motifs were exported by motifStack (v 1.30.0).
motifs <- readRDS(system.file("extdata", "PWMatrixList.rds",
package="ATACseqTFEA"))The best_curated_Human is a list of TF motifs
downloaded from TFEA github7.
There are 1279 human motifs in the data set.
motifs_human <- readRDS(system.file("extdata", "best_curated_Human.rds",
package="ATACseqTFEA"))Another list of non-redundant TF motifs are also available by downloading the data from DeepSTARR10. There are 6502 motifs in the data set.
MotifsSTARR <- readRDS(system.file("extdata", "cluster_PWMs.rds",
package="ATACseqTFEA"))To scan the binding sites along a genome, a BSgenome object is required by
the prepareBindingSites function.
# for test run, we use a subset of data within chr1:5000-100000
# for real data, use the merged peaklist as grange input.
# Drerio is the short-link of BSgenome.Drerio.UCSC.danRer10
seqlev <- "chr1"
bindingSites <-
prepareBindingSites(motifs, Drerio, seqlev,
grange=GRanges("chr1", IRanges(5000, 100000)),
p.cutoff = 5e-05)#set higher p.cutoff to get more sites.3.4 TFEA
The correct insertion site is the key to the enrichment analysis of TF binding
sites. The parameter positive and negative in the function of TFEA
are used to shift the 5’ ends of the reads to the correct insertion positions.
However, this shift does not consider the soft clip of the reads.
The best way to generate correct shifted bam files is using
ATACseqQC::shiftGAlignmentsList11 for paired-end or
shiftGAlignments for single-end of the bam file.
The samples must be at least biologically duplicated for the one-step TFEA
function.
bamExp <- system.file("extdata",
c("KD.shift.rep1.bam",
"KD.shift.rep2.bam"),
package="ATACseqTFEA")
bamCtl <- system.file("extdata",
c("WT.shift.rep1.bam",
"WT.shift.rep2.bam"),
package="ATACseqTFEA")
res <- TFEA(bamExp, bamCtl, bindingSites=bindingSites,
positive=0, negative=0) # the bam files were shifted reads3.5 View results
The results will be saved in a TFEAresults object. We will use multiple
functions to present the results.
The plotES function will return a ggplot object for single TF input and
no outfolder is defined.
The ESvolcanoplot function will provide an overview of all the TFs enrichment.
And we can borrow the factorFootprints function from ATACseqQC package
to view the footprints of one TF.
TF <- "Tal1::Gata1"
## volcanoplot
ESvolcanoplot(TFEAresults=res, TFnameToShow=TF)
### plot enrichment score for one TF
plotES(res, TF=TF, outfolder=NA)
## footprint
sigs <- factorFootprints(c(bamCtl, bamExp),
pfm = as.matrix(motifs[[TF]]),
bindingSites = getBindingSites(res, TF=TF),
seqlev = seqlev, genome = Drerio,
upstream = 100, downstream = 100,
group = rep(c("WT", "KD"), each=2))
## export the results into a csv file
write.csv(res$resultsTable, tempfile(fileext = ".csv"),
row.names=FALSE)The command-line scripts are available at extdata named as sample_scripts.R.
4 Do TFEA step by step.
The one-step TFEA is a function containing multiple steps, which include:
- count the reads in binding sites, proximal region, and distal region;
- filter the binding site not open;
- normalize the count number by the width of the count region;
- calculate the binding scores and weight the binding scores by open scores;
- differential analysis by limma for the binding score
- filter the differential results by P-value and fold change
- TF enrichment analysis
If you want to tune the parameters, it will be much better to do it step by step to avoid repeating the computation for the same step. Here are the details for each step.
4.1 Counting reads
We will count the insertion site in binding sites, proximal and distal regions
by counting the 5’ ends of the reads in a shifted bam file.
Here we suggest keeping the proximal and distal the same value.
# prepare the counting region
exbs <- expandBindingSites(bindingSites=bindingSites,
proximal=40,
distal=40,
gap=10)
## count reads by 5'ends
counts <- count5ends(bam=c(bamExp, bamCtl),
positive=0L, negative=0L,
bindingSites = bindingSites,
bindingSitesWithGap=exbs$bindingSitesWithGap,
bindingSitesWithProximal=exbs$bindingSitesWithProximal,
bindingSitesWithProximalAndGap=
exbs$bindingSitesWithProximalAndGap,
bindingSitesWithDistal=exbs$bindingSitesWithDistal)4.2 Filter the counts
We filter the binding sites by at least there is 1 reads in proximal region. Users may want to try filter the sites by more stringent criteria such as “proximalRegion>1”.
colnames(counts)## [1] "bindingSites" "proximalRegion" "distalRegion"counts <- eventsFilter(counts, "proximalRegion>0")4.3 Normalize the counts by width of count region
We will normalize the counts to count per base (CPB).
counts <- countsNormalization(counts, proximal=40, distal=40)4.4 Get weighted binding scores
Here we use the open score to weight the binding score. Users can also define
the weight for binding score via parameter weight in
the function getWeightedBindingScore.
counts <- getWeightedBindingScore(counts)4.5 Differential analysis
Here we use DBscore, which borrows the power of the limma package,
to do differential binding analysis.
design <- cbind(CTL=1, EXPvsCTL=c(1, 1, 0, 0))
counts <- DBscore(counts, design=design, coef="EXPvsCTL")4.6 Filter the DB results
We can filter the binding results to decrease the data size by
the function eventsFilter.
For the sample data, we skip this step.
4.7 TF enrichment analysis
Last, we use the function doTFEA to get the enrichment scores.
res <- doTFEA(counts)
res## This is an object of TFEAresults with
## slot enrichmentScore ( matrix: 399 x 2166 ),
## slot bindingSites ( GRanges object with 2166 ranges and 12 metadata columns ),
## slot motifID ( a list of the positions of binding sites for 399 TFs ), and
## slot resultsTable ( 399 x 5 ). Here is the top 2 rows:
## TF enrichmentScore normalizedEnrichmentScore p_value adjPval
## NRF1 NRF1 0.1923613 0.7960275 0.7253614 0.9994472
## Gfi1b Gfi1b 0.3099024 1.3769160 0.1143751 0.9585665plotES(res, TF=TF, outfolder=NA) ## will show same figure as above one
5 SessionInfo
sessionInfo()## R version 4.3.0 RC (2023-04-13 r84269)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.2 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.17-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [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: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] ATACseqQC_1.24.0 Rsamtools_2.16.0
## [3] BSgenome.Drerio.UCSC.danRer10_1.4.2 BSgenome_1.68.0
## [5] rtracklayer_1.60.0 Biostrings_2.68.0
## [7] XVector_0.40.0 GenomicRanges_1.52.0
## [9] GenomeInfoDb_1.36.0 IRanges_2.34.0
## [11] S4Vectors_0.38.0 BiocGenerics_0.46.0
## [13] ATACseqTFEA_1.2.0 BiocStyle_2.28.0
##
## loaded via a namespace (and not attached):
## [1] splines_4.3.0 later_1.3.0
## [3] BiocIO_1.10.0 bitops_1.0-7
## [5] filelock_1.0.2 tibble_3.2.1
## [7] R.oo_1.25.0 graph_1.78.0
## [9] XML_3.99-0.14 DirichletMultinomial_1.42.0
## [11] lifecycle_1.0.3 edgeR_3.42.0
## [13] lattice_0.21-8 MASS_7.3-59
## [15] magrittr_2.0.3 limma_3.56.0
## [17] sass_0.4.5 rmarkdown_2.21
## [19] jquerylib_0.1.4 yaml_2.3.7
## [21] httpuv_1.6.9 grImport2_0.2-0
## [23] DBI_1.1.3 CNEr_1.36.0
## [25] ade4_1.7-22 zlibbioc_1.46.0
## [27] R.utils_2.12.2 preseqR_4.0.0
## [29] RCurl_1.98-1.12 pracma_2.4.2
## [31] rappdirs_0.3.3 GenomeInfoDbData_1.2.10
## [33] ggrepel_0.9.3 seqLogo_1.66.0
## [35] annotate_1.78.0 codetools_0.2-19
## [37] DelayedArray_0.26.0 xml2_1.3.3
## [39] tidyselect_1.2.0 futile.logger_1.4.3
## [41] farver_2.1.1 matrixStats_0.63.0
## [43] BiocFileCache_2.8.0 base64enc_0.1-3
## [45] GenomicAlignments_1.36.0 jsonlite_1.8.4
## [47] multtest_2.56.0 motifStack_1.44.0
## [49] ellipsis_0.3.2 survival_3.5-5
## [51] motifmatchr_1.22.0 tools_4.3.0
## [53] progress_1.2.2 TFMPvalue_0.0.9
## [55] Rcpp_1.0.10 glue_1.6.2
## [57] ChIPpeakAnno_3.34.0 xfun_0.39
## [59] MatrixGenerics_1.12.0 dplyr_1.1.2
## [61] HDF5Array_1.28.0 withr_2.5.0
## [63] formatR_1.14 BiocManager_1.30.20
## [65] fastmap_1.1.1 rhdf5filters_1.12.0
## [67] fansi_1.0.4 caTools_1.18.2
## [69] digest_0.6.31 R6_2.5.1
## [71] mime_0.12 colorspace_2.1-0
## [73] GO.db_3.17.0 gtools_3.9.4
## [75] poweRlaw_0.70.6 jpeg_0.1-10
## [77] biomaRt_2.56.0 RSQLite_2.3.1
## [79] R.methodsS3_1.8.2 utf8_1.2.3
## [81] generics_0.1.3 prettyunits_1.1.1
## [83] InteractionSet_1.28.0 httr_1.4.5
## [85] htmlwidgets_1.6.2 TFBSTools_1.38.0
## [87] regioneR_1.32.0 pkgconfig_2.0.3
## [89] gtable_0.3.3 blob_1.2.4
## [91] htmltools_0.5.5 bookdown_0.33
## [93] RBGL_1.76.0 scales_1.2.1
## [95] Biobase_2.60.0 png_0.1-8
## [97] knitr_1.42 lambda.r_1.2.4
## [99] tzdb_0.3.0 reshape2_1.4.4
## [101] rjson_0.2.21 curl_5.0.0
## [103] cachem_1.0.7 rhdf5_2.44.0
## [105] stringr_1.5.0 BiocVersion_3.17.1
## [107] KernSmooth_2.23-20 parallel_4.3.0
## [109] AnnotationDbi_1.62.0 restfulr_0.0.15
## [111] pillar_1.9.0 grid_4.3.0
## [113] vctrs_0.6.2 promises_1.2.0.1
## [115] randomForest_4.7-1.1 dbplyr_2.3.2
## [117] xtable_1.8-4 evaluate_0.20
## [119] magick_2.7.4 readr_2.1.4
## [121] VennDiagram_1.7.3 GenomicFeatures_1.52.0
## [123] cli_3.6.1 locfit_1.5-9.7
## [125] compiler_4.3.0 futile.options_1.0.1
## [127] rlang_1.1.0 crayon_1.5.2
## [129] labeling_0.4.2 plyr_1.8.8
## [131] stringi_1.7.12 BiocParallel_1.34.0
## [133] munsell_0.5.0 Matrix_1.5-4
## [135] hms_1.1.3 bit64_4.0.5
## [137] ggplot2_3.4.2 Rhdf5lib_1.22.0
## [139] KEGGREST_1.40.0 shiny_1.7.4
## [141] highr_0.10 SummarizedExperiment_1.30.0
## [143] interactiveDisplayBase_1.38.0 AnnotationHub_3.8.0
## [145] GenomicScores_2.12.0 memoise_2.0.1
## [147] bslib_0.4.2 bit_4.0.5
## [149] polynom_1.4-1References
1. Bailey, T. L., Williams, N., Misleh, C. & Li, W. W. MEME: Discovering and analyzing dna and protein sequence motifs. Nucleic acids research 34, W369–W373 (2006).
2. Heinz, S. et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and b cell identities. Molecular cell 38, 576–589 (2010).
3. Chen, E. Y. et al. Enrichr: Interactive and collaborative html5 gene list enrichment analysis tool. BMC bioinformatics 14, 128 (2013).
4. Lachmann, A. et al. ChEA: Transcription factor regulation inferred from integrating genome-wide chip-x experiments. Bioinformatics 26, 2438–2444 (2010).
5. Youn, A., Marquez, E. J., Lawlor, N., Stitzel, M. L. & Ucar, D. BiFET: Sequencing bi as-free transcription factor f ootprint e nrichment t est. Nucleic acids research 47, e11–e11 (2019).
6. Berest, I. et al. Quantification of differential transcription factor activity and multiomics-based classification into activators and repressors: DiffTF. Cell Reports 29, 3147–3159 (2019).
7. Rubin, J. D. et al. Transcription factor enrichment analysis (tfea): Quantifying the activity of hundreds of transcription factors from a single experiment. Commun Biol 661 (2021).
8. Subramanian, A. et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proceedings of the National Academy of Sciences 102, 15545–15550 (2005).
9. Grant, C. E., Bailey, T. L. & Noble, W. S. FIMO: Scanning for occurrences of a given motif. Bioinformatics 27, 1017–1018 (2011).
10. Almeida, B. P. de, Reiter, F., Pagani, M. & Stark, A. DeepSTARR predicts enhancer activity from dna sequence and enables the de novo design of enhancers. bioRxiv (2021).
11. Ou, J. et al. ATACseqQC: A bioconductor package for post-alignment quality assessment of atac-seq data. BMC genomics 19, 169 (2018).