Cross-species de novo identification of cis-regulatory modules with GibbsModule: application to gene regulation in embryonic stem cells

Abstract

We introduce the GibbsModule algorithm for de novo detection of cis-regulatory motifs and modules in eukaryote genomes. GibbsModule models the coexpressed genes within one species as sharing a core cis-regulatory motif and each homologous gene group as sharing a homologous cis-regulatory module (CRM), characterized by a similar composition of motifs. Without using a predetermined alignment result, GibbsModule iteratively updates the core motif shared by coexpressed genes and traces the homologous CRMs that contain the core motif. GibbsModule achieved substantial improvements in both precision and recall as compared with peer algorithms on a number of synthetic and real data sets. Applying GibbsModule to analyze the binding regions of the Krüppel-like factor (KLF) transcription factor in embryonic stem cells (ESCs), we discovered a motif that differs from a previously published KLF motif identified by a SELEX experiment, but the new motif is consistent with mutagenesis analysis. The SOX2 motif was found to be a collaborating motif to the KLF motif in ESCs. We used quantitative chromatin immunoprecipitation (ChIP) analysis to test whether GibbsModule could distinguish functional and nonfunctional binding sites. All seven tested binding sites in GibbsModule-predicted CRMs had higher ChIP signals as compared with the other seven tested binding sites located outside of predicted CRMs. GibbsModule is available at (http://biocomp.bioen.uiuc.edu/GibbsModule).

Figure 1.

Three information sources for de novo identification of cis-regulatory motifs. The three circles represent three information sources that can be utilized for motif and CRM finding. A tool enclosed in a circle indicates that this tool utilizes that information.

Figure 2.

Motifs and CRMs in coexpressed genes and their homologous genes. A, B, C, and D represent coexpressed genes in one species. A′ and A′′ represent the homologous genes to A in two other species, and so on. (X, O, #) TFBSs for different transcription factors.

Figure 3.

GibbsModule workflow. In Step 1, a random PSWM is initialized. Steps 2-5 are the iterative steps. In Step 2, N candidate binding sites are sampled from every homologous sequence using the same PSWM. In this example, three candidate-binding sites are sampled on each sequence (N = 3). Every sampled binding site defines a candidate CRM, which includes the binding site itself and 100 bp of flanking region on each side. These candidate CRMs are marked 1, 2, 3 on the target sequence, and 1′, 2′, 3′ and 1″, 2″, 3″ on the sequences of two assisting species. In Step 3, Module-Alignment is applied to every candidate CRM on the target sequence and every CRM on the assisting sequences. In the example, the alignments are applied to CRM pairs of (1, 1′), (1, 2′), (1, 3′), (2, 1′), (2, 2′), . . . , (3, 1″), (3, 2″), and (3, 3″). In Steps 3 and 4, a most conserved CRM on the target sequence is picked up by arg ((score(n,n')) + (score(n,n′))), where n, n′, and n″ are indicators of candidate CRMs in homologous sequences SeqA, orthA1, and orthA2, respectively. (X, O, #) Other motifs close to the core motif within a CRM. In Step 5, a new PSWM is calculated from the core motifs in the most conserved CRMs.

Figure 4.

Module-Alignment. (A) An illustration of three pairs of orthologous CRMs: All three CRM pairs consist of TFBSs generated from the same motifs (squares, ellipses, and triangles). (1) Orthologous CRMs with conserved number, order, and distances. (2–4) Orthologous CRMs with different distances, order, and number of TFBSs. (B) Workflow of Module-Alignment: Module-Alignment iteratively performs local alignment and masks out conserved regions. The mutations and gaps between the alignable segments on the upper and lower sequences incur severe penalty so that Smith-Waterman can detect only the best local alignment in the first row. The conservation score from local alignment is the score of the best local alignment. However, the conservation score from Module-Alignment is the sum of the two local alignments scores from the two alignable sequence segments. Module-alignment is not designed to align any two sequences with any arbitrary lengths. Its input sequences should be potentially orthologous CRMs with lengths of several dozen to several hundred base pairs.

Figure 5.

Motifs derived from Klf ChIP-chip and mutagenesis analysis. SOX2 motif (A) and KLF motif (B) found by GibbsModule and MEME from KLF ChIP-chip data. (C) SOX2 motif in TRANSFAC.

Figure 6.

ChIP signals of predicted CRM and non-CRMs that contain putative POU5F1-binding sites.