BACKGROUND
t cells are key regulators of the adaptive immune system and have a central role in defense against pathogens and cancer as well as protection from autoimmune diseases.
cd4+ t lymphocytes can differentiate to functionally distinct effector subtypes, including t helper  <dig> , t helper  <dig>  and more recently described t helper  <dig>  cells  <cit> .
th <dig> cells secrete effector cytokine ifn-γ and regulate cell-mediated immunity and play a role in the pathogenesis of autoimmune diseases, such as multiple sclerosis.
th <dig> cells in turn produce il- <dig>  il- <dig>  and il- <dig> cytokines, and mediate immunity against extracellular pathogens and allergic reactions.
th <dig> cells, characterized by the production of a proinflammatory cytokine il- <dig>  regulate inflammatory responses on the mucosal surfaces.
for the overall health in humans and animals, the proper balance between different effector t cell types and t regulatory cells is crucial  <cit> .
aberrant activation of th <dig> and th <dig>  or th <dig> cells can trigger inflammatory autoimmune diseases as well as asthma and allergy.
previous studies utilizing genome-wide expression data and computational modeling have aimed at revealing the master regulators and regulatory networks within the differentiating th <dig> and th <dig> cells  <cit> .
however, studies in human have been less extensive than in mouse due to the difficulty in collecting sufficient amount of samples to comprehensively profile t cell differentiation over time.
in addition, lack of appropriate computational methods suitable for analyzing large-scale experimental data from multiple lineages over several time points spanning the lineage commitment process has limited the progress on revealing dynamics and molecular mechanisms underlying multiple lineage commitment.
a number of different time-series analysis approaches have been proposed to solve large-scale lineage commitment analysis problems.
the general purpose f-test  <cit>  can be used to test the difference between time-series data sets, but it does not extend to simultaneous comparison of multiple lineages and fails to take into account the correlation between the measurements at different time points.
more recent approaches to analyze time-series data, including regression, differential expression, discriminant and clustering methods, are reviewed by coffey and hinde  <cit> .
methods for differential expression analysis include e.g.
spline-based methods, generalized f-tests and hierarchical error and empirical bayes models.
spline-based edge method by storey et al.
<cit>  is relevant for our problem because it provides comparisons for multiple conditions .
although edge computes a p-value for differential expression, it does not quantify the differential expression for all lineage comparisons, such as reciprocal genes .
anova-based tanova method is based on the approach where different anova structures are defined and the optimal one is found by evaluating the effects and significancies of the factors  <cit> .
recently, stegle et al.
<cit>  proposed an approach based on gaussian processes  to determine the time interval when a gene is differentially expressed.
the methodology of stegle et al.
was limited to analyzing only two conditions.
moreover, it is often observed at transcriptional level that immediately after a treatment, such as activation of t cells by engagement of t cell receptor and cd <dig>  genes are highly dynamic for some time but activity of gene expression decreases at later time points  <cit> .
thus, an ideal computational method − that does not exist at the moment − should take into account the temporal correlation, handle a non-uniform measurement grid, cope with non-stationary processes, and be able to do a well-defined analysis of multiple conditions.
here we developed a computational methodology, ligap  which analyzes experimental data from any number of lineage commitment time-course profiles and analyzed genome-wide gene expression profiles of human umbilical cord blood t helper cells  activated through their cd <dig> and cd <dig> receptors and cultured in absence  or presence of cytokines promoting th <dig> or th <dig> differentiation.
the results give insight into differences of the three lineages in the expression landscape and provide marker genes for lineage commitment identification.
key lineage specific, that is, differentially regulated, genes discovered computationally were validated either experimentally at protein level or based on the published literature.
using a module-based analysis, we identified known and putative regulatory control mechanisms by overlaying highly coherent lineage profile clusters with genome-wide transcription factor  binding predictions and pathway information.
consistent with the previously published results on il-4/stat6-mediated control of a large fraction of genes in th <dig> program  <cit> , our analysis revealed a comparable up-regulated and down-regulated modules, which are suggested to be controlled by stat <dig> and other tfs.
interestingly, we also found that the genes which behave differently between all the lineages studied exhibit a consistent characteristic pattern, i.e., they are up-regulated in th <dig> polarizing cells, down-regulated in th <dig> polarizing cells, and in activated cells  the expression levels are between th <dig> and th <dig> cells.
in addition, our analysis revealed a large set of novel genes, which are specific for different t cell subsets in human.
all the gene expression data and differentially regulated genes as well as software implementing our computational analysis are made publicly available.
RESULTS
experimental data from primary human cd4+ t cells
we used previously published time-course gene expression measurements of activated primary human t cells  and cells polarized to differentiate to th <dig> lineage  <cit>  as well as previously unpublished data set representing th <dig> polarizing cells originating from the same naïve th precursor cells as the th <dig> and th <dig> cells.
the gene expression of th <dig> lineage was measured at time points  <dig>   <dig>   <dig>   <dig> and  <dig> hours.
the measurements from th <dig> and th <dig> samples were available at the same time points.
ligap: a computational technique to identify condition specific time-course profiles
the discovery of condition specific genes at the level of gene expression is an important first step in systems biology studies.
to capture temporal aspects of biological processes, such as cell differentiation, gene expression is commonly measured over time.
we developed a novel model-based method ligap for detecting and visualizing changes between multiple lineage commitment time-course profiles.
briefly, for each gene at a time, our method carries out all comparisons between different cell subsets.
in the case of th <dig>  th <dig> and th <dig> lineages, we assess all  <dig> alternatives;  “th <dig>  th <dig>  th <dig> time-course profiles are all similar” ,  “th <dig> and th <dig> are similar and th <dig> is different” ,  “th <dig> and th <dig> are similar and th <dig> is different” ,  “th <dig> and th <dig> are similar and th <dig> is different” , and  “th <dig>  th <dig>  and th <dig> are all different from each other” .
ligap comparisons and quantifications are illustrated in figure  <dig>  the modeling is done using gaussian processes, which provide a flexible and nonparametric approach for estimating smooth differentiation profiles.
with the help of bayesian statistics, we can quantify differences and similarities by assigning posterior probabilities for all the different profile comparisons between polarizing cell subsets.
the problem can be seen as a model selection problem, where different comparisons are thought of as different model structures  and, given experimental lineage commitment profile data d, the marginal likelihood p, j= <dig> …, <dig>  is used to score different models.
using the bayes’ theorem, the marginal likelihoods can be converted into posterior probabilities of different hypothesis.
these bayesian model scores can be used further to quantify genes, which are specific for a certain lineage.
for example, the probability of a gene being differentially regulated in th <dig> lineage, i.e., score for th <dig> is p = p + p = p + p. genes which are differentially regulated in each of the conditions can be found by quantifying the probabilities p = p or the three probabilities of differential regulation.
each score quantifies the amount of differential regulation, which refers to distinct temporal behavior from other lineages.
the methodology generalizes to any number of lineages/conditions.
our method copes with non-uniform sampling, is able to model non-stationary biological processes , can make comparisons for paired samples, and can carry out the analysis with different number of replicates and missing data.
importantly, the method affords comparison of more than two conditions of interest and is widely applicable to different experimental platforms.
ligap identifies signatures of th <dig>  th <dig> and th <dig> cell lineages
we analyzed the genome-wide gene expression time-course data from th <dig>  th <dig> and th <dig> lineages using ligap.
for all genes, the method outputs the posterior probability values for each of the five hypotheses and also computes the scores for genes being differentially regulated in the th subsets.
an overview of the differentially regulated genes is shown in figure  <dig>  where the four-dimensional data points representing the condition specificities are projected into a plane using the principle component analysis .
this demonstrates the convenience of the presented method as we are able to reduce highly complex data into a meaningful four-dimensional representation using a unified probabilistic framework.
in figure  <dig> individual points represent different genes and every gene is associated with four probabilities: p, p, p, and p. note that ifnγ has the three probabilities p, p, and p close to unity because the probability p is close to unity.
we set a criterion  for the probabilities to call the differentially regulated probe sets; this threshold is in accordance with the jeffrey’s interpretation of “strong evidence” for the bayes factor  <cit> .
in addition, we required a minimum of two-fold change between a lineage and all other lineages at some time point during the differentiation for a gene to be called as differentially regulated.
the top  <dig> and  <dig> gene symbols for th <dig> and th <dig> lineages, respectively, are listed in table  <dig>  whereas, the th <dig> list includes only  <dig> genes.
in a additional file 1: figure s <dig> are depicted two additional examples illustrating the advantage of considering temporal correlation in gene expression and thus improving the sensitivity of detecting consistent yet subtle changes.
in addition, we repeated the analysis using edge  <cit>  and tanova  <cit>  methods using the default parameter values.
tanova identified almost twice as many genes  to be differentially regulated as ligap or tanova .
a comparison of the obtained ranked lists revealed a higher correspondence between the lists produced by ligap and edge than with the list produced by tanova .
probe sets that fulfill the two-fold change criterion are marked based on the direction of the expression  in the given condition.
for example, ifnγ expression is enhanced in th <dig> compared to th <dig> and th <dig>  whereas expression of slc27a <dig> is decreased in th <dig> compared to th <dig> and th <dig>  in addition, probe sets that do not fulfill the fold change criterion are marked in the lists with “<2”.
all genes from th <dig> and th <dig> conditions as well as the top  <dig> of the th <dig> specific genes are shown.
our results of the th subset specific genes agree well with known transcriptional changes during the human t cell differentiation.
ifnγ, a hallmark molecule of th <dig> lineage, was found to be one of the most significantly up-regulated th <dig> specific transcripts .
furthermore, il18r <dig> encoding the interleukin  <dig> receptor , as well as il- <dig> receptor accessory protein  were among the top th <dig> specific genes .
expression of il18r is up-regulated specifically on th <dig> cells but not on th <dig> cells, thus, il18r can be regarded as a differentiation marker for th <dig> cells  <cit> .
in fact, il- <dig> and il- <dig> can reciprocally up-regulate expression of each other’s receptors in th <dig> cells  <cit>  and the il- <dig> - il18r system has a significant role in the synergistic effect of il- <dig> and il- <dig> in triggering efficient nf-κb signaling and enhancement of ifnγ production from human th <dig> cells  <cit> .
intriguingly, in the absence of il- <dig>  il- <dig> has also potential to induce th <dig> differentiation and cytokine response  <cit> .
the basic helix-loop-helix transcription repressor twist <dig> is also known to be expressed in th <dig> cells in il-12/stat <dig>  nf-κb and nfat dependent way and its role has been proposed to be linked to autoregulation of inflammatory cytokine production e.g.
ifnγ  <cit> .
several studies have shown that cxcr <dig> is predominantly expressed in th <dig> cells  <cit>  and, inversely, in th <dig> prone allergic conditions the expression of cxcr <dig> was reduced in allergic patients when compared to healthy individuals  <cit> .
also, an important th <dig> linked function has been observed with map3k <dig> as it acts as an upstream activator of erk via il- <dig> and tcr-dependent signaling, promotes expression of t-bet and stat <dig>  and is actually a stat <dig> target itself forming a feedback loop in the th <dig> cells  <cit> .
deficiency in map3k <dig> leads to decreased ifnγ production in t cells and in vivo impaired host defense against toxoplasma gondii <cit> .
interestingly, the retinoic acid-related orphan receptor gamma  gene encoding rorγt, the key transcription factor in the differentiation program of th <dig> cells, was also identified as a th <dig> specific gene by the ligap analysis  as its expression was up-regulated at  <dig> h time point .
in human, small numbers of t cells producing both il- <dig> and ifnγ have been characterized in peripheral blood, in lamina propria of patients with crohn’s disease as well as in patients with psoriasis  <cit> , but currently is it not known how such cells are derived from the naïve precursor cells.
other novel th <dig> specific hits identified by the ligap include two cytoskeleton associated protein-coding genes dystrophin  and palladin .
dmd encodes actin-binding cytoskeletal structure molecule, which has been mostly studied in patients with duchenne’s muscular dystrophy  <cit> .
these patients develop dystrophin specific autoreactive t cells  <cit> , however, the biological role for dystrophin or palladin in differentiating th cells is not known.
other genes novel in this context and putatively important for th <dig> cell differentiation and/or function include metrnl, , associated with rare cases of mild ring  <dig> syndrome  <cit> , glul encoding a glutamine synthetase, and associated with neuronal disorders and atherosclerotic carotid plaques  <cit> , mctp <dig> , bbs <dig> , stag <dig> , a meiotic gene, as well as pgap <dig> .
napsb coding for aspartic protease napsin b is known to be expressed in human spleen and peripheral blood leucocytes, however, it is estimated to be only a transcribed pseudogene  <cit> .
similarly, miat  is a non-protein coding gene  <cit> , and the relevance of these transcripts in t cell differentiation is not understood, yet.
the top ligap hits of th <dig> specific genes included several genes with very high probability values  and include a vast number of genes that are both specifically up-regulated and down-regulated in th <dig> conditions compared to other th subsets.
therefore, the list of th <dig> specific genes with highest probability is consistent with the previously published results obtained with other computational methods  <cit> .
importantly, gata <dig>  the well characterized master transcription regulator of th <dig> polarization  <cit>  was identified among the top th <dig> hits .
the transcriptional expression profile of gata <dig> was observed to be highly up-regulated at all time points among the cells cultured in th <dig> polarizing conditions, whereas the expression profiles in th <dig> and th <dig> cells exhibited down-regulation .
in addition to well-known subset signature molecules, the analysis identified also a number of poorly characterized molecules in relation to their function in polarized th cells.
among the highly expressed top  <dig> th <dig> hits, the specificity of these transcripts relative to th <dig>  but not to th <dig>  has already been identified at different time points with the standard limma methods  in the past  <cit> .
one of these th <dig> specific top hits was maoa, a gene encoding monoamine oxidase a, whose expression was increasingly up-regulated during the time course .
this enzyme degrades amine neurotransmitters,  and was previously found to be up-regulated in human peripheral blood monocytes after il- <dig> and il- <dig> stimulation  <cit>  as well as in th <dig> cells derived from cord blood naïve cd4+ t cells and, importantly, being indirectly controlled by stat <dig>  <cit> .
it has been hypothesized that maoa may act as an anti-inflammatory mediator by degrading serotonin which inhibits generation of tnfα from macrophages and up-regulates phagocytosis  <cit> .
the biological significance of maoa in th <dig> cells, however, remains to be elucidated.
another interesting th <dig> specific top hit was spint <dig>  encoding a transmembrane serine peptidase inhibitor kunitz type  <dig> .
spint <dig> was originally named after its homology to hepatocyte growth factor activator inhibitor  <dig> and its first isolation from human placenta  <cit> .
the kunitz inhibitory domains display potent inhibitory activity towards several trypsin-like serine proteases  <cit>  and mutations in the human spint <dig> gene cause a broad spectrum of abnormalities in organogenesis  <cit> .
in addition, spint <dig> may function as a tumor suppressor gene, as its mrna levels are down-regulated in several human cancers  and a deficiency in spint <dig> expression is linked with poor prognosis of breast cancer  <cit> .
there are no previous studies where the possible functional role of spint <dig> in human lymphocytes is unraveled, however, spint <dig> was recently found to be a stat <dig> target in human macrophages as well as in human th <dig> cells  <cit> .
we, hence, chose to experimentally validate the specificity of spint <dig> in primary human th2-polarizing cells.
we tested the specificity of spint <dig> expression at protein level on the cell surface of the th cells with flow cytometry.
at  <dig> hours after activation and induction of polarization the th <dig> cells were found to express significantly more spint <dig> than the th <dig> polarizing cells or the activated th <dig> cells .
as some of the human spint <dig> transcripts do not harbor the coding signal for the transmembrane domain  <cit> , we therefore also investigated if spint <dig> would be secreted from the th subsets.
the spint <dig> concentrations were measured from the culture supernatants by enzyme-linked immunosorbent assay  at  <dig> hours after activation and polarization, and the th <dig> cells were observed to secrete significantly more spint <dig> than th <dig> or th <dig> cells .
the th <dig> specific hits included also ppp1r14a, a phosphorylation-dependent inhibitor of smooth muscle myosin phosphatase, involved in regulation of smooth muscle contraction  <cit>  as well as dusp <dig> , responsible for dephosphorylation of erk1/ <dig>  <cit> .
recently, il- <dig> induced rna expression of signaling molecules ppp1r14a and dusp <dig> have been reported  <cit> .
as the regulation of phosphorylation of the signaling intermediates is known to be highly important in defining the cell differentiation, we wanted to experimentally validate the subset specific expression of these two signaling molecules at protein level.
we detected a clear th <dig> specific ppp1r14a and dusp <dig> protein expression at  <dig> hours time point post activation and initiation of the polarization, and very little or no expression in th <dig> and th <dig> lineages .
reciprocal regulators of lineage commitment
in context of determination of t cell subset identity, a key group of genes is the one where the expression kinetics differ between all the lineages.
the list of these significantly different genes is shown in table  <dig>  an illustrative example gene from this list is the well-known th <dig> signature cytokine gene ifng  as well as tbx <dig> encoding t-bet, a hallmark transcription factor in th <dig> differentiated cells, both of which are also known to suppress th <dig> activity  <cit> .
in addition, map3k <dig>  fas, il12rb <dig>  and il- <dig>  have been identified to play role in th <dig> polarized cells .
moreover, table  <dig> and additional file 2: table s <dig> contain numerous differentially regulated transcripts which are only poorly characterized or their role in cd4+ th cells has not been studied.
the novel th <dig> specific genes dmd and palld, encoding cytoskeletal associated proteins dystrophin and palladin, fall into the reciprocally regulated genes in the th subsets studied here.
also, th <dig> specific putative pseudogene napsb and non-coding transcript miat show reciprocal transcript profiles.
other novel genes include prr5l, which has been identified to interact with a highly conserved protein kinase tor , a central controller of cell growth and apoptosis  <cit> .
osbpl <dig> encodes oxysterol binding protein-like  <dig>  an intracellular lipid receptor that regulates cellular lipid metabolism  <cit> .
p2ry <dig>  is a membrane receptor for udp-glucose and plays a role in immune responses in human airway as well as female reproductive track epithelial cells by stimulating cytokine and chemokine production and recruitment of neutrophils  <cit> .
p2ry <dig> has also been identified to function in mouse splenic t cells as a regulator of il- <dig> induced proliferation, however, no specific link to th <dig> cells has been observed  <cit> .
also, the significance of atp9a , lpar <dig>  functioning in g-protein coupled receptor signaling, xrn <dig> , bspry , mctp <dig>  or ptpro  in th <dig> cells is yet to be studied.
recent data indicate that in b cells, ptpro dephosphorylates syk, a kinase that is critical in signal transduction of b-cell receptor  <cit> .
probe sets that fulfill the two-fold change criterion are marked with the + and - signs following the gene symbols based on the direction of the expression where + denotes up-regulation and - denotes down-regulation of th <dig> specific genes.
in addition, probe sets that do not fulfill the fold change criterion are marked in the lists with “<2”.
the known associations of genes in th cell functions or subset specific gene expressions are listed in the table.
*) annotation based on ingenuity pathway analysis ® by ingenuity systems.
nr, not reported.
the th <dig> up-regulated genes, pde7b, setbp <dig>  c9orf <dig>  tprg <dig>  igsf <dig>  or ppp1r14a have not been linked to cd4+ th cell function, although their il- <dig> mediated up-regulation has been published, and furthermore, setbp <dig>  tprg <dig> and ppp1r14a have been identified as direct targets of stat <dig>  <cit> .
interestingly, we observed that most of the genes whose expression differs between all the three lineages behave in a similar manner, i.e., they are up-regulated in th <dig> and down-regulated in th <dig>
among the reciprocally regulated genes we found  <dig> genes up-regulated in th <dig> condition and only six genes behaved in the opposite manner.
the hierarchical clustering of the kinetic profiles is depicted in figure 5a.
this suggests that there are common mechanisms that induce reverse regulatory behavior.
for example, the genes up-regulated in th <dig> condition might be controlled downstream of ifnγ.
this hypothesis is supported by the clear similarity between the profiles of ifnγ and the profiles of the clustered genes.
we prepared a similar figure showing the differences in the kinetics of all the ligap identified genes.
these results are depicted in figure 5b and they show the similarity between the th <dig> and th <dig> lineages and their dissimilarity between the th <dig> lineage.
transcription factor binding sites in th <dig> lineage
to extend our transcriptional analysis into transcriptional regulation, we decided to systematically analyze both genome-wide transcription factor  binding site predictions made in silico and comprehensive literature-derived information about target genes of selected tfs.
first, we predicted which of the transcription factors have binding sites in the refseq gene promoters  using the probtf tool  <cit>  combined with an empirical p-value computation.
we focused on genes that were identified by the previous ligap analysis and considered all transcription factors that had known binding specificities  in transfac  <cit>  .
we did not restrict our analysis only to those tfs whose transcripts are differentially expressed because, e.g., stat <dig> is not differentially expressed during the early differentiation although it is a master regulator in the early differentiation of th <dig> cells  <cit> .
an important goal is to identify master regulators of the lineage commitment processes.
recently, it was found out that most of the direct targets of stat <dig>  an important regulator of th <dig> differentiation, were up-regulated in th <dig> cells  <cit> .
here we were interested in identifying tfs whose binding sites are enriched in the promoter regions of the genes which are differentially regulated in th <dig> conditions, both among the up-regulated and down-regulated genes.
instead of looking at individual tf binding predictions that are prone to contain false positives, we used the fisher’s exact test to search for enrichment of binding sites, in comparison to randomly selected gene set.
the same analysis was carried out separately for all the differentially regulated gene sets and by taking into account the direction of regulation .
using a p-value cut-off of  <dig>  for tf binding, we identified three hits from the enrichment analysis among th <dig> specific up-regulated genes and three among the th <dig> specific down-regulated genes.
the results are depicted in figure  <dig>  the different enriched irf family motifs were combined and their targets were pooled.
in accordance to our previously published results  <cit> , the strongest hit within the th <dig> up-regulated genes was stat <dig> , followed by nkx3a , and cdp .
nkx3a is a member of the nkx family of homeobox genes that is expressed in prostate epithelium and functions as a potential prostate tumor suppressor  <cit> .
recently, in a study focusing on jurkat cells, a gata <dig> binding site on the promoter of nkx <dig> gene was identified  <cit> .
furthermore, in mouse increased expression of nkx3a was observed to be regulated by il- <dig> independently of stat <dig>  <cit> .
cdp  is highly conserved homeodomain transcription factor involved in many cellular processes, including differentiation, development and proliferation  <cit> .
interestingly, cdp has been identified as a repressive regulator of cd8α silencer region and tcrβ enhancer region and plays a role in promoting repressive chromatin modifications via association with histone deacetylase  <dig> and histone  <dig> methyltransferase  <cit> .
it is important to notice that the binding motif analysis does not give information about the possible direction of regulation, e.g., it is an open question whether cdp might up-regulate th <dig> specific genes or down-regulate the genes in th <dig> and th <dig> lineages.
the three tf hits having enriched predicted binding sites among the th <dig> down-regulated genes were the interferon regulatory factor  family of tfs , ifn-stimulated genes factor  <dig>   and stat <dig> .
irf family consists of irf <dig> to irf <dig> and has been shown to be essential in expression of type i interferon genes, ifn-stimulated genes  and other pro-inflammatory response related cytokines  <cit> .
these genes are maintained down-regulated during th <dig> proliferation and therefore, the results are in line with the th <dig> effector cells characteristics  <cit> .
moreover, ifnγ-induced expression of irf <dig> and irf <dig> has been shown to directly down-regulate il- <dig> production by repressing il- <dig> promoter sites  <cit> .
opposing to other irf family proteins, irf <dig> has been shown to directly activate il- <dig> promoter and il- <dig> regulatory elements and be essential in th <dig> cell differentiation by influencing the expression of gfi <dig>  a transcriptional repressor in th <dig> cells  <cit> .
however, the analysis relying on known tf binding specificities will not allow segregation of individual members of the irf family.
further, an essential regulator of most isgs is isgf <dig> that is composed of stat <dig>  stat <dig> and irf <dig> complex and works in conjunction with irfs  <cit> .
identification of stat <dig> as a regulator among the th <dig> down-regulated genes is well in line with our previously published results, although its effect was observed to be less profound within th <dig> down-regulated genes than among th <dig> up-regulated target genes  <cit> .
comparison analysis of the predicted stat <dig> target genes and th <dig> up-regulated and down-regulated genes gave  <dig> and  <dig> overlapping genes, respectively.
the full lists of overlapping genes are in additional file 3: table s <dig>  we further analyzed the correlation between predicted stat <dig> target promoters and experimentally observed promoter associated binding sites , and observed significant correlation  between the target sites.
the full list of predicted stat <dig> target genes and promoter associated stat <dig> binding sites identified by chip-seq as well as the overlapping genes are listed in the additional file 3: table s <dig>  the overlapping binding sites included promoters for c14orf <dig>  cish, hmmr, ino <dig>  mgat <dig>  nudcd <dig>  socs <dig>  spint <dig> and znf <dig> genes.
discussion
identification of the key t helper cell regulators provides possible targets for modulation of immune response.
to reveal t cell subset specific genes and their often subtle differences in expression, we developed a novel computational method, ligap.
traditional ways of identifying differentially expressed genes, such as the t-test, are problematic in studying time-series data since there is a need to carry out hypothesis tests on individual time points.
on the other hand, commonly used statistical tests for whole time-course, including e.g.
f-test, do not account for the inherent correlation between measurement time points.
ligap overcomes many problems that have previously prevented quantitative comparisons of multiple differentiation profiles, with or without replicates.
among several beneficial features, ligap models correlation between time points and can cope with non-stationarities and non-uniform measurement grid.
other methods, such as edge, uses splines to estimate smooth time-course profiles but does not quantify the differential expression for all lineage comparisons.
tanova uses standard regression framework and lacks explicit correlation structure between time points.
our study highlights the validity of the method by identifying known and novel differentially regulated genes and their kinetic differences during t helper cell differentiation.
in addition, the non-parametric computational analysis automatically provides informative illustrations of time-course profiles together with associated uncertainty.
ligap calculated th <dig> specific gene set contains only  <dig> genes and th <dig> specific  <dig> genes compared to  <dig> genes that are specific to th <dig> conditions.
activation of thp cells through tcr and cd <dig> results in induction of ifnγ, which in turn leads to activation of th <dig> signature genes.
addition of il- <dig>  however, results in enhanced induction of these genes and th <dig> programming.
consistent with our previous results genes differentially regulated in response to th <dig> programming are much more limited than those detected in response to initiation of th <dig> response  <cit> .
most of the th <dig> specific genes encode well-known th <dig> signature molecules.
however, also genes new in this context were discovered.
interestingly, we identified rorc as one of the th <dig> specific genes.
up-regulation of rorc in th <dig> cells and existence of th17/th <dig> cells, however, remain conflicting as the master regulator of th <dig> differentiation, t-bet, is known to inhibit transcription of rorc through runx <dig>  <cit> , and expression of il12rβ <dig> is down-regulated by il- <dig>  <cit> .
it has been suggested that the high concentration of tgfβ required for in vitro th <dig> polarization would inhibit ifnγ production  <cit> , hence, it remains an open question whether some conditions would drive the differentiation of il- <dig> and ifnγ producing cells from same naïve precursor t cell.
notably, ex vivo th <dig> cells could be induced to develop further into th17/th <dig> cells by the combined actions of ifnγ and il- <dig>  and such conditions resulted in permissive chromatin remodeling at the il12rb <dig> locus and loss of repressive histone modification at the tbx <dig> locus  <cit> .
as an example of previously uncharacterized differentially regulated genes, we validated the expression of th2-associated phosphatases dusp <dig> and ppp1r14a on protein level.
ppp1r14a was shown in human pancreatic and melanoma tumor cell lines to positively regulate ras/mapk signaling  <cit> , which are also involved in il- <dig> induced signaling cascades.
in t cells, the erks are activated though tcr stimulation and a tcr-mediated activation of ras/mapk signaling is required in differentiating murine th <dig> but not in th <dig> cells  <cit> .
furthermore, the ras/mapk cascade was shown to enhance the stability of gata <dig> protein  <cit>  as well as stat <dig> independent cd <dig> and cd <dig> induced initial il <dig> production  <cit> .
dusp <dig> on other hand is known to negatively regulate members of the mitogen-activated protein  kinase superfamily associated with cellular proliferation and differentiation  <cit> .
more specifically, dusp <dig> expression was shown to be induced by erk1/ <dig> signaling in differentiating mouse embryonic cell line and in human retinal pigment epithelial cells  <cit>  and it was hypothesized that dusp <dig> is an essential part of a negative feedback loop of erk1/ <dig> signaling  <cit> .
however, the t cell associated functions of both ppp1r14a and dusp <dig> are completely unknown.
therefore, their significance in the signaling cascades of differentiating th <dig> cells remains a highly interesting area of future research.
spint <dig> was recently identified as a direct stat <dig> target in differentiating human th <dig> cells  <cit>  and in this study we are the first to show that spint <dig> is upregulated in th <dig> cells at protein level as compared to other th cell subsets.
we found spint <dig> to be specifically expressed on th <dig> cell surface as well as secreted into the culture medium, suggesting presence of a multiple transcripts of which some may lack the anchoring transmembrane domain  <cit> .
human spint <dig>  is a physiological inhibitor of matrix cleaving proteases and decreased expression of spint <dig> has been linked to progression of several cancers  <cit> .
up-regulated expression of extracellular proteases is crucial for pro-cancerous pathways as this enables efficient remodeling of the extracellular matrix as well as cleavage and activation of growth factors and their receptors.
interestingly, a truncated and secreted spint <dig> may act as an inhibitor for the activator of hepatocyte growth factor  and hgf is prominently expressed in lung tissue and is linked to reduced expression of th <dig> cytokines and tgfβ, reduction of allergic airway inflammation, airway hyperresponsiveness and remodeling as well as reduced recruitment of eosinophils to the site of allergic inflammation in vivo <cit> .
this suggests that spint <dig> might enhance th <dig> response in allergic airway inflammation by inhibiting hgf signaling.
the ligap method elegantly identified the reciprocally regulated genes within the th <dig>  th <dig> and th <dig> conditions.
essentially, the list included genes encoding the hallmark th <dig> specific transcription factor t-bet and cytokine ifnγ as well as the transmembrane receptor for il- <dig>  this list also included few cytoskeleton associated proteins, such as dystrophin , and palladin , of which there is no current knowledge for their function in differentiating t helper cells.
the observation suggests differences in cellular structures or putatively in the interaction of apc with the th cell subsets as rearrangement of the cytoskeleton in t cells plays an important role in the organization of the immunological synapse  and th <dig> and th <dig> cells are known to form morphologically distinct iss  <cit> .
in addition to map3k <dig>  molecules that participate in phosphorylation signaling cascades e.g.
p2ry <dig>  lpar <dig>  ppp1r14a, and ptpro suggest their potential role for initiation or regulation of differentiation cascades.
importantly, the results presented here enable opportunities for further data mining and follow-up studies addressing the functions and importance of the novel th subset specific genes.
the identification of stat <dig> as the most significant tf regulating th <dig> specific enhancement of transcription by the tf binding analysis is well in line with our previous stat <dig> chip results  <cit> .
furthermore, the analysis between the predicted stat <dig> target gene promoters and experimentally observed promoter associated binding sites showed statistically significant correlation.
interestingly, the overlapping stat <dig> targets included ino <dig>  which has been identifies as a part of a chromatin remodeling complex  <cit>  and may hence, be involved in th <dig> specific epigenetical regulation of th cell differentiation.
stat <dig> specific regulation of mannosyl -glycoprotein beta- <dig> -n-acetylglucosaminyltransferase , a n-glycan-processing enzyme  <cit> , may on one hand be involved in modifying the th <dig> cell specific surface glycoprotein structures  <cit> .
the overlapping target sites included also the promoter for spint <dig>  the number of predicted stat <dig> binding sites, however, was much larger than the experimentally observed binding sites, which may reflect the typically observed high false positive rate of computational binding predictions and the cell type specific state of chromatin as well as other competing factors affecting binding in vitro.
the data created here also further suggests novel control mechanisms involving gata <dig> regulated nkx3a as well as chromatin modification associated cdp.
only less than 10% of the th <dig> down-regulated genes were reported to be direct targets of stat <dig> by elo et al.,  suggesting other major regulatory mechanisms play role among the il- <dig> induced down-regulated genes.
we found enrichment of irf family and isgf <dig> binding motifs in promoter regions of genes that are repressed in th <dig> polarizing conditions, indicating that these tfs may play a significant role in the suppressing undesired gene expression in differentiating th <dig> cells.
indeed, several irf family members have been identified as differentially expressed during th cell differentiation and necessary for both th <dig> and th <dig> polarization.
as the irf family proteins, excluding irf <dig>  share the same binding specificity model in transfac, the individual regulatory role for these factors is, however, difficult to postulate based on in silico tf binding site analysis.
CONCLUSIONS
the proposed ligap method can quantify a well-defined probabilistic specificity score for each gene and for each condition promoting a certain lineage commitment.
in addition to grouping and ranking genes based on their dynamics, ligap summarizes all time-course measurements, together with the associated uncertainty, in an illustrative summary plot for visualization and manual assessment purposes.
while here we have demonstrated the utility of ligap in analysis of gene expression dynamics, the ligap method is widely applicable to many types of datasets including quantitative time-course experiments and generalizes to any number of conditions.

 t cells are key regulators of the adaptive immune system and have a central role in defense against pathogens and cancer as well as protection from autoimmune diseases. cd4+ t lymphocytes can differentiate to functionally distinct effector subtypes, including t helper, t helper and more recently described t helper cells. th cells secrete effector cytokine ifn-γ and regulate cell-mediated immunity and play a role in the pathogenesis of autoimmune diseases, such as multiple sclerosis. th cells in turn produce il- il- and il- cytokines, and mediate immunity against extracellular pathogens and allergic reactions. th cells, characterized by the production of a proinflammatory cytokine il- regulate inflammatory responses on the mucosal surfaces. for the overall health in humans and animals, the proper balance between different effector t cell types and t regulatory cells is crucial. aberrant activation of th and th or th cells can trigger inflammatory autoimmune diseases as well as asthma and allergy. previous studies utilizing genome-wide expression data and computational modeling have aimed at revealing the master regulators and regulatory networks within the differentiating th and th cells. however, studies in human have been less extensive than in mouse due to the difficulty in collecting sufficient amount of samples to comprehensively profile t cell differentiation over time. in addition, lack of appropriate computational methods suitable for analyzing large-scale experimental data from multiple lineages over several time points spanning the lineage commitment process has limited the progress on revealing dynamics and molecular mechanisms underlying multiple lineage commitment. a number of different time-series analysis approaches have been proposed to solve large-scale lineage commitment analysis problems. the general purpose f-test can be used to test the difference between time-series data sets, but it does not extend to simultaneous comparison of multiple lineages and fails to take into account the correlation between the measurements at different time points. more recent approaches to analyze time-series data, including regression, differential expression, discriminant and clustering methods, are reviewed by coffey and hinde. methods for differential expression analysis include e.g. spline-based methods, generalized f-tests and hierarchical error and empirical bayes models. spline-based edge method by storey et al. is relevant for our problem because it provides comparisons for multiple conditions. although edge computes a p-value for differential expression, it does not quantify the differential expression for all lineage comparisons, such as reciprocal genes. anova-based tanova method is based on the approach where different anova structures are defined and the optimal one is found by evaluating the effects and significancies of the factors. recently, stegle et al. proposed an approach based on gaussian processes to determine the time interval when a gene is differentially expressed. the methodology of stegle et al. was limited to analyzing only two conditions. moreover, it is often observed at transcriptional level that immediately after a treatment, such as activation of t cells by engagement of t cell receptor and cd genes are highly dynamic for some time but activity of gene expression decreases at later time points. thus, an ideal computational method − that does not exist at the moment − should take into account the temporal correlation, handle a non-uniform measurement grid, cope with non-stationary processes, and be able to do a well-defined analysis of multiple conditions. here we developed a computational methodology, ligap which analyzes experimental data from any number of lineage commitment time-course profiles and analyzed genome-wide gene expression profiles of human umbilical cord blood t helper cells activated through their cd and cd receptors and cultured in absence or presence of cytokines promoting th or th differentiation. the results give insight into differences of the three lineages in the expression landscape and provide marker genes for lineage commitment identification. key lineage specific, that is, differentially regulated, genes discovered computationally were validated either experimentally at protein level or based on the published literature. using a module-based analysis, we identified known and putative regulatory control mechanisms by overlaying highly coherent lineage profile clusters with genome-wide transcription factor binding predictions and pathway information. consistent with the previously published results on il-4/stat6-mediated control of a large fraction of genes in th program, our analysis revealed a comparable up-regulated and down-regulated modules, which are suggested to be controlled by stat and other tfs. interestingly, we also found that the genes which behave differently between all the lineages studied exhibit a consistent characteristic pattern, i.e., they are up-regulated in th polarizing cells, down-regulated in th polarizing cells, and in activated cells the expression levels are between th and th cells. in addition, our analysis revealed a large set of novel genes, which are specific for different t cell subsets in human. all the gene expression data and differentially regulated genes as well as software implementing our computational analysis are made publicly available.  experimental data from primary human cd4+ t cells we used previously published time-course gene expression measurements of activated primary human t cells and cells polarized to differentiate to th lineage as well as previously unpublished data set representing th polarizing cells originating from the same naïve th precursor cells as the th and th cells. the gene expression of th lineage was measured at time points and hours. the measurements from th and th samples were available at the same time points. ligap: a computational technique to identify condition specific time-course profiles the discovery of condition specific genes at the level of gene expression is an important first step in systems biology studies. to capture temporal aspects of biological processes, such as cell differentiation, gene expression is commonly measured over time. we developed a novel model-based method ligap for detecting and visualizing changes between multiple lineage commitment time-course profiles. briefly, for each gene at a time, our method carries out all comparisons between different cell subsets. in the case of th th and th lineages, we assess all alternatives; “th th th time-course profiles are all similar”, “th and th are similar and th is different”, “th and th are similar and th is different”, “th and th are similar and th is different”, and “th th and th are all different from each other”. ligap comparisons and quantifications are illustrated in figure  the modeling is done using gaussian processes, which provide a flexible and nonparametric approach for estimating smooth differentiation profiles. with the help of bayesian statistics, we can quantify differences and similarities by assigning posterior probabilities for all the different profile comparisons between polarizing cell subsets. the problem can be seen as a model selection problem, where different comparisons are thought of as different model structures and, given experimental lineage commitment profile data d, the marginal likelihood p, j= …, is used to score different models. using the bayes’ theorem, the marginal likelihoods can be converted into posterior probabilities of different hypothesis. these bayesian model scores can be used further to quantify genes, which are specific for a certain lineage. for example, the probability of a gene being differentially regulated in th lineage, i.e., score for th is p = p + p = p + p. genes which are differentially regulated in each of the conditions can be found by quantifying the probabilities p = p or the three probabilities of differential regulation. each score quantifies the amount of differential regulation, which refers to distinct temporal behavior from other lineages. the methodology generalizes to any number of lineages/conditions. our method copes with non-uniform sampling, is able to model non-stationary biological processes, can make comparisons for paired samples, and can carry out the analysis with different number of replicates and missing data. importantly, the method affords comparison of more than two conditions of interest and is widely applicable to different experimental platforms. ligap identifies signatures of th th and th cell lineages we analyzed the genome-wide gene expression time-course data from th th and th lineages using ligap. for all genes, the method outputs the posterior probability values for each of the five hypotheses and also computes the scores for genes being differentially regulated in the th subsets. an overview of the differentially regulated genes is shown in figure  where the four-dimensional data points representing the condition specificities are projected into a plane using the principle component analysis. this demonstrates the convenience of the presented method as we are able to reduce highly complex data into a meaningful four-dimensional representation using a unified probabilistic framework. in figure  individual points represent different genes and every gene is associated with four probabilities: p, p, p, and p. note that ifnγ has the three probabilities p, p, and p close to unity because the probability p is close to unity. we set a criterion for the probabilities to call the differentially regulated probe sets; this threshold is in accordance with the jeffrey’s interpretation of “strong evidence” for the bayes factor. in addition, we required a minimum of two-fold change between a lineage and all other lineages at some time point during the differentiation for a gene to be called as differentially regulated. the top and gene symbols for th and th lineages, respectively, are listed in table  whereas, the th list includes only genes. in a additional file 1: figure s are depicted two additional examples illustrating the advantage of considering temporal correlation in gene expression and thus improving the sensitivity of detecting consistent yet subtle changes. in addition, we repeated the analysis using edge and tanova methods using the default parameter values. tanova identified almost twice as many genes to be differentially regulated as ligap or tanova. a comparison of the obtained ranked lists revealed a higher correspondence between the lists produced by ligap and edge than with the list produced by tanova. probe sets that fulfill the two-fold change criterion are marked based on the direction of the expression in the given condition. for example, ifnγ expression is enhanced in th compared to th and th whereas expression of slc27a is decreased in th compared to th and th in addition, probe sets that do not fulfill the fold change criterion are marked in the lists with “<2”. all genes from th and th conditions as well as the top of the th specific genes are shown. our results of the th subset specific genes agree well with known transcriptional changes during the human t cell differentiation. ifnγ, a hallmark molecule of th lineage, was found to be one of the most significantly up-regulated th specific transcripts. furthermore, il18r encoding the interleukin receptor, as well as il- receptor accessory protein were among the top th specific genes. expression of il18r is up-regulated specifically on th cells but not on th cells, thus, il18r can be regarded as a differentiation marker for th cells. in fact, il- and il- can reciprocally up-regulate expression of each other’s receptors in th cells and the il- - il18r system has a significant role in the synergistic effect of il- and il- in triggering efficient nf-κb signaling and enhancement of ifnγ production from human th cells. intriguingly, in the absence of il- il- has also potential to induce th differentiation and cytokine response. the basic helix-loop-helix transcription repressor twist is also known to be expressed in th cells in il-12/stat nf-κb and nfat dependent way and its role has been proposed to be linked to autoregulation of inflammatory cytokine production e.g. ifnγ. several studies have shown that cxcr is predominantly expressed in th cells and, inversely, in th prone allergic conditions the expression of cxcr was reduced in allergic patients when compared to healthy individuals. also, an important th linked function has been observed with map3k as it acts as an upstream activator of erk via il- and tcr-dependent signaling, promotes expression of t-bet and stat and is actually a stat target itself forming a feedback loop in the th cells. deficiency in map3k leads to decreased ifnγ production in t cells and in vivo impaired host defense against toxoplasma gondii. interestingly, the retinoic acid-related orphan receptor gamma gene encoding rorγt, the key transcription factor in the differentiation program of th cells, was also identified as a th specific gene by the ligap analysis as its expression was up-regulated at h time point. in human, small numbers of t cells producing both il- and ifnγ have been characterized in peripheral blood, in lamina propria of patients with crohn’s disease as well as in patients with psoriasis, but currently is it not known how such cells are derived from the naïve precursor cells. other novel th specific hits identified by the ligap include two cytoskeleton associated protein-coding genes dystrophin and palladin. dmd encodes actin-binding cytoskeletal structure molecule, which has been mostly studied in patients with duchenne’s muscular dystrophy. these patients develop dystrophin specific autoreactive t cells, however, the biological role for dystrophin or palladin in differentiating th cells is not known. other genes novel in this context and putatively important for th cell differentiation and/or function include metrnl,, associated with rare cases of mild ring syndrome, glul encoding a glutamine synthetase, and associated with neuronal disorders and atherosclerotic carotid plaques, mctp, bbs, stag, a meiotic gene, as well as pgap. napsb coding for aspartic protease napsin b is known to be expressed in human spleen and peripheral blood leucocytes, however, it is estimated to be only a transcribed pseudogene. similarly, miat is a non-protein coding gene, and the relevance of these transcripts in t cell differentiation is not understood, yet. the top ligap hits of th specific genes included several genes with very high probability values and include a vast number of genes that are both specifically up-regulated and down-regulated in th conditions compared to other th subsets. therefore, the list of th specific genes with highest probability is consistent with the previously published results obtained with other computational methods. importantly, gata the well characterized master transcription regulator of th polarization was identified among the top th hits. the transcriptional expression profile of gata was observed to be highly up-regulated at all time points among the cells cultured in th polarizing conditions, whereas the expression profiles in th and th cells exhibited down-regulation. in addition to well-known subset signature molecules, the analysis identified also a number of poorly characterized molecules in relation to their function in polarized th cells. among the highly expressed top th hits, the specificity of these transcripts relative to th but not to th has already been identified at different time points with the standard limma methods in the past. one of these th specific top hits was maoa, a gene encoding monoamine oxidase a, whose expression was increasingly up-regulated during the time course. this enzyme degrades amine neurotransmitters, and was previously found to be up-regulated in human peripheral blood monocytes after il- and il- stimulation as well as in th cells derived from cord blood naïve cd4+ t cells and, importantly, being indirectly controlled by stat. it has been hypothesized that maoa may act as an anti-inflammatory mediator by degrading serotonin which inhibits generation of tnfα from macrophages and up-regulates phagocytosis. the biological significance of maoa in th cells, however, remains to be elucidated. another interesting th specific top hit was spint encoding a transmembrane serine peptidase inhibitor kunitz type. spint was originally named after its homology to hepatocyte growth factor activator inhibitor and its first isolation from human placenta. the kunitz inhibitory domains display potent inhibitory activity towards several trypsin-like serine proteases and mutations in the human spint gene cause a broad spectrum of abnormalities in organogenesis. in addition, spint may function as a tumor suppressor gene, as its mrna levels are down-regulated in several human cancers and a deficiency in spint expression is linked with poor prognosis of breast cancer. there are no previous studies where the possible functional role of spint in human lymphocytes is unraveled, however, spint was recently found to be a stat target in human macrophages as well as in human th cells. we, hence, chose to experimentally validate the specificity of spint in primary human th2-polarizing cells. we tested the specificity of spint expression at protein level on the cell surface of the th cells with flow cytometry. at hours after activation and induction of polarization the th cells were found to express significantly more spint than the th polarizing cells or the activated th cells. as some of the human spint transcripts do not harbor the coding signal for the transmembrane domain, we therefore also investigated if spint would be secreted from the th subsets. the spint concentrations were measured from the culture supernatants by enzyme-linked immunosorbent assay at hours after activation and polarization, and the th cells were observed to secrete significantly more spint than th or th cells. the th specific hits included also ppp1r14a, a phosphorylation-dependent inhibitor of smooth muscle myosin phosphatase, involved in regulation of smooth muscle contraction as well as dusp, responsible for dephosphorylation of erk1/. recently, il- induced rna expression of signaling molecules ppp1r14a and dusp have been reported. as the regulation of phosphorylation of the signaling intermediates is known to be highly important in defining the cell differentiation, we wanted to experimentally validate the subset specific expression of these two signaling molecules at protein level. we detected a clear th specific ppp1r14a and dusp protein expression at hours time point post activation and initiation of the polarization, and very little or no expression in th and th lineages. reciprocal regulators of lineage commitment in context of determination of t cell subset identity, a key group of genes is the one where the expression kinetics differ between all the lineages. the list of these significantly different genes is shown in table  an illustrative example gene from this list is the well-known th signature cytokine gene ifng as well as tbx encoding t-bet, a hallmark transcription factor in th differentiated cells, both of which are also known to suppress th activity. in addition, map3k fas, il12rb and il- have been identified to play role in th polarized cells. moreover, table  and additional file 2: table s contain numerous differentially regulated transcripts which are only poorly characterized or their role in cd4+ th cells has not been studied. the novel th specific genes dmd and palld, encoding cytoskeletal associated proteins dystrophin and palladin, fall into the reciprocally regulated genes in the th subsets studied here. also, th specific putative pseudogene napsb and non-coding transcript miat show reciprocal transcript profiles. other novel genes include prr5l, which has been identified to interact with a highly conserved protein kinase tor, a central controller of cell growth and apoptosis. osbpl encodes oxysterol binding protein-like an intracellular lipid receptor that regulates cellular lipid metabolism. p2ry is a membrane receptor for udp-glucose and plays a role in immune responses in human airway as well as female reproductive track epithelial cells by stimulating cytokine and chemokine production and recruitment of neutrophils. p2ry has also been identified to function in mouse splenic t cells as a regulator of il- induced proliferation, however, no specific link to th cells has been observed. also, the significance of atp9a, lpar functioning in g-protein coupled receptor signaling, xrn, bspry, mctp or ptpro in th cells is yet to be studied. recent data indicate that in b cells, ptpro dephosphorylates syk, a kinase that is critical in signal transduction of b-cell receptor. probe sets that fulfill the two-fold change criterion are marked with the + and - signs following the gene symbols based on the direction of the expression where + denotes up-regulation and - denotes down-regulation of th specific genes. in addition, probe sets that do not fulfill the fold change criterion are marked in the lists with “<2”. the known associations of genes in th cell functions or subset specific gene expressions are listed in the table. *) annotation based on ingenuity pathway analysis ® by ingenuity systems. nr, not reported. the th up-regulated genes, pde7b, setbp c9orf tprg igsf or ppp1r14a have not been linked to cd4+ th cell function, although their il- mediated up-regulation has been published, and furthermore, setbp tprg and ppp1r14a have been identified as direct targets of stat. interestingly, we observed that most of the genes whose expression differs between all the three lineages behave in a similar manner, i.e., they are up-regulated in th and down-regulated in th among the reciprocally regulated genes we found genes up-regulated in th condition and only six genes behaved in the opposite manner. the hierarchical clustering of the kinetic profiles is depicted in figure 5a. this suggests that there are common mechanisms that induce reverse regulatory behavior. for example, the genes up-regulated in th condition might be controlled downstream of ifnγ. this hypothesis is supported by the clear similarity between the profiles of ifnγ and the profiles of the clustered genes. we prepared a similar figure showing the differences in the kinetics of all the ligap identified genes. these results are depicted in figure 5b and they show the similarity between the th and th lineages and their dissimilarity between the th lineage. transcription factor binding sites in th lineage to extend our transcriptional analysis into transcriptional regulation, we decided to systematically analyze both genome-wide transcription factor binding site predictions made in silico and comprehensive literature-derived information about target genes of selected tfs. first, we predicted which of the transcription factors have binding sites in the refseq gene promoters using the probtf tool combined with an empirical p-value computation. we focused on genes that were identified by the previous ligap analysis and considered all transcription factors that had known binding specificities in transfac. we did not restrict our analysis only to those tfs whose transcripts are differentially expressed because, e.g., stat is not differentially expressed during the early differentiation although it is a master regulator in the early differentiation of th cells. an important goal is to identify master regulators of the lineage commitment processes. recently, it was found out that most of the direct targets of stat an important regulator of th differentiation, were up-regulated in th cells. here we were interested in identifying tfs whose binding sites are enriched in the promoter regions of the genes which are differentially regulated in th conditions, both among the up-regulated and down-regulated genes. instead of looking at individual tf binding predictions that are prone to contain false positives, we used the fisher’s exact test to search for enrichment of binding sites, in comparison to randomly selected gene set. the same analysis was carried out separately for all the differentially regulated gene sets and by taking into account the direction of regulation. using a p-value cut-off of for tf binding, we identified three hits from the enrichment analysis among th specific up-regulated genes and three among the th specific down-regulated genes. the results are depicted in figure  the different enriched irf family motifs were combined and their targets were pooled. in accordance to our previously published results, the strongest hit within the th up-regulated genes was stat, followed by nkx3a, and cdp. nkx3a is a member of the nkx family of homeobox genes that is expressed in prostate epithelium and functions as a potential prostate tumor suppressor. recently, in a study focusing on jurkat cells, a gata binding site on the promoter of nkx gene was identified. furthermore, in mouse increased expression of nkx3a was observed to be regulated by il- independently of stat. cdp is highly conserved homeodomain transcription factor involved in many cellular processes, including differentiation, development and proliferation. interestingly, cdp has been identified as a repressive regulator of cd8α silencer region and tcrβ enhancer region and plays a role in promoting repressive chromatin modifications via association with histone deacetylase and histone methyltransferase. it is important to notice that the binding motif analysis does not give information about the possible direction of regulation, e.g., it is an open question whether cdp might up-regulate th specific genes or down-regulate the genes in th and th lineages. the three tf hits having enriched predicted binding sites among the th down-regulated genes were the interferon regulatory factor family of tfs, ifn-stimulated genes factor and stat. irf family consists of irf to irf and has been shown to be essential in expression of type i interferon genes, ifn-stimulated genes and other pro-inflammatory response related cytokines. these genes are maintained down-regulated during th proliferation and therefore, the results are in line with the th effector cells characteristics. moreover, ifnγ-induced expression of irf and irf has been shown to directly down-regulate il- production by repressing il- promoter sites. opposing to other irf family proteins, irf has been shown to directly activate il- promoter and il- regulatory elements and be essential in th cell differentiation by influencing the expression of gfi a transcriptional repressor in th cells. however, the analysis relying on known tf binding specificities will not allow segregation of individual members of the irf family. further, an essential regulator of most isgs is isgf that is composed of stat stat and irf complex and works in conjunction with irfs. identification of stat as a regulator among the th down-regulated genes is well in line with our previously published results, although its effect was observed to be less profound within th down-regulated genes than among th up-regulated target genes. comparison analysis of the predicted stat target genes and th up-regulated and down-regulated genes gave and overlapping genes, respectively. the full lists of overlapping genes are in additional file 3: table s we further analyzed the correlation between predicted stat target promoters and experimentally observed promoter associated binding sites, and observed significant correlation between the target sites. the full list of predicted stat target genes and promoter associated stat binding sites identified by chip-seq as well as the overlapping genes are listed in the additional file 3: table s the overlapping binding sites included promoters for c14orf cish, hmmr, ino mgat nudcd socs spint and znf genes. discussion identification of the key t helper cell regulators provides possible targets for modulation of immune response. to reveal t cell subset specific genes and their often subtle differences in expression, we developed a novel computational method, ligap. traditional ways of identifying differentially expressed genes, such as the t-test, are problematic in studying time-series data since there is a need to carry out hypothesis tests on individual time points. on the other hand, commonly used statistical tests for whole time-course, including e.g. f-test, do not account for the inherent correlation between measurement time points. ligap overcomes many problems that have previously prevented quantitative comparisons of multiple differentiation profiles, with or without replicates. among several beneficial features, ligap models correlation between time points and can cope with non-stationarities and non-uniform measurement grid. other methods, such as edge, uses splines to estimate smooth time-course profiles but does not quantify the differential expression for all lineage comparisons. tanova uses standard regression framework and lacks explicit correlation structure between time points. our study highlights the validity of the method by identifying known and novel differentially regulated genes and their kinetic differences during t helper cell differentiation. in addition, the non-parametric computational analysis automatically provides informative illustrations of time-course profiles together with associated uncertainty. ligap calculated th specific gene set contains only genes and th specific genes compared to genes that are specific to th conditions. activation of thp cells through tcr and cd results in induction of ifnγ, which in turn leads to activation of th signature genes. addition of il- however, results in enhanced induction of these genes and th programming. consistent with our previous results genes differentially regulated in response to th programming are much more limited than those detected in response to initiation of th response. most of the th specific genes encode well-known th signature molecules. however, also genes new in this context were discovered. interestingly, we identified rorc as one of the th specific genes. up-regulation of rorc in th cells and existence of th17/th cells, however, remain conflicting as the master regulator of th differentiation, t-bet, is known to inhibit transcription of rorc through runx, and expression of il12rβ is down-regulated by il-. it has been suggested that the high concentration of tgfβ required for in vitro th polarization would inhibit ifnγ production, hence, it remains an open question whether some conditions would drive the differentiation of il- and ifnγ producing cells from same naïve precursor t cell. notably, ex vivo th cells could be induced to develop further into th17/th cells by the combined actions of ifnγ and il- and such conditions resulted in permissive chromatin remodeling at the il12rb locus and loss of repressive histone modification at the tbx locus. as an example of previously uncharacterized differentially regulated genes, we validated the expression of th2-associated phosphatases dusp and ppp1r14a on protein level. ppp1r14a was shown in human pancreatic and melanoma tumor cell lines to positively regulate ras/mapk signaling, which are also involved in il- induced signaling cascades. in t cells, the erks are activated though tcr stimulation and a tcr-mediated activation of ras/mapk signaling is required in differentiating murine th but not in th cells. furthermore, the ras/mapk cascade was shown to enhance the stability of gata protein as well as stat independent cd and cd induced initial il production. dusp on other hand is known to negatively regulate members of the mitogen-activated protein kinase superfamily associated with cellular proliferation and differentiation. more specifically, dusp expression was shown to be induced by erk1/ signaling in differentiating mouse embryonic cell line and in human retinal pigment epithelial cells and it was hypothesized that dusp is an essential part of a negative feedback loop of erk1/ signaling. however, the t cell associated functions of both ppp1r14a and dusp are completely unknown. therefore, their significance in the signaling cascades of differentiating th cells remains a highly interesting area of future research. spint was recently identified as a direct stat target in differentiating human th cells and in this study we are the first to show that spint is upregulated in th cells at protein level as compared to other th cell subsets. we found spint to be specifically expressed on th cell surface as well as secreted into the culture medium, suggesting presence of a multiple transcripts of which some may lack the anchoring transmembrane domain. human spint is a physiological inhibitor of matrix cleaving proteases and decreased expression of spint has been linked to progression of several cancers. up-regulated expression of extracellular proteases is crucial for pro-cancerous pathways as this enables efficient remodeling of the extracellular matrix as well as cleavage and activation of growth factors and their receptors. interestingly, a truncated and secreted spint may act as an inhibitor for the activator of hepatocyte growth factor and hgf is prominently expressed in lung tissue and is linked to reduced expression of th cytokines and tgfβ, reduction of allergic airway inflammation, airway hyperresponsiveness and remodeling as well as reduced recruitment of eosinophils to the site of allergic inflammation in vivo. this suggests that spint might enhance th response in allergic airway inflammation by inhibiting hgf signaling. the ligap method elegantly identified the reciprocally regulated genes within the th th and th conditions. essentially, the list included genes encoding the hallmark th specific transcription factor t-bet and cytokine ifnγ as well as the transmembrane receptor for il- this list also included few cytoskeleton associated proteins, such as dystrophin, and palladin, of which there is no current knowledge for their function in differentiating t helper cells. the observation suggests differences in cellular structures or putatively in the interaction of apc with the th cell subsets as rearrangement of the cytoskeleton in t cells plays an important role in the organization of the immunological synapse and th and th cells are known to form morphologically distinct iss. in addition to map3k molecules that participate in phosphorylation signaling cascades e.g. p2ry lpar ppp1r14a, and ptpro suggest their potential role for initiation or regulation of differentiation cascades. importantly, the results presented here enable opportunities for further data mining and follow-up studies addressing the functions and importance of the novel th subset specific genes. the identification of stat as the most significant tf regulating th specific enhancement of transcription by the tf binding analysis is well in line with our previous stat chip results. furthermore, the analysis between the predicted stat target gene promoters and experimentally observed promoter associated binding sites showed statistically significant correlation. interestingly, the overlapping stat targets included ino which has been identifies as a part of a chromatin remodeling complex and may hence, be involved in th specific epigenetical regulation of th cell differentiation. stat specific regulation of mannosyl -glycoprotein beta- -n-acetylglucosaminyltransferase, a n-glycan-processing enzyme, may on one hand be involved in modifying the th cell specific surface glycoprotein structures. the overlapping target sites included also the promoter for spint the number of predicted stat binding sites, however, was much larger than the experimentally observed binding sites, which may reflect the typically observed high false positive rate of computational binding predictions and the cell type specific state of chromatin as well as other competing factors affecting binding in vitro. the data created here also further suggests novel control mechanisms involving gata regulated nkx3a as well as chromatin modification associated cdp. only less than 10% of the th down-regulated genes were reported to be direct targets of stat by elo et al., suggesting other major regulatory mechanisms play role among the il- induced down-regulated genes. we found enrichment of irf family and isgf binding motifs in promoter regions of genes that are repressed in th polarizing conditions, indicating that these tfs may play a significant role in the suppressing undesired gene expression in differentiating th cells. indeed, several irf family members have been identified as differentially expressed during th cell differentiation and necessary for both th and th polarization. as the irf family proteins, excluding irf share the same binding specificity model in transfac, the individual regulatory role for these factors is, however, difficult to postulate based on in silico tf binding site analysis.  the proposed ligap method can quantify a well-defined probabilistic specificity score for each gene and for each condition promoting a certain lineage commitment. in addition to grouping and ranking genes based on their dynamics, ligap summarizes all time-course measurements, together with the associated uncertainty, in an illustrative summary plot for visualization and manual assessment purposes. while here we have demonstrated the utility of ligap in analysis of gene expression dynamics, the ligap method is widely applicable to many types of datasets including quantitative time-course experiments and generalizes to any number of conditions. 
combining differentially regulated transcriptional profiles with transcription factor binding site and pathway information, we identified previously known and new putative transcriptional mechanisms involved in th cell subset differentiation. lineage commitmentnon-parametric analysisth cell differentiationtime-course transcriptomicstranscription factor binding we studied transcriptional dynamics of human umbilical cord blood t helper cells cultured in absence and presence of cytokines promoting th or th differentiation. a proper balance between different t helper cell subsets is necessary for normal functioning of the adaptive immune system. it summarizes all the time-course measurements together with the associated uncertainty for visualization and manual assessment purposes. revealing key genes and pathways driving the differentiation to distinct th cell lineages provides important insight into underlying molecular mechanisms and new opportunities for modulating the immune response. the ligap method is widely applicable to quantify differential time-course dynamics of many types of datasets and generalizes to any number of conditions. to identify genes that exhibit distinct lineage commitment dynamics and are specific for initiating differentiation to different th cell subsets, we developed a novel computational methodology allowing integrative analysis and visualization of multiple lineages over whole time-course profiles. applying ligap to time-course data from multiple th cell lineages, we identified and experimentally validated several differentially regulated th cell subset specific genes as well as reciprocally regulated genes. here we identified novel human th subset specific transcripts as well as regulatory mechanisms important for the initiation of the th cell subset differentiation.
 a proper balance between different t helper cell subsets is necessary for normal functioning of the adaptive immune system. revealing key genes and pathways driving the differentiation to distinct th cell lineages provides important insight into underlying molecular mechanisms and new opportunities for modulating the immune response. previous computational methods to quantify and visualize kinetic differential expression data of three or more lineages to identify reciprocally regulated genes have relied on clustering approaches and regression methods which have time as a factor, but have lacked methods which explicitly model temporal behavior.  we studied transcriptional dynamics of human umbilical cord blood t helper cells cultured in absence and presence of cytokines promoting th or th differentiation. to identify genes that exhibit distinct lineage commitment dynamics and are specific for initiating differentiation to different th cell subsets, we developed a novel computational methodology allowing integrative analysis and visualization of multiple lineages over whole time-course profiles. applying ligap to time-course data from multiple th cell lineages, we identified and experimentally validated several differentially regulated th cell subset specific genes as well as reciprocally regulated genes. combining differentially regulated transcriptional profiles with transcription factor binding site and pathway information, we identified previously known and new putative transcriptional mechanisms involved in th cell subset differentiation. all differentially regulated genes among the lineages together with an implementation of ligap are provided as an open-source resource.  the ligap method is widely applicable to quantify differential time-course dynamics of many types of datasets and generalizes to any number of conditions. it summarizes all the time-course measurements together with the associated uncertainty for visualization and manual assessment purposes. here we identified novel human th subset specific transcripts as well as regulatory mechanisms important for the initiation of the th cell subset differentiation. lineage commitmentnon-parametric analysisth cell differentiationtime-course transcriptomicstranscription factor binding 