Investigating trait variability of gene co-expression network architecture in brain by controlling for genomic risk of schizophrenia

Abstract

The effect of schizophrenia (SCZ) genetic risk on gene expression in brain remains elusive. A popular approach to this problem has been the application of gene co-expression network algorithms (e.g., WGCNA). To improve reliability with this method it is critical to remove unwanted sources of variance while also preserving biological signals of interest. In this WCGNA study of RNA-Seq data from postmortem prefrontal cortex (78 neurotypical donors, EUR ancestry), we tested the effects of SCZ genetic risk on co-expression networks. Specifically, we implemented a novel design in which gene expression was adjusted by linear regression models to preserve or remove variance explained by biological signal of interest (GWAS genomic scores for SCZ risk-(GS-SCZ), and genomic scores- GS of height (GS-Ht) as a negative control), while removing variance explained by covariates of non-interest. We calculated co-expression networks from adjusted expression (GS-SCZ and GS-Ht preserved or removed), and consensus between them (representative of a "background" network free of genomic scores effects). We then tested the overlap between GS-SCZ preserved modules and background networks reasoning that modules with reduced overlap would be most affected by GS-SCZ biology. Additionally, we tested these modules for convergence of SCZ risk (i.e., enrichment in PGC3 SCZ GWAS priority genes, enrichment in SCZ risk heritability and relevant biological ontologies. Our results highlight key aspects of GS-SCZ effects on brain co-expression networks, specifically: 1) preserving/removing SCZ genetic risk alters the co-expression modules; 2) biological pathways enriched in modules affected by GS-SCZ implicate processes of transcription, translation and metabolism that converge to influence synaptic transmission; 3) priority PGC3 SCZ GWAS genes and SCZ risk heritability are enriched in modules associated with GS-SCZ effects. Overall, our results indicate that gene co-expression networks that selectively integrate information about genetic risk can reveal novel combinations of biological pathways involved in schizophrenia.

Fig 1. Analytical pipeline.

Step 1: Preliminary processing of RNA-Seq data and calculation of variables for gene expression adjustment (genomic score of schizophrenia risk- GS-SCZ, genomic score of height- GS-Ht, cell type proportions by bulk RNA-Seq deconvolution, ancestry principal components- genomic PCs, observed and hidden technical artifacts) (part of the Step1 panel was created with templates from bioRender (https://app.biorender.com/). Step 2: Selection of covariates of interest and of non-interest for gene expression adjustment by using explained variance in gene expression and multicollinearity criteria. Step 3: Calculating eight linear regression models with gene expression as dependent variable, GS-SCZ and GS-Ht as covariates of interest alternatively preserved or removed, and other covariates of non-interest (age, sex, cell type proportions, genomic PCs and technical artifacts). Step 4: Calculation of gene co-expression networks with Weighted Gene Co-expression Analysis (WGCNA) routines from eight expression inputs (residuals extracted from step 3 linear models). Step 5: Calculation of background co-expression networks (networks neutral to genomic scores effects; schematically represented as C1-C3); identification of differential GS-SCZ and GS-Ht preserved modules (modules with weak or no correspondence in background, annotated as M1-M3); legend: M1 = example of a module with strong correspondence in background, M2 = module with weak correspondence and M3 = module with no correspondence in background. Steps 6–7: Further evidence of biological significance and specificity of modules that preserve GS-SCZ effects by convergence of SCZ genetic risk by internal and external validation.

Fig 2. Kernel density maps- the distribution of overlap between GS-SCZ, GS-Ht preserved modules and background modules.

Annotated modules on the x axis represent fragmented modules with weaker conservation in their background modules.

Fig 3. Gene overlap of modules that concentrate genetic risk of SCZ risk in all co-expression networks of preserved GS3-SCZ-GS5-SCZ and GS3-GS5 height.

Fig 4. Functional profiling of GS-SCZ preserved modules with cumulative evidence for genomic risk effect on co-expression.

Transcription, translation and metabolic ontologies are enriched in gene sets originated from modules with MEs positively correlated with GS-SCZs SCZ; nervous system development and functionality ontologies are enriched in gene sets originated from modules with MEs negatively correlated with GS-SCZs SCZ (highlighted by red rectangles). Legend: BG_GS3-SCZ_ / BG_GS5-SCZ_x = fragments from SCZ risk GS3-SCZ or GS5-SCZ preserved modules overlapped with background (BG) modules.

Fig 5. Functional profiling of GS height preserved modules with cumulative evidence for genomic risk effect on co-expression.

Similar functional divergence in biological processes enrichment like GS-SCZ preserved modules: transcription, translation and metabolic ontologies are enriched in gene sets originated from modules with MEs positively correlated with GS-SCZs SCZ; nervous system development and functionality ontologies are enriched in gene sets originated from modules with MEs negatively correlated with GS-SCZs SCZ (highlighted by red rectangles). In addition, one set from GS3-Ht preserved lightgreen module negatively correlated with genomic scores of height was enriched in synapse related GO:BP. Legend: BG_GS3_ / BG_GS5_x = fragments from height GS3 or GS5 preserved modules overlapped with background (BG) modules.

Fig 6. Gene overlap of modules that concentrate genetic risk of SCZ risk in co-expression networks of preserved GS3-SCZ-GS5-SCZ and GS3-GS5 height from LIBD dataset and the independent replication sample (PITT sample).