Gene communities in co-expression networks across different tissues

Abstract

With the recent availability of tissue-specific gene expression data, e.g., provided by the GTEx Consortium, there is interest in comparing gene co-expression patterns across tissues. One promising approach to this problem is to use a multilayer network analysis framework and perform multilayer community detection. Communities in gene co-expression networks reveal groups of genes similarly expressed across individuals, potentially involved in related biological processes responding to specific environmental stimuli or sharing common regulatory variations. We construct a multilayer network in which each of the four layers is an exocrine gland tissue-specific gene co-expression network. We develop methods for multilayer community detection with correlation matrix input and an appropriate null model. Our correlation matrix input method identifies five groups of genes that are similarly co-expressed in multiple tissues (a community that spans multiple layers, which we call a generalist community) and two groups of genes that are co-expressed in just one tissue (a community that lies primarily within just one layer, which we call a specialist community). We further found gene co-expression communities where the genes physically cluster across the genome significantly more than expected by chance (on chromosomes 1 and 11). This clustering hints at underlying regulatory elements determining similar expression patterns across individuals and cell types. We suggest that KRTAP3-1, KRTAP3-3, and KRTAP3-5 share regulatory elements in skin and pancreas. Furthermore, we find that CELA3A and CELA3B share associated expression quantitative trait loci in the pancreas. The results indicate that our multilayer community detection method for correlation matrix input extracts biologically interesting communities of genes.

Fig 1. Schematic of a multilayer gene co-expression network.

The intralayer edges, shown by the solid lines, represent co-expression. The interlayer edges, shown by the dashed lines, connect the same gene across layers.

Fig 2. Determination of the resolution parameter value by CHAMP.

(a) Multilayer correlation matrix. (b) Multilayer network obtained by graphical lasso. The convex hull of the lines in the (γ, Q) plane, each of which corresponds to a node partitioning, is a piecewise linear curve with the transition values indicated by a cross and change in the line color. Each line segment corresponds to the optimal node partitioning in the corresponding range of γ.

Fig 3. Composition of each community by layer, i.e., tissue.

(a) Multilayer correlation matrix, γ = 1. (b) Multilayer correlation matrix, γ = 3. (c) Unweighted multilayer network obtained by graphical lasso, γ = 1. (d) Unweighted multilayer network obtained by graphical lasso, γ = 3. (e) Weighted multilayer network obtained by graphical lasso, γ = 1. (f) Weighted multilayer network obtained by graphical lasso, γ = 3. The darker shades indicate nodes corresponding to genes that only appear in one layer in the given community. The lighter shades indicate nodes corresponding to genes that appear in multiple layers in the community.

Fig 4. Location of genes on chromosomes, colored by community.

There is a colored circle for each associated community pointing to each gene. Note that a gene can belong to more than one community, denoted by multiple colored circles next to each other horizontally pointing to the same gene. This figure allows us to visually see clusters of genes on specific chromosomes and their associated community.

Fig 5. Expression and co-expression analysis of a cluster of genes in community 5 on chromosome 17.

The co-expression matrices for these genes in (A) skin and in (B) pancreas are shown. The average expression for each gene in these tissues is shown in the bar graphs. The location of these genes on chromosome 17 is shown in (C), with arrows (colored according to the associated tissue) pointing from putative regulatory elements to highly co-expressed genes. (D) The panel shows different measures of the regulatory potential of this genome section. From top to bottom: 1. H3K27AC modification to histone H3 within the region, which often correlates with activation of transcription and is associated with active enhancers in a given tissue available through ENCODE database [119]. 2. DNAse1 hypersensitivity sites. They are sections of the genome that are cut by DNAse1 enzyme. Given that the chromatin has to be “open” for the DNAse to access the sequence, the sequences that are cut by DNAse indicate open chromatin, which is in turn associated with regulatory activity. Data are available through ENCODE database [119]. 3. Enhancer/promoters. These are sequences that are predicted as enhancers (gray) and promoters (red) from the GeneHancer database [120]. 4. Established interactions between regulatory regions and genes as documented by GeneHancer database [120]. These data sets combined with our co-expression analysis provide a novel outlook into potential topologically associated domains that may be regulated by specific sequences in a tissue-specific manner.

Fig 6. A schematic of SNPs in an enhancer region (gray box) that affect the expression of CELA3A (blue box) and CELA3B (orange box) in the pancreas and are associated with blood phosphate concentration.

(A) Expression levels of CELA3A and CELA3B, and blood phosphate concentration when the derived alleles for the putatively causal SNPs are absent. (B) The presence of the derived alleles for the putatively causal SNPs decreases the expression level of CELA3A, increases the expression level of CELA3B, and decreases the blood phosphate concentration.