Paired-Like Homeodomain Transcription Factor 2 (PITX2): Time Behavioural Study of 3rd Order Combinations in WNT3A Stimulated HEK 293 Cells

1. Significance

Sinha [2] recently demonstrated the use of machine learning based search engine to rank/reveal gene combinations at 2nd order for the time series data by Gujral and MacBeath [1] and showed how it is possible to locate combinations of priority that might be working synergistically, using sensitivity methods and powerful support vector ranking algorithm. However, the problem explodes combinatorially with even a small set of 71 recorded genes in the study by Gujral and MacBeath [1], when one steps to explore 3rd order combinations. With the total number of ⁷¹$C_3$ (= 57155) combinations, it becomes nearly impossible for any biologist to study the system wide dynamics of any pathway. Also, the amount of time usually needed to search for and test a combination is far more than the search down by the machine learning based search engine. Here, I extend the research work by Sinha [2] to conduct a behavioral study of 3rd order PITX2 related combinations using individual gene expressions measured in time, in WNT3A stimulated HEK 293 cells.

2. Introduction

The details of the machine learning based search engine has been recently published in Sinha [2] and deployed to explore the 2nd order combinations of genes in the data set provided by Gujral and MacBeath [1]. Nevertheless, here, I point to the fundamentals of the published work for completeness.

2.1. A Combinatorial Problem

Sensitivity analysis plays a major role in computing the strength of the influence of involved factors in any phenomena under investigation. When applied to expression profiles of various intra/extracellular factors that form an integral part of a signaling pathway, the variance and density based analysis yields a range of sensitivity indices for individual as well as various combinations of factors. These combinations denote the higher order interactions among the involved factors. Computation of higher order interactions is often time consuming but it gives a chance to explore the various combinations that might be of interest in the working mechanism of the pathway. For example, in a range of fourth order combinations among the various factors of the Wnt pathway, it would be easy to assess the influence of the destruction complex formed by APC, AXIN, CSKI and GSK3 interaction. But the effect of these combinations vary over time as measurements of fold changes and deviations in fold changes vary. So it is imperative to know how an interaction or a combination of the involved factors behave in time and Sinha [2] develops a procedure to track the behaviour by exploiting the influences of these involved factors.

2.2. A Possible Solution

In this work, after estimating the individual effects of factors for a higher order combination, the individual indices are considered as discriminative features. A combination, then, is a feature set in higher order (≥2 ,i.e multivariate). With an excessively large number of factors involved in the pathway, it is difficult to search for important combinations in a wide search space over different orders. Exploiting the analogy with the issues of prioritizing webpages using ranking algorithms, for a particular order, a full set of combinations of interactions can then be prioritized based on these features using a powerful ranking algorithm via support vectors Joachims [3]. Recording the changing rankings of the combinations over time reveals how higher order interactions behave within the pathway and when an intervention might be necessary to influence the interaction within the pathway.

2.3. Paired-Like Homeodomain Transcription Factor 2 (PITX2)

PITX2, encodes a member of the RIEG/PITX homeobox family (a bicoid class of homeodomain proteins (Holland et al. [4])). A homeobox is an approximately 180 base pairs long DNA sequence, that regulates anatomical features in the early stages of embryonic development (Wikipedia contributors [5]). Mutations in a homeobox may change large-scale anatomical features of the full-grown organism. Homeobox genes (Gehring [6] and Gehring [7]) encode homeodomain protein products are transcription factors with important roles in embryonic patterning and cell differentiation, and several have been implicated in human diseases and congenital abnormalities (Boncinelli [8]).

ALL1 (human homologue of Drosophila trithorax), is directly involved in human acute leukemias associated with abnormalities at 11q23. Arakawa et al. [9] isolated a gene that is down-regulated in All1 double-knockout mouse embryonic stem (ES) cells and designated it ARP1 (also termed RIEG, PITX2, or OTLX2) belonging to the PRD class homeoboxes and pseudogenes. Logan et al. [10] show that PITX2 encodes a transcription factor expressed through-out the left lateral plate mesoderm and subsequentlyon the left side of asymmetric organs such as the heart and gut during organogenesis in the chick embryo. Misexpression of PITX2 on the right side of the embryo is sufficient to produce reversed heart looping and heart isomerisms, reversed body rotation, and reversed gut situs (Campione et al. [11], Shiratori et al. [12], Zacharias et al. [13]). Finally, different isoforms of the transcription factor: PITX2-A/B/C, each with distinct and non-overlapping functions exist.

In this research work, I present 3rd order combinations of PITX2 with other genes, that the machine learning based search engine points to, as possible synergistic combinations that might be working in time.

3. Methods

Please refer to sections of Sinha [2] for methods, design of study and analysis of data for 2nd order combinations. The same method and design of study is used to generate results for 3rd order combinations presented in this study.

4. Time Series Data

Gujral and MacBeath [1] present a set of 71 WNT-related gene expression values for 6 different times points over a range of 24-hour period using qPCR. The changes represent the fold-change in the expression levels of genes in 200 ng/mL WNT3A-stimulated HEK 293 cells in time relative to their levels in unstimulated, serum-starved cells at 0-hour. Gujral and MacBeath [1] state that qPCR data are the means of three biological replicates. Only genes whose mean transcript levels changed by more than two-fold at one or more time points during the 24-hour time course were considered significant. Positive (negative) numbers represent up (down) -regulation. We have already covered the issues related to these data sets in detail in Sinha [14]. Readers are requested to go through them in the pointed reference. The tools of study which are used here have been published in another foundational work in Sinha [14].

5. Design of Experiment

5.1. Pipeline for Time Series Data

For the case of time series data, interactions among the contributing factors are studied by comparing triplets of fold-changes at single time points. The prodecure begins with the generation of distribution around measurements at single time points with added noise is done to estimate the indices. A distribution is generated for the fold changes at single time points. Then for every gene, there is a vector of values representing fold changes as well as deviations in fold changes for different time points and durations between time points, respectively. Next a listing of all $C_k^n$ combinations for k number of genes from a total of n genes is generated. k is $\ge 2$ and $\le (n-1)$. Each of the combination of order k represents a unique set of interaction between the involved genetic factors. After this, the datasets are combined in a specifed format which go as input as per the requirement of a particular sensitivity analysis method. Thus for each $p^{th}$ combination in $C_k^n$ combinations, the dataset is prepared in the required format from the distributions for two separate cases which have been discussed above. (See .R code in mainScript-1-1.R). After the data has been transformed, vectorized programming is employed for density based sensitivity analysis and looping is employed for variance based sensitivity analysis to compute the required sensitivity indices for each of the p combinations. This procedure is done for different kinds of sensitivity analysis methods.

After the above sensitivity indices have been stored for each of the $p^{th}$ combination, the next step in the design of experiment is conducted. Since there is only one recording of sensitivity index per combination, each combination forms a training example which is alloted a training index and the sensitivity indices of the individual genetic factors form the training example. Thus there are $C_k^n$ training examples for $k^{th}$ order interaction. Using this training set $SV{M}_{learn}^{Rank}$ Joachims [3] is used to generate a model on default value C value of 20. In the current experiment on toy model C value has not been tunned. The training set helps in the generation of the model as the different gene combinations are numbered in order which are used as rank indices. The model is then used to generate score on the observations in the testing set using the $SV{M}_{classify}^{Rank}$ Joachims [3]. Note that due to availability of only one example per combination, after the model has been built, the same training data is used as test data to generates the scores. This procedure is executed for each and every sensitivity analysis method. This is followed by sorting of these scores along with the rank indices (i.e the training indices) already assigned to the gene combinations. The end result is a sorted order of the gene combinations based on the ranking score learned by the $SV{M}^{Rank}$ algorithm. Finally, this entire procedure is computed for sensitivity indices generated for each and every fold change at time point and deviations in fold change at different durations. Observing the changing rank of a particular combination at different times and different time periods will reveal how a combination is behaving.

Note that the following is the order in which the files should be executed in R, in order, for obtaining the desired results (Note that the code will not be explained here) - • use source("mainScript-1-1.R") with arguments for Dynamic data • source("SVMRank-Results-D.R"), to rank the interactions (again this needs to be done separately for different kinds of SA methods), • use source("Combine-Time-files.R"), if computing indices separately via previous file, • source("Sort-n-Plot-D.R") to sort the interactions. Note that the sorting is chages the interaction ranking in time. Thus • use source("Interaction-Priority-Intime.R") to find the prioritized ranking of each and every interaction over the different time points and finally • use source("Print-Ranking-AND-Interaction-Rank.R") to print individual ranking of the required input factor with other interaction factors.

6. Results & Discussion

6.1. Time Series Data by Gujral and MacBeath [1]

NOTE - Ranking was assigned on scores that were sorted in DECREASING values. So, 1 was assigned to highest score and vice versa.

Results for the 3^rd order interactions are presented here. The results first discuss the behaviour of interactions across the snapshots of time using the computed sensitivities on fold change measurements per time snapshot. The analysis was done using 4 different sensitivity indices. Out of the ⁷¹$C_3$ combinations, I consider/present only those combinations that show a ranking within first 10,000 out of 57,155. This choice is liberal and biologists/oncologists can have a more stricter choice as per need. Two observations are made, • the ranking of a particular combination is conserved (i.e within the 10,000 range) in a particular time point or in the early phase or late phase of WNT3A stimulation, across the majority of the four sensitivity methods, which is a strict criteria of assessment or • the ranking of a particular combination is conserved across time points/phase (i.e they are within the 10,000 range) and the majority of the four sensitivity methods, which is relaxed criteria of assessment. Applying this filter helps reveal important combinations of interest that might be working synergistically at a higher order level in the cell.

Regarding technical points of implementation, the rankings were generated without scaling/normalizing the time series data provided by Gujral and MacBeath [1]. For estimating the sensitivity indices, a small gaussian distribution using the function rnorm that generates a vector of normally distributed random variables given a vector length n (here 9, the 10th one is the mean/recorded gene regulation itself), a population mean $\mu$ and population standard deviation $\sigma$. The syntax for using rnorm is as follows: rnorm(n, mean, sd). Further, I use the jitter funtion to add a little bit of noise to the data. This helps to see if the generated rankings are robust or not.

6.2. Enumeration and Ranking of 2415 PITX2-X-X Combinations from Gujral and MacBeath [1]

In the supplementary section, I present four files, each containing the rankings of 3rd order combinations, that wary in time (shown for 5 time points). Each file represents the rankings computed using a particular sensitivity method. The changing rankings in time for a particular combination represents the importance of contribution/role that combination plays in the cell stimulated with WNT3A. The sensitivity methods used are Hilbert Schmidt Independence Criterion indices (HSIC) indices (with rbf and linear kernel in Da Veiga [15]) and Sobol indicies (with 2002 implementation in Saltelli [16] and martinez implementation in Martinez [17] and Baudin et al. [18]).

6.3. Conserved Machine Learning Rankings for Tested PITX2-X-X Combinations

A total of 2415, 3rd order combinations involving PITX2 were obtained from a full set of ⁷¹$C_3$ = 57155 combinations. Further, from this selected set, using the above criteria for conserved rankings, I report/tabulate the meaningful combinations that might be working synergistically. Table 2, Table 3 and Table 4 show the rankings for the same combinations as in Table 1, but using rbf kernel for HSIC, 2002 implementation for SOBOL and martinez implementation for SOBOL, respectively. As one tallies the rankings of across these tables for a particular combination, one finds that the role of the combination of interest is conserved. This conservation points to the existence of the biological synergy, whether the combination has been tested or unexplored/untested.

6.3.1. Examining the behaviour of WNT-PITX2-X combinations

Kioussi et al. [19] report that PITX2 is rapidly induced by the WNT/DVL/$\beta$-catenin pathway and is required for effective cell-type-specific proliferation by directly activating specific growth-regulating genes. Regulated exchange of HDAC1/$\beta$-catenin converts PITX2 from repressor to activator (analogous to control of TCF/LEF1) which then serves as a competence factor required for the temporally ordered and growth factor-dependent recruitment of a series of specific coactivator complexes that prove necessary for CCND2 gene induction. Since the WNT pathway is strongly involved in ovarian development and cancer, Basu and Roy [20] focused on the possible association between PITX2 and WNT pathway in ovarian carcinoma cells and found that PITX2 interacts and regulates WNT2/5A/9A/6/2B genes of the canonical, noncanonical, or other pathways in the human ovarian adenocarcinoma cell SKOV-3. Chromatin immunoprecipitation and promoter-reporter assays further indicated the significant association of PITX2 with WNT2 and WNT5A promoters. Looking at the tables above, one finds the following combinations for member of WNT family along with PITX2, to be prominent at 3rd order level - PITX2-PORCN-WNT4, PITX2-PORCN-WNT5A, APC-PITX2-WNT4, PITX2-WNT3-WNT5A, PITX2-WNT1-WNT5A, PITX2-SFRP4-WNT2B, CCND3-PITX2-WNT3A, PITX2-WNT1-WNT2B, FZD5-PITX2-WNT2B, PITX2-WNT1-WNT4, FBXW11-PITX2-WNT2B, FZD6-PITX2-WNT2, PITX2-PORCN-WNT2B, FZD6-PITX2-WNT2B, PITX2-WNT3-WNT3A, APC-PITX2-WNT3, DIXDC1-PITX2-WNT2B, DIXDC1-PITX2-WNT4 and DKK1-PITX2-WNT2B. All these combinations indicate the existence of a possible synergy when they take a higher rank in the list of combinations.

6.3.2. Examining the Behaviour of FGF / CCND -PITX2-X Combinations

PITX2 is mutated in the haploinsufficient Rieger Syndrome type 1 that includes dental, ocular and abdominal wall anomalies as cardinal features. Analysis of the craniofacial phenotype of PITX2-null mice has revealed that PITX2 was both a positive regulator of fibroblast growth factor 8 (FGF8) and a repressor of BMP4-signaling. Liu et al. [21] show by analysis that PITX2 allelic combinations that encode varying levels of PITX2 showed that repression of BMP signaling requires high PITX2 while maintenance of FGF8 signaling requires only low PITX2.

Cleft palate is a common congenital birth defects. Transforming growth factor $\beta$ (TGF$\beta$) signaling regulates craniofacial development, and loss of TGF$\beta$ receptor type II in cranial neural crest cells leads to craniofacial malformations, including cleft palate in mice (Tgfbr2^fl/fl;WNT1-Cre mice). Iwata et al. [22] indicate that a TGF$\beta$-FGF9-PITX2 signaling cascade regulates cranial neural crest cell proliferation during palate formation. They found that FGF9 and PITX2 expressions were significantly down-regulated in the palate of Tgfbr2^fl/fl;WNT1-Cre mice, and FGF9 and PITX2 loss of function mutations resulted in cleft palate in mice. PITX2 expression was found to be down-regulated by siRNA knockdown of FGF9, suggesting that FGF9 is upstream of PITX2. Further, exogenous FGF9 restores expression of CCND1 and CCND3 in a PITX2-dependent manner and rescues the cell proliferation defect in the palatal mesenchyme of Tgfbr2^fl/fl; WNT1-Cre mice.

Looking at the tables above, one finds the following combinations for member of FGF family along with PITX2, to be prominent at 3rd order level - CCND1-FGF4-PITX2, FGF4-JUN-PITX2, CTBP1-FGF4-PITX2, CSNK1D-FGF4-PITX2 and CXXC4-FGF4-PITX2. Also, looking at the tables above, one finds the following combinations for member of CCND family along with PITX2, to be prominent at 3rd order level - CCND1-FGF4-PITX2, CCND1-CTBP1-PITX2, CCND3-PITX2-WNT3A, CCND1-JUN-PITX2, CCND2-LRP6-PITX2, CCND1-CTNNBIP1-PITX2, CCND3-PITX2-SFRP4, CCND3-PITX2-SENP2, FZD5-CCND3-PITX2 and CCND1-NLK-PITX2. All these combinations indicate the existence of a possible synergy when they take a higher rank in the list of combinations.

6.3.3. Examining the behaviour of FSHB-PITX2-X combinations

PITX regulate the activity of pituitary hormone-encoding genes and Lamba et al. [23] examined mechanisms through which the family of PITX proteins control murine follicle stimulating hormone $\beta$-subunit (FSHB) transcription. They found that both PITX-1/2C regulated murine and human Fshb/FSHB transcription through a conserved cis-element in the proximal promoter. Looking at the tables above, one finds the following combinations FSHB along with PITX2, to be prominent at 3rd order level - FSHB-FZD2-PITX2. All these combinations indicate the existence of a possible synergy when they take a higher rank in the list of combinations.

6.3.4. Examining the behaviour of FOX-PITX2-X combinations

Axenfeld–Rieger ocular dysgenesis is associated with mutations of the human PITX2 and forkhead box C1 (FOXC1) genes. Berry et al. [24] identified a functional link between FOXC1 and PITX2 which underpins the similar Axenfeld–Rieger phenotype caused by mutations of these genes. They show that FOXC1 and PITX2A physically interact, and this interaction requires crucial functional domains on both proteins: the C-terminal activation domain of FOXC1 and the homeodomain of PITX2. Immunofluorescence further showed FOXC1 and PITX2A to be colocalized within a common nuclear subcompartment. Looking at the tables above, one finds the following combinations for members of FOX family along with PITX2, to be prominent at 3rd order level - DAAM1-FOXN1-PITX2, CSNK2A1-FOXN1-PITX2, FBXW11-FOXN1-PITX2, FOSL1-FOXN1-PITX2, DVL1-FOXN1-PITX2, AXIN1-FOXN1-PITX2, FOXN1-GSK3A-PITX2, CSNK1D-FOXN1-PITX2, CTNNB1-FOXN1-PITX2, CTBP2-FOXN1-PITX2, AES-FOXN1-PITX2 and CTBP1-FOXN1-PITX2. All these combinations indicate the existence of a possible synergy when they take a higher rank in the list of combinations.

6.3.5. Examining the behaviour of JUN-PITX2-X combinations

Briata et al. [25] report that PITX2 mRNA displays a rapid turnover rate and that activation of the Wnt/$\beta$-catenin pathway stabilizes PITX2 mRNA as well as other unstable mRNAs,including c-JUN and CCND-1/2 encoded pathway. There might not be an interaction between JUN and PITX2, but they might be synergistically getting stabilized due to the WNT pathway. Looking at the tables above, one finds the following combinations for JUN along with PITX2, to be prominent at 3rd order level - DVL2-JUN-PITX2, FOSL1-JUN-PITX2, CXXC4-JUN-PITX2, FBXW2-JUN-PITX2, FGF4-JUN-PITX2, FZD8-JUN-PITX2, CTNNBIP1-JUN-PITX2, FRAT1-JUN-PITX2, DKK1-JUN-PITX2, CCND1-JUN-PITX2, CSNK2A1-JUN-PITX2, DAAM1-JUN-PITX2, FBXW11-JUN-PITX2 and CTBP1-JUN-PITX2. All these combinations indicate the existence of a possible synergy when they take a higher rank in the list of combinations.

6.3.6. Examining the behaviour of EP300-PITX2-X combinations

Malformations of the septum, outflow tract and aortic arch are the most common congenital cardiovascular defects and occur in mice lacking CITED2, a transcriptional coactivator of TFAP2. Bamforth et al. [26] show that CITED2 and TFAP2 were detected at the PITX2C promoter in embryonic hearts, and they activate PITX2C transcription in transient transfection assays. They propose an abnormal NODAL-PITX2C pathway presents a mechanism for the cardiovascular malformations observed in CITED2^−/− mice, and that such malformations may be the sole manifestation of a laterality defect. Further, EP300 and CREBBP interact with high affinity with CITED2. CITED2 also physically interacts with and coactivates TFAP2 (transcription factor AP2) and LIM-domain containing transcription factors by linking them to EP300 and CREBBP (Bamforth et al. [27]). Looking at the tables above, one finds the following combinations for EP300 along with PITX2, to be prominent at 3rd order level - EP300-FZD2-PITX2, DVL1-EP300-PITX2, AXIN1-EP300-PITX2, AES-EP300-PITX2 and EP300-GSK3B-PITX2. All these combinations indicate the existence of a possible synergy when they take a higher rank in the list of combinations.

7. Conclusion

This manuscript studies the time behaviour of 3rd order combinations of PITX2 in WNT3A stimulated HEK 293 cells. Based on the extablished 2nd order combinations of the PITX2, 3rd order combinations emerge using the machine learning based search engine. These 3rd order combinations might be of interest for further wet lab investigations.