INTRODUCTION {#SEC1}
============

Nowadays, DNA methylation is recognised as a key component in numerous biological mechanisms such as gene regulation, RNA splicing, nucleosome positioning and transcription factors (TF) binding dynamics ([@B1]).  Until recently, it has been well accepted that the writing of a methyl group into a cytosine can prevent TF binding events.  Indeed, methylation within the binding site of the insulator-binding protein CTCF abrogates its DNA binding ability, resulting in imprinted expression of H19 and Igf2 ([@B2]).  Similarly, AP-1 and SP1 interactions with DNA can be inhibited by methylating the CpGs adjacent to their respective DNA binding sites ([@B3], [@B4]).  Nonetheless, this rule is not universal.  Accumulating evidence during the last few years has shown that numerous TFs require methylated DNA to bind to their motifs, whereas other TFs bind to DNA regardless the methylation status ([@B5]).  For example, Hu *et* al.  characterised 41 methylated-DNA binding proteins among 1321 tested TFs using a protein microarray-based approach ([@B8]).  More recently, by using methylation-sensitive systematic evolution of ligands by exponential enrichment (SELEX), Yin *et*al.  segregated 519 TFs into five groups based on their binding strength to methylated DNA sequences ([@B9]).  The first group includes 175 TFs whose binding activities are enhanced by methylated DNA; the second group contains 117 TFs whose binding abilities are compromised by DNA methylation; the third group is composed of 33 TFs that bind to DNA regardless the methylation status; the fourth group encompasses 25 TFs whose binding activities are differently affected by DNA methylation depending on different sub-motifs or different CpG sites in the same motif; and finally, the last group of TFs (169) does not contain any CpG site within their motifs.  These systematic *in vitro* studies have unveiled the disparate mechanisms involved in TFs/DNA interactions with DNA methylation being a key leverage in this process.  However, it has been shown that only a small fraction of TF binding sites (TFBS) found *in vitro* can be validated *in vivo* ([@B10]), due to a more complex and crowded nuclear environment involving the chromatin structure and non-DNA molecules such as RNAs, histones and non-histone proteins ([@B11]).  Therefore, further efforts are required to characterise the methylation status of TFBS in the cell context. 

Most of the current TFBS databases such as JASPAR ([@B12]), HOCOMOCO ([@B13]) and GTRD ([@B14]) provide TFBS motifs derived from *in vivo* ChIP-seq datasets.  Nevertheless, none of these databases have combined TFBS motifs with a precise depiction of their respective DNA methylation profiles yet.  This is the reason why we have developed MethMotif (Methylation in Motif), a TFBS database that combines both position weight matrices (PWM), and DNA methylation profiles, in a cell specific manner. 

MATERIALS AND METHODS {#SEC2}
=====================

Data collection and pre-processing {#SEC2-1}
----------------------------------

Raw ChIP-seq datasets (fastq files) in different cell types were downloaded from ENCODE consortium and GEO databases (Figure [1A](#F1){ref-type="fig"}, Table [1](#tbl1){ref-type="table"}).  After read quality checking with FastQC (<https://www.bioinformatics.babraham.ac.uk/projects/fastqc/>) and read trimming when required, using Trimmomatic ([@B15]), short reads were aligned to the human genome assembly hg38, using STAR ([@B16]) with the spliced function blocked, for its speedy and accurate aligning performance (BioRxiv: <https://doi.org/10.1101/053686>).  Relying upon the Irreproducibility Discovery Rate (IDR) pipeline ([@B17]), TF genomic occupancy regions with high consistency in replicates were called by MACS2 ([@B18]) using default parameters. 

![ Construction of MethMotif database.  Integrative analysis of ChIP-seq and WGBS datasets, downloaded from ENCODE and GEO databases, allows the profiling of the DNA methylation landscapes surrounding transcription factor binding regions in various cell types (**A**).  Firstly, the methylation levels (beta scores or methylation scores) of CpGs located within all ChIP-ed protein peak regions are captured over the 200 bp surrounding peak summits.  The distribution of the corresponding CpG methylation scores is then profiled in a heatmap for each cell type, to present the DNA methylation levels surrounding peak regions across all ChIP-ed proteins (methylation score less than 10% is defined as homogenously unmethylated, while methylation score more than 90% is regarded as homogenously hypermethylated).  These heatmaps are accessible from the 'Explore' section of the MethMotif website (**B**).  Finally, the direct binding motifs of sequence-specific TFs are identified.  The DNA methylation level within TFBS is captured and shown in a MethMotif logo.  MethMotif logos are collected in 'MethMotif database' section of the MethMotif website (**C**).  The methylation levels within each logo can be displayed according to three states: (i) all methylation levels compiled, (ii) methylated only and, (iii) unmethylated only (**D**).](gky1005fig1){#F1}

###### 

ChIP-seq datasets used in MethMotif database

  Cell ID   Organism   Cell type/tissue   Number of ChIP-seq experiments   Source              Download date
  --------- ---------- ------------------ -------------------------------- ------------------- ----------------
  HeLa-S3   Human      Cervix             57                               ENCODE              20 August 2017
  HEK293    Human      Kidney             199                              ENCODE              12 March 2018
  IMR-90    Human      Lung               15                               ENCODE              26 March 2018
  SK-N-SH   Human      Brain              30                               ENCODE              26 March 2018
  A549      Human      Lung               42                               ENCODE              26 March 2018
  K562      Human      Blood              279                              ENCODE              2 April 2018
  HepG2     Human      Liver              138                              137-ENCODE, 1-GEO   7 April 2018
  GM12878   Human      Blood              143                              ENCODE              12 April 2018
  MCF-7     Human      Breast             99                               ENCODE              13 April 2018
  H1-hESC   Human      Stem cell          65                               ENCODE              16 April 2018
  HCT116    Human      Colorectum         22                               ENCODE              16 April 2018

We took advantage of Whole Genome Bisulfite Sequencing (WGBS) datasets to measure the genome-wide DNA methylation levels in each cell type at the CpG resolution.  WGBS datasets were obtained from ENCODE consortium and GEO databases (Table [2](#tbl2){ref-type="table"}).  Raw datasets underwent read quality check using FastQC, as well as quality and adapter trimming using Trim Galore.  High quality short reads were aligned to hg38 using Bowtie2 with Bismark (v0.16.3) and default parameters ([@B19],[@B20]).  Methylation states in CpG context were extracted with the Bismark methylation extractor module.  Whenever their correlation coefficients were confirmed over 0.8, as recommended by ENCODE standards, WGBS datasets with biological replicates were merged to increase the sequencing coverage using methylKit ([@B21]). 

###### 

WGBS datasets used in MethMotif database

  Cell ID   Organism   Cell type/tissue   Source ID     Release date
  --------- ---------- ------------------ ------------- ------------------
  HeLa-S3   Human      Cervix             GSM2175341    30 January 2017
  HEK293    Human      Kidney             GSM1254259    11 December 2015
  IMR-90    Human      Lung               ENCSR888FON   31 July 2013
  SK-N-SH   Human      Brain              ENCSR145HNT   13 December 2017
  A549      Human      Lung               ENCSR481JIW   4 December 2017
  K562      Human      Blood              ENCSR765JPC   22 March 2016
  HepG2     Human      Liver              ENCSR881XOU   13 October 2015
  GM12878   Human      Blood              ENCSR890UQO   23 February 2016
  MCF-7     Human      Breast             GSM1328112    3 July 2014
  H1-hESC   Human      Stem cell          ENCSR617FKV   13 October 2015
  HCT116    Human      Colorectum         GSM3317488    10 August 2018

In house WGBS data generation {#SEC2-2}
-----------------------------

We complemented the existing published datasets from ENCODE and GEO by generating a WGBS of the HCT116 cell line (human colorectal carcinoma).  HCT116 cells were grown in RPMI 1640 medium with 10% FBS.  The DNA was extracted using Qiagen DNeasy Blood & Tissue Kits.  After sample DNA testing, negative control DNAs were added into the initial DNA pool, and then fragmented into a 200--300 bp range using Covaris S220.  Next, terminal repairing, A-ligation, methylation sequencing adapters ligation were performed on DNA fragments.  Then, DNA fragments underwent a bisulfite treatment using the EZ DNA Methylation Gold Kit from Zymo Research.  Library concentration was firstly quantified by Qubit 2.0, before a dilution to 1 ng/ul followed by insert size verification using the Agilent Bioanalyzer 2100 and qPCR.  Finally, the library was paired-end sequenced using the Illumina HiSeq 4000.  Raw and processed data are available through the GEO portal under the accession number GSE118030.  We will use this GEO SuperSeries ID for sharing our future in-house datasets. 

Profiling the DNA methylation landscapes surrounding ENCODE ChIP-ed protein binding regions {#SEC2-3}
-------------------------------------------------------------------------------------------

Integrative analyses of ChIP-seq and WGBS datasets were performed to profile DNA methylation landscapes surrounding the genomic occupancy regions across all ENCODE ChIP-ed proteins.  For each DNA binding protein, methylation scores (beta scores) ([@B22]) of CpGs within ±100 bp surrounding peak summits were collected using the intersectBed module of bedtools ([@B23]).  For each cell type, the distributions of the collected CpG methylation scores across all ChIP-ed proteins were shown in a heatmap, and proteins were clustered into different groups based on their DNA methylation landscapes surrounding their binding regions, using hierarchical clustering with the Euclidean distance.  This representation allows us to classify ENCODE ChIP-ed proteins according to the DNA methylation levels of their respective binding loci (Figure [1B](#F1){ref-type="fig"}). 

Characterisation of DNA methylation levels within TFBS {#SEC2-4}
------------------------------------------------------

We employed the MEME-ChIP package with the default parameters to identify TF motifs in ±100 bp surrounding TF peak summits ([@B24]).  Since the actual binding sites for the ChIP-ed TFs are prone to central enrichment and maximum probability in occurrence at peak centres ([@B25]), each TF motif was determined based on its centrally enriched propensity profiled by CentriMo from the MEME-ChIP package.  Then, the exact TF binding sites were localised using the FIMO module, while the methylation states of CpGs within the binding motif were assessed based on the WGBS dataset using the intersectBed module from bedtools ([@B23]). 

In order to intuitively exhibit the DNA methylation levels within all TFBSs, we adopted a novel MethMotif logo combining the classical motif logo commonly used, with a cumulative bar chart describing the methylation level of each CpG present in the motif (Figure [1C](#F1){ref-type="fig"}, [D](#F1){ref-type="fig"}).  CpG methylation scores were segregated into three interval groups: (i) methylation scores \<10% (i.e. homogenously hypomethylated); (ii) methylation scores \>90% (i.e. homogenously hypermethylated) and (iii) methylation scores ranging from 10% to 90% (i.e. heterogeneously methylated).  We utilised WebLogo3 ([@B26]) to generate the motif logo, while the bar chart above the motif logo describing the methylation score was generated using a custom R script.  While the usage of stringent methylation score thresholds prevents mis-scoring DNA methylation in TFBS due to cell heterogeneity, a more direct genome-wide investigation by sequential ChIP-bisulfite sequencing (ChIP-BS-seq) ([@B27]) will be necessary to validate these methyl-TFBS characteristics. 

Web interface {#SEC2-5}
-------------

A web-based interface has been established at <https://bioinfo-csi.nus.edu.sg/methmotif> allowing an intuitive exploration of our datasets and facilitating queries to the database.  Users can freely access MethMotif *via* these three modes: '*Motif database direct query*'.  Like many TFBS webservers, MethMotif includes a search engine that allows users to query any DNA binding protein present in our database using either its official gene name or its MethMotif ID (see an example of MethMotif ID in Figure [1C](#F1){ref-type="fig"}).  Queries can be refined by selecting a given cell line or tissue of interest.  Query results include the MethMotif logo, information about the cell type, ChIP-seq and WGBS assays, as well as motif distributions, motif locations and downloadable processed files (PWM and methylation matrices), along with links to JASPAR 2018 and HOCOMOCO (Figure [1C](#F1){ref-type="fig"} and [D](#F1){ref-type="fig"}).'*Explore*'.  Users can intuitively explore all DNA binding proteins available in the released cell types in our database *via* dynamic heatmaps that classify each protein according to its CpG methylation pattern in a ±100 bp window surrounding ChIP-seq peak summits.  Details about any DNA binding protein of interest can be easily obtained by hovering over the corresponding heatmap\'s row (Figure [1B](#F1){ref-type="fig"}).  To date, we have made available 11 dynamic heatmaps corresponding to the 11 analysed cell types.  The heatmap dynamic interface takes advantage of InCHlib, an open source interactive JavaScript library ([@B28]).'*Batch query*'.  MethMotif allows users to analyse the occurrences of TFBSs along with their respective methylation states across a given list of genome loci.  This functionality is particularly useful to characterise binding sites for co-factors of a protein of interest, along with their respective DNA methylation levels.  Here, MethMotif generates *de novo* DNA binding motifs based only on the sequences that overlap with the input genomic regions.  TFBS methylation profiles are represented as beeswarm boxplots, where each dot represents the methylation level of a CpG from a given binding site in the genome (Figure [2](#F2){ref-type="fig"}). 

![ Workflow of MethMotif Batch Query.  MethMotif Batch Query is available via the MethMotif website, which allows users to study the occurrences of TFBSs along with DNA methylation states in a given list of genomic loci.  Users can upload the coordinates of the regions of their interest in BED format (**A**).  MethMotif database will be queried (**B**) and the presence of TFBSs together with respective DNA methylation levels is then analysed using the MethMotif Enrichment Tool (**C**) in the given regions.  If loci from the input regions overlap with any TFBS present in the MethMotif database, the Batch Query tool will output these overlapped loci along with their respective TFBS DNA methylation information *via* a beeswarm boxplot (where each dot represents the methylation level of a CpG site within the TFBS) and a MethMotif logo generated *de novo*(**D**).](gky1005fig2){#F2}

Batch query module {#SEC2-6}
------------------

The functionality of batch query module implemented in MethMotif website is achieved by mapping the given list of genome loci to all TFBS genomic coordinates, together with their respective CpG methylation scores, in MethMotif database using in-house python scripts.  The visible results, namely beeswarm boxplots and *de novo* generated MethMotif logo, are produced using WebLogo3 ([@B26]) and in-house R scripts with the grImport package ([@B29]). 

Database maintenance {#SEC2-7}
--------------------

Established ChIP-seq working standards and guidelines in ENCODE consortium foster reliable ChIP-seq data ([@B17]).  Since the beginning of the project, MethMotif database has been regularly updated with ChIP-seq datasets from ENCODE.  We plan to further enrich the MethMotif database by expanding cell types using newly generated WGBS datasets from ENCODE, GEO and in-house experiments.  Furthermore, we intend to take advantage of the Gene Transcription Regulation Database (GTRD), which is a comprehensive compendium of ChIP-seq experiments covering the existing multiple ChIP-seq databases ([@B14]), to expand the coverage of TFs in MethMotif. 

RESULTS {#SEC3}
=======

MethMotif classifies DNA-interacting proteins according to DNA methylation profiles {#SEC3-1}
-----------------------------------------------------------------------------------

The integration of WGBS with ChIP-seq datasets allowed us to profile the DNA methylation landscapes surrounding binding loci of DNA-bound proteins *in vivo*, at the genome scale.  We employed this strategy to capture the CpG methylation patterns in a ±100 bp window surrounding the peak summits of ENCODE ChIP-ed proteins with high reproducibility in different cell types (Figure [1A](#F1){ref-type="fig"}).  As described previously ([@B9]), DNA-binding proteins exhibit distinct DNA methylation levels around their occupancy regions and accordingly can be classified into different sub-groups (Figure [1B](#F1){ref-type="fig"}, [Supplementary Figure S1](#sup1){ref-type="supplementary-material"}).  In myelogenous leukaemia K562 cells, we were able to classify 250 ChIP-ed proteins based on DNA methylation landscapes around binding regions into three groups (I, high; II, medium and III, low methylated DNA binding affinity) and one outlier (very high methylated binding affinity --- ZBTB33) ([Supplementary Figure S1](#sup1){ref-type="supplementary-material"}).  Notably, our ZBTB33 epigenetic binding profile corroborates previous *in vivo* studies performed with a limited number of loci ([@B30],[@B31]), and generalises the ability of ZBTB33 to bind methylated DNA at the genome scale.  Furthermore, group III includes TFs known to bind specifically to non-methylated DNA such as RUNX1 ([@B32]) and SP1 ([@B4]).  As shown in the MethMotif website, this distinct binding propensity to methyl-CpGs among different DNA-binding proteins is not a unique observation in K562, but a widespread phenomenon across different cell types. 

MethMotif introduces a new dimension to Position Weight Matrix (PWM) {#SEC3-2}
--------------------------------------------------------------------

A key feature of TFs is their ability to bind DNA by the recognition of a specific DNA motif.  However, ChIP-seq peak regions are not necessarily associated with direct binding sites ([@B25]) and DNA methylation profiling around all aforementioned peaks may not reflect the influences of DNA methylation on the recruitment of TFs onto DNA sequences.  Therefore, we inferred TFBSs by selecting genomic loci that contained a centrally enriched DNA motif surrounding ChIP-seq centres ([@B25]).  TFBS was characterised using the MEME-ChIP package ([@B24]) in ±100 bp regions surrounding peak centres discovered by the IDR pipeline.  This strategy allowed us to capture the DNA methylation landscape within each TFBS motif.  In order to intuitively illustrate TFBS coupled with DNA methylation, we introduced a novel graphical representation combining a motif-logo ([@B26]) with cumulative bar charts displaying three different DNA methylation levels: (i) homogenously low (less than 10% sequenced CpGs of a given site are methylated); (ii) homogenously high (more than 90% of sequenced CpGs of a given site are methylated) and (iii) heterogeneously methylated (10--90% sequenced CpGs of a given site are methylated) (Figure [1C](#F1){ref-type="fig"} and [D](#F1){ref-type="fig"}).  The MethMotif logo represents these three DNA methylation levels as a cumulative bar chart where the blue, orange and green colours denote the proportions of CpGs falling in the group 1, 2 and 3, respectively.  Additionally, we display the number of total CpGs covered by WGBS on the top of each bar.  We believe that adding an epigenetic dimension to the sequence TFBS information will improve our understanding of TF binding dynamics. 

Numerous TFs exhibit distinct binding site sequences and/or distinct binding site methylation profiles across different cell types {#SEC3-3}
----------------------------------------------------------------------------------------------------------------------------------

In contrast to pre-existing TFBS databases, MethMotif records TFBSs computed in a cell specific manner.  Our meta-analysis shows that numerous TFs exhibit either distinct TFBSs, or distinct DNA methylation profiles, or both, depending on cell types.  Notably, when comparing ChIP-seq performed in K562 and GM12878 cells, we found that the basic leucine zipper factor NFE2 has the ability to bind distinct DNA motifs (Figure [3A](#F3){ref-type="fig"}).  Indeed, in K562, NFE2 recognises an extended AP-1-like motif, (T/C/A)GCTGA(C/G)TCA(T/C), while the motif identified in GM12878 is similar to the USF binding site ([Supplementary Figure S2A](#sup1){ref-type="supplementary-material"}).  Interestingly, both motifs are supported by the literature.  Specifically, NFE2 has been shown to bind an AP-1-like motif as a heterodimer with MAF family members, such as MAFF and MAFK ([@B33]); and, in another context, NFE2 was characterised to belong to a complex including USF ([@B34]).  We were able to corroborate these facts by using the MethMotif batch query module with all NFE2 binding sites in K562 and GM12878 cells ([Supplementary Figure S2B](#sup1){ref-type="supplementary-material"}). 

![ Distinct motifs and DNA methylation levels in motifs across cell types.  Examples of MethMotif logos from the MethMotif database show that TF binding propensities can be altered in terms of motifs (**A** and **B**), DNA methylation levels inside the motifs (**C**), or even both (**D**), across different cell types.](gky1005fig3){#F3}

In addition, we found that the length of spacers of leucine zipper dimers, such as JUN, varies with cell types.  For example, in hepatocyte carcinoma HepG2 and embryonic stem cell line H1-hESC, the main motif spacer of JUN is a CG dinucleotide, while it amounts to only a single base, either C or G, in breast cancer MCF-7, lung adenocarcinoma A549 and myelogenous leukaemia K562 cells (Figure [3B](#F3){ref-type="fig"}).  Interestingly, the UniProbe TFBS database, which records *in vitro* characterised DNA binding sites, contains another AP-1 leucine zipper TF, JUNDM2 (Jun dimerisation protein 2), that recognises motifs similar to JUN with a spacer that can vary from 1 to 2 nucleotides (UniProbe accession number: UP00103) ([@B35]).  Such variations in spacers\' size can be explained by different dimerisations of JUN or JUNDM2 with other TFs ([@B36]), or by post-translational modifications and modified secondary structure. 

We further noticed that some TFs have a conserved DNA motif with distinct DNA methylation profiles across cell types.  This is the case for the leucine zipper MAFF.  In K562, CpGs sites composing the MAFF motif are prone to hypomethylation, while in HeLa-S3 and HepG2, we observed a significant number of motifs with methylated CpGs (Figure [3C](#F3){ref-type="fig"}).  To the best of our knowledge, MAFF was not previously described to play different roles across these cell lines.  However, ChIP-seq results clearly show that MAFF targets distinct loci, and consequently regulates distinct biological functions ([Supplementary Figure S3](#sup1){ref-type="supplementary-material"}). 

Finally, a group of TFs such as the zinc finger REST revealed distinct binding profiles at the DNA sequence and DNA methylation levels.  Remarkably, in this specific example, changes in DNA methylation levels correlate with changes in DNA sequence.  Indeed, REST is known to bind a sequence of 6 and 7 conserved nucleotides linked by a 2-nucleotide spacer: CTGTCCNNGGTGCTG (cf.  JASPAR accession number: MA0138.1; HOCOMOCO accession number: REST_HUMAN.H11MO.0.A).  Although this motif is confirmed in A549, GM12878, K562 and H1-hESC cells, where REST binds unmethylated DNA, we observed that REST requires only the first part of its conserved motif (CTGTCC) to bind methylated DNA, as shown in HEK293, HeLa-S3, HCT116, MCF-7 and HepG2 cells (Figure [3D](#F3){ref-type="fig"}).  It is interesting to note that, in another context, non-CpG methylation facilitates REST binding events ([@B37]). 

MethMotif batch query: a case study {#SEC3-4}
-----------------------------------

MethMotif portal includes a batch query tool which allows users to study the occurrences of TFBSs and their respective DNA methylation states within a given list of genomic loci.  As a proof of concept, we ran MethMotif batch query using a list of CEBPB binding sites characterised by a ChIP-seq experiment in HCT116 cells ([Supplementary Data 1](#sup1){ref-type="supplementary-material"}).  As expected, the most enriched TFBS found by MethMotif is CEBPB (4846 overlapped occurrences, Figure [2D](#F2){ref-type="fig"} and [Supplementary Figure S4](#sup1){ref-type="supplementary-material"}).  Furthermore, MethMotif batch query brought to light 12 co-factors that overlap our input list from 1 to 300 co-occurrences ([Supplementary Figure S4](#sup1){ref-type="supplementary-material"}).  Results include known CEBPB co-factors such as ATF3 (300 co-occurrences) ([@B38]) and JUND (74 co-occurrences) ([@B39]).  Interestingly, only CEBPB binding loci harbour methylated CpG sites ([Supplementary Figure S4](#sup1){ref-type="supplementary-material"}).  This observation is in line with its pioneer factor function ([@B40]).  Indeed, here we can hypothesise that CEBPB is firstly recruited within methylated DNA across repressive chromatin, and then, once the DNA is demethylated and more accessible, co-factors are recruited. 

DISCUSSION {#SEC4}
==========

In contrast to the previously widely accepted assumption, DNA methylation is not necessarily associated with the inhibition of TF binding activities ([@B6],[@B8],[@B9]).  Interestingly, in some contexts, the recruitment of methyl-CpG specific binding TF can be abrogated by erasing DNA methylation using drugs such as 5-aza-2′-deoxycytidine (5-aza-dC) ([@B41]).  Therefore, DNA methylation becomes a key event for the recruitment of certain TFs.  As DNA methylation is also portrayed as an epigenetic mark of cellular memory ([@B42]), we can speculate that this memory requires the recruitment and/or interaction of specific factors that either inhibit or prime the chromatin.  This is the reason why the investigation of methyl-TFBS has to be done *in vivo* and in a cell specific manner.  MethMotif brought to light that, for given factors, TFBSs can vary according to the DNA methylation patterns. 

The DNA binding ability regardless of DNA methylation is a reminiscence of pioneer factors, which form a specific class of TFs required and sufficient to trigger enhancer competency ([@B43]).  MethMotif is able to detect this capacity for some pioneer factors such as CEBPA ([@B44]) and CEBPB ([@B40]).  Moreover, numerous zinc finger proteins have been characterised with methyl-DNA recognition ability ([@B45],[@B46]).  With MethMotif, we can systematically reveal that occupancy at methylated cytosine by zinc finger proteins is not a sporadic event but a widespread phenomenon.  Interestingly, the position weight matrices encompassing methylated CpGs are inclined to be specifically enriched in TpG dinucleotide at the same position, corroborating the fact that a methylated cytosine base structurally resembles a thymine ([@B47]) (Supplement Figure S5). 

Identification of methyl-cytosine binding TFs by MethMotif is the initial but crucial step to study the crosstalk between TFs and DNA methylation.  While we provide compelling evidence that numerous TFs bind differentially to methylated and non-methylated DNA, ChIP-BS-seq experiments ([@B27]) will be required to confirm these interactions.  Nonetheless, the current version of MethMotif remains a valuable resource for guiding researchers in the usage of ChIP-BS-seq as an alternative to ChIP-seq for studying methyl-binding TFs. 

In conclusion, MethMotif complements current TFBS databases by integrating DNA methylation as a new epigenetic dimension with TFBS, thereby providing a useful resource for the study of DNA methylation and TF recruitment dynamics. 

DATA AVAILABILITY {#SEC5}
=================

MethMotif is an academic and open source database under GPL licence and available at <https://bioinfo-csi.nus.edu.sg/methmotif>. 

Supplementary Material
======================

###### 

Click here for additional data file. 

The authors thank Morgane Thomas-Chollier and Daniel G. Tenen for useful discussions and comments during the preparation of this manuscript, as well as Mathanapriya Naidu for her proof-reading. 

*Authors contributions:* Q.X.X.L. and T.B. designed the platform; Q.X.X.L. programmed the website and performed all analysis; O.A. participated in the graphical user interface implementation; S.S. generated the HCT116 WGBS dataset; Q.X.X.L., D.T., S.J. and T.B. interpreted the data; Q.X.X.L. and T.B. wrote the manuscript.  T.B. directed the study. 

SUPPLEMENTARY DATA {#SEC6}
==================

[Supplementary Data](https://academic.oup.com/nar/article-lookup/doi/10.1093/nar/gky1005#supplementary-data) are available at NAR Online. 

FUNDING {#SEC7}
=======

Work in the T.B. laboratory is supported by the National Research Foundation, the Singapore Ministry of Education under its Centres of Excellence initiative, the National Medical Research Council of Singapore \[NMRC/BNIG/2035/2015\]; RNA Biology Center at the Cancer Science Institute of Singapore, NUS, as part of funding under the Singapore Ministry of Education\'s AcRF Tier 3 grants \[MOE2014-T3-1-006\].  T.B. and D.T. collaboration is supported by the Institut Français à Singapour \[Merlion Project grant number 6.10.14\].  S.J. was supported by grants from the National Research Foundation Singapore and the Singapore Ministry of Education under its Research Centers of Excellence initiative to the Cancer Science Institute of Singapore \[R-713-006-014-271\]; National Medical Research Council (NMRC) \[CBRG-NIG BNIG11nov001\]; Ministry of Education Academic Research Fund \[MOE AcRF Tier 1 T1-2012 Oct-04 and T1-2016 Apr-01\]; RNA Biology Center at CSI Singapore, NUS, from funding by the Singapore Ministry of Education\'s AcRF Tier 3 grants \[MOE2014-T3-1-006\].  Funding for open access charge: National Research Foundation, the Singapore Ministry of Education under its Centres of Excellence initiative and the Singapore Ministry of Education Academic Research Fund Tier 3 \[MOE2014-T3-1-006\]. 

*Conflict of interest statement*.  None declared.