T2-DAG: a powerful test for differentially expressed gene pathways via graph-informed structural equation modeling

With the continuous development of high-throughput sequencing technologies, there is an increasing need for powerful statistical methods for large-scale quantitative comparison between populations/conditions. For example, in gene expression data analysis, a major task is to identify genes from a large gene list that express differentially between different conditions, such as different disease statuses. Various approaches have been proposed for this task, providing critical insights into the molecular basis of many complex human traits/diseases; see Robinson et al. (2010), Love et al. (2014) and the references therein.

Besides separately analyzing individual genes, recently there is a growing interest in detecting differentially expressed gene pathways. Here, a gene pathway refers to a set of correlated genes with known interactions collectively contributing to specific biological functions. Many biologically meaningful gene pathways have been identified via experiments or statistical modeling, along with identified gene interactions, such as activation/expression or inhibition/repression of one gene by another. Public databases containing such pathways include Kyoto Encyclopedia of Genes and Genomes (KEGG, Kanehisa and Goto, 2000; Kanehisa et al., 2019), Gene Ontology (GO, Ashburner et al., 2000; Consortium, 2019), BioCarta databases (Nishimura, 2001), etc. Fig. 1 shows two directed acyclic graphs (DAGs) illustrating two KEGG gene pathways, which will be revisited in our data analysis in Section 4. $a$ $b$ and $a$ $b$ both indicate activation from gene $a$ to gene $b$; $a\dashrightarrow b$ and $a$ $b$ indicate expression and inhibition, respectively, from gene $a$ to gene $b$. A common and notable feature of the two pathways is the high sparsity of the edges in the corresponding DAGs; that is, only a small proportion of the gene pairs are connected. Such sparsity pattern is also present in most identified gene pathways (Wille et al., 2004; Li et al., 2012; Chun et al., 2015).

There are two main reasons for investigating differentially expressed gene pathways over individual genes. First, marginal tests for individual genes tend to have low statistical power due to limited sample sizes of typical gene expression datasets. On the other hand, gene pathway tests can achieve enhanced power by aggregating concerted signals from multiple related genes in the same pathway. Second, marginal analyses of individual genes often ignore interactions among genes, which potentially lead to spurious conclusions inconsistent with the underlying biological processes. Since gene interactions can provide critical insights into the molecular mechanism of the progression of traits/diseases, efficiently leveraging the auxiliary information on gene interactions can facilitate the detection of disease- and trait-specific genes and yield biologically more meaningful results.

To formalize the problem, let $\mathcal{G}=(\mathcal{V},\mathcal{D})$ denote the pathway of interest, where $\mathcal{V}=\{1,\ldots,p\}$ and $\mathcal{D}$ denote the gene set and edge set between genes, respectively. We assume that $\mathcal{G}$ is a DAG; i.e., all edges in $\mathcal{D}$ are directed edges with no loops. A DAG $\mathcal{G}$ can be summarized by an adjacency matrix $A=(a_{jk})\in\mathbb{R}^{p\times p}$, where $a_{jk}=1$ if the edge $j\rightarrow k$ is present in $\mathcal{D}$, in which case we call gene $j$ a parent of gene $k$ and gene $k$ a child of gene $j$, and $a_{jk}=0$ otherwise. Although some identified gene pathways cannot be described by DAGs, they tend to have a small number of loops and/or undirected edges (see Table S1 for a summary of the KEGG pathways included in our data analysis). Thus, for such gene pathways, analyses based on their largest directed acyclic subgraphs are still feasible and informative. Approaches to directly handling non-DAG pathways are feasible but will require techniques different from the ones in this study, which we will briefly discuss in Section 5.

Assume we have observed expression levels of the genes in $\mathcal{V}$ for individuals in two different conditions. Specifically, let $x^{(l)}_{ij}$ denote the expression level of gene $j\in\mathcal{V}$ for individual $i\in\{1,\ldots,n_{l}\}$ in group $l\in\{1,2\}$. We assume that the data are appropriately transformed so that ${\bf x}^{(l)}_{i}=(x^{(l)}_{i1},\ldots,x^{(l)}_{ip})^{\intercal}$, $i=1,\ldots,n_{l}$, are independent and identically distributed (i.i.d.) random vectors from a multivariate Gaussian distribution with mean ${\mbox{\boldmath${\mu}$}}_{l}$ ($l=1,2$) and covariance $\Sigma$. Besides observed gene expression data, there is auxiliary information on the presence/absence of directed relationships among genes (contained in $\mathcal{D}$), such as activation/inhibition of one gene by another. Our goal is then to test for the inequality of the two population mean, i.e., $H_{0}\mathrel{\mathop{\mathchar 58\relax}}{\mbox{\boldmath${\mu}$}}_{1}={\mbox{\boldmath${\mu}$}}_{2}\text{ versus }H_{a}\mathrel{\mathop{\mathchar 58\relax}}{\mbox{\boldmath${\mu}$}}_{1}\neq{\mbox{\boldmath${\mu}$}}_{2}$, based on ${\bf x}^{(l)}_{i}$s, while efficiently leveraging the auxiliary information on gene interactions.

Conventional hypothesis testing for the inequality of two mean vectors has been thoroughly investigated. A classical approach is the Hotelling’s $T^{2}$ test (Hotelling, 1931) with the test statistic

$$T^{2}=\frac{n_{1}n_{2}}{n_{1}+n_{2}}(\overline{\bf x}^{(1)}-\overline{\bf x}^{(2)})^{\intercal}\widehat{\Sigma}^{-1}(\overline{\bf x}^{(1)}-\overline{\bf x}^{(2)}),$$

(1)

where $\overline{\bf x}^{(l)}=n_{l}^{-1}\sum_{i=1}^{n_{l}}{\bf x}^{(l)}_{i}$ and $\widehat{\Sigma}=(n_{1}+n_{2}-1)^{-1}\sum_{l=1}^{2}\sum_{i=1}^{n_{l}}({\bf x}^{(l)}_{i}-\overline{\bf x}^{(l)})({\bf x}^{(l)}_{i}-\overline{\bf x}^{(l)})^{\intercal}.$ When the number of genes $(p)$ is smaller than the total sample size $(n_{1}+n_{2})$, $T^{2}$ is well-defined; under $H_{0}$, $T^{2}\left(n_{1}+n_{2}-p-1\right)/\left(p\left(n_{1}+n_{2}-2\right)\right)$ follows an $F$-distribution with parameters $p$ and $n_{1}+n_{2}-1-p$. The $T^{2}$ is a sum-of-squares type of statistic powerful against “dense alternatives", where the signals of differential expression are scattered among a large number of genes. This can often be true for gene pathway analyses, as genes in the same pathway tend to have concerted signals of expression that together form certain biological functions.

While Hotelling’s $T^{2}$ is well known as the uniformly most powerful invariant (UMPI) test, it is not applicable when $p$ is larger than $(n_{1}+n_{2})$, referred to as the high-dimensional setting hereafter. This is simply because in the high-dimensional setting the sample covariance matrix $\widehat{\Sigma}$ becomes singular, making $T^{2}$ ill-posed. Research has been conducted extensively on extending $T^{2}$ to the high-dimensional setting. For example, Bai and Saranadasa (1996) and Srivastava and Du (2008) directly modified the ill-posed $T^{2}$ by replacing $\widehat{\Sigma}$ in (1) with the identity matrix $I$ and the diagonal of $\widehat{\Sigma}$, respectively. Lopes et al. (2011) proposed to first project the high-dimensional data onto a randomly constructed low-dimensional space and then perform the Hotelling’s $T^{2}$ test. Chen et al. (2011) replaced $\widehat{\Sigma}$ with $(\widehat{\Sigma}+\lambda I)$ for some regularization parameter $\lambda>0$ and proposed a nonparametric testing procedure. Recently, Li et al. (2020) extended the work of Chen et al. (2011) and proposed an Adaptive Regularized Hotelling’s $T^{2}$ (ARHT) test that combines the test statistics corresponding to a set of $\lambda$s optimally chosen by a data-adaptive selection mechanism. Methods that are not direct modifications of Hotelling’s $T^{2}$ are also extensive. For example, Cai et al. (2014) proposed the CLX test with a supremum-type test statistic powerful against “sparse” alternatives. Gregory et al. (2015) proposed the GCT test that avoids the estimation of $\widehat{\Sigma}^{-1}$ by assuming an ordering of the variables so that the dependence between variables can be modeled according to their displacement. Xu et al. (2016) proposed the aSPU test, a robust adaptive test powerful against a wide range of alternatives.

While these existing methods can, to varying degrees, handle potentially high-dimensional gene expression data, they cannot incorporate the auxiliary information on gene interactions available from pathway databases. Incorporating such auxiliary information could substantially improve the test power, especially considering the inadequate sample size of a typical gene expression dataset. Previous efforts to address this issue include the NetGSA test (Shojaie and Michailidis, 2009 and Shojaie and Michailidis, 2010), which assumes that the magnitude of gene interactions is priorly known and incorporates such information using a linear mixed model. This approach has limited applicability in real applications because the magnitude of gene interactions is often unknown or only partially known and needs to be estimated using external data. Another approach proposed by Jacob et al. (2012), referred to as the Graph T2 test hereafter, projects data onto a low-dimensional space spanned by top eigenvectors of the graph Laplacian of the gene pathway and then performs the dimension reduction test in Lopes et al. (2011). Unlike the NetGSA test, Graph T2 can simply incorporate information on the presence/absence of gene interactions, but it may have compromised power if the signal is not coherent with the selected eigenvectors of the corresponding graph Laplacian.

To address these limitations of the existing methods, we propose T2-DAG, a novel test for identifying differentially expressed gene pathways, which efficiently leverages auxiliary information on gene interactions. Unlike the NetGSA test that requires additional information on the magnitude of gene interactions, the T2-DAG test only requires information on the presence/absence of gene interactions, which is more easily accessible in practice. Although the Gragh T2 test also requires only the information on the presence/absence of gene interactions, it simply incorporates a few eigenvectors of the corresponding graph Laplacian. On the contrary, T2-DAG directly incorporates all gene interactions into the analysis with minimal information loss. Our main idea is to link the auxiliary information on gene interactions with (partial) correlations among genes through a linear structural equation model (SEM). Based on the SEM, we construct a novel pathway-informed estimator of the precision matrix of the genes via the LDL decomposition (Krishnamoorthy and Menon, 2013), where the ordering of the genes is naturally determined by the topological ordering of the pathway. We then construct a test statistic by replacing the ill-posed $\widehat{\Sigma}^{-1}$ in (1) with the novel estimator of the precision matrix followed by a $Z$-transformation, which is denoted by T2-DAG $Z$. Under mild conditions, we prove that the proposed test statistic converges in distribution to the standard normal distribution. An alternative T2-DAG $\chi^{2}$ test can be constructed without the $Z$-transformation, which converges in distribution to $\chi^{2}_{p}$ under more stringent conditions. We also derive the asymptotic power of the tests under local alternatives that characterize the true signals of differential expression in the gene pathway between two populations/conditions. Under various simulated data scenarios, the proposed T2-DAG tests show well-controlled type I error rates and substantially higher powers compared to multiple existing methods, even with partially missing or incorrect information on gene interactions, unadjusted confounding effects, or non-Gaussian error terms. Although theoretically limited, the T2-DAG $\chi^{2}$ test compares favorably to the T2-DAG $Z$ test in the simulated scenarios, where it has slightly higher but still well-controlled type I error rates and higher powers. We also demonstrate the efficiency of the T2-DAG tests in a lung cancer-related application to detect differentially expressed gene pathways between different cancer stages.

Recall that our goal is to test for inequality of the mean expression levels of a gene pathway in two populations, $H_{0}\mathrel{\mathop{\mathchar 58\relax}}{\mbox{\boldmath${\mu}$}}_{1}={\mbox{\boldmath${\mu}$}}_{2}\text{ versus }H_{a}\mathrel{\mathop{\mathchar 58\relax}}{\mbox{\boldmath${\mu}$}}_{1}\neq{\mbox{\boldmath${\mu}$}}_{2}$, based on observed gene expression data, ${\bf x}^{(l)}_{i}=(x^{(l)}_{i1},\ldots,x^{(l)}_{ip})^{\intercal}$, $i=1,\ldots,n_{l}$, $l=1,2$, while leveraging available auxiliary information on gene interactions described by a DAG $\mathcal{G}=(\mathcal{V},\mathcal{D})$. Throughout the article, we use $\xi_{1}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}\xi_{2}$ to denote the independence between two random variables $\xi_{1}$ and $\xi_{2}$. We write $g(n)=O(f(n))$ if and only if $\lim_{n\rightarrow\infty}g(n)/f(n)<\infty$, and use $g(n)=o(f(n))$ to denote $\lim_{n\rightarrow\infty}g(n)/f(n)=0$. We also write $f(n)\asymp g(n)$ if and only if $\lim_{n\rightarrow\infty}|f(n)/g(n)|=C$ for some constant $C>0$. Let $I(\mathcal{A})$ denote the indicator function of the event $\mathcal{A}$; i.e., $I(\mathcal{A})=1$ if $\mathcal{A}$ is true, and $I(\mathcal{A})=0$ otherwise. We use lowercase letters, bold lowercase letters, and uppercase letters to denote scalars, vectors, and matrices, respectively.

We first define a new random variable, ${\bf z}_{i}=(z_{i1},\ldots,z_{ip})^{\intercal}$, based on the original gene expression data, ${\bf x}^{(l)}_{i}$, $i=1,\ldots,n_{l}$, $l=1,2$:

$${\bf z}_{i}=\begin{cases}{\bf x}^{(1)}_{i}-{\mbox{\boldmath${\mu}$}}_{1}&i=1,2,\ldots,n_{1}\\ {\bf x}^{(2)}_{i-n_{1}}-{\mbox{\boldmath${\mu}$}}_{2}&i=n_{1}+1,\ldots,n_{1}+n_{2}.\end{cases}$$

Given that ${\bf x}^{(l)}_{i}$ are i.i.d. random vectors that follow a multivariate normal distribution $\mathcal{N}({\mbox{\boldmath${\mu}$}}_{l},\Sigma)$, we can easily check that ${\bf z}_{i}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathcal{N}({\mbox{\boldmath${0}$}},\Sigma)$, $i=1,2,\ldots,n_{1}+n_{2}$. We propose to model each ${\bf z}_{i}$ through the following SEM with a vector of Gaussian noises:

$${\bf z}_{i}=Q^{\intercal}{\bf z}_{i}+{\mbox{\boldmath${\epsilon}$}}_{i},$$

(2)

where $Q=(q_{jk})\in\mathbb{R}^{p\times p}$ is an unknown coefficient matrix with $q_{jk}$ characterizing the magnitude of the effect of gene $j$ on gene $k$, and ${\mbox{\boldmath${\epsilon}$}}_{i}=(\epsilon_{i1},\ldots,\epsilon_{ip})^{\intercal}$ are $i.i.d.$ random vectors following $\mathcal{N}({\mbox{\boldmath${0}$}},R)$. We make the following assumption on $R$:

Assumption 1 states that $\epsilon_{ij_{1}}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}\epsilon_{ij_{2}}$ for all $j_{1}\neq j_{2}$ and $i=1,\ldots,n_{1}+n_{2}$. This amounts to the assumption that there is no latent confounders, which is common in the existing SEM literature; see Loh and Bühlmann (2014), Peters and Bühlmann (2014), and the references therein. We next link the coefficient matrix $Q$ in model (2) with the DAG $\mathcal{G}$ through the following assumption.

Assumption 2 is only required for theoretical analysis. Simulation studies in Section 3 will allow the auxiliary DAG to have missing or mis-specified edges. Since $\mathcal{G}$ is a DAG, we can label the variables in $\mathcal{V}$ according to the topological ordering of $\mathcal{G}$ so that $Q$ is a strictly upper triangular matrix, i.e., an upper triangular matrix with the diagonal entries being 0 (Loh and Bühlmann, 2014). As a result, model (2) can be rewritten as the following $p$ linear models:

$$z_{ij}=f_{j}(\epsilon_{i1},\ldots,\epsilon_{ij})\mbox{, }j=1,\ldots,p\mbox{, }i=1,\ldots,n_{1}+n_{2},$$

where $\{f_{j}(\cdot)\}_{j=1}^{p}$ are linear functions that are determined by $Q$. Thus, for $1\leq j_{2}<j_{1}\leq p$ and each $i$, one can check that $\mbox{Cov}(\epsilon_{ij_{1}},z_{ij_{2}})=0$, indicating that $\epsilon_{ij_{1}}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}z_{ij_{2}}$ under Gaussianity. This property will play an important role in establishing the asymptotic null distribution of our test statistic.

We next construct DAG-informed estimators of $\Sigma$ and $\Sigma^{-1}$. From model (2), we can obtain ${\mbox{\boldmath${\epsilon}$}}_{i}=(I-Q^{\intercal}){\bf z}_{i}$, $i=1,\ldots,n_{1}+n_{2}$. Applying the covariance function to both sides of this equation yields $R=(I-Q^{\intercal})\Sigma(I-Q).$ Since $Q$ is strictly upper triangular, we know that $I-Q$ is invertible, leading to

$$\Sigma=(I-Q^{\intercal})^{-1}R(I-Q)^{-1}.$$

(3)

Taking the inverse of both sides of (3) gives

$$\Sigma^{-1}=(I-Q){R}^{-1}(I-Q)^{\intercal}.$$

(4)

Equations (3) and (4) serve as the basis for the construction of the DAG-informed estimators of $\Sigma$ and $\Sigma^{-1}$, which are critical to the proposed test statistic. Note that (3) and (4) are, respectively, LDL-type decompositions of $\Sigma$ and $\Sigma^{-1}$; similar decompositions have been employed in the literature for estimating high-dimensional covariance/precision matrices (Wu and Pourahmadi, 2003; Huang et al., 2006). The rationale behind the LDL decomposition is that there exists some natural ordering of the variables, which, in the case of gene pathways, is determined by the topological ordering of the corresponding graph $\mathcal{G}$.

We observe from equations (3) and (4) that $\Sigma$ and $\Sigma^{-1}$ are functions of $Q$ and $R$. We then construct estimators for $Q$ and $R$. Letting ${\bf q}_{j}$ denote the $j$-th column of $Q$, we decompose model (2) into $p$ sub-models,

$$z_{ij}={\bf z}_{i}^{\intercal}{\bf q}_{j}+\epsilon_{ij},\mbox{ }j=1,\ldots,p.$$

(5)

To estimate the unknown parameters in (5), it is natural to consider the following square loss function for each sub-model $j$:

		$\displaystyle~{}\sum_{l=1}^{2}\sum_{i=1}^{n_{l}}\left\{x^{(l)}_{ij}-\mu_{lj}-\left({\bf x}^{(l)}_{i}-{{\mbox{\boldmath${\mu}$}}}_{l}\right)^{\intercal}{\bf q}_{j}\right\}^{2}$
	$\displaystyle=$	$\displaystyle~{}\sum_{l=1}^{2}\sum_{i=1}^{n_{l}}\left\{x^{(l)}_{ij}-\gamma_{lj}-{\bf x}^{(l)\intercal}_{i}{\bf q}_{j}\right\}^{2},$		(6)

where $\gamma_{lj}=\mu_{lj}+{{\mbox{\boldmath${\mu}$}}}_{l}^{\intercal}{\bf q}_{j}$. Letting $S_{j}=\{k\mathrel{\mathop{\mathchar 58\relax}}a_{kj}=1\}$, one can see that under Assumption 2, ${\bf q}_{j}$ has $|S_{j}|$ nonzero entries with $|S_{j}|\leq j-1$ for $j>1$ and that ${\bf q}_{1}={{\mbox{\boldmath${0}$}}}$. Since real gene pathways are far less dense than fully connected networks (Wille et al., 2004; Li et al., 2012; Chun et al., 2015), we make the following assumption on $S_{j}$:

This assumption appears to be reasonable for all 206 pathways considered in our data analysis in Section 4 (see Table S1 for details). Under Assumptions 1-3, we can estimate $\gamma_{lj}$ and the nonzeros entries of ${\bf q}_{j}$ in (2) using the ordinary least squares (OLS). Specifically, denote by ${\bf 1}_{n_{1}+n_{2}}$ the $(n_{1}+n_{2})$-dimensional vector of all ones and let ${\bf g}=(g_{1},\ldots,g_{n_{1}+n_{2}})^{\intercal}$ denote the group indicator; that is, $g_{i}$ equals 1 for $i=1,\ldots,n_{1}$ and equals 0 for $i=n_{1}+1,\ldots,n_{1}+n_{2}$. Then, letting $X=\left[{\bf x}_{1}^{(1)}~{}\cdots~{}{\bf x}_{n_{1}}^{(1)}~{}{\bf x}_{1}^{(2)}~{}\cdots~{}{\bf x}_{n_{2}}^{(2)}\right]^{\intercal}$, we rewrite (2) in the following matrix form

$$\left\|X_{j}-\theta_{1j}{\bf 1}_{n_{1}+n_{2}}-\theta_{2j}{\bf g}-X_{S_{j}}{\bf q}_{j,S_{j}}\right\|_{2}^{2},$$

(7)

where $X_{j}$ denotes the $j$-th column of $X$, $X_{S_{j}}$ denotes the sub-matrix of $X$ with columns indexed by $S_{j}$, and ${\bf q}_{j,S_{j}}$ denotes the sub-vector of ${\bf q}_{j}$ with entries indexed by $S_{j}$ that are non-zero according to the DAG and need to be estimated. Comparing (7) with (2), one can see that $\gamma_{1j}=\theta_{1j}+\theta_{2j}$ and $\gamma_{2j}=\theta_{1j}$. Thus, letting $W=[{\bf 1}_{n_{1}+n_{2}}~{}{\bf g}]$, we can obtain the OLS estimator of ${\bf q}_{j,S_{j}}$ and $(\theta_{1j},\theta_{2j})^{\intercal}$:

	$\displaystyle\widehat{\bf q}_{j,S_{j}}=$	$\displaystyle~{}\left(X_{S_{j}}^{\intercal}(I_{n_{1}+n_{2}}-\mathcal{P}_{W})X_{S_{j}}\right)^{-1}X_{S_{j}}^{\intercal}(I_{n_{1}+n_{2}}-\mathcal{P}_{W})X_{j},$
	$\displaystyle(\widehat{\theta}_{1j},\widehat{\theta}_{2j})^{\intercal}=$	$\displaystyle~{}(W^{\intercal}W)^{-1}W^{\intercal}\left(X_{j}-X_{S_{j}}\widehat{\bf q}_{j,S_{j}}\right),$		(8)

where $\mathcal{P}_{W}=W(W^{\intercal}W)^{-1}W^{\intercal}$ is the orthogonal projection matrix onto the column space of $W$. For ease of notation, we also denote the OLS estimator of $q_{kj}$ by $\widehat{q}_{kj}$ which equals 0 if $k\notin S_{j}$, and define the estimator of $Q$ as $\widehat{Q}=(\widehat{q}_{jk})_{j,k=1,\ldots,p}$. Given $\widehat{Q}$, we then construct a “plug-in” estimator of $R$ denoted by $\widehat{R}=\mbox{diag}(\widehat{r}_{1},\ldots,\widehat{r}_{p})$, with

$$\widehat{r}_{j}=\frac{1}{n_{1}+n_{2}-|S_{j}|-4}\left\|X_{j}-\widehat{\theta}_{1j}{\bf 1}_{n_{1}+n_{2}}-\widehat{\theta}_{2j}{\bf g}-X_{S_{j}}\widehat{\bf q}_{j,S_{j}}\right\|_{2}^{2}.$$

The denominator, $(n_{1}+n_{2}-|S_{j}|-4)$, makes $\widehat{r}_{j}$ slightly different from the classical OLS estimator of $r_{j}$. Although it may seem non-intuitive, we show in Lemma 1.3 in the supplementary material that $\mathbb{E}[\widehat{r}_{j}^{-1}]=r_{j}^{-1}$; this property is critical for establishing the asymptotic null distribution of our test statistic.

By replacing $Q$ and $R$ in (3) and (4) with $\widehat{Q}$ and $\widehat{R}$, we have the following DAG-informed empirical estimators of $\Sigma$ and $\Sigma^{-1}$:

$$\widehat{\Sigma}_{\text{DAG}}=(I-\widehat{Q}^{\intercal})^{-1}\widehat{R}(I-\widehat{Q})^{-1}~{}~{}\mbox{and}~{}~{}\widehat{\Sigma}_{\text{DAG}}^{-1}=(I-\widehat{Q})\widehat{R}^{-1}(I-\widehat{Q})^{\intercal}.$$

Finally, we propose our DAG-informed $Z$ test statistic:

$$T_{\text{DAG-}Z}=(2p)^{-1/2}\left(T^{2}_{\text{DAG-}\chi^{2}}-p\right),$$

(9)

where $T^{2}_{\text{DAG-}\chi^{2}}$ is a DAG-informed Hotelling’s $T^{2}$-type test statistic,

$$T^{2}_{\text{DAG-}\chi^{2}}=N(\overline{\bf x}^{(1)}-\overline{\bf x}^{(2)})^{\intercal}\widehat{\Sigma}^{-1}_{\text{DAG}}(\overline{\bf x}^{(1)}-\overline{\bf x}^{(2)}),$$

(10)

with $N=n_{1}n_{2}/(n_{1}+n_{2})$, which is an alternative test statistic for testing $H_{0}$ that will be discussed in Remark 2. As such, we call the proposed testing procedure T2-DAG hereafter.

Let $N_{e}=\sum_{j=1}^{p}|S_{j}|$ and $p_{0}=\sum_{j=1}^{p}I(|S_{j}|>0)$ denote the number of edges in $\mathcal{D}$ and the number of genes with at least one parent gene, respectively. The following result characterizes the asymptotic null distribution of $T_{\text{DAG-}Z}$.

Theorem 1 directly yields an asymptotically valid two-sided p-value for $T_{\text{DAG-}Z}$ for testing $H_{0}\mathrel{\mathop{\mathchar 58\relax}}{\mbox{\boldmath${\mu}$}}_{1}={\mbox{\boldmath${\mu}$}}_{2}$, that is, $p_{\text{DAG-}Z}=\mathbb{P}(|z|>|T_{\text{DAG-}Z}|)$, where $z$ is a random variable following the standard normal distribution. The assumptions in Theorem 1 often hold for typical gene pathways. In particular, $N_{e}=o\left(\sqrt{p}(n_{1}+n_{2})\right)$ implies that the total number of edges in the gene pathway is small compared to the square root of the number of genes multiplied by the total sample size, which is a weak assumption in real applications. See Table S1 for information on the real pathways considered in our data analysis.

To examine the asymptotic power of the proposed test, we first define the following two distance measures between the mean vectors of the two groups:

1. Kullback-Leibler divergence: $D_{\text{KL}}({\mbox{\boldmath${\mu}$}}_{1},{\mbox{\boldmath${\mu}$}}_{2})=({\mbox{\boldmath${\mu}$}}_{1}-{\mbox{\boldmath${\mu}$}}_{2})^{\intercal}(I-Q)R^{-1}(I-Q)^{\intercal}({\mbox{\boldmath${\mu}$}}_{1}-{\mbox{\boldmath${\mu}$}}_{2})$.

2. Pathway-specific divergence: $D_{\text{DAG}}({\mbox{\boldmath${\mu}$}}_{1},{\mbox{\boldmath${\mu}$}}_{2})=\sum_{j\mathrel{\mathop{\mathchar 58\relax}}|S_{j}|>0}({\mbox{\boldmath${\mu}$}}_{1}-{\mbox{\boldmath${\mu}$}}_{2})_{S_{j}}^{\intercal}\{\Sigma_{S_{j},S_{j}}\}^{-1}({\mbox{\boldmath${\mu}$}}_{1}-{\mbox{\boldmath${\mu}$}}_{2})_{S_{j}}$, where $({\mbox{\boldmath${\mu}$}}_{1}-{\mbox{\boldmath${\mu}$}}_{2})_{S_{j}}$ denotes the sub-vector of ${\mbox{\boldmath${\mu}$}}_{1}-{\mbox{\boldmath${\mu}$}}_{2}$ with entries indexed by $S_{j}$, and $\Sigma_{S_{j},S_{j}}$ denotes the sub-matrix of $\Sigma$ with rows and columns both indexed by $S_{j}$.

The following result characterizes the asymptotic power of the proposed test under local alternatives regarding $D_{\text{KL}}({\mbox{\boldmath${\mu}$}}_{1},{\mbox{\boldmath${\mu}$}}_{2})$ and $D_{\text{DAG}}({\mbox{\boldmath${\mu}$}}_{1},{\mbox{\boldmath${\mu}$}}_{2})$.

In lung cancer diagnosis and prognosis, extensive research has been conducted utilizing gene expression profiling to identify lung cancer-related gene pathways for therapeutic research (Weiss and Kingsley, 2008; Shao et al., 2010; Cai and Jiang, 2014; Dancik and Theodorescu, 2014; Chang et al., 2015; Long et al., 2019; Venugopal et al., 2019; Kong et al., 2020). In this application, we aim to identify differentially expressed gene pathways between different stages of lung cancer, the findings of which could provide insights into the underlying genetic mechanisms of lung cancer progression.

We applied T2-DAG $Z$ and T2-DAG $\chi^{2}$, as well as the existing tests considered in the simulation studies, to a gene expression dataset for lung cancer patients in The Cancer Genome Atlas (TCGA) Program (Network et al., 2014). The dataset was obtained from the Lung Cancer Explorer (LCE) (Quantitative Biomeical Research Center, 2019) with standard quality control, where the expression values of zero were set to the overall minimum value and all data were then $\log_{2}$-transformed so that the gene expression data are approximately normal for all lung cancer stages (Cai et al., 2014, 2019). Expression levels of 20429 genes were measured for lung cancer tissue samples collected from 513 patients, among whom $N_{1}=278$, $N_{2}=124$, $N_{3}=84$ and $N_{4}=27$ were diagnosed with stage I, II, III and IV lung cancer, respectively. Additionally, $N_{0}=59$ normal tissue samples were collected from 59 of the 513 patients. When comparing cancer tissues with normal tissues, we excluded the tumor samples of these 59 patients to avoid dependent samples.

We considered $H=206$ common human pathways related to signaling, metabolic and human diseases from the KEGG database (Kanehisa and Goto, 2000; Kanehisa et al., 2019). Pathway information on previously identified gene interactions is categorized into activation/expression or inhibition/repression of one gene by another (Kanehisa and Goto, 2000; Kanehisa et al., 2019); see Fig. 1 for an illustrative example. Detailed characteristics of the pathways are summarized in Table S1. Among the 206 pathways, 34 have circles (including self loops) in the corresponding graphs, among which 18 has 1 circle only, 8 have 2 circles, 4 have 3 circles, and 4 have more than 3 circles each. We removed these circles when implementing the proposed T2-DAG test but used the complete graph information when implementing the Graph T2 test. Other than this, all the tests were implemented in the same way as in Section 3. We further applied Bonferroni correction (Bonferroni, 1936) with a significance level $\alpha=\alpha_{0}/H$ for each test to control the family-wise error rate (FWER) at $\alpha_{0}=0.05$.

All tests except aSPU conclude that all 206 gene pathways are differentially expressed between normal tissues and tissues of each cancer stage (Tables S7 - S10). Table 1 summarizes the number of differentially expressed gene pathways detected by each test between stage I and stage III lung cancer and between stage II and stage III lung cancer. For the comparison between stage I and stage III lung cancer, the sample sizes ($N_{1}=278$, $N_{2}=124$) are large relative to the corresponding $p$, $p_{0}$ and $d$ for all pathways. The Hotelling’s $T^{2}$ test and Graph T2 test identify 8 and 10 significant pathways, respectively, while CH-Q, aSPU, GCT, SK and ARHT identify much more pathways (i.e., 96, 95, 82, 59 and 37, respectively). On the other hand, the CLX test, which is powerful for detecting sparse signals, only identifies 1 significant pathway. This may indicate that the signals of differential gene expression in real pathways tend to be dense, i.e., scatter among many correlated genes. The existing tests in total identify 129 unique pathways, while the T2-DAG $\chi^{2}$ test and $Z$ test identify 127 and 144 significant pathways, respectively, which cover 115 and 122, respectively, of the 129 pathways identified by the existing tests. Furthermore, the $-{\log}_{10}$-transformed p-values ($-{\log}_{10}p$) of the T2-DAG tests are positively correlated with those of all the other tests; in particular, results of T2-DAG tests are in line with those of CH-Q, SK, CLX, GCT, aSPU and ARHT, with correlations higher than 0.5. The two T2-DAG tests yield highly correlated results (the correlation of $-{\log}_{10}p$ equals 0.98). T2-DAG $Z$ detects a slightly larger number of significant pathways than T2-DAG $\chi^{2}$, and for the pathways detected by both tests, T2-DAG $Z$ tends to give smaller p-values than T2-DAG $\chi^{2}$. This is different from what was observed from simulations where T2-DAG $Z$ tends to have lower rejection rates than T2-DAG $\chi^{2}$, which could be due to the difference in data structures between the simulated data and real data. This inconsistency in relative performance is observed not only for the two T2-DAG tests, but also for other tests such as CH-Q and ARHT, which perform similarly on the simulated data but yield quite different results on the lung cancer data where they detect 96 and 37 differentially expressed pathways, respectively, between stage I and stage III lung cancer.

For the comparison between stage II and stage III lung cancer, all tests except GCT and aSPU detect no differentially expressed gene pathway. This may be due to the limited sample sizes ($N_{3}=124$, $N_{4}=84$) and/or relatively small differences between stage II and stage III lung cancer. On the other hand, aSPU identifies 1 significant pathway, while GCT identifies 21 significant pathways that do not include the one identified by aSPU. These 21 pathways are possibly false positives because of two reasons. First, none of the pathways is detected by any other test. Second, in our simulation studies, GCT could have substantially inflated type I error rates under small $p$ settings (Fig. 2, Fig. 4, Fig. S1 and Fig. S7).

To take a closer look at these results, we randomly selected 20 of the 206 pathways and report the corresponding $-{\log}_{10}(p)$s (Fig. 6: stage I versus stage III; Fig. S18: stage II versus stage III). Overall, the p-values for comparison between stage I and stage III lung cancer show consistency across different tests. Specifically, for pathways on which all tests yield insignificant results with large p-values, the T2-DAG tests give similar and sometimes even larger p-values compared to the other tests; see the Sulfur relay system pathway and the Circadian rhythm pathway in Fig. 6. For pathways on which T2-DAG yields highly significant results with the smallest p-values, many of the other tests also yield significant results; see the Calcium signaling pathway and the Type I diabetes mellitus in Fig. 6. On the contrary, for the pathways identified by GCT to be significantly differentially expressed between stage II and stage III cancer, all the other tests yield insignificant results with large $p$-values (Fig. 7).

We also compared other pairs of cancer stages and summarized the results in Section 3 of the supplementary material. Results for the comparison of stage I versus II, stage I versus IV, and stage II versus IV cancer have similar patterns as the comparison between stage I and III cancer. Similar to the comparison between stage II and III cancer, when comparing stage III to stage IV cancer, none of the tests identifies any significant pathway except for GCT which identifies 5 significant pathways. Detailed results for all the tests are summarized in Tables S1 - S10.

The proposed T2-DAG $Z$ test and the alternative T2-DAG $\chi^{2}$ test have illustrated their advantages over the existing tests through our simulation studies and application to a lung cancer dataset. When establishing asymptotic properties, the T2-DAG $Z$ test is favored over the alternative T2-DAG $\chi^{2}$ test since it requires less stringent conditions and thus is presented as our main result. Although facing theoretical limitations, the alternative T2-DAG $\chi^{2}$ test has excellent finite-sample performance where it has slightly higher but still well-controlled type I error rates and higher powers than T2-DAG $Z$, as shown by simulations. In the lung cancer application, however, T2-DAG $\chi^{2}$ appears to be slightly more conservative than T2-DAG $Z$ but still identifies a much larger number of differentially expressed pathways than the existing tests for stage I versus stage II, III, or IV cancer. In applications, the choice between T2-DAG $Z$ and T2-DAG $\chi^{2}$ can be made based on the theoretical properties and simulated finite-sample results under the specific data scenarios of sample size ($n_{1}$ and $n_{2}$), dimension ($p$), and sparsity of gene interactions ($p_{0}$ and $d$) inferred from the auxiliary pathway information.

The proposed T2-DAG test was motivated by gene pathway analyses, but it is generally applicable to problems in many other areas where similar auxiliary graph information is available. For example, in human microbiome studies for identifying host-microbe associations, the auxiliary graph information is provided by the phylogenetic tree that characterizes the evolutionary relationships among the microbes.

We have illustrated the robustness of T2-DAG with respective to model mis-specifications due to missing/inaccurate edge information or unadjusted confounding effects. While our simulation results are encouraging in these scenarios, the current estimation procedure may no longer be optimal. Thus, we consider the following two extensions of T2-DAG as our future work. First, while T2-DAG requires that the auxiliary pathway information can be represented by a DAG, it is feasible to extend it to non-DAG pathways with more complex gene interactions. For example, the recently proposed linear SEM for cyclic mixed graphs (Améndola et al., 2020) shows a possible way to account for feedback loops and undirected edges in the gene pathway. Second, the proposed test assumes a diagonal covariance matrix $R$ for the error term ${\epsilon}$ in (2). As discussed in Section 3.4, this assumption may fail due to latent confounders. In this case, we may consider a non-diagonal but sparse $R$, i.e., only a small number of off-diagonal entries of $R$ are non-zero, and estimate $R$ via graphical lasso (Friedman et al., 2008) or other penalized estimation approaches for high-dimensional precision matrix (Zhang and Zou, 2014; Kuismin et al., 2017; Fan et al., 2019). Furthermore, when measurements are available for some potential confounders, we can directly incorporate them into the SEM as covariates. We have in fact observed promising performance of this strategy through additional simulation studies, but the theoretical properties of the strategy require further investigation.

While the current T2-DAG test only incorporates information on the presence/absence of gene interactions, it is possible to further incorporate the signs of the interactions, such as activation ($+$) or inhibition ($-$) of one gene by another. One possible strategy is to conduct modeling in a Bayesian framework, where the signs of gene interactions can be accounted for by inducing truncated normal priors on the corresponding parameters in the coefficient matrix $Q$. The asymptotic properties of the test, however, need to be re-investigated.

Finally, we note that the key component of the proposed T2-DAG test is the novel DAG-informed estimator for the precision matrix (see Equations (3) and (4)), which could provide substantially improved power through incorporating external pathway information that describes gene interactions. Besides being incorporated into our T2-DAG test, this novel estimator has plenty of other potential usage. For example, as illustrated in Section 4, T2-DAG, with a Hotelling’s $T^{2}$-type test statistic, is powerful against “dense” alternatives which are usually true for differentially expressed gene pathways. In applications where the signals are sparse, one may instead consider combining the proposed DAG-informed estimator with tests that are powerful against “sparse” alternatives, such as the CLX test (Cai et al., 2014), to further boost power. Another example is that the proposed DAG-informed estimator can be used for high-dimensional linear discriminant analysis (LDA) (Bouveyron et al., 2007) to improve the classification accuracy. Specifically, consider the Gaussian case where one wants to classify a random vector $\bf x$ drawn with equal probability from one of the two Gaussian distributions, $\mathcal{N}({\mbox{\boldmath${\mu}$}}_{1},\Sigma)$ (class 1) and $\mathcal{N}({\mbox{\boldmath${\mu}$}}_{2},\Sigma)$ (class 2). When ${\mbox{\boldmath${\mu}$}}_{1},{\mbox{\boldmath${\mu}$}}_{2}$ and $\Sigma$ are known, Fisher’s LDA rule, given by

$$C(\bm{x})=\Bigg{\{}\begin{array}[]{ll}1,&~{}~{}~{}\bm{\delta}^{\top}\Omega\left(\bm{x}-\frac{\bm{\mu}_{1}+\bm{\mu}_{2}}{2}\right)<0\\ 2,&~{}~{}~{}\bm{\delta}^{\top}\Omega\left(\bm{x}-\frac{\bm{\mu}_{1}+\bm{\mu}_{2}}{2}\right)\geq 0,\end{array}$$

where ${\mbox{\boldmath${\delta}$}}={\mbox{\boldmath${\mu}$}}_{2}-{\mbox{\boldmath${\mu}$}}_{1}$, and $\Omega=\Sigma^{-1}$ denotes the precision matrix, is known to be optimal. In practice, one needs to estimate ${\mbox{\boldmath${\mu}$}}_{1},{\mbox{\boldmath${\mu}$}}_{2}$, and $\Omega$. The proposed DAG-informed estimator of $\Omega$ can thus be applied here to facilitate the classification.

The authors declare no conflict of interest.

The data underlying this article are available at https://github.com/Jin93/T2DAG.

Supplementary material available at Bioinformatics online includes technical results and proofs, detailed information on the pathways and gene expression data used in the data analysis, detailed simulation settings in the scenario of non-Gaussian error terms, as well as additional results on the simulated and real datasets.

The R (R Development Core Team, 2021) package T2DAG which implements the proposed T2-DAG test is available on Github at https://github.com/Jin93/T2DAG.