Introduction

Embryonic stem (ES) cells have limitless self-renewal capability as well as the ability to differentiate into multiple cell types. This makes them a good model system to study not only for understanding normal development, but for potential clinical applications. They facilitate in vitro study of early developmental events as well as differentiation into crucial cell types, such as hematopoietic or neuronal cells (Keller, 2005). Most importantly they have proven to be a highly valuable resource in regenerative medicine, where ES cells have shown a great potential to be used as a treatment against disease or in tissue repair after injury (Balber, 2011). This therapeutic application has strengthened the push towards under-standing how stem cells are programmed during self-renewal and differentiation.

Transcription factors (TFs) are key players in driving cellular reprogramming (Takahashi and Yamanaka, 2006). Many TFs have so far been identified for their crucial role in pluripotency and differentiation (Young, 2011), and genome-wide putative targets for them have been characterized using ChIP-sequencing or ChIP-on-chip (Xu et al., 2013). On the other hand, it is also becoming feasible to obtain binding sequence preferences for TFs in mam-malian cells for hundreds of TFs (Jolma et al., 2013). To utilize the complementary information from both of these resources for further biological discovery, we have developed a web tool called ‘TRES’ to link gene sets to likely upstream regulators in ES cells.

Databases and protocols behind this web tool

TRES predicts transcriptional control mechanisms that may regulate the coordinated expression of a given set of genes using the four databases outlined below. The gene list is interrogated against each database separately to find enriched transcription factor associations. TRES generates two outputs: a single summary integrated output as well as individual predictions from each dataset separately. For the integrated summary, as the P values of motif-based and ChIP-based approaches are not directly comparable, we rank the TFs predicted from each database and then present them in the order of aggregated rank. For a more informative output, the p-values for the transcription factors from each database are then collected and presented to the user, ordered from the lowest to highest. The output from TRES shows which transcription factors are enriched, the database used, the gene (for ChIP data) or binding site (for motif data) overlap and the Bonferroni corrected p-value. Details of how these p-values are calculated are given below.

1) ChIP sequencing datasets: Using a murine blood ChIP sequencing compendium, we have previously shown that compared to approaches based on comparing transcription factor target gene list only, counting the number of binding events for each TF performs much better at predicting the transcriptional control behind a gene set of interest (Joshi et al., 2012). We collected publicly available ChIP sequencing data for mouse and human embryonic stem cells to form two compendia of 97 murine and 49 human samples. The peaks were used if provided by the authors otherwise they were calculated using the MACS (Zhang et al., 2008) program default settings. This resulted in 97 murine and 49 human peak files in the mm9 and hg19 genome assemblies, respectively. Binding events in the promoter and gene body were associated with the corresponding gene. On the other hand, inter-genic peaks were associated with the nearest gene on either side within 50 kb, such that each peak is assigned to at most two genes.

For a set of user-defined genes the following procedure is followed. We calculate the number of genomic regions mapped to the genes (n) for all 190k murine and 170k human genomic regions bound by at least one TF (N). For each ChIP-Seq dataset, the number of peaks (m) near user defined genes (k) is calculated. The p value is calculated using a hypergeometric test (the Fischer exact test). We use Bonferroni correction to adjust the p values for multiple hypotheses testing.

The longer the gene locus, the more likely the number of binding sites is higher. In order to account for this bias, we normalised the weights with gene length. We calculated the correlation coefficient between the number of binding events in a gene locus and the gene length, which resulted in a value of 0.27 (p value < 1e-16). On the other hand, the log of the gene length shows a greater correlation (with coefficient 0.46). We have, therefore, chosen to normalise the weights with respect to the logarithm of the gene length.

2) ChIP-on-chip dataset: Traditionally, the overlap of ChIP-on-chip targets and a given gene set is calculated using a non-weighted approach. We calculate the number of binding events in each gene locus using the ChIP sequencing compendium (1) above, and the weight is assigned to each gene as the log of the number of binding events normalised by the gene length. The overlap of genes with ChIP-on-chip targets is calculated taking weights into account. To determine the statistical significance, the p value is calculated using a hypergeometric test (the Fischer exact test). Again, we use Bonferroni correction to adjust the p values for multiple hypotheses testing.

We further show that this weighted approach, where each gene is weighted according to the number of binding sites in its locus (gene length normalised), gives better performance. The ChIP-on-chip data for 41 transcription factors was obtained from (Xu et al., 2013). 11 of 16 deletion datasets (Nishiyama et al., 2009a) and 5 over expression data sets (Xu et al., 2013) were true positives (i.e. the perturbed transcription factor was given by TRES when no weights were assigned to the gene lists) whereas 14 of 16 deletion data sets and 7 overexpression datasets were true positives when weights were assigned to gene lists, clearly showing the benefits of the weighted approach.

3) JASPAR sequence motif dataset:We downloaded 684 mammalian sequence motifs from the JASPAR database (Bryne et al., 2008). Traditionally enriched sequence motifs are searched for in the promoter regions of genes of interest. In order to search through a more exhaustive list of regulatory regions, which includes putative enhancers with promoters, we used the ChIP-sequencing compendium (1) above. We searched for the enrichment of sequence motifs not only in the gene promoters, but also in all putative regulatory regions in the gene locus. The putative regulatory regions were defined as the regions bounded by at least one transcription factor in the ChIP-sequencing compendium (1). The Centrimo tool (Bailey and Machanick, 2012), which is available as part of the MEME suit, was used to find enriched motifs.

4) Jolma 2013 sequence motif dataset: Jolma et al. (2013) analyzed the sequence-specific binding of human TFs using high-throughput SELEX and ChIP sequencing, and collected binding specificities for 550 TFs. The motif enrichment was calculated using the same method as in (3).

Citing this webtool

"TRES predicts transcription control in embryonic stem cells." Christopher Pooley, David Luau, Patrick Lombart, Berthold Gottgens and Anagha Joshi. Bioinformatics (submitted).

References

Bailey,T.L. and Machanick,P. (2012) Inferring direct DNA binding from ChIP-seq. Nucleic Acids Res., 40, e128–e128.

Bryne,J.C. et al. (2008) JASPAR, the open access database of transcription factor-binding profiles: new content and tools in the 2008 update. Nucleic Acids Res., 36, D102–106.

Jolma,A. et al. (2013) DNA-Binding Specificities of Human Transcription Factors. Cell, 152, 327–339.

Joshi,A. et al. (2012) Gene Set Control Analysis (GSCA) predicts haematopoietic control mechanisms from genome-wide transcription factor binding data. Exp. Hematol.

Lachmann,A. et al. (2010) ChEA: transcription factor regulation inferred from integrating genome-wide ChIP-X experiments. Bioinforma. Oxf. Engl., 26, 2438–2444.

Liberzon,A. et al. (2011) Molecular signatures database (MSigDB) 3.0. Bioinforma. Oxf. Engl., 27, 1739–1740.

Martello,G. et al. (2012) Esrrb is a pivotal target of the Gsk3/Tcf3 axis regulating embryonic stem cell self-renewal. Cell Stem Cell, 11, 491–504.

Nishiyama,A. et al. (2009) Uncovering early response of gene regulatory networks in ESCs by systematic induction of transcription factors. Cell Stem Cell, 5, 420–433.

Xu,H. et al. (2013) ESCAPE: database for integrating high-content published data collected from human and mouse embryonic stem cells. Database J. Biol. Databases Curation, 2013, bat045.