Porcupine O-Acyltransferase (PORCN): 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₃ (= 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 PORCN 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. Porcupine (PORCN)

The Drosophila segment polarity gene product Porcupine (Porc) was first identified as being necessary for processing Wingless (Wg), a Drosophila Wnt (Wnt) family member. Tanaka et al. [4] identified Mouse (Mporc) and Xenopus (Xporc) homologs of porc and found that they encode endoplasmic reticulum (ER) proteins with multiple transmembrane domains. Further, Mporc mRNA was differentially expressed during embryogenesis and in various adult tissues, demonstrating that the alternative splicing is regulated to synthesize the specific types of Mporc. In transfected mammalian cells, they found all types of Mporc affected the processing of mouse WNT1, WNT3A, WNT4, WNT6, and WNT7B but not WNT5A. Lastly, they also found that all types of Mporc co-immunoprecipitated with various WNT proteins. Their results suggested that Mporc may function as a chaperone-like molecule for WNT.

Caricasole et al. [5] report that the human Porcupine locus (MG61/PORC) spans 15 exons over approximately 12 kb of genomic sequence on Xp11.23. Like its mouse and Xenopus homologues, MG61/PORC encodes four protein isoforms (A–D) generated through alternative splicing and expressed in a tissue-specific fashion. They present evidence indicating that MG61/PORC can influence the activity of a human WNT7A expression construct in a T-cell factor-responsive reporter assay.

Liu et al. [6] indicate that post-translational modification of WNTs includes lipid modification and glycosylation. The former is performed by PORCN. PORCN is a membrane-bound O-acyltransferase located in the endoplasmic reticulum and can add palmitoleate groups to WNT proteins that is necessary for WNT ligand secretion, and it is a member of the membrane-bound O-acyltransferases (MBOATs). Lipid modification is necessary for Wnt activity, and the opposite is true for glycosylation as observed by Willert et al. [7]. Liu et al. [8] developed a screen for small molecules that blocked WNT secretion and discovered LGK974, a potent and specific small-molecule PORCN inhibitor. They show that LGK974 inhibits WNT signaling, including reduction of the WNT-dependent LRP6 phosphorylation and the expression of WNT target genes, like as AXIN2. The inhibitor is effective in multiple tumor models at well-tolerated doses. Together, their findings provide a strategy and and a tool for targeting WNT-driven cancers through the inhibition of PORCN. Further down the line, Madan et al. [9] developed a novel potent, orally available PORCN inhibitor, ETC-1922159 that blocked the secretion and activity of all WNTs. ETC-1922159 is remarkably effective in treating RSPO-translocation bearing colorectal cancer (CRC) patient-derived xenografts. This is the first example of effective targeted therapy for this subset of CRC. By this demonstration they show that inhibition of WNT signaling by PORCN inhibition holds promise as differentiation therapy in genetically defined human cancers.

Sutton [10], conver a range of PORCN related developmental disorders. Till now, most research work has focused on the role of PORCN along with WNTs. In this research work, I present 3rd order combinations of PORCN with other genes, apart from the WNTs, 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 [11]. 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 [11].

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₃ 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 PORCN-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 [12]) and Sobol indicies (with 2002 implementation in Saltelli [13] and martinez implementation in Martinez [14] and Baudin et al. [15]).

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

A total of 2415, 3rd order combinations involving PORCN were obtained from a full set of ⁷¹C₃ = 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 CXXC4-PORCN-X Combinations

Ekici et al. [16] identified a homozygous splice site mutation in the CCDC88C gene as a novel cause of a complex hydrocephalic brain malformation, via positional cloning in a consanguineous family with autosomal recessive hydrocephalus. CCDC88C encodes DAPLE (HkRP2), a Hook-related protein with a binding domain for the central WNT signalling pathway protein DVL and DVL is inhibited by CXXC4 and CCDC88C, mediates Wnt signalling via inhibition of GSK3. Looking at the tables above, one finds the following combinations for CXXC4 along with PORCN, to be prominent at 3rd order level - CXXC4-PORCN-SENP2, CXXC4-PORCN-WNT2B, CXXC4-PORCN-TCF7, CXXC4-PORCN-WNT4, CXXC4-PORCN-SFRP4, CXXC4-PORCN-FBXW4, CXXC4-PORCN-PPP2R1A, CXXC4-PORCN-PPP2CA and CXXC4-PORCN-WNT5A. 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 FOSL1-PORCN-X Combinations

The ability of the right ventricle (RV) to adapt to an increased pressure afterload determines survival in patients with pulmonary arterial hypertension. The WNT pathway plays an important role in the development of the RV and may also be implicated in adult cardiac remodeling. Nayakanti et al. [17] show that • WNT/$\beta$-catenin signaling molecules are upregulated in RV of patients with pulmonary arterial hypertension and animal models of RV overload (pulmonary artery banding-induced and monocrotaline rat models); • Activation of WNT/$\beta$-catenin signaling led to RV remodeling via transcriptional activation of FOSL1 and FOSL2; • Further, the mRNA expression profiling of WNT signaling molecules and immunofluorescence staining of aCTNNB1 (active $\beta$-catenin) and PORCN in RV myocardial tissues from human donors (control) and patients with idiopathic PAH (IPAH) demonstrated a significant upregulation of WNT/$\beta$-catenin signaling and • finally, pharmacological inhibition of WNT signaling using inhibitor of PORCN, LGK-974 attenuated fibrosis and cardiac hypertrophy leading to improvement in RV function in both, pulmonary artery banding- and monocrotaline-induced RV overload and the blockade using LGK-974 in WNT3A-stimulated cardiac fibroblasts resulted in significant downregulation of mRNA and protein expression of Wnt/$\beta$-catenin signaling and its target genes ($\beta$-catenin, PORCN, MYC), transcription factors (FOSL1, FOSL2), proliferation marker (CCND1), and fibroblast-to-myofibroblasts transdifferentiation markers (POSTN, CTGF). Looking at the tables above, one finds the following combinations for FOSL1 along with PORCN, to be prominent at 3rd order level - FOSL1-PORCN-SFRP4, FOSL1-PORCN-SENP2, FOSL1-PORCN-WNT4, FOSL1-PORCN-WNT2B, FOSL1-PORCN-TLE2, FOSL1-PORCN-RHOU, FOSL1-PORCN-SFRP1 and FOSL1-PORCN-FBXW4. 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 PITX2-PORCN-X Combinations

PITX2 is activated by canonical WNT signaling, and can also regulate expression of several WNT ligands. Kioussi et al. [18] report that the bicoid-related transcription factor 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. Basu and Roy [19] found that PITX2 interacts and regulates WNT2/5A/9A/6/2B genes of the canonical, noncanonical, or other pathways in the human ovarian cancer cell SKOV-3. Looking at the tables above, one finds the following combinations for PITX2 along with PORCN, to be prominent at 3rd order level - PITX2-PORCN-WNT2B, PITX2-PORCN-TLE1, PITX2-PORCN-SFRP4, PITX2-PORCN-SENP2, PITX2-PORCN-WNT5A, PITX2-PORCN-FBXW4, PITX2-PORCN-TCF7L1, PITX2-PORCN-PPP2R1A, PITX2-PORCN-RHOU, and PITX2-PORCN-SFRP1. 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 CSNK1D-PORCN-X Combinations

Zhu et al. [20] show that the circadian gene CSNK1D was significantly upregulated in hepatocellular carcinoma and enhanced malignant behaviors via WNT signaling pathway by stabilizing DVL3. Looking at the tables above, one finds the following combinations for CSNK1D along with PORCN, to be prominent at 3rd order level - CSNK1D-PORCN-SENP2, CSNK1D-PORCN-WNT2B, CSNK1D-PORCN-TLE2, CSNK1D-PORCN-TCF7, CSNK1D-PORCN-SFRP1 and CSNK1D-PORCN-SLC9A3R1. 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 RHOU-PORCN-X Combinations

Having demonstrated that the expression levels of WNT1, PORCN and RSPO2 are downregulated in human Alzheimer’s disease (AD) brains, Macyczko et al. [21] examined the potential associations among the expression levels of these genes. It was found that the PORCN mRNA levels were positively correlated not only with the WNT1 mRNA levels, but also with the RSPO2 mRNA levels in human AD brains. Tao et al. [22] demonstrate via three lines of evidence evidence demonstrate that WNT1 responsive Cdc42 homolog (WRCH1 or RHOU) induction is correlated or associated with the expression of WNT1, i.e • After retroviral infection of cells, WRCH1/RHOU is selectively expressed in WNT1 expressing cells but not in control cells infected with an empty viral vector; • WNT1 can induce WRCH1 mRNA levels in cell coculture; and (3) WRCH1 is expressed in WNT1 induced mouse mammary tumors but not in wild-type mouse mammary glands. Also, Schiavone et al. [23] show that WNT1 could induce RHOU promoter transcription at similar levels in the presence or absence of signal transducer and activator of transcription 3 (STAT3). Looking at the tables above, one finds the following combinations for RHOU along with PORCN, to be prominent at 3rd order level - FRAT1-PORCN-RHOU, BCL9-PORCN-RHOU, CSNK1G1-PORCN-RHOU, FZD1-PORCN-RHOU, PITX2-PORCN-RHOU, FOSL1-PORCN-RHOU and LRP5-PORCN-RHOU. 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 FBXW-PORCN-X Combinations

In HEK-293 cells, Xiao et al. [24] identified FBXW7 as the direct downstream gene of miR-367, which consequently released the LIN28 dependent inhibition of suppressive LET7. Through their informatics analysis, miR-367 was predicated to function through WNT signaling, and decreased LET7 played an important role to maintain TCF-4/WNT pathway activity. The reintroduction of FBXW7 abolished the oncogenic stimulation of miR-367 on TCF-4 activity, with WNT signaling factors depression. In conclusion, they demonstrated that a FBXW7-centric signaling that miR-367 function through, and the downstream WNT activation was responsible for miR-367 oncogenic role of self-renewal, in a LIN28B dependent manner. ETC-1922159 (Madan et al. [9]) potently suppressed the FBXW7–wild-type HPAF-II tumor growth and Zhong and Virshup [25] show that in the FBXW7–wild-type HPAF-II xenografts, PORCN inhibition led to decreased phosphorylation of Rb protein and reduced total Rb abundance. Looking at the tables above, one finds the following combinations for members of FBXW family along with PORCN, to be prominent at 3rd order level - FBXW2-PORCN-SFRP4, FZD7-PORCN-FBXW4, FBXW2-PORCN-WNT4, FBXW11-PORCN-WNT4, FBXW11-PORCN-SFRP4, CCND3-PORCN-FBXW4, CXXC4-PORCN-FBXW4, PITX2-PORCN-FBXW4, FBXW11-PORCN-WNT2B, FZD1-PORCN-FBXW4 and FOSL1-PORCN-FBXW4. All these combinations indicate the existence of a possible synergy when they take a higher rank in the list of combinations.

6.3.7. Examining the Behaviour of KREMEN1-PORCN-X Combinations

Multiciliated cells contain hundreds of cilia whose directional movement powers the mucociliary clearance of the airways, a vital host defense mechanism and their specification requires canonical WNT signaling, which then must be turned off. Further, ciliogenesis and polarized ciliary orientation are regulated by noncanonical WNT/planar cell polarity (WNT/PCP) signaling. Cooney et al. [26] show that DKK3 and WNT4, act together to facilitate a canonical to noncanonical WNT signaling switch during multiciliated cell formation and demonstrate that DKK3 is contemporaneously expressed with WNT4 by the same population of basal cells, and suggest that it binds to the multiciliated cell surface via KREMEN1. Looking at the tables above, one finds the following combinations for KREMEN1 along with PORCN, to be prominent at 3rd order level - KREMEN1-PORCN-TCF7L1, KREMEN1-PORCN-WNT4, KREMEN1-PORCN-SENP2, KREMEN1-PORCN-WNT3A and KREMEN1-PORCN-WNT2B. 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 PORCN in WNT3A stimulated HEK 293 cells. Based on the extablished 2nd order combinations of the PORCN, 3rd order combinations emerge using the machine learning based search engine. These 3rd order combinations might be of interest for further wet lab investigations.