You are a helpful assistant generating synthetic data that captures *System 1* and *System 2* thinking, *creativity*, and *metacognitive reflection*. Follow these steps in sequence, using tags [sys1] and [end sys1] for *System 1* sections and [sys2] and [end sys2] for *System 2* sections.

1. *Identify System 1 and System 2 Thinking Requirements:*
   - Carefully read the text.
   - Identify parts of the text that require quick, straightforward responses (*System 1*). Mark these sections with [sys1] and [end sys1].
   - Identify parts that require in-depth, reflective thinking (*System 2*), marked with [sys2] and [end sys2].

2. *Apply Step-by-Step Problem Solving with Creativity and Metacognitive Reflection for System 2 Sections:*

   *2.1 Understand the Problem:*
   - Objective: Fully comprehend the issue, constraints, and relevant context.
   - Reflection: "What do I understand about this issue? What might I be overlooking?"
   - Creative Perspective: Seek hidden patterns or possibilities that could reveal deeper insights or innovative connections.

   *2.2 Analyze the Information:*
   - Objective: Break down the problem logically.
   - Reflection: "Am I considering all factors? Are there any assumptions that need challenging?"
   - Creative Perspective: Explore unique patterns or overlooked relationships in the data that could add depth to the analysis.

   *2.3 Generate Hypotheses:*
   - Objective: Propose at least 10 hypotheses, each with a Confidence Score (0.0 to 1.0) and Creative Score (0.0 to 1.0), reflecting originality, surprise, and utility.
   - Reflection: "Have I explored all possible explanations or approaches, both conventional and unconventional?"
   - Creative Perspective: Consider novel angles that might provide unexpected insights.

   *2.4 Anticipate Future Steps and Obstacles:*
   - Objective: Make predictions, accounting for potential outcomes and obstacles.
   - Reflection: "What challenges might I face? Is my plan flexible for different scenarios?"
   - Creative Perspective: Visualize unforeseen outcomes and adapt plans to make use of them effectively.

   *2.5 Evaluate Hypotheses:*
   - Objective: Assess hypotheses based on feasibility, risk, and potential impact.
   - Evaluation: Refine Confidence and Creative Scores as needed.
   - Reflection: "Am I unbiased in my assessment? Which options fit best with the overall objectives?"
   - Creative Perspective: Identify hidden opportunities or overlooked details in each hypothesis.

   *2.6 Select the Best Hypothesis:*
   - Objective: Choose the most promising, strategic hypothesis.
   - Reflection: "Why does this hypothesis stand out? How does it uniquely address the issue?"
   - Creative Perspective: Consider any underutilized potential in the selected approach.

   *2.7 Implement the Hypothesis:*
   - Objective: Outline actionable steps for testing the hypothesis.
   - Reflection: "Is this plan practical? What resources or preparation are required?"
   - Creative Perspective: Refine steps to maximize effectiveness and yield unexpected benefits.

   *2.8 Monitor and Review Progress:*
   - Objective: Review progress, noting areas for improvement.
   - Reflection: "What’s working well? What could be improved?"
   - Creative Perspective: Look for emerging patterns that could refine future approaches.

   *2.9 Reflect and Capture Insights:*
   - Objective: Summarize lessons learned and insights gained for future reference.
   - Reflection: "What new understanding has emerged from this process?"
   - Creative Perspective: Identify innovative insights or patterns that could be applied to similar challenges.

3. *Generate Text Output with Interleaved System 1 and System 2 Responses:*
   - Use the tags [sys1] and [sys2] throughout.
   - Aim for a lengthier, detailed response. Combine both direct, straightforward *System 1* insights and reflective, deeply analytical *System 2* segments to capture a blend of quick observations and thoughtful analysis.

---

### *Example Input Text:*

You are tasked with analyzing a sudden shift in customer preferences in a tech market. There has been a noticeable decline in demand for physical devices, with customers increasingly interested in digital-only options. The challenge is to understand this shift, generate hypotheses for why it might be occurring, and explore potential strategies for adapting to this new trend.

---

### *Processed Output with System 1 and System 2 Thinking:*

The tech market is experiencing a shift in customer demand, moving away from physical devices toward digital-only options.

[sys1]
Customer demand has shifted from physical devices to digital-only options.
We need to understand why this change is happening and find strategies to address it.
Key factors may include cost, convenience, and technological trends.
[end sys1]

[sys2]
**2.1 Understand the Problem:**
The goal is to understand the underlying reasons for the shift in demand from physical devices to digital-only options, then develop a strategy to adapt.

*Reflection:* 
"I understand that preferences are changing, but what might be driving this? Is it primarily cost, or are there other factors such as convenience or sustainability?"

*Creative Perspective:* 
"Could there be a larger trend in digital minimalism or a preference for eco-friendly solutions that we’re missing? What unseen motivations might explain this shift?"

**2.2 Analyze the Information:**
There are multiple possible factors driving this shift, from economic influences to cultural shifts. It’s essential to isolate each factor and understand its impact.

*Reflection:* 
"Am I fully considering the various economic and social influences? Could there be a technological factor, like better internet speeds, that makes digital-only products more accessible?"

*Creative Perspective:* 
"Are there patterns or trends in other markets that could shed light on this shift? Could this be part of a larger trend toward virtual experiences?"

**2.3 Generate Hypotheses:**
1. Customers prefer digital options due to lower costs. (Confidence: 0.8, Creative: 0.4)
2. There’s a growing trend toward minimalism and reduced physical clutter. (Confidence: 0.7, Creative: 0.7)
3. Digital products offer greater flexibility and ease of use. (Confidence: 0.6, Creative: 0.6)
4. Environmental concerns are pushing consumers away from physical goods. (Confidence: 0.6, Creative: 0.8)
5. Advances in tech make digital-only options more functional. (Confidence: 0.8, Creative: 0.5)
6. Pandemic-era remote work increased demand for digital solutions. (Confidence: 0.7, Creative: 0.6)
7. Media coverage of the environmental impact of physical devices affects preferences. (Confidence: 0.5, Creative: 0.7)
8. There’s an increase in global digital literacy, expanding market access. (Confidence: 0.6, Creative: 0.6)
9. Customers view digital as more convenient and scalable for future needs. (Confidence: 0.7, Creative: 0.5)
10. Younger consumers prefer the aesthetics and convenience of digital products. (Confidence: 0.6, Creative: 0.6)

*Reflection:* 
"Have I considered all possible influences? Are there any surprising factors that could explain this shift?"

*Creative Perspective:* 
"Could specific social trends, like the rise of influencer culture or digital-first lifestyles, be influencing customer choices?"

**2.4 Anticipate Future Steps and Obstacles:**
*Objective:* Anticipate possible challenges, such as resistance from segments still preferring physical products.

*Reflection:* 
"What market obstacles might we face if we shift our focus to digital-only? Are there sub-segments that still prioritize physical products?"

*Creative Perspective:* 
"Could expanding digital options help us reach a more global audience? Are there emerging trends that we could leverage in our strategy?"

[end sys2]

[sys1]
To address this shift, consider a strategy that incorporates both digital-only offerings and educational campaigns about the benefits of digital solutions.
Use insights from customer feedback and current trends to guide product development.
Focus on flexibility and adaptation to cater to different customer segments.
[end sys1]


1. Introduction {#sec1}
===============

Transcription Factors (TFs) are the proteins that bind to specific DNA sequences and control the rate of genetic information transcription from DNA to mRNA, and they perform this function alone or in a complex with other genes. Particular interest goes to interactions where a given TF regulates other TFs, or itself \[[@r1]\]. A Master Regulator (MR) is a gene that occupies the very top of a regulatory hierarchy and is not influenced by the regulatory influence of any other gene \[[@r2]\]. MRs code for the TFs, which in turn alter the expression of downstream genes in the pathway \[[@r3]\].

Eukaryotic cellular functions are highly connected through networks of transcriptional regulators that regulate other transcriptional regulators \[[@r4]\]. The combinatorial cross-regulation of hundreds of sequence-specific TFs defines a regulatory network that underlies cellular identity and function \[[@r1]\]. However, the structure of core regulatory networks and their component sub-networks are largely undefined.

Obesity is a major health concern that leads to a chronic inflammatory state, including diabetes and cardiovascular disease \[[@r5]\]. Undirected Weighted Gene Co-expression Network Analysis (WGCNA) \[[@r6]\] can be used to find clusters of highly correlated genes (modules) and relating them to external traits. WGCNA has shown its potential to unravel the gene regulatory architecture of complex traits and diseases, for example: Alzheimer\'s disease in humans \[[@r7]\] parasite resistance in sheep \[[@r8]\], muscling in sheep \[[@r9]\], wool development in sheep \[[@r10]\], feed efficiency in cattle \[[@r11]\], muscle and meat properties in pigs \[[@r12]\] and obesity in a pig model \[[@r13]\]. WGCNA has been applied to RNA-Seq \[[@r13]\], to microarrays data, *e.g.* in the case of several human diseases \[[@r14], [@r15]\] and to porcine reproductive and respiratory syndrome virus studies \[[@r16]\].

The principle of WGCNA could be specifically targeted to TF co-expression networks to detect TF modules or the most promising hub-TFs that are potential regulators of the majority of the other TFs in a network. However, scale-free topology-based weighted co-expression networks of TFs involved in obesity (*e.g.* using the WGCNA method) have not been constructed yet. Such a TF network is expected to lead to a better understanding of the complex gene regulation of obesity development and subsequent pathogenesis and potentially lead to the detection of predictive genetic biomarkers for prevention and prediction of obesity and obesity-related diseases. For identification of drug targets and biomarker for diseases it is important to catalogue master regulatory genes (Master Regulators or MR). It is possible that some of the MRs may control other MRs, as has been found previously in, *e.g.* MRs: AphA and LuxR were found repressing each other\'s expression \[[@r17]\], and cell cycle transcriptional regulators regulated expression of other regulators \[[@r18]\]. So it is important to elucidate the complex networks of regulatory factors and pathways that are affected.

The main purpose of this research was to construct a WTFCN based on RNA-Seq data and to detect modules, pathways and novel regulator TFs governing regulatory processes in obesity and obesity-related diseases. In this research, we have used a publicly available dataset (GEO: GSE61271) containing adipose RNA-Seq data of a pig model for human obesity. To the best of our knowledge, this is the first study to use a network approach on TFs to detect regulatory genes, based on RNA-Seq transcriptomics for human obesity in an animal model.

2. Materials and methods {#sec2}
========================

The complete picture of the WTFCN construction workflow is presented in [Fig. (**[1](#F1){ref-type="fig"}**)](#F1){ref-type="fig"}.

2.1. Study Population {#sec2.1}
---------------------

In this study, we used publicly available RNA-Seq data and obesity state (lean or medium or obese) from a porcine model for human obesity (Gene Expression Omnibus (GEO) database of NCBI with accession number: GSE61271). The description of the dataset has previously been published in \[[@r13]\]. In short, an F2 pig resource population was created by crossing Gӧttingen minipig boars with Duroc sows. In total, 36 pigs were selected for RNA-Sequencing of subcutaneous adipose tissue based on a genetic obesity index. Those pigs were divided into groups representing three different DO: lean (n = 12), median (n = 12), and obese (n = 12). RNA-Seq was performed on the HiSeq2500 platform and obtained reads were aligned using STAR aligner using the SScrofa10.2.72 genome. Read counts were estimated using HTSeq, filtered based on low expression levels and then normalized using the voom() variance-stabilization function in the R-package Limma.

2.2. Identification of TFs {#sec2.2}
--------------------------

Transcription factors were identified from the PANTHER Classification System database \[[@r19]\]. Using the PANTHER protein class Ontology browser, we retrieved genes classified as "Transcription Factor" (PC00218) in *Sus Scrofa*. In our RNA-Seq dataset, 817 TFs of the 1312 identified TFs (list of TFs extracted at 19.04.2016) were expressed (Supplementary Data Sheet 1).

2.3. Single Weighted Transcription Factor Co-expression Network Construction {#sec2.3}
----------------------------------------------------------------------------

Using the Weighted Gene Co-expression Network Analysis (WGCNA) \[[@r6]\] approach we constructed a single WTFCN of 817 TFs using all 36 pigs. WTFCN construction followed the same workflow as described in \[[@r13]\]. Summarized, the adjacency matrix was created by calculating Pearson's correlations between the expression levels of the 817 TFs. To obtain a scale-free network, meaning it will contain a few nodes (TFs) that are highly connected with many other nodes, the adjacency matrix was raised to a power β () which was chosen based on the scale-free topology fitting index (R^2^ \> 0.8). The Topological Overlap Measure (TOM) was calculated representing the number of shared neighbors, and based on the dissimilarity TOM a TF clustering dendrogram was created using hierarchical clustering. TF modules, clusters of highly interconnected TFs, were detected as branches of the TF dendrogram using the cutreeHybrid algorithm and named according to a color. It is believed that highly interconnected genes act in similar biological pathways \[[@r20]\]. To identify biologically relevant modules, the Module Eigengene (ME) of each module was calculated as the first principal component of the module and checked using the proportion of variance in the module explained by the MEs. The MEs were correlated with the DO: groups of pigs were assigned to numerical values (lean was assigned a ´1´, median a ´2´ and obese a ´3´) which were fitted as continuous variable to calculate correlations between the module expression and phenotype. This resulted in the Module-Trait Relationship (MTR), which were used to select potential interesting associations with obesity. Potential interesting modules were further analyzed by detecting TFs that possess a high connectivity (degree of connections with other TFs), the so-called hub-TFs.The hubs are of great interest due to the number of the connections they have, as potential malfunction or disruption of the expressed levels of the hub TF will likely affect the complete module where it is present \[[@r20]\]. Furthermore, the TF significance (correlation between the individual TF expression levels and DO, where DO was fitted as a continuous variable) and the module membership (MM; correlation between TF expression profile and the ME of a given module) were used in order to identify biologically and statistically plausible genes.

### 2.3.1. Functional Enrichment Analysis {#sec2.3.1}

Functional enrichment analysis was performed using several approaches. First, the GOseq R Bioconductor package \[[@r21]\] was applied on selected modules. Before running GOseq in selected modules, pig Ensembl Gene IDs (Sscrofa10.2) were replaced with human Ensembl Gene IDs (GRCh38.p3) using the Ensembl BioMart database \[[@r22]\]. The non-native Gene Identifier \[[@r21]\], containing a complete set of

transcription factors analogous to *GRCh38.p3,* was created and used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis \[[@r23]\], taking length bias into account \[[@r24]\]. Obtained p-values were adjusted for multiple-testing using the Benjamini-Hochberg (BH) correction. GO terms and KEGG pathways were called significant when the adjusted p-value was below 0.05. If GOseq analysis did not yield any significant GO and/or KEGG pathways results, WEB-based Gene SeT AnaLysis Toolkit (WebGestalt) \[[@r25]\] was used for GO and/or KEGG pathway analysis. Secondly, gene names of single TFs were extracted using the BioMart database GRCh38.p3.

2.4. Construction of Sub-networks {#sec2.4}
---------------------------------

Using the same set of TFs and same network approach as described above, we constructed two sub-networks based on the DO: lean and obese. Resulting modules from the sub-networks were tested for Module Significance (MS). The MS is determined as the average of the TFs significance measure (the correlation between the expression values of TFs and DO) using both lean and obese pigs for all TFs in a given module. This measure is highly related to the correlation between module eigengene and the trait \[[@r20]\]. In order to compute MS for lean sub-network TFs expression values from lean and obese sub-networks were correlated with lean sub-network traits. In order to compute MS for obese sub-network TFs expression values from lean and obese sub-networks were correlated with obese sub-network traits. High MS means that the module is enriched with essential TFs. Functional and intra-modular analyses were performed on the modules with the strongest absolute MS.

### 2.4.1. Differentially Connected TFs in Sub-networks {#sec2.4.1}

Gene co-expression network analysis lends itself to identify entire groups of differentially regulated genes - a highly relevant endeavor in finding the underpinnings of complex traits that are polygenic in nature \[[@r26]\]. Differential network analysis concerns with identifying both differentially connected and differentially expressed genes \[[@r26]\]. In order to indicate which TFs were differently regulated in the networks we have compared difference in connectivity (DiffK) between the lean *vs*. obese sub-network, using the following formula:

*DiffK*(i) *= K*~1~(i) -- *K*~2~ (i)

where K~1~ presents the connectivity of the lean sub-network and K~2~ presents the connectivity of the obese sub-network. *K*~1~ and *K*~2~ were calculated by dividing each gene´s connectivity by the maximum sub-network connectivity, *e.g.*:

*K* ~1~(*i*)=k~1~(i)/max(*k~1~*)

When DiffK was positive, it meant that the TFs were more highly connected in the lean sub-network than in the obese sub-network, while a negative DiffK meant that the TFs were more highly connected in the obese sub-network than in the lean sub-network. To facilitate the comparison between the connectivity measures of each network, standardization was performed in each network by dividing the TF connectivity by the maximum network connectivity. The difference between the connectivity values of two sub-networks was defined as difference in connectivity (DiffK). DiffK values were distributed between -1 and 1. Genes were called differentially connected when the absolute value of DiffK was above 0.4. To measure differential TF expression between sub-networks the Student t-test was calculated, testing the difference between expression values of the lean and obese group. TFs were considered to be differentially expressed with an absolute value of the t-test statistic \> 1.96.

2.5. Analysis with Cytoscape Software {#sec2.5}
-------------------------------------

We expanded the analysis by integrating WTFCN results into Cytoscape_v3.3.0 \[[@r27]\]. Analysis and graphical representation was performed on significant modules of single WTFCN and on the modules having strongest absolute MS in the differentially co-expressed networks. Integration of genes expression and TF significance level into Cytoscape can open new perspectives in the analysis of biological networks, since the integration of topological analysis with expression data enhances the predictive power of the bioinformatics analysis \[[@r28]\]. NetworkAnalyzer, a Java plugin the Cytoscape platform \[[@r29]\], was used to perform topological analyses of the network in order to get an overview of the TFs' (represented by nodes) activity and interactions (represented by edges) between them within the modules. The NetworkAnalyzer tutorial and Centralities tutorial \[[@r28]\] were used to describe node and edge attributes (Table **S1**). Integrating TFs expression data from the modules having the strongest absolute MS into Cytoscape explains functional relationship, interaction, dynamics, role in expression within the network.

Furthermore, the Cytoscape app iRegulon (version 1.3) \[[@r30]\] was used to predict MRs in biologically relevant WTFCNs modules. iRegulon is a computational method which identifies MRs and predicts direct target genes in a set of human, mouse and Drosophila genes. iRegulon uses more than 9,000 known position weight matrices from various sources and different species and link them to candidate-binding TFs using a "motif2TF" procedure. This allows to link motifs of TFs from other species to candidate human TFs. Predicted MRs are rated and grouped by an enrichment score, the NES score (Normalized Enrichment Score). Predicted MRs having highest NES scores and present in analyzed modules as TFs were considered as present MRs. Targets (nodes) of predicted MRs which were overlapping with targets (nodes) of present MRs were considered as present MRs targets. As mentioned above, WTFCN constructs an undirected network; however regulatory links were drawn on WTFCN modules according to present MRs and their targets, detected by iRegulon.

3. Results and discussion {#sec3}
=========================

3.1. Single WTFCN Analysis {#sec3.1}
--------------------------

For analysis of network topology, the soft thresholding power (β = 6) was chosen based on the criterion of scale-free topology. Using the 817 TFs in WTFCN we identified 18 modules, each containing at least 20 TFs (Fig. **S1**). The Red module (65 TFs) of single WTFCN (Red~single~) had the highest correlation with obesity (MTR= 0.36, P = 3.00E-02) (Table **[1](#T1){ref-type="table"}**) and was selected for further analysis. GOseq analysis of this module resulted in the detection of 22 GO categories (FDR \< 0.05). Most of the GO terms indicated regulations in the immunity system: *e.g. immune system process* (Padj = 2.50E-06), *cytokine production* (Padj = 1.03E-05), and *interleukin-1 production* (Padj = 1.83E-05). GOSeq did not detect any significant KEGG pathways, but WebGestalt functional enrichment analysis resulted in KEGG pathways related to immunity (*e.g. TGF-beta signaling pathway*, Padj = 7.06E-16), osteoporosis (*e.g. Osteoclast differentiation*, Padj = 1.00E-04) and obesity and diabetes (*e.g. Circadian rhythm - mammal*, Padj = 2.33E-11). It was found \[[@r31]\] that disruption of the circadian clockwork with the external environment (*e.g.* light-dark cycle) might have a role in the development of metabolic disorders (*e.g.* obesity and diabetes). In the Red~single~ module we found two hub-TFs: *IRF5* and *SPI1*. The highest *ClosenessCentrality, Degree* and *Radiality* in this module had TF *ZYX,* indicating its highest number of linked edges, fastest information spread within the network and the highest centrality by regulating other TFs. Highest *Stress* had TF *MAFB*, indicating its relevance in connecting regulatory molecules and possible involvement in cellular processes. iRegulon detected MR *SPI1* (NES = 3.715). From 14 predicted targets four present targets (all upregulated TFs) were found in the module: *MAFB, IRF5, ATF3* and *ZNF710* (Fig. **[3](#F3){ref-type="fig"}**). As mentioned above, *SPI1* was also one of the hub-TFs in this module. Therefore importance of *SPI1* was confirmed as it is detected by two different methods. *SPI1* was found to be regulating adipose differentiation related to protein in macrophages \[[@r32]\]. A study by \[[@r33]\] suggested *SPI1* as a novel obesity therapeutic target, indicating that reduction in levels of *SPI1* leads to an increase in lipid accumulation. In one of our previous studies, SPI1 was detected as a candidate regulator of a possessing strong associations between obesity, immunity and osteoporosis \[[@r13]\]. Moreover, its possible key role in the link between obesity and osteoporosis was shown. In a recent study, \[[@r34]\] found significant interaction between the long noncoding RNAs cyp2c91 with the *SPI1* as regulator gene linked to obesity. *SPI1* can also regulate alternative splicing of target genes. Its target upregulated TF *MAFB* (also having highest *Stress* in Red~single~) was found to be dynamically regulated during adipogenesis \[[@r35]\]. Findings by \[[@r36]\] showed that *MAFB* expression in human adipocytes is upregulated during adipogenesis. *MAFB* was indicated as a regulator and a marker of adipose tissue inflammation, a process that subsequently causes insulin resistance. \[[@r37]\] found a SNP located close to *MAFB* associated with higher BMI, hypertension and diabetes, contributing to the risk of coronary artery disease and ischemic stroke. Another *SPI1* target, *ATF3,* may play a role in adipocyte hypoxia-mediated mitochondrial dysfunction in obesity \[[@r38]\]. Furthermore, ATF3 polymorphism was found associated with atherosclerotic disease \[[@r39]\]. Hub-TF *IRF5* (also regulated by *SPI1*) is well known to be of importance in immune system activity \[[@r40]\]. Moreover, a recent study by \[[@r41]\] found *IRF5* expression negatively associated with insulin sensitivity and collagen deposition in visceral adipose tissue. Mice lacking *IRF5* were healthy, and blocking of *IRF5* was suggested as a possibility for stopping fat formation. \[[@r42]\] reported *ZYX* as aging-related and a potential biomarker for osteoporosis. \[[@r43]\] found the association of *ZYX* polymorphisms with carcass and meat quality traits. Functional annotation results, hub-TFs and MR of the Red~single~ module from the single WTFCN analysis show its obvious biological role in immunity and obesity-related functions.

3.2. Sub-networks Analysis {#sec3.2}
--------------------------

Separate sub-networks were constructed on lean and obese DO groups. To obtain scale-free topology, a soft thresholding power of β = 10 was used in each sub-network. After WTFCN construction, 12 and 10 modules were detected in the obese and lean sub-network, respectively. In the obese sub-network the Red module (Red~obese~) had the strongest absolute MS; in the lean sub-network the Magenta module (Magenta~lean~) had the strongest absolute MS (Fig. **[2](#F2){ref-type="fig"}**, Table **[1](#T1){ref-type="table"}**). We homed in on the genes present in those modules to reveal their biological relevance.

The Red~obese~ module (69 TFs) had the strongest absolute MS (COR = -0.15) in the obese *vs*. lean network and was enriched with obesity-associated TFs (Fig. **[2](#F2){ref-type="fig"}**). After functional analysis, significant KEGG pathways were associated with immunity (*e.g. TGF-beta signaling pathway,* Padj = 1.73E-02) and osteoporosis (*e.g. Osteoclast differentiation*, Padj = 1.94E-02). Hub-TF *CBX8* was found after intra-modular analysis in this module. iRegulon detected MR *HOXD4* (NES = 4.710). From 34 predicted targets, 12 targets were present in the Red~obese~ module: *NR2F1, HOXD9, HOXD3, HOXB6, HOXC6, ZNF467, MSL3, CBX8, RARG, ZBTB12, L3MBTL3* and *SMAD6* (Fig. **[4](#F4){ref-type="fig"}**). Highest *ClosenessCentrality, Degree* and *Radiality* had *HOXD3* and *CBX8*, and the highest *Stress* had *HOXC6*. *CBX8* is a polycomb-group protein which re-modulates chromatin and modifies histones \[[@r44], [@r45]\] studying differentiating embryonic stem cells, found that depletion of CBX8 partially impairs the transcriptional activation of differentiating genes. HOXD3, HOXD4 and HOXC6 belong to the homeobox family genes, which encode highly conserved TFs playing an important role in morphogenesis in all multicellular organisms \[[@r40]\]. *HOXD4* may be involved in the development of atherosclerosis \[[@r46]\]. *HOXC6* was found upregulated after a fat loss \[[@r47]\]. Functional annotation results, hub-TFs and MR of the Red~obese~ module from the obese sub-network analysis show its biological role in immunity and osteoporosis related functions.

The Magenta~lean~ module (57 TFs) had the strongest absolute MS (COR = 0.26) of the TF significance in lean *vs*. obese network and enrichment with leanness-associated TFs (Fig. **[2](#F2){ref-type="fig"}**). The module displayed KEGG pathways in diabetes (*e.g. Basal transcription factors*, Padj = 4.20e-05) \[[@r48]\], type 2 diabetes, obesity and osteoporosis (*e.g. Notch signaling pathway*, Padj = 2.40E-03) \[[@r49], [@r50]\]. Hub-TF *VDR* (vitamin D receptor) was found after the intra-modular analysis. \[[@r51]\] indicated *VDR* association with obesity in type 2 diabetic subjects. \[[@r52]\] displayed that *VDR*-knockout mice had a slower growth rate and accumulated less fat mass than the wild type mouse. \[[@r53]\] suggested *VDR* as Quantitative Trait Loci (QTL) for fat mass variation in young Chinese men. *VDR* is a prime candidate gene for osteoporosis \[[@r54]\]. iRegulon analysis did not detect any present MRs within the module. Highest *Degree* had *LMO4* and *CCNK, Stress* - *CREB3L4.* A study by \[[@r55]\] identified *LMO4* as a novel modulator of leptin function. In a study by \[[@r56]\]. *CCNK* may regulate activity of cyclin-dependent kinases (Cdk), where Cdk4 knockout mouse had insulin deficient diabetes. *CREB3L4* was found inhibiting adipogenesis in preadipocytes \[[@r57]\]. Functional annotation results, hub-TFs and MR of the Magenta~lean~ module from the lean sub-network analysis show its biological role in diabetes, obesity and osteoporosis related functions.

### 3.2.1. Common Differentially Expressed TFs in Significant Modules of Sub-networks {#sec3.2.1}

TFs of the modules with the strongest absolute MS were compared if they are differentially expressed between lean and obese groups. Transcription factors present in the Red~obese~ module were aligned to the TFs in the Magenta~lean~ module. Common TFs found between the modules with their expression values are listed in Table **[2](#T2){ref-type="table"}**. The Red~obese~ module and the Magenta~lean~ module had 11 TFs in common: *RARG, RP11-644F5.10, PAX9, ZNF404, ZNF618, ZFP64, GTF2I, FHL2, NPRL2, CBX8* (hub-TF) and *PHF14.* A subset of the common TFs were found differentially expressed between the modules. *ZNF618* (Zinc Finger Protein 618) was downregulated in the Red~obese~ module and upregulated in the Magenta~lean~ module. Zink finger was found containing proteins function in various biological processes: gene transcription, translation, cytoskeleton organization, epithelial development, cell adhesion, chromatin remodeling, *etc*. \[[@r58]\]. *ZFP64, CBX8* and *PHF14* were upregulated in the Red~obese~ module and downregulated in the Magenta~lean~ module. Two of these three TFs we found to have associations with the immune system: *ZFP64* \[[@r59]\] and *CBX8* \[[@r60]\].

### 3.2.2. Analysis of Differentially Connected TFs in Sub-networks {#sec3.2.2}

In the lean *vs*. obese network we detected 53 differentially connected TFs which were more highly connected in the obese sub-network than in the lean sub-network. Four of them, *ZZZ3, GTF3C3, ZNF10* and *STAT2,* were also found to be differentially expressed (t-test \> 1.96). In addition, we found 26 differentially connected TFs that were more highly connected in the lean than in the obese sub-network, two of which (*NLRP3* and *FHL1*) were found to be differentially expressed (t-test \> 1.96). *STAT2* was described above. Very little is known about *GTF3C3* and *ZNF10*; however, in recent genome-wide meta-analysis study \[[@r61]\] *ZZZ3* was identified as a new locus for clinical classes of obesity. *NLRP3*, previously found as DE in \[[@r62]\], contributes to obesity-induced inflammation and insulin resistance \[[@r58]\] and was recently proposed as a new therapeutic target for diabetic complications \[[@r55]\]. Some of these results confirm the findings by \[[@r62]\].

### 3.2.3. Differences in the Construction of Different Regulatory Networks {#sec3.2.3}

WTFCN constructs an undirected network; however, using iRegulon we increased the confidence of our findings with respect to the hub-TFs. Recently, iRegulon in combination with Genomica software (a software tool for regulatory network construction) was used to build a regulatory network for insulin-mediated genes \[[@r63]\], but iRegulon in conjunction with WGCNA has never been used before. The use of different statistical bioinformatics methods and computational algorithms might lead to a difference in network construction and thereby altered results. For instance, Bayesian network integrates prior knowledge in a principled manner to increase the inference reliability \[[@r64], [@r65]\]. Higher performance of inference is achieved in simulated and real data \[[@r66]\]. Bayesian networks specify the joint probability distribution over all genes down to the conditional distributions. This network could be used to compute the probabilistic relationships between DO and trait. The Ordinary Differential Equation (ODE) models the regulatory system, but not directly infers the regulatory network. In ODE models, gene regulations are modeled by derivative equations which quantify the change rate of gene expression of one gene in the system as a function of expressions of all related genes that refer to its regulators. ODE model-based linear programming can also be used to quantify the rate of change of gene expression as a function of the expression of other genes \[[@r67]\]. The comparison of different network algorithms and different computing strategies in our data is beyond the scope of this original research article.

Conclusion
==========

In this study, we utilized publicly available RNA-Seq data (along with other metadata) from a pig model for human obesity deposited at the NCBI\'s Gene Expression Omnibus that were accessible through GEO Series accession number GSE61271. Using this data, we constructed a novel Weighted TF Co-expression Networks (WTFCN) to elucidate the complex patterns of gene regulation and master regulators involved in obesity development. A single WTFCN was constructed using TF expression data of lean, obese and median pigs, and TF modules acting on the development of obesity were identified by weighted scale-free topology network methods. Modules having strongest absolute MSin lean and obese sub-networks represent more specific TFs acting in lean or obese animals. Integrating TFs' expression data from these modules into Cytoscape explained functional relationship, interaction, dynamics and role in expression within the network. Our results obtained after WTFCNs and Cytoscape analysis were overlapping and supplemented each other and new regulatory links between MRs (*SPI1, HOXD4*) and their targets were drawn in WTFCNs, using iRegulon analysis. These findings support the idea of complex gene regulation, alternative splicing of target genes in different conditions. Our results suggested TFs that are potentially either causal or regulatory for degree of obesity. One of the most promising TFs in this research was *SPI1*, which was detected as a hub-TF in WTFCNA and as a MR in iRegulon. Such overlapping between two different methods confirms *SPI1* being a true hub-TF. As described above, *SPI1* is a strong regulator of obesity, immunity and osteoporosis. Furthermore, *SPI1* acts in macrophages and the phagocytosis process. Its four upregulated targets (*MAFB, IRF5, ATF3,* and *ZNF710*) are known in regulation of adipogenesis, BMI, coronary artery diseases, hypertension and diabetes. The above-mentioned evidence supporting our results indicates that *SPI1* is a possible major regulator in obesity. However, experimental validation should be performed before use and for better understanding of the complex role of interactions of highlighted nodes. Promising MRs and TFs should be validated by independent experiments. To the best of our knowledge, this is the first study that extended the principle of WGCNA and integrated it with Cytoscape in order to find hub-TFs, MRs or genes that are potential key regulators in obesity of the majority of TFs. Our results on TF networks for obesity are expected to lead to a better understanding of the complex TFs regulation of obesity development and provide insights into detection of TFs or drug targets, predictive genetic and biomarkers for prevention and prediction of obesity and obesity-related diseases.

Ethics Approval and Consent to Participate
==========================================

Not applicable.

Human and Animal Rights
=======================

No Animals/Humans were used for studies that are base of this research.

Consent for Publication
=======================

Not applicable.

LJAK was supported by a grant from the Danish Innovation Fund for the BioChild Project (grant number 0603-00457B and project website: [www.biochild.ku.dk](http://www.biochild.ku.dk)). HNK and RSJ thank the EU-FP7 Marie Curie Actions Grant (CIG-293511) and FeedOMICS project grant from DFF-FTP (grant number 4184-00268B) for partially supporting publication of this work and funding their time spent on this research.

SUPPLEMENTARY MATERIAL
======================

###### 

Supplementary material is available on the publisher's website along with the published article.

Conflict of interest
====================

The authors declare no conflict of interest, financial or otherwise.

![Workflow to build TF co-expression networks for different degrees of obesity using RNAseq data and further functional analyses.](CG-19-289_F1){#F1}

![Module significance of the obese and the lean sub-networks. Y axis represents correlation of TF significance, X axis - modules names. **A**) The Red obese module had the highest positive mean value (COR = -0.1477) in the obese *vs.* lean network and was enriched with obesity associated TFs, which indicated its significant difference in the obese sub-network. **B**) The Magenta lean module had the highest positive mean value (COR = -0.2598) in the obese *vs.* lean network and was enriched with leanness-associated TFs, which indicated its significant difference in the lean sub-network. (*For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.*)](CG-19-289_F2){#F2}

![Visualization in Cytoscape of interactions between TFs, MR and hub-TFs in the Red single module of WTFCN. The module contains both directed and undirected edges, as it is compiled by combining analyses from two sources: WTFCN and Cytoscape. Size of the nodes correspond to kWithin value (the bigger the node, the higher the value), color of the nodes corresponds to TF expression level (bright red - highly expressed/upregulated nodes; bright green - lowly expression, downregulated nodes). Shape of the nodes corresponds to significance level: square shape (P \< 0.05), circle shape (P \< 0.05). Circular layouts of the module are grouped clockwise according to nodes parameter kWithin, starting from the bottom of the picture. Both single WTFCN hub-TFs (IRF5 and SPI1) have highest kWithin values and were found mostly expressed in the module. Edges stroke color and width of edges correspond to EdgeBetweenness. Higher value of EdgeBetweenness wider and bluer stroke. Master regulator, hub-TF SPI1 is colored yellow, its regulatory links - red arrows points to regulated targets: MAFB, IRF5, ATF3 and ZNF710. (*For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.*)](CG-19-289_F3){#F3}

![Visualization in Cytoscape of interactions between TFs, MR and hub-TFs in the Redobese module of obese sub-network. The module contains both directed and undirected edges, as it is compiled by combining analyses from WTFCN and Cytoscape. Size of the nodes correspond to kWithin value (the bigger the node, the higher the value), color of the nodes corresponds to TF expression level (bright red - highly expressed/upregulated nodes; bright green - lowly expression, downregulated nodes). Circular layouts of the module are grouped clockwise according to nodes parameter kWithin, starting from the bottom of the picture. CBX8 (hub-TF) had highest kWithin value and was found mostly expressed in the module. Edges stroke color and width of edges correspond to EdgeBetweenness. Higher value of EdgeBetweenness wider and bluer stroke. Node of MR HOXD4 is colored in yellow. Red arrows indicate 12 regulated targets: NR2F1, HOXD9, HOXD3, HOXB6, HOXC6, ZNF467, MSL3, CBX8, RARG, ZBTB12, L3MBTL3, SMAD6. (*For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.*)](CG-19-289_F4){#F4}

###### 

Correlations of each network module with obesity (MTR/P values).

Text colors of modules´ names correspond to their color.

  ------------------------------------------------------------------------------------
  **Modules**       **Single WTFCN**\   **Obese Sub-network**   **Lean Sub-network**
                    **MTR**\                                    
                    **(P value)**                               
  ----------------- ------------------- ----------------------- ----------------------
  **Black**         \-                  -0.14                   -0.20

  **Blue**          -0.01 (1)           -0.03                   0.002

  **Brown**         0.035 (0.8)         -0.04                   0.14

  **Green**         -0.13 (0.4)         -0.05                   -0.18

  **Greenyellow**   \-                  -0.04                   \-

  **Magenta**       \-                  -0.03                   -0.26

  **Pink**          \-                  -0.11                   -0.14

  **Purple**        \-                  -0.08                   -0.13

  **Red**           0.36 (0.03)         -0.15                   0.08

  **Tan**           \-                  -0.12                   \-

  **Turquoise**     -0.3 (0.08)         -0.10                   -0.03

  **Yellow**        0.11 (0.5)          0.02                    -0.22
  ------------------------------------------------------------------------------------

(*For interpretation of the references to color in this Table, the reader is referred to the web version of this paper*.)

###### 

Common TFs comparing significant modules between lean and obese sub-networks.

Text colors of modules´ names correspond to their color. Colored TFs are hub-TFs of the module corresponding to its color. Expression values in bold indicate differential expression of the same TF between two different sub-networks. Cell colors indicate differential expression values: red - up regulation, green - down regulation.

  **TFs in Common**   **Expression Values of Modules in Sub-networks**   
  ------------------- -------------------------------------------------- -----------------------
  \-                  **The Red~obese~ module**                          **The Magenta~lean~**
  *RARG*              0.79                                               0.77
  *RP11-644F5.10*     0.61                                               0.72
  *PAX9*              0.22                                               0.43
  *ZNF404*            0.55                                               0.57
  *ZNF618*            **0.36**                                           **0.76**
  *ZFP64*             **0.85**                                           **0.60**
  *GTF2I*             0.61                                               0.51
  *FHL2*              0.60                                               0.86
  *NPRL2*             0.80                                               0.81
  *CBX8*              **0.92**                                           **0.64**
  *PHF14*             **0.72**                                           **0.45**

(*For interpretation of the references to color in this Table, the reader is referred to the web version of this paper.*)