Covariance-based sample selection for heterogeneous data: Applications to gene expression and autism risk gene detection

TL;DR

This paper introduces COBS, a covariance-based sample selection method that enhances gene risk detection in heterogeneous brain tissue data by selecting more homogeneous samples, improving the power of downstream analysis.

Contribution

The paper presents COBS, a novel method for selecting homogeneous samples based on covariance, improving gene detection in heterogeneous datasets for autism research.

Findings

COBS improves DAWN's prediction of autism risk genes.

Sample selection based on covariance enhances detection power.

Method validated using data from two different time points.

Abstract

Risk for autism can be influenced by genetic mutations in hundreds of genes. Based on findings showing that genes with highly correlated gene expressions are functionally interrelated, "guilt by association" methods such as DAWN have been developed to identify these autism risk genes. Previous research analyzes the BrainSpan dataset, which contains gene expression of brain tissues from varying regions and developmental periods. Since the spatiotemporal properties of brain tissue is known to affect the gene expression's covariance, previous research have focused only on a specific subset of samples to avoid the issue of heterogeneity. This leads to a potential loss of power when detecting risk genes. In this article, we develop a new method called COBS (COvariance-Based sample Selection) to find a larger and more homogeneous subset of samples that share the same population covariance matrix for the downstream DAWN analysis. To demonstrate COBS's effectiveness, we utilize genetic risk scores from two sequential data freezes obtained in 2014 and 2019. We show COBS improves DAWN's ability to predict risk genes detected in the newer data freeze when utilizing the risk scores of the older data freeze as input.

Figures

Equations

$$I^{(p)}$$

$${\mathbf{X}}_1^{(p)},\dots ,{\mathbf{X}}_{n_p}^{(p)}$$

$$\mathcal{P}=\Big{\{}p~{}:~{}I^{(p)}=1\Big{\}}.$$

$$\mathcal{X}=\bigcup_{p\in\widehat{\mathcal{P}}^{(1)}}\big{\{}\bm{X}^{(p)}_{1},\ldots,\bm{X}_{n_{p}}^{(p)}\big{\}}.$$

$$\widehat{T}=\max_{ij}\big{(}\widehat{t}_{ij}\big{)}\quad\text{where }\widehat{t}_{ij}=\frac{\big{(}\widehat{\sigma}_{X,ij}-\widehat{\sigma}_{Y,ij}\big{)}^{2}}{\widehat{s}_{X,ij}/n_{1}+\widehat{s}_{Y,ij}/n_{2}},\quad i,j\in 1,\ldots,d,$$

$$\widehat{T}^{(b)}=\max_{ij}\big{(}\widehat{t}^{(b)}_{ij}\big{)}\quad\text{where }\widehat{t}^{(b)}_{ij}=\frac{\big{(}\widehat{\sigma}^{(b)}_{X,ij}-\widehat{\sigma}^{(b)}_{Y,ij}\big{)}^{2}}{\widehat{s}_{X,ij}/n_{1}+\widehat{s}_{Y,ij}/n_{2}},\quad i,j\in 1,\ldots,d,$$

$$\text{p-value}=\frac{\big{|}\{b:|\widehat{T}^{(b)}|\geq|\widehat{T}|\}\big{|}}{B}.$$

$$\widehat{T}^{(b)}=\max\Big{\{}\widehat{T}^{(b)}_{\ell}:\ell\in\mathcal{L}\Big{\}}.$$

$$\widehat{T}_{\ell}\geq\text{quantile}\Big{(}\{\widehat{T}^{(1)},\ldots,\widehat{T}^{(b)}\};1-\alpha\Big{)}.$$

$$\mathbb{P}\Big{(}\text{no true null hypothesis among }\mathcal{H}\text{ null hypotheses are rejected}\Big{)}\geq 1-\alpha+o(1)$$

$$V=\big{\{}1,\ldots,r\big{\}},\quad E=\big{\{}(i,j):H_{0,(i,j)}\text{ is accepted by the Stepdown method}\big{\}}.$$

$$G^{\prime }\subseteq G\Rightarrow |\mathcal{A}(G^{\prime })|\le |\mathcal{A}(G)|,\text{for any two graphs }G,G^{\prime }.$$

$$\alpha \ge {\alpha }^{\prime }\Rightarrow G\subseteq G^{\prime },$$

$$\displaystyle\text{True positive rate (TPR) for hypotheses}=\frac{\big{|}\big{\{}(i,j)\in\mathcal{L}:i\leq r_{1}\text{ and }j\leq r_{1}\big{\}}\big{|}}{{r_{1}\choose 2}},$$

$$\displaystyle\text{False positive rate (FDR) for hypotheses}=\frac{\big{|}\big{\{}(i,j)\in\mathcal{L}:i>r_{1}\text{ or }j>r_{1}\big{\}}\big{|}}{{r\choose 2}-{r_{1}\choose 2}}.$$

$$\displaystyle=\frac{\big{|}\big{\{}p\in\widehat{\mathcal{P}}:p\leq r_{1}\big{\}}\big{|}}{r_{1}},$$

$$\displaystyle=\frac{\big{|}\big{\{}p\in\widehat{\mathcal{P}}:p>r_{1}\big{\}}\big{|}}{r_{2}+r_{3}}.$$

$${\hat{T}}_{(i,k)}^{(b)}\le {\hat{T}}_{(i,j)}^{(b)}+{\hat{T}}_{(j,k)}^{(b)}.$$

$$\widehat{T}=\max_{ij}\big{(}\widehat{t}_{ij}\big{)}\quad\text{where }\widehat{t}_{ij}=\big{|}\widehat{\sigma}_{X,ij}-\widehat{\sigma}_{Y,ij}\big{|},\quad i,j\in 1,\ldots,d,$$

$$p(\bm{x})=\frac{1}{(2\pi)^{d/2}|\bm{\Sigma}|^{1/2}}\exp\Big{\{}-\frac{1}{2}\big{(}f(\bm{x})-\bm{\mu}\big{)}^{\top}\bm{\Sigma}\big{(}f(\bm{x})-\bm{\mu}\big{)}\Big{\}}\prod_{j=1}^{d}\big{|}f^{\prime}_{j}(x_{j})\big{|},$$

$${\Sigma }_{ij}^{(1)}=\{\begin{array}{ll}1& \text{if }i=j,\\ a& \text{if }i\ne j,i\text{ is in the same cluster as }j,\\ b& \text{if }i\ne j,i\text{ is not in the same cluster as }j,\end{array}$$

$$a=0.9-\beta \cdot 0.4,\text{and}b=0.1+\beta \cdot \mathrm{0.4.}$$

$$\Phi\big{(}t;N(0,\Sigma_{jj})\big{)}=\Phi\big{(}f^{-1}_{j}(t);\widehat{p}_{g_{j}}\big{)},\quad\forall t\in\mathbb{R}.$$

$$\bm{x}=\big{(}x_{1},\ldots,x_{d}\big{)}=\big{(}f^{-1}_{1}(z_{1}),\ldots,f^{-1}_{d}(z_{d})\big{)}$$

$$Z_j\sim \mathbb{P}(I_j=0)N(0,{\sigma }^2)+\mathbb{P}(I_j=1)N(\mu ,{\sigma }^2),$$

$$\mathbb{P}(\bm{I}=\bm{\eta})\propto\exp\Big{(}b\cdot\sum_{j=1}^{d}\eta_{j}+c\cdot\bm{\eta}^{T}\widehat{\bm{A}}\bm{\eta}\Big{)},$$

Full text

Keywords: Bootstrap covariance test, Microarray, Multiple testing with dependence

1 Introduction

The genetic cause of autism spectrum disorder (ASD), a neurodevelopmental disorder that affects roughly 1-2% individuals in the United States, remains an open problem despite decades of research (Autism and Investigators,, 2014). ASD is characterized primarily by impaired social functions and repetitive behavior (Kanner et al.,, 1943; Rutter,, 1978). To better understand this disorder, scientists identify specific genes that are liable for increasing the chance of developing ASD when damaged or mutated (Sanders et al.,, 2015). These are genes are called risk genes. While breakthroughs in genomic technologies and the availability of large ASD cohorts have led to the discovery of dozens of risk genes, preliminary studies suggest there are hundreds of risk genes still unidentified (Buxbaum et al.,, 2012). In this work, we build upon the current statistical methodologies to further improve our ability to identify risk genes.

We focus on statistical methods that use gene co-expression networks to help identify risk genes. These networks are estimated from brain tissue’s gene expression data. Since these gene co-expression networks provide insight into genes that regulate normal biological mechanisms in fetal and early brain development, it was hypothesized that risk genes that alter these mechanisms should be clustered in these networks (Šestan et al.,, 2012). Early findings confirmed this hypothesis (Parikshak et al.,, 2013; Willsey et al.,, 2013). These results led to the development of the Detection Association With Networks (DAWN) algorithm which uses a “guilt by association” strategy – implicating new risk genes based on their connectivity to previously identified risk genes (Liu et al.,, 2014, 2015). However, the previous DAWN analyses suffer from statistical limitations that we will investigate and resolve in this article.

We challenge previous analyses’ assumptions regarding the homogeneity of the covariance matrix in gene expression data. Previous DAWN analyses assume that gene expression samples from the same brain tissue type share the same covariance matrix. This assumption was influenced by the findings in Kang et al., (2011) and Willsey et al., (2013), which showed that gene co-expression patterns differ among different brain regions and developmental periods on average. Statistically, this means that the covariance matrix among the gene expressions may differ with respect to the spatiotemporal properties of the brain tissue. Hence, previous DAWN analyses (Liu et al.,, 2014, 2015) use only samples from a particular brain tissue type chosen by the findings in Willsey et al., (2013). However, no further statistical analysis is performed to check for homogeneity of this specific subset of samples. In addition, since previous DAWN analyses limit themselves to a subset of gene expression samples, many other samples assumed to be heterogeneous are excluded. This leads to a potential loss of power when identifying risk genes.

To overcome these limitations, we develop a method called COBS (COvariance-Based sample Selection), a two-staged procedure in order to select a subset of gene expression samples in a data-driven way that is more homogeneous and larger in sample size than the fixed subset used previously. In the first stage, we take advantage of the recent developments in high-dimensional covariance testing (Cai et al.,, 2013; Chang et al.,, 2017) to determine whether if the gene expression from two different brain tissues share the same population covariance matrix. We combine this with a multiple-testing method called Stepdown that accounts for the dependencies among many hypothesis tests (Romano and Wolf,, 2005; Chernozhukov et al.,, 2013). In the second stage, after determining which pairs of brain tissues have statistically indistinguishable covariance matrices, we develop a clique-based procedure to select which brain tissues to use in the downstream DAWN analysis. We show that COBS selects brain tissues within the BrainSpan dataset that align with current scientific knowledge and also leads to an improved gene network estimate for implicating risk genes. This article addresses the numerous algorithmic challenges needed to implement this idea.

In Section 2, we describe the data and statistical model for heterogeneity in the covariance matrix. In Section 3, we provide a visual diagnostic to investigate the homogeneity assumptions of previous DAWN analyses. In Section 4, we describe the different stages of COBS to find a subset of homogeneous samples within a dataset. In Section 5, we illustrate the properties of COBS on synthetic datasets. In Section 6, we apply our procedure on gene expression data to show that, when combined with DAWN, we have an improved gene network that can better implicate risk genes. Section 7 provides an overall summary and discussion.

2 Data and model background

Due to the challenge of obtaining and preserving brain tissue, datasets recording the gene expression patterns of brain tissue are rare. The BrainSpan project contributes one of the largest microarray expression datasets available (the “BrainSpan dataset” henceforth), sampling tissues from 57 postmortem brains with no signs of large-scale genomic abnormalities (Kang et al.,, 2011). Many studies have favored this dataset because its 1294 microarray samples capture the spatial and temporal changes in gene expression that occur in the brain during the entirety of development (De Rubeis et al.,, 2014; Dong et al.,, 2014; Cotney et al.,, 2015). While our paper focuses on this particular microarray expression dataset, our method would apply to other gene expression datasets such as RNA sequencing data.

The heterogeneity of gene expression due to the spatiotemporal differences in brain tissues presents statistical challenges. As documented in Kang et al., (2011), the region and developmental period of the originating brain tissue contribute more to the heterogeneity than other variables such as sex and ethnicity. To understand this heterogeneity, we use the following schema to model the BrainSpan dataset. Each of the 1294 microarray samples is categorized into one of 16 spatiotemporal windows, or windows for short, depending on which brain region and developmental period the brain tissue is derived from. Within each window, all microarray samples originating from the same brain are further categorized into the same partition. There are 212 partitions in total. Figure 1 summarizes how many partitions and microarray samples belong in each window in the BrainSpan dataset. This schema allows us to model the microarray samples more realistically since the gene co-expression patterns vary greatly on average from window to window (Willsey et al.,, 2013). Additionally, Willsey et al., (2013) find that among all the windows, the known risk genes in Window 1B are most tightly co-expressed. Window 1B is highlighted in Figure 1 and contains the 107 microarray samples from the prefrontal cortex and primary motor-somatosensory cortex from 10 to 19 post-conceptual weeks. Due to this finding, previous DAWN analyses focus on all 107 samples from 10 partitions, assuming that these samples were homogeneous without further statistical investigation, and discard the remaining 1187 samples, (Liu et al.,, 2014, 2015). We seek to improve upon this heuristical sample selection procedure, first by formalizing a statistical model.

2.1 Statistical model

We use a mixture model that assumes that microarry samples from the same partition are homogeneous while samples from different partitions could be heterogeneous. For the $p$th partition, let $\bm{X}^{(p)}_{1},\ldots\bm{X}^{(p)}_{n_{p}}\in\mathbb{R}^{d}$ denote $n_{p}$ i.i.d. samples, and let $w(p)$ denote the window that partition $p$ resides in. These $n_{p}$ samples are drawn from either a distribution with covariance $\bm{\Sigma}$, or another distribution with a different covariance matrix $\bm{\Sigma}_{p}$. Our notation emphasizes that the distributions in consideration are not necessarily Gaussian, and $\bm{\Sigma}$ is the covariance matrix shared among all partitions, while $\bm{\Sigma}_{p}$ may vary from partition to partition. A fixed but unknown parameter $\gamma_{w(p)}\in[0,1]$ controls how frequently the partitions in window $w$ are drawn from these two distributions, meaning it controls the amount of heterogeneity. For each partition $p$, this mixture model is succinctly described as,

[TABLE]

where $D(\bm{\Sigma})$ denotes an arbitrary distribution with covariance matrix $\bm{\Sigma}$, and $I^{(p)}$ is the latent variable that determines whether or not the samples in partition $p$ have covariance $\bm{\Sigma}$ or $\bm{\Sigma}_{p}$. With this model setup, our task is to determine the set of partitions that originate from the covariance matrix $\bm{\Sigma}$, which we will call

[TABLE]

The findings of Kang et al., (2011) and Willsey et al., (2013) inform us on how much heterogeneity to expect within a window via $\gamma_{w(p)}$. While analyses such as Liu et al., (2015) assume that all the samples in Window 1B are homogeneous, it is noted in Kang et al., (2011) that sampling variability in brain dissection and in the proportion of white and gray matter in different brain tissues can cause variability in the gene co-expression patterns. This means that scientifically, we do not expect all the partitions in Window 1B to be homogeneous (i.e., $\gamma_{w(p)}=1$). Furthermore, Willsey et al., (2013) find a hierarchical clustering among the four brain regions. This is illustrated in Figure 1, where the gene co-expression patterns in the brain regions represented in first row are most similar to those in the second row and least similar to those in the fourth row. The authors also find a smooth continuum of gene expression patterns across different developmental periods, represented as the columns of the table in Figure 1. Hence, we expect $\gamma_{w(p)}$ to decrease smoothly as the window $w$ becomes more dissimilar to Window 1B, in both the spatial and temporal direction.

2.2 Connections to other work

Other work use models similar to (2.1) on microarray expression data to tackle the different co-expression patterns among different tissues and subjects, but their methods differ from ours. One direction is to directly cluster the covariance matrices of each partition (Ieva et al.,, 2016). However, this approach does not account for the variability in the empirical covariance matrix, unlike our hypothesis-testing based method. Another approach is to explicitly model the population covariance matrix for each partition as the summation of a shared component and a partition-specific heterogeneous component. This is commonly used in batch-correction procedures where the analysis removes the heterogeneous component from each partition (Leek and Storey,, 2007). However, we feel such an additive model is too restrictive for analyzing the BrainSpan dataset, as we do not believe there is a shared covariance matrix across all windows of the brain. Instead, our approach will find specific set of partitions with statistically indistinguishable covariance matrices.

3 Elementary analysis

In this section, we develop a visual diagnostic to investigate if the 10 partitions in Window 1B used in previous work (Liu et al.,, 2014, 2015) are as homogeneous as these previous analyses assume. Using a hypothesis test for equal covariances, our diagnostic leverages the following idea: we divide the partitions into two groups and apply a hypothesis test to the samples between both groups. If all the partitions were truly drawn from distributions with equal covariances, then over many possible divisions, the empirical distribution of the resulting p-values should be roughly uniform. We can visualize this distribution by using a QQ-plot. The less uniform the p-values look, the less we are inclined to interpret our partitions to be all drawn from distributions with equal covariances. The following algorithm summarizes this diagnostic.

Algorithm 1: Covariance homogeneity diagnostic

Loop over trials $t=1,2,\ldots,T$:

1. (a)

Randomly divide the selected partitions in the set $\widehat{\mathcal{P}}$ into two sets, $\widehat{\mathcal{P}}^{(1)}$ and $\widehat{\mathcal{P}}^{(2)}$, such that $\widehat{\mathcal{P}}^{(1)}\cup\widehat{\mathcal{P}}^{(2)}=\widehat{\mathcal{P}}$ and $\widehat{\mathcal{P}}^{(1)}\cap\widehat{\mathcal{P}}^{(2)}=\emptyset$. 2. (b)

For each partition $p\in\widehat{\mathcal{P}}^{(1)}$, center the samples $\bm{X}^{(p)}_{1},\ldots,\bm{X}^{(p)}_{n_{p}}$. Then aggregate all samples in $\widehat{\mathcal{P}}^{(1)}$ to form the set of samples

[TABLE]

Similarly, form the set of samples $\mathcal{Y}$ from the set of partitions $\widehat{\mathcal{P}}^{(2)}$. 3. (c)

Compute the p-value for a hypothesis test that tests whether or not the samples in $\mathcal{X}$ and $\mathcal{Y}$ have the same covariance matrix. 2. 2.

Produce a QQ-plot of the resulting $T$ p-values to see if empirical distribution of the p-values is close to a uniform distribution.

We remind the reader that the above procedure is a diagnostic. This is not necessarily a recipe for a goodness-of-fit test since the $T$ p-values are not independent, which makes it difficult to analyze its theoretical properties without a carefully designed global null test. However, as we will demonstrate in later sections of this article, this diagnostic is nonetheless able to display large-scale patterns in our dataset.

3.1 Specification of covariance hypothesis test

To complete the above diagnostic’s description, we describe the procedure to test for equality of covariance matrices. Following the model (2.1), let $\mathcal{X}=\{\bm{X}_{1},\ldots,\bm{X}_{n_{1}}\}$ and $\mathcal{Y}=\{\bm{Y}_{1},\ldots,\bm{Y}_{n_{2}}\}$ be $n_{1}$ and $n_{2}$ i.i.d. samples from $d$-dimensional distribution with covariance $\bm{\Sigma}_{X}$ and $\bm{\Sigma}_{Y}$ respectively, both with an empirical mean of $\mathbf{0}$. We define $\mathbb{X}\in\mathbb{R}^{n_{1}\times d}$ and $\mathbb{Y}\in\mathbb{R}^{n_{2}\times d}$ as the matrices formed by concatenating these samples row-wise. Define the empirical covariance matrices as $\widehat{\bm{\Sigma}}_{X}=\mathbb{X}^{\top}\mathbb{X}/n_{1}$, and $\widehat{\bm{\Sigma}}_{Y}=\mathbb{Y}^{\top}\mathbb{Y}/n_{2}$, where we denote the individual elements of these matrices as $\widehat{\bm{\Sigma}}_{X}=[\widehat{\sigma}_{X,ij}]_{1\leq i,j\leq d}$ and likewise for $\widehat{\bm{\Sigma}}_{Y}$.

We now discuss the hypothesis test for equal covariance, $H_{0}:\bm{\Sigma}_{X}=\bm{\Sigma}_{Y}$, that we will consider in this article based on the test statistic defined in Chang et al., (2017) which extends Cai et al., (2013). In these works, the authors note that if $\bm{\Sigma}_{X}=\bm{\Sigma}_{Y}$, then the maximum element-wise difference between $\bm{\Sigma}_{X}$ and $\bm{\Sigma}_{Y}$ is 0. Hence, Chang et al., (2017) defines the test statistic $\widehat{T}$ as the maximum of squared element-wise differences between $\widehat{\bm{\Sigma}}_{X}$ and $\widehat{\bm{\Sigma}}_{Y}$, normalized by its variance. Specifically,

[TABLE]

where $\widehat{s}_{X,ij}=\sum_{m=1}^{n_{1}}(\mathbb{X}_{mi}\mathbb{X}_{mj}-\widehat{\sigma}_{X,ij})^{2}/n_{1}$ is the empirical variance of the variance-estimator $\widehat{\sigma}_{X,ij}$, and $\widehat{s}_{Y,ij}$ is defined similarly.

Then, Chang et al., (2017) constructs an empirical null distribution of $\widehat{T}$ under $H_{0}:\bm{\Sigma}_{X}=\bm{\Sigma}_{Y}$ using the multiplier bootstrap (Chernozhukov et al.,, 2013). On each of the $b\in\{1,\ldots,B\}$ trials, the multiplier bootstrap computes a bootstrapped test statistic $\widehat{T}^{(b)}$ by weighting each of the $n_{1}+n_{2}$ observations by a standard Gaussian random variable drawn independently of all other variables, denoted collectively as $(g^{(b)}_{1},\ldots,g_{n_{1}}^{(b)},g_{n_{1}+1}^{(b)},\ldots,g_{n_{1}+n_{2}}^{(b)})$. Specifically, we construct the bootstrap statistic for the $b$th trial as

[TABLE]

where $\widehat{\sigma}^{(b)}_{X,ij}=\sum_{m=1}^{n_{1}}g_{m}^{(b)}(\mathbb{X}_{mi}\mathbb{X}_{mj}-\widehat{\sigma}_{X,ij})/n_{1}$ and $\widehat{\sigma}^{(b)}_{Y,ij}=\sum_{m=1}^{n_{2}}g_{n_{1}+m}^{(b)}(\mathbb{Y}_{mi}\mathbb{Y}_{mj}-\widehat{\sigma}_{Y,ij})/n_{2}$. We compute the p-value by counting the proportion of bootstrap statistics that are larger than the test statistic,

[TABLE]

Chang et al., (2017) prove that this test has asymptotically $1-\alpha$ coverage under the null hypothesis as long as the all distributions in the distribution family $D$ in (2.1) have sub-Gaussian and sub-exponential tails, even in the high-dimensional regime where $d\gg\max(n_{1},n_{2})$.

3.2 Application to BrainSpan

Equipped with a complete description of the diagnostic, we apply it to the BrainSpan dataset. Among the 10 partitions in Window 1B, we divide the partitions into two groups uniformly at random 250 times, and compute a p-value using Method 1 (with normalization) for each division using 200 bootstrap trials. The QQ-plot of the resulting p-values are shown in Figure 2A, where we see that the p-values are biased towards 0. This implies the 10 partitions in Window 1B are heterogeneous since they do not seem to all share the same covariance matrix. Furthermore, we apply this diagnostic to all partitions in the BrainSpan dataset with 5 or more samples. This results in using only 125 of the 212 partitions shown in Figure 1. The resulting p-values become more biased towards 0 (Figure 2B), implying the dataset as a whole is more heterogeneous than the partitions in Window 1B. In the next section, we develop a method to resolve this issue by finding the largest subset of partitions possible among the 125 partitions in the BrainSpan dataset that share the same covariance matrix.

4 Methods: COBS (Covariance-based sample selection)

While we have discussed a method to test for equivalent covariance matrices between any two partitions in Section 3, we cannot directly apply this method to select a large number of homogeneous partitions in the BrainSpan dataset without suffering a loss of power due to multiple testing. Since there are $r=125$ partitions with more than 5 samples, applying the hypothesis test to each pair of partitions results in ${r\choose 2}=7750$ dependent p-values. These p-values are dependent since each of the $r$ partitions is involved in $r-1$ hypothesis tests. Hence, standard techniques such as a Bonferroni correction are too conservative when accounting for these dependencies, likely leading to a loss of power.

To properly account for this dependency, we introduce our new method called COBS, which comprises of two parts. First, we use a Stepdown method in Subsection 4.1 that simultaneously tests all ${r\choose 2}$ hypothesis tests for equal covariance matrices, which builds upon the bootstrap test introduced previously in Section 3. After determining which of the ${r\choose 2}$ pairs of partitions do not have statistically significant differences in their covariance matrices, we develop a clique-based method in Subsection 4.2 to select a specific set of partitions $\widehat{\mathcal{P}}$.

4.1 Stepdown method: multiple testing with dependence

We use a Stepdown method developed in Chernozhukov et al., (2013) to control the family-wise error rate (FWER). We tailor the bootstrap-based test in Subsection 3.1 to our specific setting in the algorithm below. We denote $\widehat{T}_{(i,j)}$ as the test statistic formed using (3.1) to test if the covariance of samples between partition $i$ and partition $j$ are equal. Similarly, let $\widehat{T}^{(b)}_{(i,j)}$ denote the corresponding bootstrap statistics on the $b$th bootstrap trial. Here, $\text{quantile}(\{x_{1},\ldots,x_{n}\};1-\alpha)$ represents the empirical $(1-\alpha)\cdot 100$% quantile of the vector $(x_{1},\ldots,x_{n})$.

Algorithm 2: Stepdown method

Initialize the list enumerating all ${r\choose 2}$ null hypotheses corresponding to the set of partition pairs, \mathcal{L}=\Big{\{}(1,2),\ldots,(r-1,r)\Big{\}}. 2. 2.

Calculate $\widehat{T}_{\ell}$ for each $\ell\in\mathcal{L}$, as stated in (3.1). 3. 3.

Loop over steps $t=1,2,\ldots$:

1. (a)

For each bootstrap trial $b=1,\ldots,B$:

1. i.

Generate $N=\sum_{p}n_{p}$ i.i.d. standard Gaussian random variables, one for each sample in each partition, and compute $\widehat{T}^{(b)}_{\ell}$ for all $\ell\in\mathcal{L}$, as stated in (3.2). 2. ii.

Compute

[TABLE] 2. (b)

Remove any $\ell\in\mathcal{L}$ if

[TABLE]

If no elements are removed from $\mathcal{L}$, return the null hypotheses corresponding to $\mathcal{L}$. Otherwise, continue to step $t+1$.

Using techniques in Romano and Wolf, (2005) and Chernozhukov et al., (2013), it can be proven that this method has the following asymptotic FWER guarantee,

[TABLE]

under the same assumptions posed in Chang et al., (2017). The reason LABEL:label:alg2 is able to control the FWER without a Bonferroni correction is because the null distribution in the Stepdown method is properly calibrated to account for the joint dependence among the ${r\choose 2}$ tests. Specifically, when ${r\choose 2}$ tests are individually performed as in Subsection 3.1, the test statistics (3.1) are dependent, but the bootstrapped null distributions do not account for this dependence. Hence, accounting for the dependence via a Bonferroni correction after-the-fact can lead to a substantial loss in power. However, in the Stepdown procedure, the bootstrapped null distributions retain the dependencies jointly since they are generated from the same $N$ Gaussian random variables in each trial. See Chernozhukov et al., (2013) (Comment 5.2) for a further discussion.

Robustness concerns.

In practice, due to the maximum function in the test statistic $\widehat{T}_{\ell}$ displayed in (3.1), the Stepdown method could possibly erroneously reject a hypothesis due to the presence of outliers. One way to circumvent this problem to purposely shrink the value of the test statistic $\widehat{T}_{\ell}$ while leaving the bootstrapped statistics $\widehat{T}^{(b)}_{\ell}$ in (3.2) the same. Specifically, we can replace $\max_{ij}(\widehat{t}_{ij})$ in (3.1) with the $\text{quantile}(\{\widehat{t}_{ij}\}_{ij};1-\epsilon)$, where $\epsilon$ is a positive number extremely close to 0. This has the desired effect of “discarding” the large values in $\{\widehat{t}_{ij}\}_{ij}$. Observe that this procedure would potentially lead to a slight loss in power, but the inferential guarantee in (4.2) still holds since there can only be strictly less rejections.

Computational concerns.

While we use the test statistics (3.1) when describing the Stepdown method, we note that this method applies to a broader family of test statistics. In Appendix C, we discuss in detail one alternative to the test statistic in (3.1) that can dramatically reduce up the computation complexity of the Stepdown method. However, we defer this to the appendix because in our specific problem setting of testing equality of covariances, it does not seem to perform well empirically.

4.2 Largest quasi-clique: selecting partitions based on testing results

After applying the covariance testing with the Stepdown method described in the previous subsection, we have a subset of null hypotheses from $\mathcal{H}$ that we accepted. In this subsection, we develop a clique-based method to estimate $\mathcal{P}$, the subset of partitions that share the same covariance matrix defined in (2.2), from our accepted null hypotheses.

We conceptualize the task of selecting partitions as selecting vertices from a graph that form a dense subgraph. Let $H_{0,(i,j)}$ denote the null hypothesis that the population covariance matrices for partition $i$ and $j$ are equal. Let $G=(V,E)$ be an undirected graph with vertices $V$ and edge set $E$ such that

[TABLE]

Since each of the ${|\mathcal{P}|\choose 2}$ pairwise tests among the partitions in $\mathcal{P}$ satisfies the null hypotheses, the vertices corresponding to $\mathcal{P}$ would ideally form the largest clique in graph $G$. However, this ideal situation is unlikely to happen. Instead, due to the probabilistic nature of our theoretical guarantee in (4.2), there are likely to be a few missing edges in $G$ among the vertices corresponding to $\mathcal{P}$. Hence, a natural task is to instead find the largest quasi-clique, a task that has been well-studied by the computer science community (see Tsourakakis, (2014) and its references within). We say a set of $k$ vertices form a $\gamma$-quasi-clique if there are at least $\gamma\cdot{k\choose 2}$ edges among these $k$ vertices for some $\gamma\in[0,1]$. The largest $\gamma$-quasi-clique is the largest vertex set that forms a $\gamma$-quasi-clique . We justify the choice to search for this $\gamma$-quasi-clique since, by our model (2.1), the prevalent covariance matrix among the $r$ partitions is the desired covariance matrix $\bm{\Sigma}$ we wish to estimate. Here, $\gamma$ is an additional tuning parameter, but we set $\gamma=0.95$ by default throughout this entire paper.

Unfortunately, many algorithms that could be used to find the largest $\gamma$-quasi-clique do not satisfy a certain monotone property in practice, which hinders their usability. Specifically, consider an algorithm $\mathcal{A}$ that takes in a graph $G$ and outputs a vertex set, denoted by $\mathcal{A}(G)$, and for two graphs $G^{\prime}$ and $G$, let $G^{\prime}\subseteq G$ denote that $G^{\prime}$ is a subgraph of $G$. We say that algorithm $\mathcal{A}$ has the monotone property if

[TABLE]

We are not aware of such a property being important in the quasi-clique literature, but it is a natural property to inherit from the multiple testing community. That is, a multiple testing procedure has the monotone property if increasing the signal-to-noise ratio (i.e., decreasing the p-values) yields more rejections (see (Hahn,, 2018) and references within). Similarly in the quasi-clique setting, it is natural to expect that increasing the signal-to-noise ratio (i.e., removing edges in $G$) yields less partitions selected. The monotone property is crucial in practice since it can be shown that the chosen FWER level $\alpha$ and the graph $G$ defined in (4.3) have the following relation,

[TABLE]

where $G$ and $G^{\prime}$ are the graphs formed by FWER level $\alpha$ and $\alpha^{\prime}$ respectively. Hence, an algorithm that does not exhibit the property in (4.4) will be fragile – using a smaller $\alpha$ to accept more null hypotheses might counterintuitively result in less partitions being selected. As we will demonstrate in Section 5 through simulations, many existing algorithms to find the largest quasi-clique do not satisfy the monotone property empirically. Therefore, we develop the following new algorithm to remedy this.

We describe the algorithm below. It starts by constructing a list containing all maximal cliques in the graph based on (4.3). A maximal clique is a vertex set that forms a clique but is not subset of a larger clique. The algorithm then proceeds by determining if the union of any two vertex sets forms a $\gamma$-quasi-clique. If so, this union of vertices is added to the list of vertex sets. The algorithm returns the largest vertex set in the list when all pairs of vertex sets are tried and no new $\gamma$-quasi-clique is found. We demonstrate in Section 5 that this algorithm exhibits the monotone property (4.4) empirically.

Algorithm 4: Clique-based selection

Form graph $G$ based on (4.3). 2. 2.

Form $\mathcal{Q}$, the set of all vertex sets that form a maximal clique in $G$. Each vertex set is initialized with a child set equal to itself. 3. 3.

While there are vertex sets $A,B\in\mathcal{Q}$ the algorithm has not tried yet:

1. (a)

Determine if $C=A\cup B$ forms a $\gamma$-quasi-clique in $G$. If so, add $C$ as a new vertex set into $\mathcal{Q}$, with $A$ and $B$ as its two children sets. 4. 4.

Return the largest vertex set in $\mathcal{Q}$.

A naive implementation of the above algorithm would require checking if an exponential number of vertex set unions $C=A\cup B$ forms a $\gamma$-quasi-clique, and each check requires $O(r^{2})$ operations. However, we are able to dramatically reduce the number of checks required by using the following heuristic: we only check whether the union of $A$ and $B$ forms a $\gamma$-quasi-clique if the union of two children sets, one from each $A$ and $B$, forms a $\gamma$-quasi-clique. This heuristic allows us to exploit previous calculations and reduce computational costs. We implement this idea by using one hash table to record which vertex sets are children of other vertex sets, and another hash table table to record if the union of two vertex sets forms a $\gamma$-quasi-clique. This idea is illustrated in Figure 4. Additional details on how to initialize and optionally post-process LABEL:label:alg4 are given in Appendix D.

5 Simulation study

We perform empirical studies to show that COBS has more power and yields a better estimation of the desired covariance matrix $\bm{\Sigma}$ over conventional methods as the samples among different partitions are drawn from increasingly dissimilar distributions.

Setup: We generate synthetic data in $r=25$ partitions, where the data in each partition has $n=15$ samples and $d=1000$ dimensions drawn from a non-Gaussian distribution. Among these $r$ partitions, the first group of $r_{1}=15$ partitions, second group of $r_{2}=5$ partitions and third group of $r_{3}=5$ partitions are drawn from three different nonparanormal distributions respectively (Liu et al.,, 2009). The goal in this simulation suite is to detect these $r_{1}$ partitions with the same covariance structure. The nonparanormal distribution is a natural candidate to model genomic data with heavier tails and multiple modes (Liu et al., (2012) and Xue and Zou, (2012)), and serves to demonstrate that our methods in Section 4 does not rely on the Gaussian assumption. Formally, a random vector $\bm{X}=(X_{1},\ldots,X_{d})\in\mathbb{R}^{d}$ is drawn from a nonparanormal distribution if there exists $d$ monotonic and differentiable functions $f_{1},\ldots,f_{d}$ such that when applied marginally, $\bm{Z}=(f_{1}(X_{1}),\ldots,f_{d}(X_{d}))\sim N(\bm{\mu},\bm{\Sigma})$, a Gaussian distribution with proxy mean vector $\bm{\mu}$ and proxy covariance matrix111We emphasize “proxy” covariance matrix, for example, since the covariance of $\bm{X}$, the random variable we sample, is not $\bm{\Sigma}$. $\bm{\Sigma}$. We provide the details of how we generate the three nonparanormal distributions in Appendix E, but we highlight the key features regarding $\bm{\Sigma}$ below.

We construct three different proxy covariance matrices $\bm{\Sigma}^{(1)},\bm{\Sigma}^{(2)},$ and $\bm{\Sigma}^{(3)}$ in such a way that for a given parameter $\beta\in[0,1]$, we construct $\bm{\Sigma}^{(2)}$ and $\bm{\Sigma}^{(3)}$ to be more dissimilar from $\bm{\Sigma}^{(1)}$ as $\beta$ increases. We highlight the key features of our constructed proxy covariance matrices here. All three proxy covariance matrices are all based on a stochastic block model (SBM), a common model used to model gene networks (Liu et al.,, 2018; Funke and Becker,, 2019). The first $r_{1}$ partitions are generated using proxy covariance matrix $\bm{\Sigma}^{(1)}$, which is an SBM with two equally-sized clusters where the within-cluster covariance is $a=0.9$ and the between-cluster covariance is $b=0.1$. The second $r_{2}$ partitions are generated using proxy covariance matrix $\bm{\Sigma}^{(2)}$, which is similar to $\bm{\Sigma}^{(1)}$ except $a$ and $b$ are shrunk towards $0.5$ depending on the magnitude of $\beta$. The last $r_{2}$ partitions are generated using proxy covariance matrix $\bm{\Sigma}^{(3)}$, which is similar to $\bm{\Sigma}^{(1)}$ except an equal fraction of variables from both clusters break off to form a third cluster, depending on the magnitude of $\beta$. By generating $\bm{\Sigma}^{(1)},\bm{\Sigma}^{(2)}$, and $\bm{\Sigma}^{(3)}$ in this fashion, the parameter $\beta$ can control the difficulty of the simulation setting – a larger $\beta$ means COBS would ideally have more power in distinguishing among the first $r_{1}$ partitions from the other partitions. Figure 5 visualizes the resulting covariance matrices for the three nonparanormal distribution we generate in this fashion for $\beta=0.3$ and $\beta=1$.

Multiple testing: We use the Stepdown method described in Subsection 4.1 on our simulated data where $\beta=\{0,0.3,0.6,1\}$ to see how the true positive rates and false positive rates vary with $\beta$. Let $\mathcal{L}=\{(i_{1},j_{1}),(i_{2},j_{2}),\ldots\}$ denote the returned set of partition pairs that correspond to the accepted null hypothesis. Since our goal is to find the first $r_{1}$ partitions, we define the true positive rate and false positive rate for individual hypotheses to be

[TABLE]

We plot the RoC curves visualizing the TPR and FPR in Figure 6. Each curve traces out the mean true and false positive rate over 25 simulations as $\alpha$ ranges from [math] (top-right of each plot) to $1$ (bottom-left of each plot), where we use 200 bootstrap trials per simulation. Figure 6A shows the naive analysis where we compute all ${r\choose 2}$ p-values, one for each hypothesis test comparing two partitions, and accept hypotheses for varying levels of $\alpha$ after using a Bonferroni correction. Figure 6B shows the Stepdown method. In both plots, we see that as $\beta$ increases, each method has more power. However, as we mentioned in Subsection 4.1, there is a considerable loss of power when comparing the Bonferroni correction to the Stepdown method. This is because the Bonferroni correction is too conservative when accounting for dependencies.

Partition selection: After using Stepdown, we proceed to select the partitions as in Subsection 4.2 to understand the monotone property and see how the true and false positive rates for partitions vary with $\beta$.

Figure 7 shows that three methods currently in the literature that can be used to find the largest quasi-clique in (4.3) fail the monotone property (4.4), whereas COBS succeeds. In Figure 7A, we compare our clique-based selection method, described in Subsection 4.2, against spectral clustering, a method used in network analyses designed to find highly connected vertices (Lei and Rinaldo,, 2015), whereas in Figure 7B, two methods recently developed in the computer science community are compared (Chen and Saad, (2010) and Tsourakakis et al., (2013)). These three methods are detailed in Appendix D, and all the methods receive the same set of accepted null hypotheses as the FWER level $\alpha$ varies. Recall that since the Stepdown method accepts more hypotheses as $\alpha$ decreases, the graph formed by (4.3) becomes denser as $\alpha$ increases. However, as we see in Figure 7, the number of partitions selected by all but our method sometimes decreases as number of accepted null hypotheses increases, hence violating the desired monotone property.

Figure 8A shows the RoC curves for varying $\beta$ as the FWER level $\alpha$ varies. This figure is closely related to Figure 6B. We use our clique-based selection method to find the largest $\gamma$-quasi-clique for $\gamma=0.95$. Let $\widehat{\mathcal{P}}$ denote the selected set of partitions. Similar to before, we define the TPR and FPR in this setting as

[TABLE]

We see that the power of the COBS increases as $\beta$ increases, as expected.

Covariance estimation: Finally, we show that COBS is able to improve the downstream covariance estimation compared to other approaches. To do this, we use four different methods to select partitions and compute the empirical covariance matrix among the samples in those partitions. The first three methods resemble analyses that could be performed on the BrainSpan dataset in practice. The first method uses the COBS. The second method always selects all the partitions, which resembles using all the partitions in the BrainSpan dataset. The third method always selects the same 5 partitions – 3 partitions contain samples drawn from the nonparanormal distribution with proxy covariance $\bm{\Sigma}^{(1)}$, while the other 2 partitions contain samples from each of the remaining two distributions. This resembles previous work (Liu et al.,, 2015) that consider only partitions in Window 1B. For comparison, the last method resembles an oracle that selects exactly the $r_{1}$ partitions containing samples drawn the nonparanormal distribution with proxy covariance $\bm{\Sigma}^{(1)}$.

Figure 8B shows that our partition selection procedure performs almost as well as the oracle method over varying $\beta$ level. Notice that for low $\beta$, COBS and the method using all partitions yield a smaller spectral error than the oracle method. This is because for low $\beta$, the covariance matrices $\bm{\Sigma}^{(1)}$, $\bm{\Sigma}^{(2)}$, and $\bm{\Sigma}^{(3)}$ are almost indistinguishable. However, as $\beta$ increases, the dissimilarities among $\bm{\Sigma}^{(1)}$, $\bm{\Sigma}^{(2)}$, and $\bm{\Sigma}^{(3)}$ grow. This means methods that do not adaptively choose which partitions to select become increasingly worse. However, our procedure remains competitive, performing almost as if it knew which partitions contain samples drawn the nonparanormal distribution with proxy covariance $\bm{\Sigma}^{(1)}$. Additional simulations that go beyond the results in this section are deferred to Appendix F.

6 Application on BrainSpan study

We demonstrate the utility of COBS by applying it within the DAWN framework established in Liu et al., (2015). Specifically, in this section, we ask two questions. First, does COBS select reasonable partitions within the BrainSpan data, given our current scientific understanding outlined in Section 2? Second, does using COBS within the DAWN framework lead to a more meaningful gene co-expression network that can implicate genes using a “guilt-by-association” strategy?

Here, we discuss the different datasets relevant to the analysis in this section. DAWN relies on two types of data to identify risk genes: gene expression data to estimate a gene co-expression network and genetic risk scores to implicate genes associated with ASD. For the former, we use the BrainSpan microarray dataset (Kang et al.,, 2011), which has been the primary focus of this article so far. For the latter, we use the TADA scores published in De Rubeis et al., (2014) which are p-values, one for each gene, resulting from a test for marginal associations with ASD based on rare genetic variations and mutations.222TADA stands for Transmission and De novo association (He et al.,, 2013). For enrichment analysis, we use a third dataset consisting of TADA scores from Satterstrom et al., (2020). We use this third dataset only to assess the quality of our findings, and these TADA scores are derived as in De Rubeis et al., (2014), but include additional data assimilated since 2014. Relying on a later “data freeze,” this 2019 study has greater power to detect risk genes compared to the 2014 study: the two studies report 102 and 33 risk genes, respectively, with FDR cutoff of $10\%$. Additional details of our analysis in this section can be found in Appendix G.

6.1 Gene screening

We first preprocess the BrainSpan data by determining which genes to include in our analysis. This is necessary since there are over 13,939 genes in the BrainSpan dataset, most of which are probably not correlated with any likely risk genes. Including such genes increases the computationally cost and is not informative for our purposes. Hence, we adopt a similar screening procedure as in Liu et al., (2015), which involves first selecting genes with high TADA scores based on De Rubeis et al., (2014), and then selecting all genes with a high Pearson correlation in magnitude with any of the aforementioned genes within the BrainSpan dataset. We select a total of 3,500 genes to be used throughout the remainder of this analysis.

6.2 Partition selection

Motivated by the findings in Willsey et al., (2013), we analyze the BrainSpan dataset using COBS to find many partitions that are homogeneous with most partitions in Window 1B (Figure 1). We use the Stepdown method with 200 bootstrap trials and FWER level $\alpha=0.1$. This simultaneously finds which null hypotheses are accepted among the ${125\choose 2}$ hypotheses tested. Based on these results, we select the partitions that form the largest $\gamma$-quasi-clique for $\gamma=0.95$.

We visualize the results of the Stepdown method in Figure 9, illustrating that COBS finds 24 partitions which have statistically indistinguishable covariance matrices, 7 of which are in Window 1B. We form the graph $G$ based on the accepted null hypotheses, as described in (4.3). Figure 9A shows the full graph with all 125 nodes, while Figure 9B shows the connected component of $G$ as an adjacency matrix. We can see that the 24 partitions we select, which contain 272 microarray samples, correspond to 24 nodes in $G$ that form a dense quasi-clique.

We visualize the proportion of selected partitions per window in the BrainSpan dataset in Figure 10A to demonstrate that our findings are consistent with the findings in Willsey et al., (2013). As mentioned in Section 2, Willsey et al., (2013) find that partitions in Window 1B are mostly homogeneous and are enriched for tightly clustered risk genes. The authors also found that, on average, gene expression varies smoothly across developmental periods, meaning there is greater correlation between the gene expressions belonging to adjacent developmental windows. The authors also estimate a hierarchical clustering among the four brain regions. Indeed, our results match these finding. We select a large proportion of partitions in Window 1B, and the proportion of selected partitions smoothly decreases as the window representing older developmental periods as well as brain regions become more dissimilar to Window 1B.

Lastly, we apply the same diagnostic as in Section 3 to show in Figure 10B that the 272 samples within our 24 selected partitions are much more homogeneous than the 107 samples among the 10 partitions in Window 1B. The p-values we obtain after 250 divisions are much closer to uniform that those shown in Figure 2.

6.3 Overview of DAWN framework

As alluded to in Section 1, DAWN estimates a gene co-expression network using the microarray partitions to boost the power of the TADA scores using a “guilt-by-association” strategy. Figure 11 illustrates this procedure as a flowchart. The first step uses COBS to select 24 partitions from the BrainSpan dataset, as stated in the previous subsection. In the second step, DAWN estimates a Gaussian graphical model via neighborhood selection (Meinshausen and Bühlmann,, 2006) from the 272 samples in these partitions to represent the gene co-expression network. In the third step, DAWN implicates risk genes via a Hidden Markov random field (HMRF) model by combining the Gaussian graphical model with the TADA scores. The details are in Liu et al., (2015), but in short, this procedure assumes a mixture model of the TADA scores between risk genes and non-risk genes, and the probability that a gene is a risk gene depends on the graph structure. An EM algorithm is used to estimate the parameters of this HMRF model, after which a Bayesian FDR procedure (Muller et al.,, 2006) is used on the estimated posterior probabilities of being a risk gene to output the final set of estimated risk genes. The methodology in the second and third step are the same as those in Liu et al., (2015), as we wish to compare only different ways to perform the first step.

6.4 Investigation on gene network and risk genes

In this subsection, we compare the DAWN analysis using the 24 partitions selected by COBS (i.e., the “COBS analysis”) to using the 10 partitions in Window 1B originally used in Liu et al., (2015) (i.e., the “Window 1B analysis”) to show how COBS improves the estimated gene network.

Closeness of genes within co-expression network.

We demonstrate that the 102 genes detected by the newer TADA scores (Satterstrom et al.,, 2020) are roughly 10%-30% closer to the 33 genes detected by the older TADA scores (De Rubeis et al.,, 2014) in the gene network estimated in the COBS analysis than in the Window 1B analysis. This suggests that the COBS analysis estimates a more useful gene network, because when future TADA scores are published after Satterstrom et al., (2020), the next wave of detected risk genes are more likely to also be closer to the original risk genes detected in De Rubeis et al., (2014). We defer the details to Appendix G, but highlight the procedure to derive this result here. Effectively comparing the distances between genes in a network is a difficult problem since the estimated gene networks in the COBS and Window 1B analyses have different number of edges. In addition, current research has suggested that natural candidates such as the shortest-path distance or the commute distance do not accurately capture the graph topology (Alamgir and Von Luxburg, (2012) and Von Luxburg et al., (2014)). Hence, we use two different distance metrics to measure the closeness of two sets of genes that potentially overcome this problem. The first is using the path distance via the minimum spanning tree, and the second is using the Euclidean distance via the graph root embedding (Lei,, 2018). Using either of these metrics lead to the same conclusion.

Enrichment analysis.

We demonstrate that COBS improves DAWN’s ability to predict risk genes based on the newer TADA scores (Satterstrom et al.,, 2020) when utilizing the older TADA scores (De Rubeis et al.,, 2014) as input. Specifically, the COBS analysis and the Window 1B analysis implicate 209 and 249 risk genes respectively an FDR cutoff of 10%, respectively. The risk genes implicated in the COBS analysis have a better enrichment for the 102 genes detected using the newer TADA scores (Satterstrom et al.,, 2020): 18.8% (COBS analysis) versus 14.6% (Window 1B analysis). We note that genes implicated by DAWN but not by the TADA scores are not considered false positives. In fact, He et al., (2013) suggests that there are upwards of 500 to 1000 genes that increase risk for ASD. Hence, we are unlikely to detect all of the true risk genes based on tests that rely on rare genetic variation alone.

Robustness to $\gamma$.

We additionally verify the robustness of the above enrichment results to the parameter $\gamma$. Recall that $\gamma$ controls the density of the edges in the quasi-clique, as introduced in Subsection 4.2, and we typically set $\gamma=0.95$ by default. When we re-run the entire analysis with different values of $\gamma$ varying from 0.85 to 0.97 at intervals of 0.01, we obtain 13 different sets of estimated risk genes. We stop at $\gamma=0.97$ since larger values result in no partitions selected outside of Window 1B. When we intersect all 13 sets of risk genes together, we find that 144 risk genes are implicated regardless of the value of $\gamma$, of which 22.9% are in the list of 102 risk genes found using only the newer TADA scores (Satterstrom et al.,, 2020). This is a promising result, as it demonstrates that the implicated risk genes in the COBS analysis are more enriched than those in the Window 1B analysis for a wide range of $\gamma$.

7 Conclusion and discussions

In this article, we develop COBS to select many partitions with statistically indistinguishable covariance matrices in order to better estimate graphical models for ASD risk gene detection. Our procedure first applies a Stepdown method to simultaneously test all ${r\choose 2}$ hypotheses, each testing whether or not a pair of partitions share the same population covariance matrix. The Stepdown method is critical since it can account for the dependencies among all ${r\choose 2}$ hypotheses via bootstrapping the joint null distribution. Then, our procedure uses a clique-based selection method to select the partitions based on the accepted null hypotheses. The novelty in this latter method is its ability to preserve monotonicity, a property stating that less partitions should be selected as the number of accepted null hypotheses is smaller. We demonstrate empirically that the COBS achieves this property while common methods such as spectral clustering do not. When we apply COBS to the BrainSpan dataset, we find scientifically meaningful partitions based on the results in Willsey et al., (2013). We also find that COBS aids in clustering the risk genes detected in Satterstrom et al., (2020) closer to the risk genes detected in (De Rubeis et al.,, 2014) within the estimated gene co-expression network and in getting a better enrichment in implicated risk genes via the DAWN analysis.

The theoretical role of the FWER level $\alpha$ is not well understood mathematically. Specifically, while (4.2) provides a theoretical guarantee about the set of null hypothesis accepted, we would like to prove a theoretical guarantee about the set of selected partitions $\widehat{\mathcal{P}}$. Towards this end, we suspect that with some modification to COBS, closed testing offers a promising theoretical framework (see Dobriban, (2018) and references within). This will be investigated in future work.

COBS is applied directly to help implicate risk genes for ASD, but this line of work has broader implications in genetics. Due to the improvement of high throughput technologies, it has become increasingly accessible to gather large amounts of gene expression data. This includes both microarray and RNA sequencing data. However, as we have seen in this article, gene expression patterns can vary wildly among different tissues. Hence, it is challenging to select samples that are relevant for specific scientific tasks. Beyond analyzing brain tissues, Greene et al., (2015) develop procedures to select relevant samples amongst a corpus of microarray expression data to estimate gene co-expression networks for different tissue types. While Greene et al., (2015) does not motivate their method from a statistical model, our work provides a possible statistical direction for this research field to move towards.

Acknowledgments: We thank Bernie Devlin and Lambertus Klei for the insightful discussions about our analysis and results. We thank Li Liu and Ercument Cicek for providing the code used in Liu et al., (2015) to build off of. We also thank the anonymous reviewers for helpful suggestions on how to restructure the simulations and analyses.

Appendix A Code and dataset

The R code for replicating all analyses and figures in this article are hosted on GitHub in the repository https://github.com/linnylin92/covarianceSelection. The three major datasets used in this article are also included in the repository. The first dataset is the BrainSpan microarray samples collected by (Kang et al.,, 2011). While the original dataset is publicly available on GEO (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE25219), we provide a locally preprocessed dataset, which was created to be amendable for our analysis in R. The second dataset is the older TADA scores (De Rubeis et al.,, 2014). The third dataset is the list of 102 risk genes detected using the newer TADA scores (Satterstrom et al.,, 2020).

Appendix B Brain region details

There are four primary brain regions, each containing smaller subregions.

• •

PFC-MSC: The prefrontal cortex (PFC) and primary motor-somatosensory cortex (MSC) consist of six smaller regions: primary motor cortex, primary somatosensory cortex, ventral prefrontal cortex, medial prefrontal cortex, dorsal prefrontal cortex and orbital prefrontal cortex.

• •

V1C, ITC, IPC, A1C, STC: A region consisting of the primary visual cortex (V1C), inferior temporal cortex (ITC), primary auditory cortex (A1C), and superior temporal cortex (STC).

• •

STR, HIP, AMY: A region consisting of the stratum (STR), hippocampal anlage or hippocampus (HIP) and amygdala (AMY).

• •

MD, CBC: A region consisting of the mediodorsal nucleus of the thalamus (MD) and the cerebellar cortex (CD).

Appendix C Extension to the Stepdown method

One of the largest drawbacks of the Stepdown method lies in its intensive computational cost. For $r$ partitions, at most ${r\choose 2}$ bootstrap statistics need to be computed in each bootstrap trial, each requiring a computational cost of $O(d^{2}\cdot n_{p})$. In this section, we develop a computational extension to the Stepdown method that yields a more computationally efficient algorithm as long as the test statistic $\widehat{T}$ satisfies the triangle inequality and the number of variables $d$ is large. That is, for any bootstrap trial $b$ and for any partitions $i$, $j$ and $k$, we require that the bootstrap statistics satisfy

[TABLE]

This property can potentially save expensive calculations when calculating (4.1) in LABEL:label:alg2 by reducing the number of bootstrap statistics we need to explicitly calculate. Since we only care about the maximum bootstrap statistic $\widehat{T}^{(b)}$ in each trial, the triangle inequality gives an upper bound on the bootstrap statistic $\widehat{T}^{(b)}_{(i,k)}$ between partitions $i$ and $k$, leveraging bootstrap statistics already calculated within a specific bootstrap trial. As we sequentially iterate through all pairs of partitions $(i,k)$, if the upper bound for $\widehat{T}^{(b)}_{(i,k)}$ is smaller than the current maximum bootstrap statistic within a specific bootstrap trial $b$, we do not need to explicitly compute $\widehat{T}^{(b)}_{(i,k)}$.

Unfortunately, the test statistic (3.1) described in Subsection 3.1 originally from Chang et al., (2017) does not satisfy the triangle inequality (C.1). Hence, we consider a new test statistic defined as

[TABLE]

and we make a similar modification for its bootstrap counterpart, $\widehat{T}^{(b)}$. It can easily be shown that the above bootstrap statistics satisfies the desired triangle inequality. Additionally, using the techniques in Chernozhukov et al., (2013), it can be proven that this test statistic will still yield a hypothesis test with asymptotic $1-\alpha$ coverage under the null, analogous to (4.2). We will call the Stepdown procedure that uses (C.2) the “Accelerated Stepdown” procedure.

To formalize how to take advantage of this triangle inequality property, we describe a subroutine that leverages this property to compute $\widehat{T}^{(b)}$ in (4.1) by representing the individual bootstrap statistics $\widehat{T}^{(b)}_{(i,j)}$ as weighted edges in a graph. The algorithm uses Dijsktra’s algorithm to find the shortest path between vertices. This implicitly computes the upper bound in the bootstrap statistic between two partitions using the triangle inequality. This algorithm can provide substantial improvement in computational speed by leveraging the fact that determining the shortest path on a fully-dense graph has a computational complexity of $O(r^{2})$, whereas computing $T^{(b)}_{(i,j)}$ has a computational complexity of $O(d^{2}\cdot n_{p})$.

Algorithm 3: Distance metric-based procedure to compute $\widehat{T}^{(b)}$

Form graph $G=(V,E)$ with $r$ nodes and all ${r\choose 2}$ edges, and initialize each edge to have a weight of infinity. 2. 2.

Arbitrarily construct a spanning tree $\mathcal{T}$ and compute all $\widehat{T}^{(b)}_{(i,j)}$’s corresponding to edges $(i,j)\in\mathcal{T}$. Record $z=\max_{(i,j)\in\mathcal{T}}\widehat{T}^{(b)}_{(i,j)}$. 3. 3.

Construct a set of edges $\mathcal{S}=\mathcal{L}\backslash\mathcal{T}$ which represents the bootstrap statistics between specific pairs of partitions that have yet to be computed. 4. 4.

While $\mathcal{S}$ is not empty:

1. (a)

Arbitrarily select an edge $(i,j)\in\mathcal{S}$ and remove it from $\mathcal{S}$. Compute the shortest-path distance from vertex $i$ to $j$ in $G$. 2. (b)

If the shortest-path distance is larger than $z$, update the edge $(i,j)$ to have weight $\widehat{T}^{(b)}_{(i,j)}$, and update $z$ to be $\max(z,\widehat{T}^{(b)}_{(i,j)})$. 5. 5.

Return $z$.

As we will see in Section F, while (C.2) can take advantage of this computational speedup, it yields a much less powerful test when compared to test using (3.1). This is intuitive, as (C.2) does not normalize by the sum of the empirical variances, unlike (3.1). Hence, we do not use the Accelerated Stepdown procedure within COBS when analyzing the BrainSpan dataset in this paper. However, we believe there are potentially other settings outside of covariance testing where this computational speedup idea can be utilized more effectively. We leave this as direction for future work.

Appendix D Details of algorithms to find quasi-cliques

The first subsection remarks on possible extensions to the clique-based selection method described in LABEL:label:alg4. The second subsection describes the three other algorithm used in Section 5 for us to compare against. Throughout this section, for a generic graph $G$, we use $V$ to denote the set of vertices in $G$, $G_{S}$ to denote a subgraph formed by a vertex set $S\subseteq V$, and $E(G)$ to denote the number of edges in $G$.

D.1 Extensions to clique-based selection method

We mention two extensions to clique-based selection method (LABEL:label:alg4) that can be useful in practice.

• •

Initializing algorithm around a desired set of vertices: In certain cases, the user would want the $\gamma$-quasi-clique to be initialized around a desired subset of vertices in $G=(V,E)$. For instance, in our setting, since Liu et al., (2015) applies DAWN to the 10 partitions in Window 1B, it is natural for us to encourage COBS to select as many partitions in Window 1B as possible to enable a meaningful comparison.

To resolve this, first, we run LABEL:label:alg4 at the desired level $\gamma$ on $G_{S}$. This would output a subset of vertices $S_{\text{core}}\subseteq S$ that form the largest $\gamma$-quasi-clique in $G_{S}$. Then, we run LABEL:label:alg4 at the same level $\gamma$ on the full graph $G$ but perform an additional operation after (2.): after $\mathcal{Q}$ is initialized with all maximal cliques in $G$, we check each vertex set $A\in\mathcal{Q}$ if $A\cup S_{\text{core}}$ forms a $\gamma$-quasi-clique. If yes, we replace $A$ with $A\cup S_{\text{core}}$ in $\mathcal{Q}$. If not, we remove $A$ from $\mathcal{Q}$. The algorithm then proceeds to (3.) as usual. By applying this simple change, we are ensured the returned vertex set by LABEL:label:alg4 contains $S_{\text{core}}$.

• •

Post-processing the returned vertex set: In certain cases, the returned vertex set of LABEL:label:alg4 has a few vertices with a very low degree when compared to the other vertices. To resolve this, we post-process this vertex set by removing vertices that are connected to less than half the other vertices in the returned set.

In this paper, we use the initialization extension only when analyzing the BrainSpan dataset in Section 6, where we initialize the largest quasi-clique around the 10 partitions in Window 1B.

D.2 Overview of other algorithms

We overview the three algorithms introduced in Section 5 that are designed to find large quasi-cliques.

• •

Chen and Saad, (2010): This algorithm recursively splits a graph $G$ into two in a hierarchical-clustering type approach with respect to a carefully constructed weight matrix. This forms a tree-type data structure, and then the algorithm scans the tree in a breath-first-search type fashion for the largest subgraph with an edge density larger than $\gamma$.

• •

Tsourakakis et al., (2013): This algorithm performs a local search by adding vertices mypoically and then removing vertices occasionally until no more myopic improvements can be made. Specifically, it first initializes the set $S$ of vertices to contain a vertex that maximizes the ratio between the number of triangles and the degree, and includes all of the neighbors of said vertex. Then algorithm iteratively tries to incrementally improve the $f_{\gamma}(S)=E(G_{S})-\gamma{|S|\choose 2}$ as much as possible by adding neighbors of $S$. When it is no longer able to improve $f_{\gamma}(S)$, the algorithm tries removing a vertex from $S$ to improve $f_{\gamma}(S)$. The algorithm then iterates between such adding and removing vertices from $S$ for a fixed number of iterations.

• •

Spectral clustering: While many different community detection methods for random graphs now exist (for example, see Abbe, (2017) and Athreya et al., (2017) and the references within), we choose spectral clustering as described in Lei and Rinaldo, (2015) as a prototypical example of how many of such methods fail to demonstrate the monotone property as described in Subsection 4.2. Specifically, this method applies K-means clustering to the top $K$ eigenvectors of the adjacency matrix, where $K$ is a tuning parameter to specify. To find large quasi-cliques, we iteratively try spectral clustering for a range of $K$’s (i.e., $K=2,\ldots,5$), and for each detected cluster in any of the estimated clusterings, we compute if the corresponding vertices of said cluster forms a $\gamma$-quasi-clique. If any $\gamma$-quasi-clique is found, we return the largest $\gamma$-quasi-clique discovered in this fashion.

Appendix E Formal description of simulation setup

We say a multivariate vector $\bm{X}\in\mathbb{R}^{d}$ is distributed based a nonparanormal distribution with proxy mean vector $\bm{\mu}$, proxy covariance matrix $\bm{\Sigma}$, and monotonic and differentiable functions $f_{1},\ldots,f_{d}$ if the density of $\bm{X}$ is

[TABLE]

where $f(\bm{x})=(f_{1}(x_{1}),\ldots,f_{d}(x_{d}))$. This is defined in Liu et al., (2009). In our simulation suite, we set $\bm{\mu}=0$. Let this distribution be denoted as $\text{NPN}(\mathbf{0},\bm{\Sigma},f)$. In the next two subsections, we formalize the details of $\bm{\Sigma}$ and $f_{1},\ldots,f_{d}$.

E.1 Details on proxy covariance matrices $\mathbf{\Sigma}$

The following three bullet points detail the construction of $\bm{\Sigma}^{(1)}$, $\bm{\Sigma}^{(2)}$ and $\bm{\Sigma}^{(3)}$ respectively. As mentioned in Section 5, $\beta\in[0,1]$ is a user-defined parameter that controls the dissimilarity among these three matrices.

• •

Construction of $\bm{\Sigma}^{(1)}$: As mentioned in Section 5, $\bm{\Sigma}^{(1)}\in\mathbb{R}_{+}^{d\times d}$ follows an SBM with two equally-sized clusters. Specifically, the first cluster contains indices $1,\ldots,\lfloor d/2\rfloor$ and the second cluster contains indices $\lfloor d/2\rfloor+1,\ldots,d$. Then, we construct $\bm{\Sigma}^{(1)}$ where

[TABLE]

for all $i,j\in 1,\ldots,d$ and $a=0.9$ and $b=0.1$.

• •

Construction of $\bm{\Sigma}^{(2)}$: $\bm{\Sigma}^{(2)}$ is constructed the same as $\bm{\Sigma}^{(1)}$, except

[TABLE]

When $\beta=1$, this means that $\bm{\Sigma}^{(2)}$ is a matrix with 0.5 everywhere along the off-diagonal.

• •

Construction of $\bm{\Sigma}^{(3)}$: $\bm{\Sigma}^{(3)}$ is constructed in a similar way to $\bm{\Sigma}^{(1)}$, except there are three clusters. The first cluster contains indices $1,\ldots,\lfloor\beta\cdot d/6\rfloor,\lfloor d/2\rfloor+1,\ldots,\lfloor d/2+\beta\cdot d/6\rfloor$. The second cluster contains indices $\lfloor\beta\cdot d/6\rfloor+1,\ldots,\lfloor d/2\rfloor$. The third cluster contains indices $\lfloor d/2+\beta\cdot d/6\rfloor+1,\ldots,d$. Observe that this partitions $1,\ldots,d$, and when $\beta=1$, this results in three clusters of the roughly the same size. We then construct $\bm{\Sigma}^{(3)}$ using (E.2) but using these three clusters.

E.2 Details on functions $f_{1},\ldots,f_{d}$

At a high-level, these functions $f_{1},\ldots,f_{d}$ ensure that these marginal distributions of our sampled nonparanormal random variables are similar to the marginal distributions of the BrainSpan data. These marginal distributions are constructed in the following way. We first randomly sample $d$ variables (i.e., genes) uniformly from the BrainSpan dataset, $\{g_{1},\ldots,g_{d}\}\subseteq\{1,\ldots,n\}$. Next, for each $j$, let $\widehat{p}_{g_{j}}$ denote the kernel density estimate of variable $g_{j}$ in the BrainSpan dataset, using the default bandwidth selection used by the stats::density function in R.

We now formalize how to construct $f_{1},\ldots,f_{d}$. As described in Liu et al., (2009), we actually construct the inverse of these functions $f^{-1}_{1},\ldots,f^{-1}_{d}$ as they are more amendable for sampling, which must exist since $f_{1},\ldots,f_{d}$ are monotonic and differentiable. Recall that $\bm{\mu}=\mathbf{0}$. We first sample a vector $\bm{z}=(z_{1},\ldots,z_{d})$ from a Gaussian distribution $N(\mathbf{0},\bm{\Sigma})$. Let $\Phi(t;P)$ denote the cumulative distribution function evaluated at $t$ for a univariate density $P$. For any $j\in 1,\ldots,d$, we construct $f^{-1}_{j}$ such that

[TABLE]

That is, we construct $f^{-1}_{j}$ so that $z_{j}$ is at the same quantile with respect to $N(0,\Sigma_{jj})$ as $f^{-1}_{j}(z_{j})$ is with respect to the kernel density estimate $\widehat{p}_{g_{j}}$. Notice that by constructing $f^{-1}_{1},\ldots,f^{-1}_{d}$ in this fashion, each function is monotone and differentiable. We then set

[TABLE]

as one sample from the nonparanormal distribution $\text{NPN}(\mathbf{0},\bm{\Sigma},f)$.

Notice that by introducing non-Gaussianity into our simulation suite in this fashion, we ensure that the marginal distribution of all $r$ partitions resemble the BrainSpan dataset, and also ensure that the first $r_{1}$ partitions still are drawn from the same population covariance matrix. Also, by generating data in this fashion, we are able to obtain complicated dependencies between the mean and variance, as well as observe multi-modal distributions and heavier-tailed distributions compared to the Gaussian. See Liu et al., (2009) for a more detailed discussion.

E.3 Example of sampled nonparanormal distribution

We provide a visual illustration of what the sampled nonparanormal distribution could look like. We sample 375 samples from $\text{NPN}(\mathbf{0},\bm{\Sigma}^{(1)},f)$ when $\beta=0$, and plot two of the resulting pairwise scatterplots in Figure 12. We can think of the 375 samples as equivalent to aggregating all $r=25$ partitions together, each having $n=15$ samples. These two scatterplots show that the nonparanormal can display multiple modes marginally or heavier tails.

Appendix F Additional simulation results

F.1 Covariance homogeneity diagnostic in simulation

In this section, we apply the diagnostic developed in Section 3 to the simulation suite described in Section 5. Our goal is to determine how the QQ-plots behave as the selected partitions become less homogeneous. As done in Section 5, we consider four partition selection strategies: COBS (using $\alpha=0.1$ and $\gamma=0.95$), Base (which selects 3 partitions contain samples drawn from the nonparanormal distribution with proxy covariance $\bm{\Sigma}^{(1)}$, while the other 2 partitions contain samples from each of the remaining two distributions), All (which selects all $r$ partitions) and Oracle (which selects exactly the $r_{1}$ partitions containing samples drawn the nonparanormal distribution with proxy covariance $\bm{\Sigma}^{(1)}$).

We see in Figure 13 and Figure 14 that the QQ-plot is a reasonable diagnostic in this simulation suite. Between these two figures, we vary $\beta$ among 0, 0.3, 0.6 and 1. We notice that as $\beta$ increases, the QQ-plot derived from COBS remains relative uniform, similar to that of the Oracle. When $\beta=1$, COBS selects one erroneous partition in this particular trial shown, which results in the QQ-plot showing a deviation away from uniform. The QQ-plots derived from the Base procedure looks relative uniform when $\beta=0$ (which is to be expected, as all $r$ partitions share the same covariance matrix when $\beta=0$), but quickly has QQ-plots that deviate from uniform as $\beta$ increases. Note that the since the Base procedure selects only 5 partitions, there are a limited number of ways to split the partitions into two groups, which yields a limited number of points in the QQ-plot. The QQ-plots derived from the All procedure follow a similar trend as the Base procedure, but not as severe. These plots match the findings shown in Figure 8.

F.2 Simulation under Gaussian setting

While the simulations in Section 5 use nonparanormal distributions, we demonstrate that similar results hold for Gaussian distributions. This demonstrates that there is nothing particularly special about the nonparanormal or the Gaussian distribution that enable COBS to work well, and suggests COBS can work in much more general settings. Specifically, in this simuation suite, everything is the same as in Section 5, except all the functions $f_{1},\ldots,f_{d}$ are set to be the identity function. Hence, this means that the first $r_{1}$ partitions are drawn from Gaussian distributions with covariance $\bm{\Sigma}^{(1)}$, the next $r_{2}$ partitions are drawn from Gaussian distributions with covariance $\bm{\Sigma}^{(2)}$, and so on.

When we use Bonferroni or the Stepdown method in this Gaussian setting, we observe ROC curves for the individual hypotheses that strongly resemble Figure 6. This is shown in Figure 15.

Similarly, when we use COBS to select partitions, the ROC curves as well as the spectral error curves strongly resemble Figure 8A and B. This is shown in Figure 16.

F.3 Simulation using Accelerated Stepdown

In this subsection, we apply the Accelerated Stepdown procedure described in Section C within the COBS procedure in the simulation setting described in Section 5. Specifically, we use the test statistic (C.2) and analogous bootstrap statistics, but keep all other parts of the simulation suite the same.

When we plot the ROC curve for the individual hypotheses in Figure 17, we already notice a dramatic loss of power when compared to its original counterpart using the test statistic (3.1) shown in Figure 6. In fact, it seems like the Bonferroni procedure has almost no power at all, even when $\beta=1$.

Due to the loss of power for the individual hypotheses, we observe a loss of power for the selected partitions as well (Figure 18A) and spectral errors that strongly resemble selecting all the partitions (Figure 18B).

Appendix G Additional details on BrainSpan analysis

The first subsection describes the analysis pipeline we used throughout Section 6 in more detail. The second subsection describes the two distance metrics used in Subsection 6.4. The third subsection describes additional results alluded to in Subsection 6.4.

G.1 Description of analysis pipeline

We now summarize the pipeline used in Section 6 for clarity.

Screening of genes: This is the step described in Subsection 6.1, derived from Liu et al., (2015). We first select all genes whose p-value in the older TADA dataset (De Rubeis et al.,, 2014) is less than 0.01. Then, we rank all remaining genes by their maximum Pearson correlation in magnitude with any of the formerly selected genes in decreasing order based on the BrainSpan partitions within Window 1B aggregated. We select genes based on this ranking in order until we have selected a combined total of $d=3500$ genes. We analyze all 125 partitions using only these $d$ genes for the remainder of the analysis. 2. 2.

Applying COBS: This is the two-staged procedure we developed in this paper, detailed in Section 4. In the first stage, we apply the Stepdown procedure using $\alpha=0.1$. In the second stage, we select the clique-based selection method where $\gamma=0.95$, as well as using both extensions discussed in Subsection D.1. This results in 24 selected partitions within the BrainSpan dataset, as detailed in Subsection 6.2. We then combine all the 24 selected partitions to form a dataset $\mathbb{X}\in\mathbb{R}^{n\times d}$ with $n=272$ microarray samples and $d$ genes to be used for the remainder of the analysis. 3. 3.

Estimating the Gaussian graphical model: This step is described in Subsection 6.3 and is the same as in Liu et al., (2015). We fit a Gaussian graphical model using neighborhood selection (Meinshausen and Bühlmann,, 2006) based on $\mathbb{X}$, where the tuning parameter $\lambda$ (which controls the sparsity of the graphical model) is chosen such that the resulting graph has high scale-free index as well as a comparable number of edges to the estimated graph when COBS is not used. This choice of $\lambda$ is detailed at the end of this subsection. We defer the remaining estimation and computation details to Liu et al., (2015). We denote the adjacency matrix of the estimated graphical model as $\widehat{\bm{A}}\in\{0,1\}^{d\times d}$. 4. 4.

Estimating the HMRF: This step is also described in Subsection 6.3 and is the same as in Liu et al., (2015). We briefly summarize this step here, as it is less common in the statistical literature. Let $\bm{Z}\in\mathbb{R}^{d}$ denote the Z-scores for the selected genes, derived from the TADA scores in De Rubeis et al., (2014). We model $\bm{Z}$ using a HMRF, where for each gene in $j=1,\ldots,d$, $Z_{j}$ is an i.i.d. random variable drawn from a mixture of two Gaussians,

[TABLE]

where $I_{j}\in\{0,1\}$ is an unobserved Bernoulli random variable and $\mu\in\mathbb{R}$ and $\sigma\in\mathbb{R}_{+}$ are two unknown scalars to be estimated. The first Gaussian distribution represents the Z-scores for genes that are not associated with ASD, and the second Gaussian distribution represents the Z-scores for risk genes. The distribution of entire vector $\bm{I}\in\{0,1\}^{d}$ follows an Ising model with probability mass function,

[TABLE]

for any $\bm{\eta}\in\{0,1\}^{d}$ and two unknown scalars $b,c\in\mathbb{R}$ to be estimated. An EM algorithm is used to fit this HMRF model, and we obtain the estimated posterior probability $\widehat{p}_{j}=\mathbb{P}(I_{j}=0|\bm{Z})$, representing the probability gene $j$ is not a risk gene given the risk scores. We defer the estimation and computation details to Liu et al., (2015). 5. 5.

Applying Bayesian FDR: This step is also described in Subsection 6.3 and is the same as in Liu et al., (2015). We apply a procedure (Muller et al.,, 2006) to $\widehat{\bm{p}}$ to select a set of genes where the Bayesian FDR is controlled at level $10\%$. We defer the computation details to Liu et al., (2015). This results in the set of 209 detected risk genes detailed in Subsection 6.4.

Usage of De Rubeis et al., (2014).

We note that the older risk scores dataset (De Rubeis et al.,, 2014) is used twice, once in the screening stage (Step 1 above) and again to estimate the parameters of the HMRF (Step 4 above). As argued by Liu et al., (2015), it is important for this dataset to be the same in both steps, as the goal of DAWN is to boost the power of the risk scores by a “guilt-by-association” strategy. Hence, it is important to ensure the genes with low TADA scores remain in the analysis after screening, so they can implicate genes with TADA scores that are not as low.

Choice of $\lambda$.

We use the following procedure to tune $\lambda$ when estimating the Gaussian graphical model using only the 10 partitions from Window 1B as well as when using the 24 partitions selected by COBS. We tune $\lambda$ on a grid between 0.05 and 0.1, equally spaced into 15 values, for both graphical models. Our criteria for selecting $\lambda$ within this grid is inspired by Liu et al., (2015), who use a scale-free index, a number between 0 and 1 that measures how well the graph follows a power law. Specifically, we ensure the scale-free indices from both graphical models are approximately comparable as well as that both estimated graphical models have about 10,000 edges. Our focus on this number of edges comes from Liu et al., (2015), which estimated a graphical model with 10,065 edges. By ensuring both of our estimated graphical models have around 10,000 edges, we are able to ensure that both graphical models pass roughly the same amount of information into the HMRF stage of DAWN.

Using this procedure, we set $\lambda=0.05$ when estimating the graphical model using only the 10 partitions from Window 1B (for 9990 edges and a scale-free index of 0.77) and $\lambda=0.064$ when estimating the graphical model using the 24 partitions selected by COBS (for 9142 edges and scale-free index of 0.83).

G.2 Methods to measure distance of two nodes in a graph

As alluded to in Subsection 6.4, the shortest path distance and the commute distance do not seem like appropriate candidates to measure the distance between two genes (i.e., vertices) in a gene co-expression network (i.e., graph) due to the fact that the network estimated in the Window 1B analysis has more edges than in the COBS analysis (9990 and 9142 edges respectively). Hence, both of these distance metrics would naturally favor the denser graph.

To overcome this problem, we use two distance metrics that we believe enable a more fair comparison.

• •

Minimal spanning tree (MST) distance: This is a natural alternative to measure the distance between two vertices. Given a graph $G=(V,E)$, we first find the MST $G_{\text{(MST)}}\subseteq G$, and then compute the path distance between the two vertices in $G_{\text{(MST)}}$.

• •

Graph root embedding distance: A more statistically motivated way to measure the distance between two vertices is to first embed all vertices $V$ into a latent space. As shown in Lei, (2018), the graph root embedding is a natural candidate to do this, as it can theoretically represent a wide range of random graphs. This is essentially a more sophisticated spectral embedding. We first represent the graph $G$ as an adjacency matrix $\bm{A}$, and compute the top-$k$ eigenvectors (corresponding both the largest $k$ positive eigenvalues and largest $k$ negative eigenvalues in magnitude). Each vertex is then represented as a latent vector of length $2k$. The distance between two vertices is then defined as the Euclidean distance between their corresponding latent vectors. We defer the remaining details to Lei, (2018).

It is important to use both positive and negative eigenvalues since a scree plot reveals there are almost the same number of positive and negative eigenvalues for the adjacency matrices estimated in both the COBS and Window 1B analyses.

G.3 Additional results about closeness of genes

We provide more details that the 102 genes detected by the newer TADA scores (Satterstrom et al.,, 2020) are roughly 10%-30% closer to the 33 genes detected in the older TADA scores (De Rubeis et al.,, 2014) in the gene network estimated in the COBS analysis than in the Window 1B analysis. We call the 33 genes detected in De Rubeis et al., (2014) as the De Rubeis genes, and the 102 genes detected in Satterstrom et al., (2020) that are not part of the former 33 genes as the Satterstrom genes.

We use the MST distance defined above to ask: how far away are the closest $k$ De Rubeis genes from any Satterstrom gene on average (mean). Figure 19A plots this average distance against $k$. We use the graph root embedding distance to ask: how close is the nearest De Rubeis gene from any Satterstrom gene on average (mean) when using an embedding of latent dimension $2k$. Figure 19B plots this average distance against $k$. In both instances, regardless of how the parameter $k$ is chosen, the plot shows that the Satterstrom genes are closer to the De Rubeis genes on average. Both metrics show that the red curve is roughly 10%-30% lower than the pale curve across all values of $k$, hence giving our stated result.