cisRED: a database system for genome-scale computational discovery of regulatory elements.

We describe cisRED, a database for conserved regulatory elements that are identified and ranked by a genome-scale computational system (www.cisred.org). The database and high-throughput predictive pipeline are designed to address diverse target genomes in the context of rapidly evolving data resources and tools. Motifs are predicted in promoter regions using multiple discovery methods applied to sequence sets that include corresponding sequence regions from vertebrates. We estimate motif significance by applying discovery and post-processing methods to randomized sequence sets that are adaptively derived from target sequence sets, retain motifs with p-values below a threshold and identify groups of similar motifs and co-occurring motif patterns. The database offers information on atomic motifs, motif groups and patterns. It is web-accessible, and can be queried directly, downloaded or installed locally.

INTRODUCTION

Approaches for identifying transcriptional regulatory motifs computationally have been reviewed previously (1–3). Recently, progress has been made towards identifying a ‘comprehensive catalog’ of mammalian elements (4,5). As genome resources, data types and tools evolve, predictive approaches developed for such work can also be directed at a number of closely connected issues. Some of these include taking advantage in motif discovery of increasing numbers of genomes, including low-coverage sequences; quantifying the contributions to motif discovery of different genomes or sets of genomes; improving the predictive reliability of motifs by genome-scale clustering and co-occurrence; determining a best minimal set of motif discovery methods, probably in a discovery approach that uses multiple methods (6); and using coexpression and other functional data types.

In this report, we describe a new cisRED database that contains predictions for whole-genome discovery of regulatory elements in mammals and other eukaryotes. We also describe the predictive system behind the database, which uses genome-scale approaches to predict deeply conserved ab initio motifs and identifies groups of similar motifs and co-occurring patterns of motifs. Results are available in a web-accessible and downloadable MySQL database. The system is designed to be readily maintained and extended in the context of rapidly evolving resources, data types and tools.

DATA SOURCES AND PROCESSING

The system is outlined in Figure 1, and is described in detail at . The upstream section of the pipeline loads the database with significant discovered atomic motifs, and the section downstream of the database identifies groups of similar atomic motifs and co-occurring patterns of motifs. An atomic motif consists of a set of sequences, typically with a common length between 6 and 12 bp, members of which are present in a sequence region on the target species and in corresponding regions on other genomes.

Figure 1

Data processing system for high-throughput motif discovery, clustering, co-occurrence, annotation and performance assessment.

The system uses genome resources that are a combination of directly downloaded and processed sequences, annotations and relationships (e.g. orthology, coexpression and interactions). These are stored in an automatically updated local resource that holds a wide range of public and commercial databases.

Motif discovery was carried out in a search region based on a single major transcript for each gene. An input target sequence set consisted of a sequence from the target species, and corresponding sequences from homologous vertebrate genes. The following rules were used to assemble target sequence sets for cisRED 1.2e. We identified homologous genes by combining data from Compara (7), HomoloGene (8), Inparanoid (9) and KEGG (10). For each target human gene and each of its orthologues, we took the major Ensembl (7) transcript if its protein sequence was N-terminal complete. We further required that the human gene had an annotated Ensembl 5′-untranslated region (5′-UTR) that was at least 10 bp long. For each target gene's orthologue set, we required at least one of Dog, Mouse and Rat, and at least one of Chicken, Frog, Fugu, Tetraodon and Zebrafish. If the human 5′-UTR was <500 bp, we applied no UTR requirements to orthologues; however, if the human 5′-UTR was >500 bp, we required that all orthologues had annotated Ensembl 5′-UTRs that were at least 10 bp long. After this filtering, for any orthologous region that was missing from an input sequence set, we added a sequence region from the current UCSC (11,12) multiple sequence alignment. For a small subset of the target genes, corresponding regions from ENCODE-specific species were then added from the current ENCODE multiple sequence alignment (13). Finally, regions were added from an internally processed version of the unannotated Macaca mullata genome from . These rules delivered ∼7500 human target genes with sequence sets that contained a minimum of three and an average of six orthologous sequences from other genomes. Sequence regions extended 1.5 kb upstream and 100 bp downstream of transcription start sites (TSS), net of all types of repeats except LTR/ERV1, LTR/ERVL, LTR/MaLR and of coding sequences, which were masked.

In addition to a target sequence set, we supplied each discovery program with a genomic ‘background’ input file. The background input for each species consisted of 1000 concatenated search regions that were randomly selected from the genome's entire set of search regions.

We used multiple discovery methods in parallel, running each method with a range of parameter settings; typically target motif width and motif occurrence model were varied. We used a compact but diverse base set of discovery methods that consisted of CONSENSUS (14), MEME (15) and MotifSampler (16). Raw discovered motifs were post-processed, for example, to remove identical motifs reported by the same algorithm, and to merge strongly overlapped motifs.

We assigned p-values to motifs discovered by multiple methods across large sets of target genes whose sequence sets varied in species composition, as follows. First, we identified a representative subset of target sequence sets that sampled the range of species compositions of the target sets. Then, for each representative target sequence set, we created random sequence sets by retaining the original sequence from the target species and replacing each orthologous sequence with a synthetic sequence. The random synthetic sequence was generated from the target sequence by a tool we developed to simulate neutral evolution using published substitution rates and indel rates and lengths. Each random sequence set was generated 30 times in order to avoid statistical bias. We then submitted each target sequence set and all random sequence sets to identical motif discovery and post-processing procedures. We assigned method-independent (MI) scores to all motifs discovered in target and random sequence sets, using a trainable function that contained four non-negative parameters: where coefficients took the following possible values: D, W ∈ R, B ∈ {0, 1} and C ∈ [0, 1]. In the above equation, D characterized the number of site sequences in a motif; W and B characterized the shape of a motif's information content profile; and C characterized motif sequence conservation. The score increased when the motif was conserved for species that are evolutionarily remote from the target, and decreased when the motif was not conserved for species that are close to the target (primates, in the case of human). We used the distribution of MI scores for motifs from a target gene's random sequences to transform the MI scores for the target gene's motifs into p-values. Finally, we loaded the database with motifs whose p-values were below a threshold (which, for cisRED 1.2e, was 0.05).

A library of known transcription factor binding sites, split into mutually exclusive training and testing fractions, was used to optimize the scoring function and to characterize the performance of the system. The library contained ∼1000 sites for ∼300 human genes from TRANSFAC v9.1 (17), and ∼250 binding sites that we curated from the literature. We optimized the scoring function by simulated annealing, using a training fraction of known sites from randomly selected genes and two objective functions: the area under a receiver operating characteristic (ROC) curve and the number of experimentally known motifs that were not predicted. We assessed the system's predictive performance with a test fraction of known sites, using observables like sensitivity, specificity and positive predictive value (PPV) (e.g. ). Although comparisons of method performance are constrained by many factors, current system performance compared favourably to results in a recent study (6).

To identify groups of similar motifs, we defined two pairwise motif similarity metrics. For the first metric, we used a version of the Levenstein edit distance between two sequences that was modified to permit no internal gaps (18). For each motif, M, and its reverse complement, RC(M), we scanned motif pairs relative to each other, and reported the overall minimum average mutual edit distance to motif K, i.e. min(d(M,K), d(RC(M),K)). The second metric was based on the maximum information content shared between position frequency matrices derived for each motif, and also treated a motif and its reverse complement as equivalent. We hierarchically clustered pairwise dissimilarity matrices with the local density-based OPTICS algorithm (19). We extracted clusters from OPTICS' reachability output by applying an automatic cluster recognition method that identifies cluster boundaries as inflection points in the reachability plot (20), then traversing these hierarchical segmentation results with an algorithm that traced a deepest available path, constrained by a maximum preset depth.

The large size of mammalian genomes makes it challenging to organize the computational hardware and software infrastructure required to address such issues routinely. We did the large-scale discovery, similarity and co-occurrence calculations on a Beowulf-style, ∼400 CPU (Pentium III, Xeon, Opteron) OSCAR cluster () running Red Hat Linux 9; and remotely on the Beowulf ∼1700 CPU ‘glacier’ cluster at WestGrid (). We clustered motifs on a 12 dual-core CPU (UltraSPARC IV) SMP server with 96 GB of RAM running Solaris 9.

CISRED DATABASE CONTENTS, STRUCTURE AND ACCESS

Figure 2 shows a schematic diagram of the database design from a user perspective. The database infrastructure is designed to evolve to hold results from a range of mammals, as well as results from model organisms. Because promoter regions are enriched sources of regulatory elements, motif discovery and the cisRED design were both based on regions around TSSs. cisRED human v1.2 contains motifs from promoter regions of ∼7500 human genes using Ensembl v30 (NCBI 35) data, as well as a pilot result set of ∼250 mouse genes.

Figure 2

Database contents and high-level links, from a web user perspective, as of cisRED human v1.2. ‘Atomic motifs’ are motifs discovered in target sequence sets. ‘Groups’ are clusters of similar motifs that are identified by large-scale OPTICS (19) clustering. ‘Patterns’ are co-occurring sets of group-labelled motifs.

The current database design makes three types or levels of information available for regulatory elements: (i) atomic motifs, which are discovered independently in each target region sequence set; (ii) groups of similar motifs, each of which is a putative model for the binding site of a single transcription factor; and (iii) co-occurring patterns of models, which are putative regulatory modules.

Predicted regulatory elements can be viewed directly in cisRED's web user interface. From this interface, motifs can be viewed in the UCSC genome browser, in the Ensembl genome browser via a DAS server, or in the Sockeye comparative genomics workspace (21). A user can filter the displayed motifs by criteria like the p-value threshold, the orthologous species present in a motif and the discovery method.

The database contains a table of high-confidence globally coexpressed genes (22); genes coexpressed with each cisRED target gene are listed on that gene's page. Coexpression resources are available for download at . The database also contains a table of single nucleotide polymorphisms (SNPs) from dbSNP (8) that occur in predicted motifs. When a predicted human motif contains a SNP, the cisRED atomic motif page highlights this variation on the target sequence of a motif site sequence set, and the highlighted base hyperlinks to the SNP's primary source information in dbSNP.

The database can be accessed in several ways. A web user interface is available at . A current schema diagram is available from the ‘Databases & Methods’ page, and direct SQL queries can be run on the MySQL databases at . A user can download the data, with SQL files, as well as a compressed file that contains all input FASTA sequence sets. As genome resources and our data processing evolve, we update database contents. Because older versions of the database may not be compatible with the current user interface, only recent versions for each species can be accessed through the web or via a MySQL client. However, historical releases of cisRED databases are archived and can be downloaded.

Certain parts of the system's software are available from a tab on the cisRED home page. Sockeye (21) permits, for example, a user to assess details of conserved regions relative to genomic annotations in multiple sequence alignments. The HitPlotter visualizer displays large sets of discovered motifs from multiple-method discovery runs, and is available on request as a beta release. The database infrastructure is designed to facilitate installing cisRED locally.

FUTURE WORK

Database contents will be extended to include large-scale results for human, mouse, rat, Caenorhabditis elegans and Drosophila melanogaster. A new schema design supports these species. For example, cisRED human v2 will contain ∼18 000 human genes. The input sequence sets for this database were based on human TSSs that were identified by considering Ensembl (7) and RefSeq (23) annotations. To address the limitations of gene and transcript annotations for non-human species, corresponding vertebrate search regions were taken from UCSC (11,12) and ENCODE multiple sequence alignments (13). We are extending our ability to take advantage of unannotated and low-coverage genome sequence data.

We are continuing to improve motif post-processing and scoring. Optimizing the scoring function depends on having a large library of known motif sites; we anticipate that a newly created web database for submitting and curating binding sites from the literature may help us to enlarge and improve this resource (; S. B. Montgomery and O. L. Griffith, manuscript in preparation).

Given a scalable clustering method, we continue to assess motif similarity metrics and how best to use the hierarchical information output from OPTICS. Given groups of similar motifs, we are applying group labels to atomic motifs, then identifying overrepresented co-occurring motif patterns using hypergeometric statistics, imposing separation constraints on neighbouring motifs, and searching in two stages for patterns larger than pairs (24,25). We are implementing methods for annotating discovered motifs as known or novel against known site resources. Given annotated motifs, we will annotate overrepresented co-occurring patterns of human motif pairs as known versus novel using the TRANSCompel resource (26). We are applying genome-scale motif clustering and co-occurrence as filters for predicted motifs that may improve the predictive reliability and the resulting catalogue of conserved regulatory elements.

We have assembled a large multi-species coexpression resource that contains public microarray and SAGE data from diverse sources (22) (Table 1). We have shown that combining global coexpression data from multiple platforms improves confidence in coexpression predictions when assessed against the Gene Ontology (GO) (27). From this, we established GO-based Pearson correlation thresholds that identified high-confidence globally coexpressed gene pairs. The coexpression database makes results available from this global analysis and from two other recent analyses (28,29). Although coexpressed genes can have similar regulatory elements, the system's predictive performance improved only marginally when inputs included coexpressed genes in addition to orthologous genes (data not shown). Given these results, currently we are assessing an approach that includes no coexpressed genes in motif discovery inputs, but uses coexpression information to assess groups and co-occurring patterns identified in genome-scale sets of atomic motifs.

Table 1

Contents of coexpression database (22)

Species	Platform	Experiments	Unique genes
H.sapiens	SAGE (short)	272	20 312
	Oligo microarray	1640	12 452
	cDNA microarray	2852	13 111
M.musculus	SAGE (short)	85	12 715
	Oligo microarray	1802	8 164
	cDNA microarray	366	8 102
C.elegans	SAGE (long/short)	26	15 685
	cDNA microarray	1059	15 956
Total		8102	54 434

We will extend the database user interface to offer more complex user filtering, as well as motif searches based on consensus strings or matrices. For work with classes of regulatory elements that are defined by wet lab data types based on, for example, ChIP or DNase I hypersensitivity (e.g. see ENCODE tracks within the UCSC genome browser) (11,13), we have designed a new schema that is based on search regions rather than on genes, and will extend the user interface to support this.

We will continue to assess the contributions to regulatory element predictions of different genomes and sets of genomes. We will integrate and assess new motif discovery methods, and will identify a best minimal set of methods on an ongoing basis.