ACDC is designed to detect dependencies between gene-gene co-expression (or connectivity) and a set of external features that can be either exposures or responses. That is, ACDC is applied to test for evidence of association between measures of co-expression and measures of external features. Notably, the external features are not constrained to be categorical, the typical requirement (2), but could be continuous or ordinal.

The concept of covariance can be used to quantify the dependence between two random variables and thus to quantify gene-gene co-expression. It is possible for the covariance of a pair of genes to depend on external features. For example, suppose in a biological pathway, two genes tend to be co-regulated and thus co-expressed, resulting in positive covariance. A perturbation to the pathway could alter that relationship, resulting in a change in co-expression and thus a change in covariance. If candidate perturbagens and the expression of genes in the pathway are measured, ACDC may be applied to detect these types of effects simultaneously for the multiple genes and perturbagens. Using a similar rationale, ACDC could be applied to detect downstream results of pathway perturbations if the affected phenotypes are measured.

Suppose all individuals have measurements for all M gene expression features in the set, referred to here as a “module”, and all P external features. Assume the vector of P external features are distributed as multivariate normal,

with x_p representing each external feature. Though we describe x as multivariate normal here, we can relax this assumption in practice and allow other distributions and variable types, as in a design matrix.

Suppose the M gene expression features are also distributed as multivariate normal with the covariance matrix depending on x,

where g_j denotes the expression of gene j. The covariance matrix can be represented by,

The off-diagonal elements of Σ_g can be considered measures of co-expression, and (for given values of x) estimated in the conventional way,

Note that this is essentially an average over individuals. Letting s denote an individual, each contribution is,

These individual components have approximately the same expectation as the scaled sum, therefore they can also be described as estimators for σ_{j, k}. We leverage this property to test for dependencies between the covariances and the external features.

We can denote the co-expression profile for a given module as,

We are interested in dependencies that may exist between the external features, x, and the gene-pair covariances, the off diagonals of Σ_g. If we have a single external feature or a single pair of genes, conventional general linear modeling (GLM) approaches could be used to relate x to C. For multiple gene pairs and external features, CCA can be applied, or sparse CCA for high dimensional settings. CCA finds min [G, P] linear combinations, a∈ℝ^P, b∈ℝ^G, of C and x, respectively, that maximize the correlation,

for example, for the first pair of canonical variables. Note that CCA can be applied even if G and/or P is equal to one (4). Wilks–Lambda can be used to conduct a joint hypothesis test of whether the correlation coefficients found by CCA are significantly different from zero,

A rejected test implies dependent co-expression, i.e., that there are linear combinations of gene-gene covariances associated with linear combinations of external features.

False discovery rates (FDR) can be computed using the Benjamini–Hochberg (BH) (5) method when multiple modules are tested and parametric assumptions apply. If severe departures from the assumed distributions may be present, permutation-based approaches such as the Millstein and Volfson (MV) FDR (6) method can be used.

and is also available from the CRAN repository, https://cran.r-project.org/web/packages/modACDC/index.html

To identify modules of correlated genes, we applied the Partition data reduction method (10, 11), an agglomerative approach that requires the user to specify an acceptable proportion of information loss when collapsing all features to a single measure such as the mean. Selection of the information loss threshold was guided by the aim to maximize information explained in the ACT score while minimizing noise. Further explanation is provided in the Supplementary material (12). We used an information loss constraint of 0.35 which corresponds to a minimum of 65% information from the non-reduced data captured by each new feature, as assessed by the intraclass correlation coefficient (ICC). This reduction threshold resulted in roughly 50% reduction in features when compared to the full dataset (Supplementary Figure 1).

In analyses of individual genes within modules of interest, gene-ACT score relationships were modeled using ordinal logistic regression, adjusting for patient age, race, sex, data collection site, and the first three principal components (PCs), which may capture global dependencies due to cell-type composition and technical artifacts, of each gene expression data set. Associations were identified at the 0.05 FDR level.

To clarify the novel attributes of ACDC, a comparative analysis was conducted in the ABRIDGE cohort using CoXpress. The ACT score was dichotomized at the median value to indicate better vs. worse asthma control to conform to the coXpress requirement of a binary phenotype. The Pearson correlation coefficient was used as the similarity measure, and for module identification the dendrogram was cut at a height of 0.35 for consistency with the Partition approach.

ACDC was performed on 65 modules identified by Partition in the ABRIDGE dataset. The results for the top five modules based on BH FDR can be found in Table 2. Evidence suggestive of differential co-expression as determined by CCA Wilks–Lambda p ≤ 0.05 was found for two modules including genes NOD-like Receptor Family Pyrin Domain Containing 12 (NLRP12), Meteorin Like, Glial Cell Differentiation Regulator (METRNL), and Ghrelin And Obestatin Prepropeptide (GHRL) in module A (BH FDR = 0.0737), and Adenosine A3 Receptor (ADORA3), Arachidonate 15-Lipoxygenase (ALOX15), and Indoleamine 2,3-Dioxygenase 1 (IDO1) in module B (BH FDR = 0.1569). We also computed the non-parametric, permutation based FDR estimate Millstein–Volfson (MV) to account for departures from the normality assumption by the ACT variable, which is ordinal. However, the results of the MV FDR test are in approximate agreement with the BH FDR results, yielding two modules with evidence of differential co-expression [(Supplementary Figure 2), FDR = 0.0554, 95% CI: (0.0054, 0.5742)].

To further explore the relationship between co-expression of genes in modules A and B and asthma control, Kruskal–Wallis tests were performed to determine whether covariance measures for all possible pairs of these genes differ across levels of the ACT score components. Eight of the total 24 tests resulted in p-values less than 0.05, with the top six coming from module B. The most significant test involved the co-expression of IDO1 and ADORA3 and the frequency of waking from wheezing in the past 6 months (p = 0.0021; Figure 2).

We performed ACDC using data from CAMP in an attempt to replicate results observed for the top five modules identified in ABRIDGE. We found evidence of differential co-expression for module B (p = 0.0315) but not module A (p = 0.6823; Table 2). Also, evidence of differential co-expression was observed for gene pairs in modules C and D, which were not significant in ABRIDGE, Interleukin 5 Receptor Subunit Alpha (IL5RA) and Peripheral Myelin Protein 22 (PMP22) in module C, and Interleukin 17 Receptor B (IL17RB) and Interleukin 6 (IL6) in module D. Note that these results have not been adjusted for multiple testing.

Kruskal–Wallis tests were also performed for the same gene-pair covariances tested in ABRIDGE. Of the 24 tests performed, there were three with p-values less than 0.05, all from module B. The most significant test compared the co-expression of IDO1 and ADORA3 across levels of rescue bronchodilator use in the past 7 days (p = 0.02) (Figure 3). Additionally, we performed Kruskal–Wallis tests for all gene-pair covariances and 7-day ACT components for the three modules with CCA Wilks–Lambda p-values below 0.05. Of the 20 tests performed, the same three pairs from module B showed evidence of differential co-expression, but no others had p-values less than 0.05.

Following the differential co-expression analysis, we performed ordinal logistic regression for each of the 13 genes in the top five modules and found increased risk of suboptimal acute asthma control for all genes in modules B and C, after adjusting for covariates (Table 3). Higher expression of ADORA3, ALOX15, and IDO1 was associated with suboptimal 6-month ACT scores (Figure 4).

Adjusted ordinal logistic regressions were performed for the same 13 genes as the ABRIDGE cohort (Section 3.3). In the CAMP cohort, the regressions also showed highly statistically significant associations for all genes in modules B and C, and non-significant associations for modules A and D (Table 3). Unlike the results from ABRIDGE, a significant protective effect was seen for NOD-like Receptor Family CARD Domain Containing 3 (NLRC3) [OR: 0.3926, 95% CI: (0.1864, 0.8268)]. Associations between these genes and 7-day ACT scores (Figure 5) also imply that increasing gene expression is associated with suboptimal acute asthma control.

The five most differentially co-expressed modules identified by the CoXpress analysis can be seen in Table 4. As a rule of thumb for identifying differentially co-expressed modules, the coXpress authors suggest pr_g1 ≤ 0.05 and pr_g2≥0.05, which implies correlations different than zero in one of the classes but not the other. None of the ABRIDGE modules met this threshold, and values of pr_g1, pr_g2 ≤ 0.05 for all of the five top modules indicate that the intra-module correlations are non-zero for patients with both better and worse asthma control. We note that genes ADORA3 and ALOX15 appear in module 1, the most differentially co-expressed module.