‘Motif2Site’: an R package to detect binding sites from ChIP-seq and recenter them
April 25, 2023
Abstract
‘Motif2Site’ is an R package to detect transcription factor binding sites from motif sets and ChIP-seq experiments. The motif sets are either user provided bed files or user provided DNA nucleotide sequence motifs. It also combines andcompares binding sites of a transcription factor across various experiments and conditions.
Package
Motif2Site 1.4.0
Contents
1 Introduction
Transcription factors often bind to specific DNA nucleotide patterns referred to as sequence motifs. Although sequence motifs are well-characterized for many transcription factors, detecting the actual binding sites is not always a straightforward task. Chromatin immunoprecipitation (ChIP) followed by DNA sequencing (ChIP-seq) is the major technology to detect binding regions of transcription factors. However, the binding regions detected by ChIP-seq may contains several or no DNA sequence motifs. Motif2Site is a novel R package which uses ChIP-seq information and detect transcription factor binding sites from a user provided DNA sequence motifs set.
Motif2Site gets two different input, motif and ChIP-seq alignment information, to detect binding sites. First input is ChIP-seq alignment short reads in the bam or bed format. For each aligned short read the center of the alignment is calculated using GenomicAlignments (Lawrence et al. 2013) and Rsamtools (Morgan et al. 2020) packages. Motif information is the second input which is provided by user either as a bed file or a DNA string with a mismatch number. In the case of DNA string input, Biostrings (Pagès et al. 2019) and BSgenome (H. Pagès 2019) packages are used to find motif locations on the genome in GenomicRanges (Lawrence et al. 2013) format.
Negative binomial distribution is used to model count data by using edgeR package (Robinson, McCarthy, and Smyth 2010). mixtools (Benaglia et al. 2009) is used to deconvolve binding intensities of closely spaced binding sites. Motif2Site also combines binding sites across different experiments. It calls differential binding sites using TMM normalization and GLM test using edgeR package (Robinson, McCarthy, and Smyth 2010).
2 Major functions of Motif2Site
First install and load the libraries needed to run the examples of this document:
library(GenomicRanges)
library(Motif2Site)
library(BSgenome.Scerevisiae.UCSC.sacCer3)
library(BSgenome.Ecoli.NCBI.20080805)The functions, implemented in Motif2Site, perform three tasks: 1. To assist users to select better sequence motif input. 2. To detect binding sites from sequence motifs and ChIP-seq datasets 3. To combine and compare binding sites across different experiments, conditions, or tissues. Each of these functions are explained in a separate section.
3 Selecting sequence motif
Motif2Site uses DNA motif information as one of its input. To facilitate choosing proper sequence motif, compareBedFiless2UserProvidedRegions and compareMotifs2UserProvidedRegions functions compare motif regions with a user provided confident regions in terms of precision and recall.
As the first example, an artificially generated “YeastSampleMotif.bed” bed file in yeast is considered as a confident binding regions set. compareMotifs2UserProvidedRegions function compares these regions with locations of ‘TGATTSCAGGANT’ 1-mismatch, ‘TGATTCCAGGANT’ 0-mismatch, ‘TGATWSCAGGANT’ 2-mismatches on the yeast genome in terms of precision and recall.
yeastExampleFile =
system.file("extdata", "YeastSampleMotif.bed", package="Motif2Site")
YeastRegionsChIPseq <- Bed2Granges(yeastExampleFile)
SequenceComparison <-
compareMotifs2UserProvidedRegions(
givenRegion=YeastRegionsChIPseq,
motifs = c("TGATTSCAGGANT", "TGATTCCAGGANT", "TGATWSCAGGANT"),
mismatchNumbers = c(1,0,2),
genome="Scerevisiae",
genomeBuild="sacCer3"
)
SequenceComparison## motifName regionCoverage motifCoverage
## 1 TGATTSCAGGANT_1mismatch 0.350877193 0.6060606
## 2 TGATTCCAGGANT_0mismatch 0.005847953 1.0000000
## 3 TGATWSCAGGANT_2mismatch 1.000000000 0.0475132In the following example, an artificially generated “YeastSampleMotif.bed” bed file in yeast is considered as a confident binding regions set. compareBedFiless2UserProvidedRegions function compares these regions with two bed files “YeastBedFile1.bed” and “YeastBedFile2.bed” in terms of precision and recall.
# Yeast artificial dataset for comparison bed files
yeastExampleFile =
system.file("extdata", "YeastSampleMotif.bed", package="Motif2Site")
YeastRegionsChIPseq <- Bed2Granges(yeastExampleFile)
bed1 <- system.file("extdata", "YeastBedFile1.bed", package="Motif2Site")
bed2 <- system.file("extdata", "YeastBedFile2.bed", package="Motif2Site")
BedFilesVector <- c(bed1, bed2)
SequenceComparison <-
compareBedFiless2UserProvidedRegions(
givenRegion=YeastRegionsChIPseq,
bedfiles = BedFilesVector,
motifnames = c("YeastBed1", "YeastBed2")
)
SequenceComparison## motifName regionCoverage motifCoverage
## 1 YeastBed1 0.3099415 1.000000
## 2 YeastBed2 0.4853801 0.2090684 Detecting binding sites
DetectBindingSitesMotif function detects binding sites from provided sequence motif information. Here, Artificial ChIP-seq data for FUR transcription factor in Escherichia coli was generated in two conditions fe2+ and dpd inspired by (Seo et al. 2014). In the following examples, artificial sequence motif locations have been and provided as a bed file called ‘FurMotifs.bed’. The alignment of ChIP-seq short reads is the other input of this function. The alignment files can be passed to the function as bam or bed files. In the following examples both IP and background alignment files have been passed as single-end bed files to this function.
# FUR candidate motifs in NC_000913 E. coli
FurMotifs = system.file("extdata", "FurMotifs.bed", package="Motif2Site")
# ChIP-seq FUR fe datasets binding sites from user provided bed file
# ChIP-seq datasets in bed single end format
IPFe <- c(system.file("extdata", "FUR_fe1.bed", package="Motif2Site"),
system.file("extdata", "FUR_fe2.bed", package="Motif2Site"))
Inputs <- c(system.file("extdata", "Input1.bed", package="Motif2Site"),
system.file("extdata", "Input2.bed", package="Motif2Site"))
FURfeBedInputStats <-
DetectBindingSitesBed(BedFile=FurMotifs,
IPfiles=IPFe,
BackgroundFiles=Inputs,
genome="Ecoli",
genomeBuild="20080805",
DB="NCBI",
expName="FUR_Fe_BedInput",
format="BEDSE"
)
FURfeBedInputStats## $FRiP
## [1] 0.3668289 0.1997165
##
## $sequencingStatitic
## nonSequenced underBinding overBinding
## 1 0.02293399 0 0.005022041
##
## $motifStatistics
## skewnessTestRejected decompositionRejected accepted
## 1 98 57 194# ChIP-seq FUR dpd datasets binding sites from user provided bed file
# ChIP-seq datasets in bed single end format
IPDpd <- c(system.file("extdata", "FUR_dpd1.bed", package="Motif2Site"),
system.file("extdata", "FUR_dpd2.bed", package="Motif2Site"))
FURdpdBedInputStats <-
DetectBindingSitesBed(BedFile=FurMotifs,
IPfiles=IPDpd,
BackgroundFiles=Inputs,
genome="Ecoli",
genomeBuild="20080805",
DB="NCBI",
expName="FUR_Dpd_BedInput",
format="BEDSE"
)
FURdpdBedInputStats## $FRiP
## [1] 0.2826046 0.2823270
##
## $sequencingStatitic
## nonSequenced underBinding overBinding
## 1 0.02254338 0 0.005022041
##
## $motifStatistics
## skewnessTestRejected decompositionRejected accepted
## 1 63 59 192DetectBindingSitesMotif function also works with DNA string motifs. In the following example, FUR binding sites in fe2+ condition are detected from ‘GWWTGAGAA’ with 1-mismatch motif. The dataset is generated only for ‘NC_000913’ build. Therefore, in this example the coordinates of this regions are provided as a ‘GivenRegion’ field. Providing this field ensures that binding sites are only detected in the given regions. This will accelerate the peak calling, and also improve the prediction accuracy.
# Granages region for motif search
NC_000913_Coordiante <- GRanges(seqnames=Rle("NC_000913"),
ranges = IRanges(1, 4639675))
# ChIP-seq FUR fe datasets binding sites from user provided string motif
# ChIP-seq datasets in bed single end format
FURfeStringInputStats <-
DetectBindingSitesMotif(motif = "GWWTGAGAA",
mismatchNumber = 1,
IPfiles=IPFe,
BackgroundFiles=Inputs,
genome="Ecoli",
genomeBuild="20080805",
DB= "NCBI",
expName="FUR_Fe_StringInput",
format="BEDSE",
GivenRegion = NC_000913_Coordiante
)
FURfeStringInputStats## $FRiP
## [1] 0.07620002 0.08057296
##
## $sequencingStatitic
## nonSequenced underBinding overBinding
## 1 0.02464332 0 0.003891051
##
## $motifStatistics
## skewnessTestRejected decompositionRejected accepted
## 1 3 2 365 Combining binding sites across experiments
recenterBindingSitesAcrossExperiments function combines the binding sites of different tissues or conditions into a single count matrix. In the FUR example, at the first step it combines fe2+ and dpd binding sites. At the next step, it recalculates the p-adjusted values. To ensure the high quality of the combined binding site set, an stringent cross-experiment FDR cutoff is applied (default 0.001). The accepted binding sites should fullfill this cutoff at least in one experiment. Another FDR cutoff value (default 0.05) is used to assign binding or non-binding labels to each binding site for each experiment.
# Combine FUR binding sites from bed input into one table
corMAT <- recenterBindingSitesAcrossExperiments(
expLocations=c("FUR_Fe_BedInput","FUR_Dpd_BedInput"),
experimentNames=c("FUR_Fe","FUR_Dpd"),
expName="combinedFUR"
)
corMAT## FUR_Fe FUR_Fe FUR_Dpd FUR_Dpd
## FUR_Fe 1.0000000 0.4062142 -0.4255622 -0.3734487
## FUR_Fe 0.4062142 1.0000000 -0.4638111 -0.4260953
## FUR_Dpd -0.4255622 -0.4638111 1.0000000 0.7113370
## FUR_Dpd -0.3734487 -0.4260953 0.7113370 1.0000000FurTable <-
read.table(file.path("combinedFUR","CombinedMatrix"),
header = TRUE,
check.names = FALSE
)
FurBindingTotal <-
GRanges(seqnames=Rle(FurTable[,1]),
ranges = IRanges(FurTable[,2], FurTable[,3])
)
FurFe <- FurBindingTotal[which((FurTable$FUR_Fe_binding =="Binding")==TRUE)]
FurDpd <- FurBindingTotal[which((FurTable$FUR_Dpd_binding =="Binding")==TRUE)]
findOverlaps(FurFe,FurDpd) ## Hits object with 108 hits and 0 metadata columns:
## queryHits subjectHits
## <integer> <integer>
## [1] 1 1
## [2] 7 5
## [3] 9 7
## [4] 10 8
## [5] 16 11
## ... ... ...
## [104] 204 189
## [105] 205 190
## [106] 206 191
## [107] 207 192
## [108] 210 194
## -------
## queryLength: 213 / subjectLength: 195pairwisDifferential function uses edgeR TMM normalization and GLM test to detect differential binding sites across two experiments. In the following example it takes combined FUR count matrix and detect differential binding sites across fe2+ and dpd.
# Differential binding sites across FUR conditions fe vs dpd
diffFUR <- pairwisDifferential(tableOfCountsDir="combinedFUR",
exp1="FUR_Fe",
exp2="FUR_Dpd",
FDRcutoff = 0.05,
logFCcuttoff = 1
)
FeUp <- diffFUR[[1]]
DpdUp <- diffFUR[[2]]
TotalComparison <- diffFUR[[3]]
head(TotalComparison)## seqnames start end logFC logCPM LR PValue
## 1 NC_000913 8346 8355 0.3979131 11.73185 0.02575615 0.8724972918
## 2 NC_000913 58862 58874 7.7314291 11.55399 13.25391337 0.0002720114
## 3 NC_000913 70522 70531 5.5034730 11.83521 9.75140375 0.0017918583
## 4 NC_000913 111527 111535 5.0683186 11.46908 7.65483916 0.0056620239
## 5 NC_000913 146492 146501 -6.8220279 10.79775 3.26337403 0.0708432210
## 6 NC_000913 167150 167158 -4.4357274 11.39174 5.16263159 0.02307790426 Session Info
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] BSgenome.Ecoli.NCBI.20080805_1.3.1000
## [2] BSgenome.Scerevisiae.UCSC.sacCer3_1.4.0
## [3] BSgenome_1.68.0
## [4] rtracklayer_1.60.0
## [5] Biostrings_2.68.0
## [6] XVector_0.40.0
## [7] Motif2Site_1.4.0
## [8] GenomicRanges_1.52.0
## [9] GenomeInfoDb_1.36.0
## [10] IRanges_2.34.0
## [11] S4Vectors_0.38.0
## [12] BiocGenerics_0.46.0
## [13] BiocStyle_2.28.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.0 viridisLite_0.4.1
## [3] dplyr_1.1.2 mixtools_2.0.0
## [5] bitops_1.0-7 fastmap_1.1.1
## [7] RCurl_1.98-1.12 lazyeval_0.2.2
## [9] GenomicAlignments_1.36.0 XML_3.99-0.14
## [11] digest_0.6.31 lifecycle_1.0.3
## [13] survival_3.5-5 kernlab_0.9-32
## [15] magrittr_2.0.3 compiler_4.3.0
## [17] rlang_1.1.0 sass_0.4.5
## [19] tools_4.3.0 utf8_1.2.3
## [21] yaml_2.3.7 data.table_1.14.8
## [23] knitr_1.42 htmlwidgets_1.6.2
## [25] DelayedArray_0.26.0 BiocParallel_1.34.0
## [27] purrr_1.0.1 grid_4.3.0
## [29] fansi_1.0.4 colorspace_2.1-0
## [31] edgeR_3.42.0 ggplot2_3.4.2
## [33] scales_1.2.1 MASS_7.3-59
## [35] SummarizedExperiment_1.30.0 cli_3.6.1
## [37] rmarkdown_2.21 crayon_1.5.2
## [39] generics_0.1.3 httr_1.4.5
## [41] rjson_0.2.21 cachem_1.0.7
## [43] splines_4.3.0 zlibbioc_1.46.0
## [45] parallel_4.3.0 BiocManager_1.30.20
## [47] restfulr_0.0.15 matrixStats_0.63.0
## [49] vctrs_0.6.2 Matrix_1.5-4
## [51] jsonlite_1.8.4 bookdown_0.33
## [53] magick_2.7.4 locfit_1.5-9.7
## [55] limma_3.56.0 plotly_4.10.1
## [57] jquerylib_0.1.4 tidyr_1.3.0
## [59] glue_1.6.2 codetools_0.2-19
## [61] gtable_0.3.3 BiocIO_1.10.0
## [63] munsell_0.5.0 tibble_3.2.1
## [65] pillar_1.9.0 htmltools_0.5.5
## [67] GenomeInfoDbData_1.2.10 R6_2.5.1
## [69] evaluate_0.20 lattice_0.21-8
## [71] Biobase_2.60.0 highr_0.10
## [73] Rsamtools_2.16.0 segmented_1.6-4
## [75] bslib_0.4.2 Rcpp_1.0.10
## [77] nlme_3.1-162 xfun_0.39
## [79] MatrixGenerics_1.12.0 pkgconfig_2.0.3References
Benaglia, T., D. Chauveau, D. R. Hunter, and D. Young. 2009. “mixtools: An R Package for Analyzing Finite Mixture Models.” Journal of Statistical Software 32 (6): 1–29. http://www.jstatsoft.org/v32/i06/.
Lawrence, M., W. Huber, H. Pagès, P. Aboyoun, M. Carlson, R. Gentleman, M. Morgan, and V. Carey. 2013. “Software for Computing and Annotating Genomic Ranges.” PLoS Computational Biology 9 (8). https://doi.org/10.1371/journal.pcbi.1003118.
Morgan, M., H. Pagès, V. Obenchain, and N. Hayden. 2020. Rsamtools: Binary Alignment (Bam), Fasta, Variant Call (Bcf), and Tabix File Import. http://bioconductor.org/packages/Rsamtools.
Pagès, H., P. Aboyoun, R. Gentleman, and S. DebRoy. 2019. Biostrings: Efficient Manipulation of Biological Strings.
Pagès, Hervé. 2019. BSgenome: Software Infrastructure for Efficient Representation of Full Genomes and Their Snps.
Robinson, M. D., D. J. McCarthy, and G. K. Smyth. 2010. “EdgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics 26 (1): 139–40. https://doi.org/10.1093/bioinformatics/btp616.
Seo, S. W., D. Kim, H. Latif, and others. 2014. “Deciphering Fur Transcriptional Regulatory Network Highlights Its Complex Role Beyond Iron Metabolism in Escherichia Coli.” Nature Communications 5 (1). https://doi.org/10.1038/ncomms5910.