The comprehensive transcriptional analysis in Caenorhabditis elegans by integrating ChIP-seq and gene expression data

Summary

The fundamental step of learning transcriptional regulation mechanism is to identify the target genes regulated by transcription factors (TFs). Despite numerous target genes identified by chromatin immunopre-cipitation followed by high-throughput sequencing technology (ChIP-seq) assays, it is not possible to infer function from binding alone in vivo. This is equally true in one of the best model systems, the nematode Caenorhabditis elegans (C. elegans), where regulation often occurs through diverse TF binding features of transcriptional networks identified in modENCODE. Here, we integrated ten ChIP-seq datasets with genome-wide expression data derived from tiling arrays, involved in six TFs (HLH-1, ELT-3, PQM-1, SKN-1, CEH-14 and LIN-11) with tissue-specific and four TFs (CEH-30, LIN-13, LIN-15B and MEP-1) with broad expression patterns. In common, TF bindings within 3 kb upstream of or within its target gene for these ten studies showed significantly elevated level of expression as opposed to that of non-target controls, indicated that these sites may be more likely to be functional through up-regulating its target genes. Intriguingly, expression of the target genes out of 5 kb upstream of their transcription start site also showed high levels, which was consistent with the results of following network component analysis. Our study has identified similar transcriptional regulation mechanisms of tissue-specific or broad expression TFs in C. elegans using ChIP-seq and gene expression data. It may also provide a novel insight into the mechanism of transcriptional regulation not only for simple organisms but also for more complex species.

1. Introduction

Nowadays, Caenorhabditis elegans (C. elegans) is regarded as a model biological system, in which dozens of transcription factors (TFs) have now been identified and characterized by a combination of genetics, molecular screens and the genome sequencing project. Sufficient information for those TFs is available to suggest the important roles in cell determination, differentiation or other functions of them. However, only a small fraction of C. elegans TFs have so far been identified and studied on the candidate target genes and their expression patterns.

As we all know, the most basic step in the study of gene regulation is to uncover the target genes regulated by TFs. Recently, with the extensive application of chromatin immunoprecipitation followed by high-throughput sequencing technology (ChIP-seq), genome-wide identification of TF binding events becomes more facilitated. Nowadays, the numerous genome-wide TF binding data are freely available from the databases of ENCODE for human and mouse as well as modENCODE for C. elegans and Drosophila melanogaster, which may provide useful biological insights for understanding the mechanisms of transcriptional regulation. Nevertheless, the binding sites alone are not sufficient to infer regulation. Lots of previous studies were focused on identifying the target genes controlled by TFs and dissecting the transcriptional regulatory networks underlying the relevant biological process through a simple peak-based method (Boyer et al., 2005; Chen et al., 2008; Kim et al., 2008). In those studies, the genes with promoters overlapping with one or more peaks were identified as the targets in most cases, but less information on the genes of other regions in genome was reported, especially for the regions out of 2 kb transcription start site (TSS) of a gene. In this case, some of the important long-range interactions were missed and it is unable to distinguish functional issues from non-functional binding sites. Furthermore, some of TF networks using ChIP-seq and gene expression data based on matrix decomposition have been constructed to unravel transcriptional regulatory programmes in complex systems (Yan et al., 2013; Guan et al., 2014).

Meanwhile, gene expression profiles of more than 30 different cells and developmental stages as well as whole-animal RNA isolated from seven different developmental stages using tiling arrays have been generated to produce a spatial and temporal map of C. elegans gene expression (Spencer et al., 2011). It provides a basis for establishing the roles of individual genes in cellular differentiation by measuring native transcripts from a broad array. The datasets of relative gene expression levels across tissue types and developmental stages are available online for further study.

It will be clear that the integration of the genome sequencing information and gene expression data in genome wide using C. elegans as a model biological system will result in uncovering the comprehensive transcriptional mechanisms. We hypothesized that a particular gene set regulated by a TF should behave in some measurably consistent manner. In this study, we tested the hypothesis by integrating TF binding and gene expression data of C. elegans in genome wide to further explore the mechanism of transcriptional regulation, which can be applied to other species.

2. Materials and methods

The summary of collected ChIP-seq and gene expression data was shown in Table 1. The workflow of gene expression and ChIP-seq data collection and integration analysis was shown in Fig. 1 and described as follows.

Fig. 1. The workflow of gene expression and ChIP-seq data collection and integration analysis.

Table 1. The summary of collected ChIP-seq and gene expression data

(i) ChIP-seq datasets collection

The ten ChIP-seq datasets of six TFs (HLH-1, ELT-3, PQM-1, SKN-1, CEH-14 and LIN-11) with tissue-specific and four TFs (CEH-30, LIN-13, LIN-15B and MEP-1) with broad expression patterns were collected from modENCODE database (Gerstein et al., 2010; Niu et al., 2011; Landt et al., 2012). Briefly, the raw sequences were aligned to C. elegans genome (WS232) using ELAND allowing up to two mismatches per read. And peaks were called using PeakSeq with the threshold of P-value as 0·05 (Rozowsky et al., 2009).

(ii) Target gene assignment and classification

To identify candidate target genes possibly regulated by the TFs involved in the ChIP experiments, 5 kb regions both up- and downstream of a peak summit were first searched for transcriptional start (TSS) and end (TES) sites in the C. elegans genome (WS232) as described previously (Zhong et al., 2010). If no TSS was annotated for a gene, the TSSs were arbitrarily set to the positions 150 bases upstream of the translational start sites. A total of 31 640 protein-coding transcripts from 20 553 genes were used for target identifications. The targets of each ChIP-seq study were assigned to unique gene and the target genes were classified into 15 classes based on their positions relative to the peak summits. Those with peak summit within the gene region were defined as the class of C0; those with peak summit located 0–5 kb upstream the TSS were defined as C1–5 for each 1 kb interval, respectively, whereas 0–5 kb downstream the TES as C11–15, respectively (Fig. 2(a)). For those peaks with no candidate genes in the 5-kb flanking region, the two nearest genes were found. The nearest gene upstream the peak was defined as C21, and the second nearest was C22. Similarly, the nearest gene downstream the peak was defined as C23, and the second nearest was C24 (Fig. 2(b)). If the same gene could be assigned to several classes due to the different peaks nearby, it was assigned into a class with the smallest number (Fig. 2(c)). For example, if one gene could belong to both C1 and C4, it was defined as C1 target finally.

Fig. 2. The target genes assignment and classification. The targets of each ChIP-seq study were assigned to unique gene and the target genes were classified into 15 classes based on their positions relative to the peak summits. Those with peak summit within the gene region were defined as the class of C0; those with peak summit located 0–5 kb upstream the TSS were defined as C1–5 for each 1 kb interval, respectively, whereas 0–5 kb downstream the TES as C11–15, respectively (a). For those peaks with no candidate genes in the 5-kb flanking region, the two nearest genes were found. The nearest gene upstream the peak was defined as C21, and the second nearest was C22. Similarly, the nearest gene downstream the peak was defined as C23, and the second nearest was C24 (b). If the same gene could be assigned to several classes due to the different peaks nearby, it was assigned into a class with the smallest number (c).

(iii) Association of binding site locations with tissue-specific or broad expression

For association of the target genes of each TF with those enriched in tissue-specific cells or broad expression of whole animal based on microarray expression assays, expression values of all genes with corresponding patterns were derived from tiling array data (Spencer et al., 2011). The summary of involved ChIP-seq and gene expression data were shown in Table 1. For each gene, the normalized log2 intensities from three replicates were averaged to give rise to its expression value. The expression values of the target genes were extracted for each gene. To plot the tissue-specific or broad expression of different classes of the target genes, those that all of whole genome genes were used as a reference and the plotting was performed with ‘vioplot’ package in R (Hintze JaN, 1998). Mann–Whitney U-test was performed to calculate the significance levels (P-value) of gene expression difference between reference and the target genes of each class with cutoff of P-value as 0·001.

(iv) Calculation of TF activity

With both expression values and initial binding relationship provided, we can estimate transcription factor activity using proper mathematical model. Here, we applied network component analysis (NCA) method to calculate the transcription factor activity in 15 different classes for each TF (Liao et al., 2003). Generally, the values for target genes were set to be 1 in initial transcription factor activity matrix, whereas all others were set to be 0. The final transcription factor activity was estimated by matrix decomposition, and the value will vary from 0 to some number bigger than 1.

3. Results and discussion

(i) Most of binding sites are located within 4 kb upstream of TSS

Based on the assignment and classification of ChIP-seq targets for all ten TFs, the target genes were classified into 15 classes, from C0 to C24 (Fig. 2, see in the Materials and methods). According to the proportion of target genes in each class of ten ChIP-seq studies, more than 60% of targets were classified from C0 to C4 (Fig. 3). It suggested that most of binding sites are located within 4 kb upstream of TSS, not only for the regulation of TFs with tissue-specific expression pattern but also for the regulation of TFs with broad expression pattern.

Fig. 3. The proportion of target genes in each class of ten ChIP-seq studies. The bar charts show the proportion of target genes in each class of ten ChIP-seq studies, including the TFs of HLH-1, ELT-3, PQM-1, SKN-1, CEH-14 and LIN-11 with tissue-specific expression patterns as well as CEH-30, LIN-13, LIN-15B and MEP-1 with broad expression patterns. X-axis represents the different ChIP-seq studies; Y-axis represents the percentage (%) of target genes in each class (C0–C24).

(ii) Integration analysis of binding site locations with tissue-specific expression data

Actually, there are several major transcriptional regulators of tissue-specific differentiation in C. elegans, which are expressed in specific tissues where they are known to bind and regulate diverse well-studied tissue-specific genes. Among the 22 TFs involved in the study of ChIP-seq as part of the modENCODE Consortium, there are some of tissue-specific TFs, including HLH-1 for body wall muscle (BWM), ELT-3 for hypodermis (Hyps), PQM-1 for intestine, SKN-1 for pharynx as well as CEH-14 and LIN-11 both for neurons (Table 1) (Niu et al., 2011). Through read mapping with ELAND and peak calling with PeakSeq as well as the peak gene annotation, 13 384 target genes for HLH-1, 7989 target genes for ELT-3, 8337target genes for PQM-1, 10 553 target genes for SKN-1, 7759 target genes for CEH-14 and 4360 target genes for LIN-11 were identified, respectively. According to the binding site locations, the target genes were classified into 15 different categories (C0–C24) for all the six TFs studies (Table 2, see in Section 2).

Table 2. The number of target genes in each class for ten TFs

Numerous target genes have been identified by the ChIP-seq assays, but we are not certain whether the bindings are functional or not in vivo as well as whether the corresponding expression behaviours of each tissue-specific TF target genes are similar or not. We addressed the question by the integration analysis of the tissue-specific expression data derived from tiling arrays (Spencer et al., 2011). The results of normalized expression levels in specific tissues for different categories of corresponding TF target genes, respectively, can be seen in Fig. 4. Moreover, the expression patterns of the mean expression distribution of each class in six tissue-specific studies were shown in Fig. 5(a). In common, all first three classes (C0–C2) of the target genes showed significantly elevated level of expression as opposed to that of reference whole genome genes for those six different TFs (Mann–Whitney U-test, P < 10⁻³), supporting that the majority of binding events are likely to be functional mostly through up-regulating its target genes within 2 kb upstream from TSS. Another class C3 was also reported as significantly elevated one in some of studies, including for ELT-3, SKN-1, CEH-14 and LIN-11. It indicated that bindings of those four TFs within 2–3 kb upstream from a gene's TSS may be also considered as one candidate group that more likely control its function by up-regulating its activities. Intriguingly, expression of all first six classes of the LIN-11 target genes representing genes with a peak inside (C0) or 0–5 kb upstream (C1–C5) from their TSS showed significantly higher than reference. In contrast, most classes of the 0–5 kb downstream (C11–C15) target genes for those six TFs showed decreased level of expression, suggesting a diverse expression behaviour of binding targets downstream from TSS that likely inhibit the TF activities.

Fig. 4. The expression patterns of tissue specific TFs target genes. The Violin plots of normalized expression levels in specific tissues for different categories of corresponding TF target genes were performed for six TFs with tissue-specific expression pattern, including HLH-1 for body wall muscle (BWM), ELT-3 for hypodermis (Hyps), PQM-1 for intestine, SKN-1 for pharynx as well as CEH-14 and LIN-11 both for neurons, respectively. X-axis represents 15 different classes of target genes and reference of genome gene expression; Y-axis represents the normalized gene expression value. Significant levels by Mann–Whitney U-test (P < 10⁻³) are indicated in the plots, dark-red for significantly up-regulated and dark-green for significantly down-regulated.

Fig. 5. The mean expression distribution and NCA results of each class in six tissue-specific studies. (a) It showed the mean expression distribution of each class in six tissue-specific studies, which is based on the integration analysis of ChIP-seq and gene expression data, including HLH-1 for body wall muscle (BWM), ELT-3 for hypodermis (Hyps), PQM-1 for intestine, SKN-1 for pharynx as well as CEH-14 and LIN-11 both for neurons, respectively. X-axis represents 15 different classes of target genes and reference of genome gene expression; Y-axis represents the mean level of normalized gene expression value. (b) It showed the transcription factor activity in 15 different classes for each of six TFs with tissue-specific expression pattern based on NCA method, including HLH-1 for BWM, ELT-3 for Hyps, PQM-1 for intestine, SKN-1 for pharynx as well as CEH-14 and LIN-11 both for neurons, respectively. X-axis represents 15 different classes of target genes and reference of genome gene expression; Y-axis represents the percentage of target genes with the corresponding A scores larger than 1 compared with the raw number of targets in ChIP-seq data.

Obviously, the mean expression level of target genes in C0 and C1 classes reached the summit as opposed to that of the other classes for all six TFs (Fig. 5(a)), which also suggested their most important roles in regulatory function. What is more interesting is that all of the class of C21 representing the nearest target genes out of 5 kb upstream from TSS showed sharply increased level of expression. However, the gene expression difference test compared with the reference is not significant, which is mainly due to the small sample size of genes count in C21.

In order to confirm our results, we have also performed NCA on the ChIP-seq and gene expression data. Compared to the traditional statistical methods of principal component analysis and independent component analysis, NCA is more comparable to determine low-dimensional representations of high-dimensional datasets for reconstruction of regulatory signals in biological systems (Liao et al., 2003). Here, we applied NCA method to calculate the transcription factor activity in 15 different classes for each TF. According to the percentage of target genes with the corresponding A scores larger than 1 compared with the raw number of targets in ChIP-seq data, there appears to be maximum plot ratio in the class of C0, C1 and C21 (Fig. 5(b)). It suggests that they are three of most essential classes of target genes in the TF transcriptional regulation issues, which is consistent with our above integration analysis.

(iii) Integration analysis of binding site locations with broad expression data

Similarly, there are also several involved TFs with the broad expression pattern (expressed in the whole animal) in TF binding sites ChIP-seq database of modENCODE project, including CEH-30, LIN-13, LIN-15B and MEP-1. There were 9469 target genes for CEH-30, 8185 target genes for LIN-13, 7445 target genes for LIN-15B and 12 013 target genes for MEP-1, respectively, identified in those four broad expression TFs ChIP-seq. We then integrated all of the ChIP-seq targets information with the gene expression data in different whole animal developmental stages derived from the same tiling arrays (Spencer et al., 2011). As expected, the bindings within 2 kb upstream of or within its target gene (C0–C2) showed significantly elevated level of expression compared to that of the reference for those four broad expression TFs (Mann–Whitney U-test, P < 10⁻³) (Fig. 6), which is the same as above results of six TFs with tissue-specific expression pattern. Moreover, similar results were observed for the bindings within 2–3 kb upstream (C3) of those TFs target genes except for MEP-1 as well as for the bindings within 3–4 kb upstream (C4) of LIN-15B. Generally, the bindings within 0–5 kb downstream (C11–C15) of those TFs showed constant or decreased expression values as opposed to reference, the significance appeared in C15 for LIN-13 and C5 to C15 for MEP-1. According to the mean expression distribution of each class in those four TF studies, the mean expression level of nearest bindings out of 5 kb upstream of target genes (C21) was also found to have a summit like C0 and C1 (Fig. 7(a)), which suggested its equally essential roles in regulatory function of TFs with broad expression pattern. The results were validated by the following NCA approach for those four TFs (Fig. 7(b)).

Fig. 6. The expression patterns of broad expressed TFs target genes. The Violin plots of normalized expression levels in whole animal for different categories of corresponding TF target genes were performed for four TFs with broad expression pattern, including CEH-30, LIN-13, LIN-15B and MEP-1, respectively. X-axis represents 15 different classes of target genes and reference of genome gene expression; Y-axis represents the normalized gene expression value. Significant levels by Mann–Whitney U-test (P < 10⁻³) are indicated in the plots, dark-red for significantly up-regulated and dark-green for significantly down-regulated.

Fig. 7. The mean expression distribution and NCA results of each class in four broad expression TFs studies. (a) It showed the mean expression distribution of each class in four broad expression TFs studies, which is based on the integration analysis of ChIP-seq and gene expression data, including CEH-30, LIN-13, LIN-15B and MEP-1, respectively. X-axis represents 15 different classes of target genes and reference of genome gene expression; Y-axis represents the mean level of normalized gene expression value. (b) It showed the transcription factor activity in 15 different classes for each of four TFs with broad expression pattern based on NCA method, including CEH-30, LIN-13, LIN-15B and MEP-1, respectively. X-axis represents 15 different classes of target genes and reference of genome gene expression; Y-axis represents the percentage of target genes with the corresponding A scores larger than 1 compared with the raw number of targets in ChIP-seq data.

(iv) Similar transcriptional regulation mechanisms of tissue-specific or broad expression transcription factors in C. elegans

Although diverse transcription factor binding features revealed by genome-wide ChIP-seq in C. elegans (Niu et al., 2011), some of similar issues in regulatory mechanisms may be found between tissue-specific and broad expression TFs. Firstly, the majority of binding events in C. elegans are likely to be functional mostly through up-regulating its target genes, particularly three regions of within 2 kb upstream or within its target genes rather than those located elsewhere from its targets. Most of transcriptional regulation issues were previously reported to occur by a TF binding in a promoter region, which was defined as 2 kb from a TSS, particularly in simple organisms such as bacteria or yeast (Capaldi et al., 2008). However, the regulation within the target needs further validation. Secondly, some of regulation issues in C. elegans were also found to occur through long-range regions like enhancers, including 2–3 kb (C3) or 3–4 kb (C4) upstream from TSS. It looks similar as more complex species, in which more regulation occurs through long-range enhancers often spanning many tens of thousands of base pairs and the enhancer interactions become the principal form of gene regulation (Arnosti & Kulkarni, 2005; Lee et al., 2005).

As we all know, transcription regulatory networks play a pivotal role in the development, function and pathology of metazoan organisms, which comprises direct or indirect protein–protein or protein–DNA interactions between TFs and their target genes. For our ten TFs, there existed synergistic effects or antagonistic effects between the TFs by co-binding with the same DNA, which have been reported in the previous studies (Lum et al., 2000; Vermeirssen et al., 2007; Cheng et al., 2011).