Multi-tissue Analysis of Co-expression Networks by Higher-Order Generalized Singular Value Decomposition Identifies Functionally Coherent Transcriptional Modules

<div><p>Recent high-throughput efforts such as ENCODE have generated a large body of genome-scale transcriptional data in multiple conditions (e.g., cell-types and disease states). Leveraging these data is especially important for network-based approaches to human disease, for instance to identify coherent transcriptional modules (subnetworks) that can inform functional disease mechanisms and pathological pathways. Yet, genome-scale network analysis across conditions is significantly hampered by the paucity of robust and computationally-efficient methods. Building on the Higher-Order Generalized Singular Value Decomposition, we introduce a new algorithmic approach for efficient, parameter-free and reproducible identification of network-modules simultaneously across multiple conditions. Our method can accommodate weighted (and unweighted) networks of any size and can similarly use co-expression or raw gene expression input data, without hinging upon the definition and stability of the correlation used to assess gene co-expression. In simulation studies, we demonstrated distinctive advantages of our method over existing methods, which was able to recover accurately both common and condition-specific network-modules without entailing <i>ad-hoc</i> input parameters as required by other approaches. We applied our method to genome-scale and multi-tissue transcriptomic datasets from rats (microarray-based) and humans (mRNA-sequencing-based) and identified several common and tissue-specific subnetworks with functional significance, which were not detected by other methods. In humans we recapitulated the crosstalk between cell-cycle progression and cell-extracellular matrix interactions processes in ventricular zones during neocortex expansion and further, we uncovered pathways related to development of later cognitive functions in the cortical plate of the developing brain which were previously unappreciated. Analyses of seven rat tissues identified a multi-tissue subnetwork of co-expressed heat shock protein (Hsp) and cardiomyopathy genes (<i>Bag3</i>, <i>Cryab</i>, <i>Kras</i>, <i>Emd</i>, <i>Plec</i>), which was significantly replicated using separate failing heart and liver gene expression datasets in humans, thus revealing a conserved functional role for Hsp genes in cardiovascular disease.</p></div

Performance comparison for C3D, WGCNA and DiffCoEx methods

<p><i>Top</i>, three cluster types (“common” “nested” and “overlapping”) were simulated in conditions where the cluster size () is reported for both the intersection and union part of the clusters. <i>Bottom</i>, for each method the average TPR and FPR () across 20 replicated datasets were calculated and reported for the simulated cluster densities. For C3D analysis (blue lines) we required each cluster to be detected with a misclassification error rate (MER) of 5% or 20% and . For WGCNA (red line) and DiffCoEx (green line) we considered two “default values” for the cut-off threshold, which were chosen according to the WGCNA guidelines (see <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1004006#pgen.1004006.s010" target="_blank">Text S1</a> for details).</p

Description of the cluster structures used in the simulation studies.

<p>We simulated three cluster types: “common” (<i>Cluster pattern 1</i>), “nested” (<i>Cluster pattern 2</i>) and “overlapping” (<i>Cluster pattern 3</i>) that are shared across three or more conditions. For <i>Cluster pattern 2</i> and <i>Cluster pattern 3</i>, the “intersection cluster” is defined by the nodes in common to all conditions (red square) whereas the “union cluster” is defined by the nodes in common to all conditions plus the nodes present in individual conditions (black square).</p

Human co-expression cluster 1.

<p><i>Top left</i>, each node in the network represents a gene and, in keeping with <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1004006#pgen.1004006-Liu1" target="_blank">[61]</a>, for each gene we highlight significant up-regulation in VZ (red) or CP (green) as compared with the other neocortex regions. Genes that are were not differentially expressed between neocortex regions are coloured in grey. Genes present in relevant KEGG pathways (p53 signaling, ECM-receptor interaction, Cell cycle and DNA replication) are extracted from the main network and highlighted. <i>Top right</i>, functional annotation for the network: top five significant GO biological processes and KEGG pathways (full list in <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1004006#pgen.1004006.s007" target="_blank">Table S3</a>). <i>Bottom left</i>, summary of cell-type enrichment analysis expressed as (Benjamini and Hochberg (BH)-adjusted <i>p</i>-value, Cten analysis). <i>Bottom right</i>, graph with the significant protein-protein interactions (PPI), including the overall significance of the directed PPI network (DAPPLE analysis, ). The colour scale on the right indicate the significance of the detected PPI.</p

Co-expression clusters identified in all rat tissues.

<p>For each rat cluster detected in all seven tissues we report the number of probe sets, the top five functional categories and their statistical significance (full list in <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1004006#pgen.1004006.s006" target="_blank">Table S2</a>), the summary of cell-type enrichment statistics expressed as (Benjamini and Hochberg (BH)-adjusted <i>p</i>-value, Cten analysis) and the graph with the significant protein-protein interactions (PPI), including the overall significance of the directed PPI network (DAPPLE analysis). The colour scale on the right indicate the significance of the detected PPI.</p

<i>Rat cluster 1</i> shows co-expression between Hsp and cardiomyopathy genes which is conserved with human heart and liver tissues

<p>(A) Network of 135 annotated rat genes identified by C3D as co-expressed in heart, aorta, liver and skeletal muscle tissues (). In each tissue we selected the top 5% of edges based on the (absolute) covariance between gene expression profiles and then calculated the average covariance across the four tissues. Edges are represented by lines connecting nodes (genes) and the thickness of the line is proportional to the average covariance value. Within the network, heat shock protein (Hsp) and cardiomyopathy genes are highlighted in blue and red, respectively. The Kendall correlations between the expression profiles of Hsp and cardiomyopathy genes are graphically represented as sub-networks separately for each tissue. Line thickness is proportional to the value of the Kendall correlation. (B) Enrichment for functional categories (, full list in <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1004006#pgen.1004006.s006" target="_blank">Table S2</a>) and for disease association (adjusted , details in <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1004006#pgen.1004006.s007" target="_blank">Table S3</a>). (C) Significant protein-protein interaction (PPI) network () where the Hsp and cardiomyopathy genes showing conserved PPI are highlighted (blue and red circles). (D) Conserved co-expression network detected in heart tissue samples from patients with advanced idiopathic or ischemic cardiomyopathy. The network includes all human orthologous genes of the genes in <i>rat cluster 1</i> that have significant edges by covariance selection (). (E) Conserved co-expression network detected in liver tissue samples from healthy volunteers. The network includes all human orthologous genes of the genes in <i>rat cluster 1</i> that have significant edges by covariance selection ().</p

Computational methods for direct cell conversion

<p>Directed cell conversion (or transdifferentiation) of one somatic cell-type to another can be achieved by ectopic expression of a set of transcription factors. Since the experimental identification of transcription factors for transdifferentiation is extremely time-consuming and expensive, there are still relatively few transdifferentiations achieved in comparison to the number of human cell-types. However, the growing volume of transcriptional data available and the recent introduction of data-driven algorithmic approaches that predict factors for transdifferentiation holds great promise for accelerating this field. Here we review those computational methods whose in-silico predictions have been experimentally validated, highlighting differences and similarities. Our analysis reveals that the factors predicted by each method tend to be different due to varying source cells used, gene expression quantification and algorithmic steps. We show these differences have an impact on the regulatory influences downstream, with some methods favoring transcription factors regulating developmental progression and others favoring factors regulating mature cell processes. These computational approaches offer a starting point to predict and test novel factors for transdifferentiation. We argue that collecting high-quality gene expression data from single-cells or pure cell-populations across a broader set of cell-types would be necessary to improve the quality and consistency of the in-silico predictions.</p

Mesangial cells contain soluble factors polarising macrophages depending on their genetic background.

<p>Supernatant transfer from P5 MCs (WKY or LEW) onto WKY or LEW bone marrow-derived macrophages (BMDMs). <b>A</b>. qRT-PCR assessment of M1 macrophage markers (<i>Nos2</i>, <i>Il1b</i>, <i>Tnfa</i>, <i>Il12b</i>) in WKY (left panel) and LEW (right panel) BMDMs following incubation with MCs supernatants from both strains. Note the increased expression of all M1 markers in WKY BMDMs following incubation with WKY MC supernatants. <b>B</b>. qRT-PCR assessment of M2 macrophage markers (<i>Mrc1</i>, <i>Il10</i>) in WKY (left panel) and LEW (right panel) BMDMs following incubation with MCs supernatants from both strains. Note the increased expression of all M2 markers in LEW BMDMs following incubation with WKY MC supernatants. *, P<0.01; **P<0.001; the results are representative of two experiments, n = 4 rats/strain used per experiment.</p

Genome-wide expression analysis by microarrays identifies two distinct transcriptomes in WKY and LEW MCs.

<p><b>A</b>. Dendogram representing the hierarchically clustering using all differentially expressed genes between WKY and LEW MCs in basal (unstimulated) and TNFα-stimulated cells. This shows that the strain (WKY or LEW) and the treatment (basal or TNFα-stimulated) effects cluster in 4 distinct groups. <b>B</b>. Validation of the markedly differentially expressed transcripts (fold change >10) by qRT-PCR. The upper panel shows the top differentially expressed candidates identified by microarray analysis with a false discovery rate (FDR) <0.01. The positive fold change (FC) values designate up-regulation in WKY MCs and negative FC values designate up-regulation in LEW MCs in basal conditions. The lower panel shows qRT-PCR validation for the transcripts showing differential expression in the microarray dataset. n = 4 rats, **, P<0.001. <b>C</b>. KEGG pathway analysis applied to differentially expressed genes between WKY and LEW MCs in basal state (WKY-LEW)_basal and (WKY-LEW)_TNFα identified the DNA replication and the vasculature development pathways as the most significant ones respectively. <b>D</b>. KEGG pathway analysis in LEW and WKY MCs treated with TNFα [(LEW)-(LEW) _TNFα and (WKY)-(WKY)_TNF] identified strain-specific pathways (shown in blue in the LEW and red in WKY) upon TNFα stimulation. Note that the DNA-replication pathway is the most significant one in the WKY MCs following TNFα stimulation.</p

Mesangial cells (MCs) inhibit splenocyte proliferation.

<p>Rat MCs were cultured until passage 5 (P5) and incubated with Con-A stimulated or non-stimulated splenocytes for 72 hours. Cell proliferation was measured by incorporation of tritiated thymidine (³H-TdR). <b>A</b>. Co-culture of WKY MCs and splenocytes resulted in a significant decrease in ConA-stimulated splenocyte proliferation in both 1∶10 and 1∶20 ratio of MC:Splenocyte. **, P<0.001 compared with splenocytes alone. The results are representative of three independent experiments. Cpm, counts per minute <b>B</b>. Co-culture of LEW MCs and splenocytes resulted in a significant decrease in Con-A-stimulated splenocyte proliferation in both 1∶10 and 1∶20 ratio of MC:Splenocyte. **, P<0.001 compared to splenocytes alone. The results are representative of three independent experiments. Cpm, counts per minute.</p