Combining assumptions and graphical network into gene expression data analysis

Analyzing gene expression data rigorously requires taking assumptions into consideration but also relies on using information about network relations that exist among genes. Combining these different elements cannot only improve statistical power, but also provide a better framework through which gene expression can be properly analyzed. We propose a novel statistical model that combines assumptions and gene network information into the analysis. Assumptions are important since every test statistic is valid only when required assumptions hold. So, we propose hybrid p-values and show that, under the null hypothesis of primary interest, these p-values are uniformly distributed. These proposed hybrid p-values take assumptions into consideration. We incorporate gene network information into the analysis because neighboring genes share biological functions. This correlation factor is taken into account via similar prior probabilities for neighboring genes. With a series of simulations our approach is compared with other approaches. Area Under the ROC Curves (AUCs) are constructed to compare the different methodologies; the AUC based on our methodology is larger than others. For regression analysis, AUC from our proposed method contains AUCs of Spearman test and of Pearson test. In addition, true negative rates (TNRs) also known as specificities are higher with our approach than with the other approaches. For two group comparison analysis, for instance, with a sample size of n=10, specificity corresponding to our proposed methodology is 0.716146 and specificities for t-test and rank sum are 0.689223 and 0.69797, respectively. Our method that combines assumptions and network information into the analysis is shown to be more powerful. These proposed procedures are introduced as a general class of methods that can incorporate procedure-selection, account for multiple-testing, and incorporate graphical network information into the analysis. We obtain very good performance in simulations, and in real data analysis.


Introduction
x Gene expression data can be analyzed in a multiple testing setting as well as many other statistical methods. The validity of each test depends on the underlying distributional assumptions of the test. A proper analysis of gene expression data requires taking assumptions, usually normality, into consideration (Pounds and Fofana 2012;Pounds and Rai 2009). In addition to incorporating distributional assumptions into the overall testing, it may also be informative to incorporate any prior knowledge of association between entities (Bowman and George 1995). Such associations are often recorded by graphical networks (Wei and Pan 2008). Combining these different elements, besides gaining statistical power, provides a framework through which analysis of gene expression data can be improved. We propose a novel statistical approach that incorporates testing for distributional assumption validity with prior information provided by gene graphical network. In particular, we use graphical networks to incorporate spatial dependence into the analysis of gene expression data. The spatial correlation is taken into account by assuming similar prior probabilities for neighboring genes.
We compare our approach with other methods through a series of simulations, and demonstrate that hybrid-network leads to an improvement on power over other approaches in most of the settings. The comparison of the different methodologies is based on specificities and/or Area under the ROC Curve (AUC). The specificity of a test is called the true negative rate; it is the proportion of samples that test negative using the test in question that are genuinely negative. An ROC curve or a receiver operating characteristic curve shows the performance of a classification model at all classification thresholds. An ROC curve is constructed by reporting sensitivities, true positive rates, on the y-axis and false positive rates on the x-axis.
The network analysis we use is the conditional autoregressive (CAR) model. CAR models are commonly used to represent spatial autocorrelation in data relating to a set of non-overlapping areal units. Those models are typically specified in a hierarchical Bayesian framework, with inference based on Markov Chain Monte Carlo (MCMC) simulation. The most widely used software to fit CAR models is WinBUGS or OpenBUGS. In our work, we use an R function BUGS(·) that helps run OpenBUGS inside R software. Another R function, CARBayes(·), is described in Lee (2013) that can be used for Bayesian spatial modeling with conditional autoregressive priors. Using CARBayes the spatial adjacency information can be specified as a neigbourhood matrix, whereas, with BUGS(·), the user has to specify an adjacency matrix.

Material and methods
Network information can be represented by directed or undirected graphs. Graphs are structures of discrete mathematics and have found applications in scientific disciplines that consider networks of interacting elements, such as genes that interact by sharing some biological resemblances. A graph consists of a set of nodes and a set of edges that connect the nodes. Usually the nodes are the entities of interest. For instance, each gene can be considered a node and the edges the relationships among the genes. A graph can be used in a practical way by developing software to translate between representations, a process sometimes referred to as "coercion".
In data analysis, graphs provide a data structure for knowledge representation, for example in the Gene Ontology (GO). Many studies incorporate gene network information Fofana et al. Journal of Statistical Distributions and Applications (2021) 8:9 Page 3 of 17

Fig. 1 Undirected Graphs
in data analysis through the GO project. Graphs provide a computational object that can easily and naturally be used to reflect physical objects and relationships of interest. Graphs are important to statistical methodology for exploratory data analysis. A knowledgerepresentation graph can be juxtaposed with observed data to guide the discovery of important phenomena in the data. In statistical inference, inferential statements about relations between genes due to significantly frequent co-citation, or relation between gene expression and protein complex can be made, (Wei and Pan 2008). A graph may be directed or undirected. A directed edge is an ordered pair of endvertices that can be represented graphically as an arrow drawn between the end-vertices. In such an ordered pair the first vertex is called initial vertex or tail and the second the terminal vertex or head. An undirected graph disregards any sense of direction and treats both head and tail identically, see Fig. 1.

Statistical models for hybrid testing
Considerthe following multiple hypothesis testings with θ g , a parameter for gene g, and G is the total number of genes, H og , the null hypothesis, and H 1g is the alternative hypothesis. Suppose two mutually exclusive test procedures, M 1 and M 2 , can be used to perform these statistical tests. When M 1 is used, suppose T 1 = {T 11 , · · · , T 1G } represents the test statistics and P 1 = {P 11 , · · · , P 1G } the corresponding set of p-values, and suppose T 2 = {T 21 , · · · , T 2G } and P 2 = {P 21 , · · · , P 2G } the corresponding quantities for procedure M 2 .

Fig. 2 Simulation Network
Suppose A g = i is an indication that the assumption for gene g holds for procedure M i for testing H og vs H 1g , i = 1, 2. For testing suppose T a = {T a1 , · · · , T aG } are the test statistics obtained from A g with the corresponding set of p-values P a = {P a1 , · · · , P aG }. And then, from this method, we define an appropriate summary statistic and denote it by P = {P 1 , · · · , P G } with P g = P 1g , ifA g = 1 P 2g , ifA g = 2 g = 1, · · · , G.
The following theorem states the distribution of P g under the null hypothesis H og of Eq. (1).
Theorem Suppose there are only two mutually exclusive procedures M 1 and M 2 that can be used to test the null hypothesis Let P 1 be the p-value obtained if method M 1 is used for testing the null hypothesis H 0 , and P 2 be the p-value if method M 2 is used instead. Let P be defined by Then P is uniformly distributed under the null hypothesis H 0 .
Proof First, we recall some probability theory basics. Let M 1 , M 2 , · · · , M n be a partition of a sample space , that is M i ∩ M j = ∅ ∀ i = j and n i M i = . Then, for any event i (E ∩ M i ) and then for any probability P, Fofana et al. Journal of Statistical Distributions and Applications (2021) 8:9 Page 5 of 17 . Also, the law of total probability holds for conditional probability, that is The question is to show that the hybrid p-value, P, follows a uniform distribution under the null hypothesis (H 0 ); that is F P (p) = P(P < p | H 0 ) = p, ∀p ∈ (0, 1), with F P the cumulative distribution function of P.
Under the null hypothesis (H 0 ) of primary interest (gene is expressed, say) and under M 1 and M 2 , respectively, both P 1 and P 2 are uniformly distributed, that is P(P 1 < p | M 1 , H 0 ) = p and P(P 2 < p | M 2 , H 0 ) = p, see (Pounds and Rai 2009) for instance. Recall that P is a random variable, since P 1 and P 2 are random variables. For the proof, we consider M 1 and M 2 as two events. The notation | H 0 means under the null hypothesis (H 0 ).
Thus P is uniformly distributed under H 0 . Now, transform the p-values by where is the cumulative distribution function of the standard normal distribution N(0, 1), and P g is the p-value corresponding to test g. The null distribution of Z g is the standard normal under H og of Eq. (1). Assume that under the alternative Z g ∼ N μ 1 , σ 2 1 , then where φ(·; μ 1 , σ 2 1 ) is the probability density function of N(μ 1 , σ 2 1 ), f is a density function.

Bayesian hierarchical models for spatial data
Conditional autoregressive (CAR) models are commonly used to represent spatial autocorrelation in data relating to a set of non-overlapping areal units. Those data are prevalent in many fields like agriculture (Besag and Higdon 1999), and epidemiology (Lee 2011). There are three different CAR priors commonly used to model spatial autoregression. Each model is a special case of a Gaussian Markov random field (GMRF) that can be written in a general form as where Q is a precision matrix that controls for the spatial autocorrelation structure of the random effects, and is based on a non-negative symmetric G × G neighborhood or weight matrix W, W = (w kj ) where w kj = 1 if genes k and j are neighboring genes and w kj = 0 otherwise, and φ = (φ 1 , · · · , φ G ), is a set of random effects. CAR priors are commonly specified as a set of G univariate fully conditional distributions , and G is the total number of genes (Lee 2013;Lee 2011). The first CAR prior proposed by Besag et al. (1991) is as The conditional expectation is the average of the random effects in neighboring genes, while the conditional variance is inversely proportional to the number of neighbors. The inverse proportionality of conditional variance is due to the fact that if random effects are spatially correlated then the more neighbors a node has the more information there is from its neighbors about the value of its random effect (subject-specific effect). This first CAR prior is used to implement the hybrid-network methodology as in Wei and Pan (2008). The second CAR prior proposed by Leroux et al. (1999) is given by while the third CAR prior proposed by Stern and Cressie (1999) is defined by where ρ is a spatial autocorrelation parameter, with ρ = 0 corresponding to independence and with ρ = 1 corresponding to a strong spatial autocorrelation. A uniform prior on the unit interval is specified for ρ, that is ρ ∼ ∪(0, 1), while the usual uniform prior on (0, M τ ) is assigned to τ 2 , with the default value being M τ = 1000. The intrinsic CAR prior by Besag et al. (1991) is obtained from the second and third CAR priors when ρ = 1, while when ρ = 0 the difference is on the denominator in the conditional variances.

Standard and spatial normal mixture model
Multipletesting is often an essential step in the analysis of high-dimensional data, such as genomic or proteomic data. The data analysis can be based on p-values, z-scores, tscores, etc. These test statistics are obtained from data reduction techniques. The hybrid p-values discussed in "Statistical models for hybrid testing" section is an example. Consider for example a test statistic Z. We can assume that across hypotheses g = 1, · · · , G the test statistic Z g follows a two-component mixture with density f as in (5). From this two-component mixture two different types of mixture models, the standard and spatial normal mixture models are considered. While spatial normal mixture models consider network information in the analysis, the standard normal mixture models do not.

Standard normal mixture model
In a standard two-component mixture model, Z g has a density function f of the form where π 0 is the proportion of genes that are not expressed (null hypothesis), f o is the distribution of Z g under the null hypothesis, and f 1 is the distribution of Z g under the alternative hypothesis.

Spatial normal mixture model
In a spatial normal mixture model, one defines gene-specific prior probabilities where T g is defined by where z g is the expression value of gene g for g = 1, · · · , G, and π g1 = 1 − π g0 . It is believed that genes on the same network, that is a group of genes with the same function, share the same prior probability of expression while different networks have possibly varying prior probabilities. The prior probabilities π gs , based on a gene network, are related to two latent Markov random fields x s = x gs ; g = 1, · · · , G , s = 0, 1 by a logistic transformation: Each of the G-dimensional latent vectors x s is distributed according to an intrinsic Gaussian conditional auto-regression model (ICAR) (Besag and Kooperberg 1999). The distribution of each spatial latent variable x gs conditional on x −gs = {x ks ; k = g} depends only on its direct neighbors. To be more specific, where δ g is the set of indices for the neighbors of gene g, and m g is the corresponding number of neighbors. The other model specifications are articulated in this way g = 1, · · · , G and s = 0, 1. Network structure is summarized in a matrix format called an adjacent matrix: Adj = (a ij ), i = 1, · · · , G; j = 1, · · · , G, where a ij = 1, if i = j and genes i and j are related 0, otherwise.

Prior distributions
In a standard normal mixture model, a beta distribution is often assumed as the prior distribution for π 0 . In a spatial normal mixture model, gene-specific prior probabilities are introduced. For the spatial normal mixture model, the prior probabilities for π gs , based on a gene network, are related to two latent Markov random fields (MRFs), as mentioned previously. From Eq. (14), we assume priors on the variance components σ 2 cs ∼ inverse-gamma(0.01, 0.01), the corresponding precision 1 σ 2 cs has gamma(0.01, 0.01) with mean 1 and variance 100. σ 2 cs acts as a smoothing parameter for the spatial field and consequently controls the degree of dependency among the prior probabilities of the genes. The size of σ 2 cs determines how similar the π gs are. The smaller the σ 2 cs are, the more similar the π gs . with θ s = μ s , σ 2 s , s = 0, 1 and Z is a gene expression test statistic. A direct approach to estimate π 0 , π 1 , θ 0 , and θ 1 is to compute the likelihood function

Maximum a posterior estimation
and the log likelihood as Obtaining MAPE's of the parameters directly is not possible. To estimate the parameters the expectation-maximization (EM) algorithm may be used. In order to use the EM algorithm, define latent variables v = (v gk , z gk ) | k = 1, · · · , n and g = 1, with G 0 (genes not expressed) and G 1 (expressed genes) are null hypothesis and alternative groups respectively, n is the sample common to all genes. If we include latent variables we get complete data, the observed z s and the unobserved v s. The maximum a posterior function for the complete data is Taking the log on Eq. (19) we get the log maximum a posterior function as The EM algorithm can be used to obtain MAPE's of π 0 , π 1 , θ 0 and θ 1 , if (z 1k , z 2k , · · · , z Gk ) are assumed to be independent. However, since there is a graphical network among genes, (z 1k , z 2k , · · · , z Gk ) are not independent. In order to take into account gene graphical network a Bayesian methodology is used. Network analysis is brought into the analysis by generating latent variables according to GMRFs as in Eq. (14). After assigning prior distributions to the parameters, posterior distributions can be found using a partial Gibbs sampler and some Metropolis Hasting algorithms. We use OpenBugs software to get the MAPE's of π 0 , π 1 , θ 0 , and θ 1 .

Statistical inference
The decision rule and acceptance of null hypotheses is based on probabilities from posterior distributions. For each gene g, the point estimate of p(H 0g | Data) is computed and compared to a threshold τ , for g = 1, · · · G. H 0g is rejected whenp(H 0g | Data), the point estimate, of p(H 0g | Data) is less than a threshold τ .
The p-values p g obtained from the hybrid method are transformed, and the transformed statistics z g = −1 (1 − p g ) are used, with −1 the standard normal quantile function. Through Bayesian modeling, network information is added to the analysis. With the Bayesian inference these posterior estimates areπ g0 =p(H 0g | Data). Inferences for the Bayesian hierarchical models are obtained using MCMC simulations, with a combination (2021) 8:9 Page 10 of 17

Simulations
To compare the hybrid-network method with other methods, we conducted simulation studies designed to mimic real data analysis. We conducted standard two-group comparison studies (treatment vs control), k-group (k > 2) comparison (ANOVA), and regression analysis. The k-group comparison is directly applicable to a genomic study comparing human ependymoma, a brain tumor that occurs in three distinct anatomic regions: Posterior Fossa (PF), Spine (SP), and Supratentorial (ST). Regression analysis is often useful to determine whether, for example, gene expression levels are related to a particular covariate such as DNA synthesis rate (INHIBO). For each of the three types of analyses conducted in the simulation studies, two different tests can be used. The first one requires the normality assumption while the second may be appropriate when the normality assumption does not hold. For the two-group comparison the hybrid-network method chooses between the standard t-test for normally distributed data and the Wilcoxon test when the normality assumption fails. For k-group (k > 2) comparison, the hybrid-network method chooses between the standard ANOVA test and the Kruskal-Wallis test. For the regression analysis, the Pearson test for linear dependency is chosen when the normal assumption holds and the Spearman test if the normality assumption does not hold. Fofana et al. Journal of Statistical Distributions and Applications (2021) 8:9 Page 11 of 17 In Eq. (12), we useπ g0 , the estimate of π g0 . And, the decision rule consists of rejecting the null hypothesis, H g0 , for gene g, ifπ g0 is less than a threshold, τ . The conclusion is that the corresponding gene g is expressed. For cancer data analysis, for instance, if a gene is expressed, health researchers will target that gene in finding cure.
The comparison of the different methodologies is mainly based on specificities (not to reject the null hypotheses when they are true, we call them sometimes true negatives). We could provide both specificities and sensitivities (reject the null hypotheses when they are not true, we call them sometimes true positives); but we have decided to compute only specificities because the simulations are computationally intensive.

K-group comparison study
In a group comparison study, gene expression data can be modeled as: where Y gij is the expression level for gene g of the j th individual in the i th group, g = 1, · · · , G, i = 1, · · · , k; j = 1, · · · n i , k is the number of groups, n i is the sample size of group i, and gij ∼ N(0, 1) or gij ∼ t(ν), or gij ∼ another distribution.
A 2-group comparison (k = 2), interest is in statistical tests of the form g = 1, · · · , G. Some gene expression levels may be normally distributed while others are not normally distributed. In the two-group comparison study, two tests are often used. The t-test is used when the normality assumption holds and the Wilcoxon test, a non parametric test, is often used when the normality assumption does not hold. For each gene g, a t-test, a Wilcoxon-Mann-Whitney rank sum test, and a Shapiro-Wilk test statistics are computed. Diagnoses for adequacy of the t-test statistics are made through residuals. We compute the residuals from the t-test statistic. We define the residuals on observation, j, in treatment, i, for gene, g, as whereŶ gij is an estimate of the corresponding observation Y gij obtained as follows: If the model is adequate, residuals should be structure-less; that is, they should contain no obvious patterns. Through an analysis of residuals, many types of model inadequacies and violations of the underlying assumptions can be discovered. We use the residuals to check for normality. A probit plot of residuals is an extremely useful procedure to test for normality. If the underlying error distribution is normal, this plot will resemble a straight line. Also outliers can be detected through residuals. Outliers show up on probability plots as being very different from the main body of the data. Plotting the residuals in time order of data collection is helpful in detecting correlation between the residuals. This is useful for checking independence assumptions on the errors. To compare the hybrid-network method with other methods, we perform a simulation study. In this setup, there are two groups of sample sizes varying from 5, 10, 25, and 50. The number of gene expressions having a normal distribution, N(μ, 1), is 30. For these gene expressions, μ = 0 for the null hypothesis and μ = 1 for the alternative. The remaining gene expressions have log-normal distribution, log-normal (μ, 1), with μ = 0 in some cases and μ = 1 in other cases. And a graphical network, Fig. 2, is built among genes with 212 number of edges. We translate this graphical network into an adjacent matrix.
The results are presented in Table 1, they show that hybrid-network procedure dominates the other methodologies in most of the settings, since the hybrid-network test specificities are higher than the specificities of the other methods. When the sample size is equal to 5, for instance, the specificity corresponding to the t-test is 0.571726, the specificity corresponding to the Wilcoxon test is 0.557244, and the specificity for the hybrid-network test is 0.575314.

Hybrid ANOVA-Kruskal Wallis study
In a k-group comparison study, a statistical model can be written as Eq. (21). For the model (21), μ g is a parameter common to all treatments for gene g called the overall mean, and τ gi is a parameter unique to the ith treatment for gene g called the ith treatment effect. Consider the following multiple hypothesis tests H g0 : μ g1 = μ g2 = · · · = μ gk vs H gA : μ gi = μ gl for at least one pair (i, l) or equivalently, by using the effects models H g0 : τ g1 = τ g2 = · · · = τ gk = 0 vs H gA : τ gi = 0 for at least one i.
The hypotheses may be tested using an ANOVA test or the Kruskal-Wallis depending on the normality assumption. If the normality assumption is valid, the ANOVA test is more powerful than the Kruskal-Wallis; and the latter may be more powerful when the normality assumption does not hold. The proposed methodology, hybrid-network, combines a test of assumptions and graphical network information into the analysis. For each gene g, an ANOVA p-value, p a g , a Kruskal-Wallis p-value, P w g , and a Shapiro-Wilk p-value, P s g are computed. We define a hybrid p-value, P h g , as for g = 1, · · · , G where α is a given threshold. The hybrid p-value P h g is transformed into a hybrid z-statistic, z h g , as follows: We use z h g to build a CAR model from the given network with the marginal distribution of z h g given by where z h g is the expression value for gene g, g = 1, · · · , G.  This shows the human ependymoma expression data: genes as gene annotation, groups (Gr1 and Gr2) as sample annotation and real values as gene expression levels.
The prior probabilities π gs , based on a gene network, are related to two latent Markov random fields x s = {x gs ; g = 1, · · · , G}, s = 0, 1 by a logistic transformation: .
The distribution of each spatial latent variable x gs conditional on x −gs = {x ks ; k = g} depends only on its direct neighbors. The proposed CAR prior distribution from (Besag and Kooperberg 1999) is used as where δ g is the set of indices for the neighbors of gene g, and m g is the corresponding number of neighbors. The hybrid-network methodology, through a series of simulations, is compared to other methods. The setup of these simulations consists of three groups of sample size varying from 5, 10, 25, and 50. The number of genes with the normal distribution N(μ, 1), μ = 0 for the null hypothesis and μ = 1 for the alternative, is 30. The number of genes with the log-normal distribution, log-normal(μ, 1), with μ = 0 in some cases and μ = 1 in other cases, is 7 and the number of genes with the Cauchy distribution, Cauchy(θ, 1), with θ = 0 in some cases and θ = 1 in other cases, is 7. A graphical network is built among genes with 212 edges. We present the simulations results in Table 2. They show that hybridnetwork procedure dominates other procedures in most of the cases. When the sample size is 25, for instance, the specificities from the ANOVA test, the Kruskal Wallis and the hybrid-network test are 0.89141, 0.918197, and 0.929054, respectively.

Regression analysis
In microarray regression analysis, a statistical model can be written as where Y gj is the gene expression level for the g th gene in the j th individual with g = 1, · · · , G, j = 1, · · · , n and some gj ∼ N(0, 1) or gj ∼ t(ν), or gj ∼ another distribution.

Fig. 5 Tumor Data 2-Group Comparison
The question is whether a response variable and a covariate are correlated. To test for correlation between gene expression with a covariate such as a phenotype, the analysis can be based on Pearson test p-values (P p ), and on Spearman test p-values (P sp ). We can use Shapiro-Wilk p-values (P s ) to test for the normality assumptions. Consider, the regression analysis in matrix format where We denote the least squares estimators of β g as b g Let the vector of the fitted valuesŶ gi be denoted asŶ g , and the vector of the residual terms e gi = Y gi −Ŷ gi be as e g . The fitted values are represented bŷ and the residuals by For each gene g, compute its Pearson p-value, P p g , compute its Spearman p-value, P sp g , and from the residuals from Pearson test, a Shapiro-Wilk test of normality is performed, and for each gene g a p-value, P s g , is calculated. Finally, a hybrid p-value, P h g is computed as where α is a given threshold. Each hybrid p-value, P h g , is transformed into a hybrid z-statistic, z h g , as follows: Using z h g , the marginal distribution of z h g is given as where z h g is the expression value of gene, g, g = 1, · · · , G. We compare the hybrid-network with the other procedures through a simulation setup. The setup consists of a sample size of 25. The number of genes with the normal distribution, N(μ, 1), is 30, μ = 0 for the null hypothesis and μ = 1 for the alternative, and the number of genes with the log-normal distribution, log-normal(μ, 1), with μ = 0 in some cases and μ = 1 in other cases, is 14. We vary the cutoff point, τ , as in Wei and Pan (2008). And a graphical network is built among genes with 212 number of neighbors. The results of the analysis are presented in Fig. 3. In order to compare the hybrid-testing with other methods, we use AUCs to judge the performance of the proposed method. A greater AUC corresponds to a better methodology. They show that the hybrid-network performs better than the other competing procedures.

Application to human ependymoma microarray
We compare the hybrid-network procedure with the t-test and the Wilcoxon test using human ependymoma data. The data consists of gene expression levels, gene annotation, sample annotation, and a gene graphical network. Figure 4 illustrates a graphical network of the genes under consideration, and Table 3 is a subset of the human ependymoma expression data. In this analysis, there are two groups, the sample sizes are n 1 = 37 for group1, n 2 = 42 for group2, with the total number of genes of 102, and the number of edges is 196. The data and the R codes can be requested from the corresponding author.
Using Shapiro-Wilk p-values, it appears that some of the expression data are normally distributed and the others are not, with Shapiro-Wilk test p-values less than α = 5% for some genes. Figure 5 shows histograms of p-values from the t-test, p-values from the rank sum test, and p-values based on the Shapiro-Wilk test of normality, respectively. The last graph of Fig. 5 presents the plot of the p-values from the t-test with respect to Fofana et al. Journal of Statistical Distributions and Applications (2021) 8:9 Page 16 of 17 the corresponding p-values from the rank sum test. Using the t-test when the normality assumption is assumed, and the Wilcoxon test otherwise. We apply the hybrid-testing procedure to analyze the data. We incorporate a graphical network to accommodate interactions between genes, as these have been noted to play a crucial role in cell functions (Shojaie and Michailidis 2009). In order to compare the hybrid-network procedure with the other procedures, we report results for the first six genes. We use box plots as visual methods of comparing groups. Under each Box plot, we report the results,π ·0 , with t representing the t-test statistic , rs for Wilcoxon test statistic, and hybN for hybrid-network statistic. We also present the p-values from Shapiro Wilk test (Shp) under each box plot. The results are reported on Fig. 6.
With a cutoff point of τ = 0.1, (τ is is a classification threshold, it is like, say α, the level of significance, see Wei and Pan (2008)), all the three methods find that genes AKT1, ATF2, and CDC25B are not expressed. Only the hybrid-network test finds that the other three genes, ARHGEF2, BDNF and BRAF are expressed. This finding is in accordance with the box plot results. The gene selection is based on R head(·) function, that selects the 6-first results. By doing so, we have tried to avoid criticism of biasness in selecting genes to analyze. First, we sort the genes and then pick the 6-first genes for comparing the different methodologies.