Identification of non-exonic regions

We classified all non-exonic genomic regions into several types based entirely on the exon boundary annotations from General Feature Format (GFF) files http://www.sanger.ac.uk/Projects/C_elegans/WORMBASE/GFF_files.shtml. We defined nine types of regions, based on the exons surrounding the region (Fig. ​(Fig.3).3). These region types were used only to label the final hits, and were not used in the algorithm to rank them or filter them. They are included merely to aid the user who wants to further inspect the results manually. We excluded exonic regions from our analysis since a) they are unlikely to contain transcription factor bindings sites and b) the high conservation of coding sequences between species would negate the utility of "phylogenetic footprinting".

Open in a separate window

Figure 3

Classification of non-exonic regions. A hypothetical gene arrangement is shown. "5' intergenic": between exon1 and exon1 of two separate genes; "3' intergenic": between the last exon of both genes. "5'/3' intergenic": between first exon of one gene and last exon of the other gene; "intronic#": between any two exons of one gene; "other": all other possible combinations. In cases where the gene flanking a segment is known to exhibit alternative splicing, the segment was prefixed with 'alt_', i.e. 'alt_intronic#', 'alt_3'intergenic', etc. Two other categories, BEGIN and END, denote regions at the beginning or ending of the chromosome, in the case of C. elegans, or of the sequencing reads in the case of C. briggsae. There were two exceptions to the procedure. The first was due to the fact that the C. briggsae genome we used was an unassembled collection of 578 individual sequence reads. 112 of these reads had no exon annotations, and were ignored in this study. Of these 112, only two were greater than 10,000 bases long, with an average length of 3679.3 nucleotides. Secondly, there were 16 C. elegans and 35 C. briggsae exon annotations one nucleotide long. By visual inspection, we determined that for C. elegans these exons were in fact longer than one nucleotide, but noncoding: in all cases the single nucleotide is 'A' and when spliced forms a TGA stop codon. They were treated as non-existent for this study, which has very little effect on the procedure except that the last true intron of the gene will be considered its 3' region. For C. briggsae, they appear to be errors in the gene annotations and fall within introns. Thus, they were treated as part of the intron in which they occur.

Scanning window scoring procedure

Starting with the input of experimentally defined binding sites of a dimeric homeodomain transcription factor complex composed of the TTX-3 and CEH-10 proteins [16](ASW and OH, submitted), we first use the hidden Markov model software package HMMER [17] with the command 'hmmbuild --null <background-frequency-file> --prior <prior-frequency-file> <output-position-weight-matrix> <input-binding-site-alignment>.' The position weight matrix is an n × 4 matrix where n is the length of the transcription factor binding site alignment input used (in our case, either 14 or 16 nucleotides). The matrix is defined as:

where the effective frequency is determined as the normalized, prior-adjusted, weighted sum of counts of each nucleotide in the alignment columns, in which the sequence weights are determined using a tree-based scheme. We use the resulting position weight matrix as input to CisOrtho. Then, CisOrtho finds approximately the N (option -t) highest scoring hits, with the restriction that no single gene has more than D (option -d) hits. The score for a given sequence window is simply the sum of the n cells of the position weight matrix corresponding to the position and nucleotide of the sequence window. For example, the sequence AACTCG would be given the score m_1A + m_2A + m_3C + m_4T + m_5C + m_6G for a 4 × 6 position weight matrix |m_ij|. This scoring scheme is simply the log-odds scoring used by the original HMMER software, adapted as a scanning window algorithm via CisOrtho. We note that known target genes of TTX-3 do not contain clustered binding sites [16](ASW and OH, submitted) which allowed us to eschew clustering of binding sites as a search criteria. It is possible that previously described searches for transcription factor targets which used the clustering criteria [2,3] may result in too many false-negative results.

The effect of options -t and -d warrants further comment. To find approximately N (option -t) highest-scoring hits in a large amount of DNA, CisOrtho must first estimate the raw minimum score cutoff which yields this many hits. Technically, this is achieved by scoring every 100^th sequence window, finding the top N/100 hits, and using the lowest score of those as the cutoff in the full search. Since the criterion for acceptance is a raw score cutoff, the actual number of hits retrieved in the whole sequence cannot be controlled precisely, but depends on the statistics of the sequence and position weight matrix. Furthermore, since no gene is allowed more than D associated hits (option -d), the total number retrieved will be lower if there is a lot of clustering and thus many genes have more than D actual hits. If the user wants to perform an exhaustive search, the special value of zero for -t indicates using no score cutoff, thus accepting all windows as hits. Likewise, the -d option can be set to an arbitrarily high value (in the downloaded software). These settings may be of interest in preliminary runs to determine the clustering behavior of the hits. More restrictive settings should then be used based on the results of the preliminary run.

Orthology-based filtering and HTML tables

First, CisOrtho retrieves the set of all C. elegans/C. briggsae ortholog pairs provided as a one-to-one mapping in the file orthologs-2.00 ftp://ftp.wormbase.org/pub/wormbase/briggsae/run25_analysis_freeze2.00/orthologues/. Then, for each ortholog pair in this list, CisOrtho finds the highest scoring hit for each of the C. elegans and C. briggsae orthologs, and stores the two hits as a 'hit-pair'. Note that there is no restriction on which region type (3' intergenic, 5' intergenic, etc.) each hit comes from. For example, the C. elegans highest-scoring hit may be 'intronic1' while the C. briggsae hit may be 3'-intergenic. Each hit in the hit-pair is described by the score, nucleotide sequence and type of region in which it occurs. In addition, for each of the directly adjacent genes, the coding strand, gene name and relative position of the hits in relation to the flanking genes are provided (Fig. ​(Fig.3).3). Finally, for the hit-pair itself, the number of mismatches between C. elegans and C. briggsae hits is given, along with the average, maximum and minimum scores from the pair of hits. Several HTML output files are generated with the information sorted into columns, and each hit-pair appearing as a grouped pair of lines. The HTML tables are each sorted primarily and secondarily by some combination of the mismatch number and average, maximum or minimum score (Fig. ​(Fig.4).4). The user has also the option to display multiples target sites located in a hit pair (this option takes into account that transcription factors often bind to multiple sites in a promoter)(Fig. ​promoter)(Fig.22).

Open in a separate window

Figure 4

Output of the program pipeline. Hits of a search with the TTX-3 consensus binding site is shown. num: number in list. mis: Number of base mismatches between first C. elegans and first C. briggsae hits. segtype: Type of non-exonic region (see Figure ​Figure3).3). str1/2: negative (N) or positive (P), strand on which the first/second of the two genes that flank the identified target site are located; offset1/2: distance of the target site to the flanking gene(s) (in relation to the start codon if the target site is 5' or located in an intron; in relation to the stop codon if the site is 3' to the gene; in the latter two cases, the number has a positive value); ID: cosmid name of the flanking genes, name: flanking gene names (if available). Gene IDs/names are linked to the WormBase gene model at http://www.wormbase.org, which contains further information about the gene. In case there are multiple target sites located in a defined inter/intragenic region, there is an option to report the n highest scoring hits for each ortholog. If this option is used, the top-scoring C. elegans or C. briggsae hit in each hit-pair will be highlighted, and the next n-1 hits will be gray. Color coding: Orthologous C. elegans/C. briggsae genes ("hit-pairs") are color coded in blue (Y39A3B.5 and CBG15122 are orthologs) and green (M01E10.2 and CBG15118 are orthologs).