{'fraction_non_alphanumeric': 0.051635014106375325, 'fraction_numerical': 0.013737148050669343, 'mean_word_length': 4.652129708904109, 'pattern_counts': {'":': 3, '<': 3, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 21, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multiple tests of association with biological annotation metadata

We propose a general and formal statistical framework for multiple tests of association between known fixed features of a genome and unknown parameters of the distribution of variable features of this genome in a population of interest. The known gene-annotation profiles, corresponding to the fixed features of the genome, may concern Gene Ontology (GO) annotation, pathway membership, regulation by particular transcription factors, nucleotide sequences, or protein sequences. The unknown gene-parameter profiles, corresponding to the variable features of the genome, may be, for example, regression coefficients relating possibly censored biological and clinical outcomes to genome-wide transcript levels, DNA copy numbers, and other covariates. A generic question of great interest in current genomic research regards the detection of associations between biological annotation metadata and genome-wide expression measures. This biological question may be translated as the test of multiple hypotheses concerning association measures between gene-annotation profiles and gene-parameter profiles. A general and rigorous formulation of the statistical inference question allows us to apply the multiple hypothesis testing methodology developed in [Multiple Testing Procedures with Applications to Genomics (2008) Springer, New York] and related articles, to control a broad class of Type I error rates, defined as generalized tail probabilities and expected values for arbitrary functions of the numbers of Type I errors and rejected hypotheses. The resampling-based single-step and stepwise multiple testing procedures of [Multiple Testing Procedures with Applications to Genomics (2008) Springer, New York] take into account the joint distribution of the test statistics and provide Type I error control in testing problems involving general data generating distributions (with arbitrary dependence structures among variables), null hypotheses, and test statistics.


Motivation
Experimental data, such as microarray gene expression measures, gain much in relevance from their association with biological annotation metadata, i.e., data on data, such as, GenBank sequences, Gene Ontology terms, KEGG pathways, and PubMed abstracts. A challenging and fascinating area of research for statisticians concerns the development of methods for relating experimental data to the wealth of metadata available publicly on the WWW. Tasks include accessing and pre-processing the data, making inference from these data, and summarizing and interpreting the results.
In this context, an important class of statistical problems involves testing for associations between known fixed features of a genome and unknown parameters of the distribution of variable features of this genome in a population of interest. Here, features of a genome are said to be fixed, if they remain constant among population units. In contrast, variable features are allowed to differ among population units. Fixed features typically consist of gene annotation metadata, that reflect current knowledge on gene properties, such as, nucleotide and protein sequences, regulation, and function. Variable features often consist of gene expression measures, that reflect cellular type/state/activity under particular conditions. The fixed and variable features define, respectively, gene-annotation profiles and gene-parameter profiles; the parameter of interest then corresponds to measures of association between known gene-annotation profiles and unknown gene-parameter profiles.
For instance, for the yeast Saccharomyces cerevisiae (in short, S. cerevisiae), one may be interested in detecting associations between the vector of mean transcript (i.e., mRNA) levels for all (approximately 6,500) genes under heat-shock conditions and Gene Ontology (GO) annotation for these genes. The reader is referred to the Gene Ontology Consortium website (www.geneontology.org) and to Section 4, below, for more information on gene ontologies, and to the Saccharomyces Genome Database (SGD) website (www.yeastgenome.org), for details on S. cerevisiae. In this example, the population of interest may consist of all heat-shocked cells from well-defined cultures of a particular strain of S. cerevisiae (e.g., strain S288C). For each of the three gene ontologies (BP, CC, and MF, as described in Section 4.1), each gene is annotated with a fixed set of GO terms (i.e., this set is constant across population units for a given version of the GO Database). Thus, for a given GO term, one may define a gene-annotation profile as a known, fixed binary vector indicating for each gene whether it is annotated or not with the particular GO term. The transcript levels, however, vary among population units and the gene-parameter profile, i.e., the vector of genome-wide mean transcript levels in the population of heat-shocked yeast cells, is unknown and may be estimated, for example, from a microarray experiment involving a sample of yeast cells from the population. The association parameter of interest, between GO annotation and transcript levels, is then a vector of association measures (e.g., two-sample t-statistics) between the known binary gene-annotation profiles and the unknown continuous gene-parameter profile.
Similar inference questions arise in many other contexts and involve a variety of definitions for the gene-annotation profiles, the gene-parameter profiles, and the association parameters of interest. For example, in cancer microarray studies, one may seek associations between GO gene-annotation profiles and a gene-parameter profile of regression coefficients relating (censored) patient survival data to genomewide transcript levels or DNA copy numbers. Furthermore, gene-annotation profiles need not be binary or even polychotomous, and may correspond to pathway membership, regulation by particular transcription factors, nucleotide sequences, and protein sequences.
Note that, for the sake of illustration, we focus on gene-level features. However, our proposed methodology is generic and may be applied to other types of features, such as those concerning gene isoforms and proteins. For instance, as in alternative splicing microarray analysis, one may collect data at the finer level of gene isoforms, where one gene may have multiple isoforms [10]. In this context, isoformparameter profiles may refer to the distribution of microarray isoform expression measures in a well-defined population, while isoform-annotation profiles may consist of exon/intron counts/lengths/nucleotide distributions. One may also consider protein-level features, where, for example, protein-parameter profiles correspond to antibody microarray expression measures and protein-annotation profiles refer to protein function, domain structure, and post-translational modification (e.g., Swiss-Prot, www.expasy.org/sprot).
Methods proposed thus far suffer from a number of limitations, related, to a large extent, to the absence of a clear and precise statement of the statistical inference question. As a result, the analyses often lack statistical rigor and tend to be ad hoc and dataset-specific.
One of our main contributions is the systematic and precise translation of a general class of biological questions into a corresponding class of multiple hypothesis testing problems. A key step in this process is the proper definition of the geneannotation profiles, gene-parameter profiles, and association parameters of interest. This general formulation then allows us to apply the multiple hypothesis testing methodology developed in [14] and related articles [8,15,16,31,32,33,34,39,40,41,42], to control a broad class of Type I error rates, defined as generalized tail probabilities (gTP), gT P (q, g) = Pr(g(V n , R n ) > q), and generalized expected values (gEV), gEV (g) = E[g(V n , R n )], for arbitrary functions g(V n , R n ) of the numbers of false positives V n and rejected hypotheses R n .
We wish to emphasize the crucial and often ignored distinction between: (i) defining a parameter of interest, measuring the association between gene-annotation and gene-parameter profiles, i.e., the statistical formulation of the biological question; (ii) making inferences, i.e., estimating and testing hypotheses concerning this parameter, based on a sample drawn from the population under consideration. Most methods proposed to date focus on (ii), without providing a clear statement of the question being answered in (i), that is, various estimation and testing procedures are proposed for an undefined parameter of interest.
Due to its general and rigorous statistical framework, our approach to multiple tests of association with biological annotation metadata differs in a number of important ways from current approaches, such as those developed for inference with Gene Ontology metadata and implemented in the software packages listed on the Gene Ontology Tools webpage (www.geneontology.org/GO.tools.shtml).
General gene-annotation profiles. Existing approaches typically consider binary gene-annotation profiles, e.g., vectors of indicators of GO term annotation. Our general definition of gene-annotation profiles allows consideration of arbitrary qualitative and quantitative fixed features of a genome, e.g., membership of genes to any number of pathways or clusters, exon/intron counts/lengths/ nucleotide distributions, mean transcript levels. General gene-parameter profiles. Existing approaches typically consider binary gene-parameter profiles, e.g., vectors of indicators of differential expression. Our general definition of gene-parameter profiles allows consideration of a much broader class of testing problems, concerning arbitrary qualitative and quantitative parameters, such as, differences in mean expression levels or regression coefficients relating expression levels to clinical outcomes. Estimated gene-parameter profiles. Existing approaches typically assume known gene-parameter profiles. For example, the list of differentially expressed genes from a microarray experiment is usually treated as known and fixed in subsequent analyses with GO, while in fact it corresponds to an unknown and estimated parameter. Distinguishing between the definition of a parameter and inference concerning this parameter, as in Section 3, provides a more rigorous and general formulation of the statistical question. General tests of association. Common approaches to tests of association with GO annotation are typically limited to tests of independence in 2×2 contingency tables (e.g., based on the hypergeometric distribution, Fisher's exact test). As in Table 2, rows correspond to gene annotation with a given GO term (fixed binary gene-annotation profile) and columns to a gene property of interest, such as differential expression (treated as a fixed binary gene-parameter profile). The approach proposed in Section 3 allows consideration of a broader class of biological testing problems, while properly accounting for the fact that gene-parameter profiles are usually unknown and replaced by a random (i.e., data-driven) estimator.

Outline
This article proposes a general and formal statistical framework for multiple tests of association with biological annotation metadata, using the multiple hypothesis testing methodology developed in [14] and related articles. Section 2 provides an introduction to multiple hypothesis testing. Section 3 presents the proposed statistical framework for multiple tests of association with biological annotation metadata and discusses in detail the main components of the inference problem, namely, the gene-annotation profiles, the gene-parameter profiles, and the association parameters. Multiple testing procedures (MTP) for tests of association between gene-annotation profiles and gene-parameter profiles are outlined. Section 4 gives an overview of the Gene Ontology (GO) and R software for accessing and analyzing GO annotation metadata (e.g., for assembling GO gene-annotation profiles). The proposed statistical and computational methods are illustrated in Section 5, using the acute lymphoblastic leukemia (ALL) microarray dataset of [13], with the aim of relating GO annotation to differential gene expression between B-cell ALL with the BCR/ABL fusion and cytogenetically normal NEG B-cell ALL. Finally, Section 6 summarizes our findings and outlines ongoing work.

Multiple hypothesis testing
This section, based on Chapter 1 of [14], introduces a general statistical framework for multiple hypothesis testing and summarizes in turn the main ingredients of a multiple testing problem, including: the data generating distribution; the parameters of interest; the null and alternative hypotheses; the test statistics; the rejection regions (i.e., cut-offs) for the test statistics; errors in multiple hypothesis testing; Type I error rate and power; adjusted p-values.
The section also provides an overview of multiple testing procedures developed in [14,, for controlling generalized tail probability and expected value error rates, including the key choices of a joint null distribution and rejection regions for the test statistics.

Null and alternative hypotheses
Hypothesis testing is concerned with using observed data to make decisions regarding properties of (i.e., hypotheses for) the unknown data generating distribution.
Let X n ≡ {X i : i = 1, . . . , n} denote a random sample of n independent and identically distributed (IID) random variables from a data generating distribution P , i.e., X i IID ∼ P , i = 1, . . . , n. Suppose that the data generating distribution P is an element of a particular statistical model M, i.e., a set of possibly non-parametric distributions. Let P n denote the empirical distribution corresponding to the sample X n , i.e., the distribution which places probability 1/n on each realization of X.
In order to cover a broad class of testing problems, define M pairs of null and alternative hypotheses in terms of a collection of M submodels, M(m) ⊆ M, m = 1, . . . , M , for the data generating distribution P . Specifically, the M null hypotheses and corresponding alternative hypotheses are defined as (1) H 0 (m) ≡ I (P ∈ M(m)) and respectively. The general submodel representation accommodates tests of means, quantiles, covariances, correlation coefficients, and regression coefficients in linear and nonlinear models (e.g., logistic, survival, time-series models).
The hypothesized null values, ψ 0 (m), are frequently zero. Let denote the set of h 0 ≡ |H 0 | true null hypotheses, where the longer notation H 0 (P ) emphasizes the dependence of this set on the data generating distribution P . Likewise, let The goal of a multiple testing procedure is to accurately estimate, i.e., reject, the set H 1 , while probabilistically controlling false positives.

Test statistics
A testing procedure is a random or data-driven rule for estimating the set of false null hypotheses H 1 = {m : H 0 (m) = 0} = {m : P / ∈ M(m)}, i.e., for deciding which null hypotheses should be rejected.
The decisions to reject or not the null hypotheses are based on an M -vector of test statistics, T n = (T n (m) : m = 1, . . . , M ), that are functions T n (m) = T (m; X n ) = T (m; P n ) of the data X n , i.e., of the empirical distribution P n . Denote the typically unknown (finite sample) joint distribution of the test statistics T n by Q n = Q n (P ).
As in [14,Chapter 1], for the test of single-parameter null hypotheses of the form H 0 (m) = I (ψ(m) ≤ ψ 0 (m)) or H 0 (m) = I (ψ(m) = ψ 0 (m)), m = 1, . . . , M , consider two main types of test statistics, difference statistics, and t-statistics (i.e., standardized differences), Here and we may adopt the shorter notation σ 2 (m) = σ(m, m) = E[IC 2 (X|P )(m)] for variances. Assume that σ 2 n (m) are consistent estimators of the IC variances σ 2 (m). The influence curve of a given estimator can be derived as its mean-zero-centered functional derivative (as a function of the empirical distribution P n for the entire sample of size n), applied to the empirical distribution for a sample of size one [19,20].
This general representation for the test statistics covers standard one-sample and two-sample t-statistics for testing hypotheses concerning mean parameters, but also test statistics for correlation coefficients and regression coefficients in linear and nonlinear models. Test statistics for other types of null hypotheses include F -statistics, χ 2 -statistics, and likelihood ratio statistics.

Rejection regions
A multiple testing procedure (MTP) provides rejection regions C n (m), i.e., sets of values for each test statistic T n (m) that lead to the decision to reject the corresponding null hypothesis H 0 (m) and declare that P / ∈ M(m), m = 1, . . . , M . In other words, a MTP produces a random (i.e., data-dependent) set R n of rejected hypotheses that estimates the set H 1 of false null hypotheses, where C n (m) = C(m; T n , Q 0n , α), m = 1, . . . , M , denote possibly random rejection regions.
The long notation R(T n , Q 0n , α) and C(m; T n , Q 0n , α) emphasizes that the MTP depends on the following three ingredients: • the data, X n = {X i : i = 1, . . . , n}, through the M -vector of test statistics, T n = (T n (m) : m = 1, . . . , M ) (Section 2.2); • an (estimated) M -variate test statistics null distribution, Q 0n , which replaces the unknown true test statistics distribution Q n = Q n (P ) for the purpose of deriving rejection regions, confidence regions, and adjusted p-values (Section 2.7); • the nominal Type I error level α, i.e., a user-supplied upper bound for a suitably defined Type I error rate (Section 2.5).
Given a proper test statistics null distribution Q 0 (or estimator thereof, Q 0n ), the main task is to specify rejection regions for each null hypothesis, so that the resulting procedure probabilistically controls Type I errors. We consider MTPs based on nested rejection regions, that is, Rejection regions are typically defined in terms of intervals, such as, C n (m) = (u n (m), +∞), C n (m) = (−∞, l n (m)), or C n (m) = (−∞, l n (m)) ∪ (u n (m), +∞), where l n (m) = l(m; T n , Q 0n , α) and u n (m) = u(m; T n , Q 0n , α) are to-be-determined lower and upper critical values, or cut-offs, computed under the null distribution Q 0n for the test statistics T n . Two-sided rejection regions of the form C n (m) = (−∞, l n (m))∪(u n (m), +∞) allow the use of asymmetric cut-offs for two-sided tests.
Unless specified otherwise, we assume that large values of the test statistic T n (m) provide evidence against the corresponding null hypothesis H 0 (m), that is, we consider one-sided rejection regions of the form C n (m) = (c n (m), +∞), where c n (m) = c(m; T n , Q 0n , α). For two-sided tests of single-parameter null hypotheses using difference or t-statistics, as in Equations (6) and (7), one could take absolute values of the test statistics.
Among the different approaches for defining rejection regions, we distinguish the following.
Marginal vs. joint multiple testing procedures. In marginal multiple testing procedures, rejection regions are based solely on the marginal distributions of the test statistics (e.g., FWER-controlling single-step Bonferroni procedure [12]). In contrast, joint procedures take into account the dependence structure of the test statistics (e.g., FWER-controlling single-step maxT Procedure 1). Joint MTPs tend to be more powerful than marginal MTPs. Note that while a procedure may be marginal, proof of Type I error control by this MTP may require certain assumptions on the dependence structure of the test statistics (e.g., FWER-controlling step-up Hochberg procedure [23]). Single-step vs. stepwise multiple testing procedures. In single-step multiple testing procedures, each null hypothesis H 0 (m) is tested using a rejection region that is independent of the results of the tests of other hypotheses and is not a function of the data X n (unless these data are used to estimate the null distribution). In contrast, in stepwise procedures, the decision to reject a particular null hypothesis depends on the outcome of the tests of other hypotheses. That is, the (single-step) testing procedure is applied to a sequence of successively smaller nested random (i.e., data-dependent) subsets of null hypotheses, defined by the ordering of the test statistics (common-cut-off MTPs) or unadjusted p-values (common-quantile MTPs). The rejection regions are therefore allowed to depend on the data X n via the test statistics T n .
Stepwise procedures are of two main types, depending on the order in which the null hypotheses are tested. In step-down procedures, the most significant null hypotheses (i.e., the null hypotheses with the largest test statistics for commoncut-off MTPs or smallest unadjusted p-values for common-quantile MTPs) are considered successively, with further tests depending on the outcome of earlier ones. As soon as one fails to reject a null hypothesis, no further hypotheses are rejected. In contrast, for step-up procedures, the least significant null hypotheses are considered successively, again with further tests depending on the outcome of earlier ones. As soon as one null hypothesis is rejected, all remaining more significant hypotheses are rejected.
Stepwise MTPs tend to be more powerful than single-step MTPs. Common-cut-off vs. common-quantile multiple testing procedures. In common-cut-off multiple testing procedures, the same cut-off c 0 is used for each test statistic (e.g., FWER-controlling single-step and step-down maxT procedures, based on maxima of test statistics). In contrast, in common-quantile procedures, the cut-offs are the δ 0 -quantiles of the marginal null distributions of the test statistics (e.g., FWER-controlling single-step and step-down minP procedures, based on minima of unadjusted p-values). The latter p-value-based procedures place the null hypotheses on an "equal footing", i.e., are more balanced than their common-cut-off counterparts, and may therefore be preferable. However, this comes at the expense of increased computational complexity.

Errors in multiple hypothesis testing
In any testing problem, two types of errors can be committed. A Type I error, or false positive, is committed by rejecting a true null hypothesis (R n ∩ H 0 ). A Type II error, or false negative, is committed by failing to reject a false null hypothesis (R c n ∩ H 1 ). The situation can be summarized as in Table 1, where the number of rejected null hypotheses is the number of Type I errors or false positives is the number of Type II errors or false negatives is and the number of true positives is Note that S n , U n , V n , and W n each depend on the unknown data generating distribution P through the unknown set of true null hypotheses H 0 = H 0 (P ). Therefore, the numbers h 0 = |H 0 | and h 1 = |H 1 | = M − h 0 of true and false null hypotheses are unknown parameters (row margins of Table 1), the number of rejected hypotheses R n is an observable random variable (column margins of Table  1), and S n , U n , V n , and W n are unobservable random variables (cells of Table 1).
Ideally, one would like to simultaneously minimize both the number of Type I errors and the number of Type II errors. Unfortunately, this is not feasible and one seeks a trade-off between the two types of errors. A standard approach is to specify an acceptable level α for a suitably defined Type I error rate and derive testing procedures (i.e., rejection regions) that aim to minimize a Type II error rate (i.e., maximize power) within the class of tests with Type I error level at most α.

Type I error rate and power
When testing multiple hypotheses, there are many possible definitions for the Type I error rate and power of a testing procedure. Accordingly, we define a Type I error rate as a parameter θ n = Θ(F Vn,Rn ) of the joint distribution F Vn,Rn of the numbers of Type I errors V n = |R n ∩ H 0 | and rejected hypotheses R n = |R n |. Likewise, we define power as a parameter ϑ n = Θ(F Un,Rn ) of the joint distribution F Un,Rn of the numbers of Type II errors U n = |R c n ∩ H 1 | and rejected hypotheses R n = |R n |. We focus primarily on mappings such that θ n ∈ [0, 1] and ϑ n ∈ [0, 1].
This parametric representation covers a broad class of Type I error rates, including generalized tail probability (gTP) error rates, (16) gT P (q, g) ≡ Pr(g(V n , R n ) > q), and generalized expected value (gEV) error rates, for arbitrary functions g(V n , R n ) of the numbers of Type I errors V n and rejected hypotheses R n and user-supplied bounds q. Generalized tail probability and expected value error rates include as special cases the following commonly-used Type I error rates.
The generalized family-wise error rate (gFWER), corresponding to g(v, r) = v and q ∈ {0, . . . , (h 0 − 1)}, is the probability of at least (q + 1) Type I errors, When q = 0, the gFWER reduces to the usual family-wise error rate (FWER), controlled by the classical Bonferroni procedure. The tail probability for the proportion of false positives (TPPFP) among the rejected hypotheses, corresponding to g(v, r) = v/r and q ∈ (0, 1), is defined as with the convention that V n /R n ≡ 0 if R n = 0. The false discovery rate (FDR), corresponding to g(v, r) = v/r, is the expected proportion of false positives among the rejected hypotheses, again with the convention that V n /R n ≡ 0 if R n = 0.
Error rates Θ(F Vn/Rn ), based on the proportion of false positives (e.g., TPPFP and FDR), are especially appealing for the large-scale testing problems encountered in genomics, compared to error rates Θ(F Vn ), based on the number of false positives (e.g., gFWER), as they do not increase exponentially with the number M of tested hypotheses. However, error rates Θ(F Vn/Rn ) tend to be more difficult to control than error rates Θ(F Vn ), as they are based on the joint distribution of V n and R n , rather than only the marginal distribution of V n .

Adjusted p-values
As in the case of single hypothesis testing, one can report the results of a multiple testing procedure in terms of the following quantities: rejection regions for the test statistics, confidence regions for the parameters of interest, and adjusted p-values.
Adjusted p-values, for the test of multiple hypotheses, are defined as straightforward extensions of unadjusted p-values, for the test of individual hypotheses. Consider any multiple testing procedure R n (α) = R(T n , Q 0 , α), with rejection regions C n (m; α) = C(m; T n , Q 0 , α). That is, the adjusted p-value P 0n (m), for null hypothesis H 0 (m), is the smallest nominal Type I error level (e.g., gFWER, TPPFP, or FDR) of the multiple hypothesis testing procedure at which one would reject H 0 (m), given T n . As in single hypothesis tests, the smaller the adjusted p-value P 0n (m), the stronger the evidence against the corresponding null hypothesis H 0 (m). Thus, one rejects H 0 (m) for small adjusted p-values P 0n (m).
For instance, the adjusted p-values for the classical FWER-controlling marginal single-step common-quantile Bonferroni procedure are P 0n (m) = min {M P 0n (m), 1}. Adjusted p-values for FWER-controlling joint single-step common-cut-off maxT Procedure 1 are given in Equation (26).
Under the nestedness assumption of Equation (10), one has two equivalent representations for a MTP, in terms of rejection regions for the test statistics and in terms of adjusted p-values. Specifically, the set of rejected null hypotheses at multiple test nominal Type I error level α is Reporting the results of a MTP in terms of adjusted p-values, as opposed to only rejection or not of the null hypotheses, offers several advantages, including the following.
• Adjusted p-values can be defined for any Type I error rate (e.g., gFWER, TPPFP, or FDR). • They reflect the strength of the evidence against each null hypothesis in terms of the Type I error rate for the entire MTP. • They are flexible summaries of a MTP, in the sense that results are supplied for all Type I error levels α, i.e., the level α need not be chosen ahead of time. • They provide convenient benchmarks to compare different MTPs, whereby smaller adjusted p-values indicate a less conservative procedure. • Plots of sorted adjusted p-values allow investigators to examine sets of rejected hypotheses associated with various Type I error rates (e.g., gFWER, TPPFP, or FDR) and nominal levels α. Such plots provide tools to decide on an appropriate combination of the number of rejected hypotheses and tolerable false positive rate for a particular experiment and available resources.

Test statistics null distribution
As detailed in Chapter 2 of [14], a key feature of our proposed multiple testing procedures is the test statistics null distribution (rather than data generating null distribution) used to obtain rejection regions for the test statistics, confidence regions for the parameters of interest, and adjusted p-values. Indeed, whether testing single or multiple hypotheses, one needs the (joint) distribution of the test statistics in order to derive a procedure that probabilistically controls Type I errors. In practice, however, the true distribution Q n = Q n (P ) of the test statistics T n is unknown and replaced by a null distribution Q 0 . The choice of a proper null distribution is crucial, in order to ensure that (finite sample or asymptotic) control of the Type I error rate under the assumed null distribution does indeed provide the desired control under the true distribution. This issue is particularly relevant for large-scale testing problems, such as those involving biological annotation metadata, which concern high-dimensional multivariate distributions, with complex and unknown dependence structures among variables.
Chapter 2 of [14] provides a general characterization of a proper test statistics null distribution in terms of null domination conditions for the joint distribution of the test statistics (T n (m) : m ∈ H 0 ) for the true null hypotheses H 0 (Section 2.2). This general characterization leads to the explicit proposal of the following two main types of test statistics null distributions: a null shift and scale-transformed test statistics null distribution, based on user-supplied upper bounds for the means and variances of the test statistics for the true null hypotheses (Section 2.3), and a null quantile-transformed test statistics null distribution, based on user-supplied marginal test statistics null distributions (Section 2.4).
In practice, the test statistics null distribution Q 0 = Q 0 (P ) is unknown, as it depends on the unknown data generating distribution P . Resampling procedures based on the bootstrap are provided to conveniently obtain consistent estimators of the null distribution and of the corresponding test statistic cut-offs, parameter confidence regions, and adjusted p-values [14 As argued in [14,Chapter 2], the following two main points distinguish our approach from existing approaches to Type I error control and the choice of a test statistics null distribution (e.g., [24] and [45]). Firstly, we are only concerned with control of the Type I error rate under the true data generating distribution P , i.e., under the joint distribution Q n = Q n (P ), implied by P , for the test statistics T n . The concepts of weak and strong control of a Type I error rate and the related restrictive assumption of subset pivotality are therefore irrelevant in our context [45, p. 9-10, 42-43]. Secondly, we propose a null distribution for the test statistics (T n ∼ Q 0 ) rather than a data generating null distribution (X ∼ P 0 ∈ ∩ M m=1 M(m)). The latter practice does not necessarily provide proper Type I error control under the true distribution P . Indeed, the test statistics assumed null distribution Q n (P 0 ) and their true distribution Q n (P ) may have different dependence structures for the true null hypotheses H 0 and, as a result, may violate the required null domination conditions for Type I error control.
We stress the generality of our proposed test statistics null distributions: Type I error control does not rely on restrictive assumptions such as subset pivotality and holds for general data generating distributions (with arbitrary dependence structures among variables), null hypotheses (defined in terms of submodels for the data generating distribution), and test statistics (e.g., t-statistics, χ 2 -statistics, F -statistics).

Null shift and scale-transformed test statistics null distribution
The first original null distribution of [16,33,41], is defined as the asymptotic distribution Q 0 = Q 0 (P ) of the M -vector Z n of null shift and scale-transformed test statistics, where λ 0 (m) and τ 0 (m) are, respectively, user-supplied upper bounds for the means and variances of the H 0 -specific test statistics. In this construction, the null shift values λ 0 (m) are chosen so that the H 0 -specific statistics (Z n (m) : m ∈ H 0 ) are asymptotically stochastically greater than the original test statistics (T n (m) : m ∈ H 0 ). The resulting null distribution therefore satisfies the required null domination conditions for Type I error control.
In contrast, the null scale values τ 0 (m) are not needed for Type I error control. The purpose of τ 0 (m) is to avoid a degenerate null distribution and infinite cut-offs for the false null hypotheses (m ∈ H 1 ), an important property for power considerations. This scaling is needed, in particular, for F -statistics that have asymptotically infinite means and variances for non-local alternative hypotheses.
For a broad class of testing problems, such as the test of single-parameter null hypotheses using t-statistics (Equation (7)), the null distribution Q 0 is an M -variate Gaussian distribution, with mean vector zero and covariance matrix σ * = Σ * (P ), that is, Q 0 = N(0, σ * ). For tests where the parameter of interest is the M -dimensional mean vector Ψ(P ) = ψ = E[X], the estimator ψ n is simply the M -vector of empirical means and σ * = Σ * (P ) = Cor[X] is the correlation matrix of X ∼ P , that is, Q 0 = N(0, Cor[X]). More generally, for an asymptotically linear estimator ψ n , Σ * (P ) is the correlation matrix of the vector influence curve. This situation covers standard one-sample and two-sample t-statistics for tests of means, but also test statistics for correlation coefficients and regression coefficients in linear and non-linear models.
For testing the equality of K population mean vectors using F -statistics, an F -statistic-specific null distribution may be defined as the joint distribution of an M -vector of quadratic forms of Gaussian random variables.

Null quantile-transformed test statistics null distribution
The second and most recent proposal of [42] is defined as the asymptotic distribution Q 0 = Q 0 (P ) of the M -vectorZ n of null quantile-transformed test statistics, where q 0,m are user-supplied marginal test statistics null distributions that satisfy marginal null domination conditions. According to the generalized quantile-quantile function transformation of [46], define Q ∆ n,m (z) ≡ ∆Q n,m (z) + (1 − ∆)Q n,m (z − ), where Q n,m are the marginal distributions of the test statistics T n (m) and the random variable ∆ is uniformly distributed on the interval [0, 1], independently of the data X n .
This latest proposal has the additional advantage that the marginal test statistics null distributions may be set to the optimal (i.e., most powerful) null distributions one would use in single hypothesis testing (e.g., permutation marginal null distributions, Gaussian or other parametric marginal null distributions).

Overview of multiple testing procedures
Having identified a suitable test statistics null distribution Q 0 (or estimator thereof, Q 0n ), there remains the main task of specifying rejection regions (i.e., cut-offs) for the test statistics, confidence regions for the parameters of interest, and adjusted p-values.
As detailed in [14,, we have developed resampling-based single-step and stepwise multiple testing procedures for controlling a broad class of Type I error rates, defined as generalized tail probabilities, gT P (q, g) = Pr(g(V n , R n ) > q), and generalized expected values, gEV (g) = E[g(V n , R n )], for arbitrary functions g(V n , R n ) of the numbers of false positives V n and rejected hypotheses R n . Our proposed procedures take into account the joint distribution of the test statistics and provide Type I error control in testing problems involving general data generating distributions (with arbitrary dependence structures among variables), null hypotheses (defined in terms of submodels for the data generating distribution), and test statistics (e.g., t-statistics, χ 2 -statistics, F -statistics).
An overview of available MTPs is provided in Chapter 3 of [14]. Core methodological Chapters 4-7 discuss the following main approaches for deriving rejection regions.
Joint single-step common-cut-off and common-quantile procedures for controlling general Type I error rates Θ(F Vn ), defined as arbitrary parameters of the distribution of the number of Type I errors V n (Chapter 4 in [14], [16,33]). Error rates of the form Θ(F Vn ) include the generalized family-wise error rate, gF W ER(q) = 1 − F Vn (q) = Pr(V n > q), i.e., the chance of at least (q + 1) Type I errors.

FWER-controlling single-step common-cut-off maxT procedure
This section focusses on control of the family-wise error rate, using the singlestep maxT procedure, a common-cut-off procedure exploiting the joint distribution of the test statistics. The method is summarized below; details are given in [14, Chapter 4] and [16].
Procedure 1 [FWER-controlling single-step common-cut-off maxT procedure]. Given an M -variate test statistics null distribution Q 0 , the single-step commoncut-off maxT procedure is based on the distribution of the maximum test statistic, The adjusted p-value p 0n (m) for null hypothesis H 0 (m) is the probability, under Q 0 , that max m Z(m) exceeds the corresponding observed test statistic t n (m), that is, Procedure 1 provides proper FWER control when based on either of the two nulltransformed test statistics null distributions Q 0 introduced in Section 2.7. Consistent estimators Q 0n of the null distribution Q 0 and corresponding single-step maxT cut-offs and adjusted p-values may be obtained using the bootstrap, as in Procedure 2.9, below, for the null shift and scale-transformed test statistics null distribution [14,Section 4.4].
1. Let P ⋆ n denote a bootstrap estimator of the data generating distribution P . For the non-parametric bootstrap, P ⋆ n is simply the empirical distribution P n , that is, samples of size n are drawn at random, with replacement from the observed data X n = {X i : i = 1, . . . , n}. For the model-based bootstrap, P ⋆ n belongs to a model M for the data generating distribution P , such as a family of multivariate Gaussian distributions. That is, compute

Obtain an
For t-statistics as in Equation (7), the null values are λ 0 (m) = 0 and τ 0 (m) = 1. 6. The bootstrap estimator Q 0n of the null shift and scale-transformed null dis- 8. For controlling the FWER at nominal level α ∈ (0, 1), the bootstrap single- . . , B} that exceed the corresponding observed test statistic t n (m), that is,

Statistical framework for multiple tests of association with biological annotation metadata
Sections 3.1-3.3 introduce the main components of our approach to multiple tests of association with biological annotation metadata, namely, the gene-annotation profiles A, the gene-parameter profiles λ, and the association measures ψ = ρ(A, λ) between gene-annotation and gene-parameter profiles. We stress that the choice of a suitable association parameter ψ is perhaps the most important and hardest aspect of the inference problem, as this parameter represents the statistical translation of the biological question of interest. Once the association parameter ψ is appropriately and precisely defined, one can rely on a variety of statistical methods to estimate and test hypotheses concerning this parameter. Section 3.4 describes how the multiple testing methodology of [14] and related articles may be used to detect associations between gene-annotation and gene-parameter profiles. Note that, for the sake of illustration, we focus on gene-level features. However, as mentioned in Section 1.1, the methodology is generic and may be applied to other types of features, such as those concerning gene isoforms and proteins.

Gene-annotation profiles
Gene-annotation profiles refer to features of a genome that are assumed to be known and constant among units in a population of interest. Such features typically consist of gene annotation metadata, that reflect current knowledge on gene properties, such as, nucleotide and protein sequences, regulation, and function.
Specifically In many applications, the element A(g, m) is a binary indicator, coding the YES/NO answer to the mth question, among a collection of M questions one may ask about gene g. For example, A(g, m) could indicate whether gene g is annotated with a particular GO term m, among M terms in one of the three ontologies (BP, CC, or MF), i.e., whether gene g is an element of the node corresponding to the mth term in the GO directed acyclic graph (DAG). Other gene-annotation profiles of interest may refer to exon/intron counts/lengths/nucleotide distributions, gene pathway membership (e.g., from the Kyoto Encyclopedia of Genes and Genomes, KEGG, www.genome.ad.jp/kegg), or gene regulation by particular transcription factors. Regarding transcription regulation, one could use data from the Transcription Factor DataBase (TRANSFAC, www.gene-regulation.com) to generate gene-annotation profiles as follows. For a given transcription factor binding motif, a binary gene-annotation profile could consist of indicators for the presence or absence of the motif in the upstream control region of each gene. A continuous gene-annotation profile could be based on the position weight matrix of the binding motif.
Note that the aforementioned features are only fixed in time for a given version/release of the corresponding database(s), i.e., such biological data are constantly evolving as our knowledge of the roles of genes and proteins is accumulating and changing. The dynamic nature of biological annotation metadata is an important issue in terms of software design (Section 4.2; [18]).
Note also that gene-annotation profiles are not restricted to be binary or even polychotomous and, in particular, could be continuous gene-parameter profiles, suitably estimated from previous studies.
The main point, regarding the formulation of the statistical inference question, is that gene-annotation profiles are known and constant among population units.

Gene-parameter profiles
Gene-parameter profiles are generally unknown and concern the distribution of variable features of a genome in a well-defined population. Gene-specific variables of interest reflect cellular type/state/activity under particular conditions and include microarray measures of transcript levels and comparative genomic hybridization (CGH) measures of DNA copy numbers.
Specifically, let X = (X(j) : j = 1, . . . , J) be a J-dimensional random vector, containing G gene-specific random variables (X(g) : g = 1, . . . , G). In addition to the G gene-specific variables, X may include various biological and clinical covariates (e.g., age, sex, treatment, timepoint) and outcomes (e.g., survival time, response to treatment, tumor class). Let P denote the data generating distribution for the random J-vector X and suppose that P belongs to a (possibly non-parametric) model M.
Let the parameter mapping Λ : is a gene-specific real-valued parameter. For example, λ(g) could be the mean expression measure E[X(g)] of gene g or a regression coefficient relating an outcome component of X to the expression measure X(g) of gene g, g = 1, . . . , G.
While gene-annotation profiles are known and fixed, gene-parameter profiles are typically unknown and need to be estimated, e.g., from a microarray experiment involving a sample of population units. The sample X n = {X i : i = 1, . . . , n} is assumed to consist of n independent and identically distributed (IID) copies of X ∼ P , corresponding to n randomly sampled population units.

Association measures for gene-annotation and gene-parameter profiles
Let the parameter mapping Ψ : M → R M specify an M -dimensional association parameter vector, defined in terms of an association measure ρ : R G×M × R G → R M , known fixed gene-annotation profiles A, and an unknown gene-parameter profile λ = Λ(P ). The choice of a suitable association parameter is subject matter-dependent and requires careful consideration. For instance, for Gene Ontology annotation, it is desirable that the association parameter reflect the structure of the GO directed acyclic graph (Section 4.1). In principle, the dimension of the association parameter vector ψ could differ from the number M of features under consideration. In addition, one could accommodate several gene-parameter profiles λ.
The various quantities in the inference problem are summarized in Figure 1; examples of association parameters are given next and in Section 5.

Univariate association measures
In the simplest case, one could define the M association parameters univariately, i.e., define ψ(m) based only on the mth gene-annotation profile A(·, m), m = 1, . . . , M . Specifically, for the mth feature, let where ρ m : R G × R G → R provides a measure of association (e.g., correlation coefficient) between the G-dimensional gene-annotation profile A(·, m) and geneparameter profile λ = Λ(P ). In many situations, the same association measure ρ m may be used for each of the M features.
Continuous gene-annotation profiles and continuous gene-parameter profiles For continuous gene-annotation and gene-parameter profiles, one may use as association measure the Pearson correlation coefficient between two G-vectors. That is, whereĀ(m) ≡ g A(g, m)/G andλ ≡ g λ(g)/G denote, respectively, the averages of the G elements of the gene-annotation profile A(·, m) and gene-parameter profile λ.
Binary gene-annotation profiles and binary gene-parameter profiles For binary gene-annotation and gene-parameter profiles, one may build 2 × 2 contingency Table 2 and use as association measure the χ 2 -statistic (or corresponding p-value) for the test of independence of rows and columns. That is, 1}. Note that in this context the χ 2 -statistic ψ(m) is a parameter, i.e., it is a function of the data generating distribution P , via the gene-parameter profile λ = Λ(P ), and is therefore unknown and constant among population units. Table 2. Binary gene-annotation and gene-parameter profiles. Given a binary gene-annotation profile A(·, m) and a binary gene-parameter profile λ, one may build a 2 × 2 contingency table, with rows corresponding to the gene-annotation profile and columns to the gene-parameter profile. Cell counts are defined as . For example, for tests of association between GO annotation and differential gene expression, g 11 (m) could correspond to the number of genes that are annotated with GO term m and differentially expressed.
Binary gene-annotation profiles For binary gene-annotation profiles, one may consider association parameter vectors of the form That is, the association parameter for the mth feature is the sum, of the parameters λ(g) for genes g that have the property of interest, i.e., such that A(g, m) = 1. Such an association parameter is considered by [37], to relate continuous microarray differential expression gene-parameter profiles to binary pathway gene-annotation profiles.
The following standardized association parameters, corresponding to association measures based on two-sample t-statistics, may also be considered, denote, respectively, the numbers, averages, and variances of annotated (k = 1) and unannotated (k = 0) gene-parameters λ(g).
In commonly-encountered combined GO annotation and microarray data analyses, a binary gene-parameter profile could indicate whether genes are differentially expressed in two populations of cells, a continuous gene-parameter profile could consist of coefficients for the regression of a (censored) clinical outcome on gene expression measures, and binary gene-annotation profiles could denote whether genes are annotated with particular GO terms (Section 5; [1,2,4,22]).

Multivariate association measures
More generally, the mth association parameter could be based on the entire geneannotation matrix A or a subset of columns thereof, that is, Ψ(P )(m) = ψ(m) ≡ ρ m (A, Λ(P )), for an association measure ρ m : Association parameters of interest include: linear combinations of association parameters for several features, partial correlation coefficients, χ 2 -statistics for higher-dimensional contingency tables (e.g., with one dimension corresponding to a gene-parameter profile λ and other dimensions to several gene-annotation profiles A(·, m)), and (contrasts of) regression coefficients of a gene-parameter profile λ on several gene-annotation profiles A(·, m).
In the case of Gene Ontology annotation, the association parameter ψ should preferably reflect the structure of the GO directed acyclic graph, by taking into account, for instance, annotation information for ancestor (i.e., less specific) or offspring (i.e., more specific) terms (Section 4.1). Specifically, let P(m) denote the set of (immediate) parents of a term m. As the genes annotated by the child term m are subsets of the genes annotated by the parent terms P(m), then A(g, m) = 1 implies A(g, p) = 1 for p ∈ P(m).
Following the causal inference literature [38,43], an association parameter of interest for GO term m is the marginal causal effect parameter, defined as where A(·, P(m)) denotes the submatrix of gene-annotation profiles for parent terms P(m) and the expected values are defined with respect to the empirical distribution of {(A(g, m), A(g, P(m)), λ(g)) : g = 1, . . . , G}.
In the special case of binary (differential expression) gene-parameter profiles, the so-called parent-child method of [22] takes into account the structure of the GO DAG by testing for associations between gene-annotation and gene-parameter profiles using hypergeometric p-values computed conditionally on the annotation status of parent terms.
One could also consider Boolean combinations of annotation indicators for multiple features, that is, a transformed gene-annotation matrix whose columns are Boolean combinations of the columns of the original gene-annotation matrix. Such an approach would be particularly relevant in the context of transcription regulation, where individual features correspond to single transcription factor binding motifs and Boolean combinations to binding modules for multiple transcription factors.

Null and alternative hypotheses
Certain biological annotation metadata analyses may involve the two-sided test of the M null hypotheses of no association between gene-annotation profiles A(·, m), m = 1, . . . , M , and a gene-parameter profile λ, i.e., the test of (38) H 0 (m) = I (ψ(m) = ψ 0 (m)) vs.
Other analyses may call for the one-sided test of One-sided tests are appropriate if, for example, one is interested in determining whether differentially expressed genes are enriched regarding GO annotation. Note that in many situations, the same association measure ρ m is used for each of the M features and one only has a single, common null value ψ 0 (m).

Test statistics
As in Section 2.2, above, and Chapter 1 of [14], consider the general situation where, given a random sample X n from the data generating distribution P , one has an asymptotically linear estimator ψ n =Ψ(P n ) of the association parameter vector ψ = Ψ(P ), with M -dimensional vector influence curve IC(X|P ). LetΣ(P n ) = σ n = (σ n (m, m ′ ) : m, m ′ = 1, . . . , M ) denote a consistent estimator of the covariance matrix Σ(P ) = σ = (σ(m, m ′ ) : m, m ′ = 1, . . . , M ) of the vector influence curve IC(X|P ). For example, σ n could be a bootstrap-based estimator of the covariance matrix σ or could be computed from an estimator IC n (X) of the influence curve IC(X|P ).
A broad range of association parameters ψ and corresponding estimators ψ n satisfy the above conditions. In particular, suppose λ n =Λ(P n ) is an asymptotically linear estimator of the gene-parameter profile λ = Λ(P ), based on a random sample X n from P . Let ψ n ≡ ρ(A, λ n ) denote the corresponding resubstitution, or plug-in, estimator of the association parameter vector ψ = ρ(A, λ). Then, if the function ρ(A, λ) is differentiable with respect to λ, the resubstitution estimator ψ n is also asymptotically linear.
One can therefore handle tests where the gene-parameter profiles λ are (functions of) means, variances, correlation coefficients, and regression coefficients, and where the association measures ρ are correlation coefficients, two-sample t-statistics, and χ 2 -statistics. Examples are provided in Section 5, in the context of tests of association between differential gene expression in ALL and GO annotation.
Each null hypothesis H 0 (m) may then be tested using a (unstandardized) difference statistic, or a (standardized) t-statistic, where we adopt the shorter notation σ 2 n (m) = σ n (m, m) for variances. Certain testing problems may call for other test statistics T n , such as, F -statistics, χ 2 -statistics, and likelihood ratio statistics.
Let Q n = Q n (P ) denote the typically unknown (finite sample) joint distribution of the M -vector of test statistics T n = (T n (m) : m = 1, . . . , M ), under the data generating distribution P .

Multiple testing procedures
As mentioned in Section 2.7, above, and detailed in [14, Chapter 2], a key feature of our proposed multiple testing procedures is the test statistics null distribution Q 0 used in place of the unknown true test statistics distribution Q n = Q n (P ), for the purpose of obtaining rejection regions for the test statistics, confidence regions for the parameters of interest, and adjusted p-values.
Given a suitable test statistics null distribution Q 0 (or estimator thereof, Q 0n ), the multiple testing procedures developed in [14] and related articles may be applied to control a broad class of Type I error rates, defined as generalized tail probabilities, gT P (q, g) = Pr(g(V n , R n ) > q), and generalized expected values, gEV (g) = E[g(V n , R n )], for arbitrary functions g(V n , R n ) of the numbers of false positives V n and rejected hypotheses R n (Section 2.8).
For the purpose of illustration, we focus, as in Section 2.9, on control of the family-wise error rate, using the single-step common-cut-off maxT procedure, based on a non-parametric bootstrap estimator of the null shift and scale-transformed test statistics null distribution (Procedures 1 and 2.9).

Overview of the Gene Ontology
The Gene Ontology (GO) Consortium (www.geneontology.org) provides ontologies, i.e., structured and controlled vocabularies, to describe gene products in terms of their associated biological processes, cellular components, and molecular functions. The ontologies specify terminologies and relationships among terms. They are organism-independent and can be applied even as our knowledge of the roles of genes and proteins is accumulating and changing.
The GO Consortium and other organizations supply mappings between GO terms and genes in various organisms.
Detailed documentation is available on the Gene Ontology Documentation webpage (www.geneontology.org/GO.contents.doc.html).

The three gene ontologies: BP, CC, and MF
The GO Consortium provides three ontologies, each consisting of a structured network of terms describing gene products. A gene product may be used in one or more biological processes, may be associated with one or more cellular components, and may have one or more molecular functions.

GO directed acyclic graphs
For each of the three gene ontologies, GO terms are organized in a directed acyclic graph (DAG), where a directed graph has one-way edges and an acyclic graph has no path starting and ending at the same vertex. Each GO term is associated with a single vertex, or node, in the DAG. The words term, node, and vertex, may therefore be used interchangeably. For a given GO term, an ancestor refers to a less specialized term; an offspring refers to a more specialized term. A parent is an immediate/direct ancestor of a term; a child is an immediate/direct offspring of a term. A root node has no parents (i.e., no incoming edges); a leaf node has no children (i.e., no outgoing edges). In a DAG, a child may have several parents.
GO terms must obey the so-called true path rule: if a (child) term describes a gene product, then all its immediate parent and more distant ancestor terms must also apply to the gene product.
The DAG structure of GO terms and corresponding true path rule are germane to the definition of a suitable association measure between gene-annotation profiles and gene-parameter profiles (Section 3.3). Furthermore, as discussed in Sections 4.2-4.5, in the context of Bioconductor annotation software, the true path rule is also relevant when assembling gene-annotation matrices.

GO software tools
Many software tools have been developed to deal with GO annotation metadata. The Gene Ontology Tools webpage (www.geneontology.org/GO.tools.shtml) provides a list of consortium and non-consortium software for searching and browsing the three gene ontologies, for annotating genes and gene products using GO, and for combined GO and gene expression microarray data analysis.
For instance, the AmiGO browser (www.godatabase.org) allows: searching for a GO term and viewing all gene products annotated with this term; searching for a gene product and viewing all its associated GO terms; and browsing the ontologies to view relationships among terms and gene products annotated with a given term.
The QuickGO browser (www.ebi.ac.uk/ego), developed by the European Bioinformatics Institute (EBI), permits searches and graphical displays of the Gene Ontology by GO term, GO term identifier (ID), gene product, and other identifiers.
Software packages developed as part of the Bioconductor Project are discussed in Sections 4.2-4.5. Example 2. GO term protein-tyrosine kinase activity. To get a sense of the information provided by the GO Consortium, consider the Molecular Function ontology and the GO term protein-tyrosine kinase activity, with GO term ID GO:0004713.
Go to the AmiGO browser (www.godatabase.org), enter the GO term ID GO:0004713 in the Search GO box, select Exact Match, select Terms, and click on the Submit Query button. There are two main options for displaying information on a GO term: a "tree view" and a "graphical view". Click on the small tree-like icon (top-left corner of the table) to display the tree view with all ancestors (i.e., less specific terms) of the GO term protein-tyrosine kinase activity. Click on the Graphical View button to display the portion of the MF DAG corresponding to the GO term. Additional information may be obtained by clicking on the hyperlinked text protein-tyrosine kinase activity.
The GO term protein-tyrosine kinase activity has one (immediate) parent, protein kinase activity (GO:0004672), which itself has two parents, kinase activity (GO:0016301) and phosphotransferase activity, alcohol group as acceptor (GO:0016773). Altogether, the term protein-tyrosine kinase activity has 7 ancestors. According to the true path rule, any gene annotated with the GO term proteintyrosine kinase activity should also be annotated with all of its less specific ancestor terms.
The portion of the MF DAG for the GO term protein-tyrosine kinase activity is displayed in Figure 2 using the QuickGO browser. Section 4.5 provides sample R code for assembling GO gene-annotation matrices using Bioconductor annotation metadata packages.

GO gene-annotation matrices
As detailed in Section 3, detecting associations between GO annotation and other interesting features of a genome may be viewed as the multiple test of the null hypotheses of no association between gene-annotation profiles A(·, m) and a geneparameter profile λ = Λ(P ). The multiple testing methodology proposed in [14] and related articles is well-suited to handle the complex and unknown dependence structure among test statistics implied by the DAG structure of GO terms. The methodology is summarized in Section 2 and illustrated in Section 5, for tests of association between differential gene expression in ALL and GO annotation.

Overview of R and Bioconductor software for GO annotation metadata analysis
As discussed in [18], the Bioconductor Project provides R packages for accessing and performing statistical inference with GO annotation metadata (www.bioconductor.org; www.r-project.org). The packages include: • a general annotation software package: annotate; • packages for graph theoretical analyses: e.g., graph, Rgraphviz; • a GO-specific metadata package for navigating the three GO DAGs: GO; • an Entrez Gene 1 -specific metadata package, providing bi-directional mappings between Entrez Gene IDs and GO term IDs: humanLLMappings (www.ncbi.nlm.nih.gov/entrez/query.fcgi?db=gene); • various Affymetrix chip-specific metadata packages, providing bi-directional mappings between Affymetrix probe 2 IDs and GO term IDs: e.g., hgu95av2, hu6800 (www.affymetrix.com); • a package for annotating and generating HTML reports for Affymetrix chip data: annaffy.
Bioconductor metadata packages are updated regularly to reflect the evolving nature of biological annotation metadata; it is therefore crucial to keep track of version numbers. For information on Bioconductor software, please consult [17] and the Documentation (www.bioconductor.org/docs) and Workshops (www.bioconductor.org/workshops) sections of the Bioconductor Project website, in addition to the standard R help facilities (e.g., help function, manuals, etc.).
The remainder of this section provides sample R code demonstrating Bioconductor software (results reported for R Release 2.2.1 and Bioconductor Release 1.7). In order to run through the examples, one needs to install and load the Bioconductor packages annotate, GO, and hgu95av2. The annotation metadata used in the examples correspond to the following package versions. Accessing and analyzing annotation metadata from databases such as GenBank (www.ncbi.nlm.nih.gov/Genbank), GO (www.geneontology.org), and PubMed (www.pubmed.gov), presupposes the ability to perform the following essential bookkeeping task: mapping between different identifiers (ID) for a given gene/probe. Bioconductor annotation metadata packages consist of environment objects that provide key-value mappings between different sets of gene/probe identifiers.

> library(annotate) > library(GO) > library(hgu95av2) > > packageDescription("annotate")$Version
For instance, in the annotation metadata package hgu95av2, for the Affymetrix chip series HG-U95Av2, the hgu95av2PMID environment provides mappings from Affymetrix probe IDs (keys) to PubMed IDs (values); similarly, the hgu95av2GO environment provides mappings from Affymetrix probe IDs (keys) to GO term IDs (values). Example 3. Affymetrix probe ID 1635 at.As of Version 1.10.0 of the hgu95av2 package, the Affymetrix probe with ID 1635 at corresponds to the gene with symbol ABL1 and long name v-abl Abelson murine leukemia viral oncogene homolog 1, located on the long arm of chromosome 9. This probe maps to one GenBank accession number, one Entrez Gene ID, 14 distinct GO term IDs, and 160 distinct PubMed IDs.

The annotation metadata package GO
The GO package provides environment objects containing key-value pairs for mappings between GO term IDs, GO terms, GO term ancestors, GO term parents, GO term children, GO  For information on any of the GO environments, use the help function, e.g., help(GOTERM) or ?GOBPPARENTS. For instance, the environment GOTERM provides mappings from GO term IDs (keys) to GO terms (values); the environments GOBPPARENTS, GOCCPARENTS, and GOMFPARENTS, provide ontology-specific mappings from GO term IDs (keys) to GO term parent IDs (values). The environments GOALLLOCUSID, GOLOCUSID2GO, and GOLOCUSID, provide mappings between GO term IDs and Entrez Gene IDs and are used in Section 4.5, below, to assemble an Entrez Gene ID-by-GO term ID gene-annotation matrix for the MF gene ontology. Example 4. GO term ID GO:0004713. Let us use the GO package to obtain information on (all) ancestors, (immediate) parents, (immediate) children, and (all) offspring of the term corresponding to the GO term ID GO:0004713.

Affymetrix chip-specific annotation metadata packages: The hgu95av2 package
The Bioconductor Project provides Affymetrix chip-specific annotation metadata packages for the main chip series for the human, mouse, rat, and other genomes (e.g., HG-U133, HG-U95, HU-6800, MG-U74, and RG-U34 series). These packages, built using the infrastructure package AnnBuilder, contain environment objects for mappings between Affymetrix probe IDs and other types of gene/probe identifiers. Note that analogous packages are not supplied for two-color spotted microarrays, as there is no standard microarray design for this type of platform and specialized annotation metadata packages may have to be created for each microarray facility (e.g., using AnnBuilder). Once annotation metadata packages are available to provide mappings between different sets of gene/probe identifiers, the tools in annotate and related packages may be used in a similar manner for any type of microarray platform.
Consider the hgu95av2 package, for the Affymetrix chip series HG-U95Av2. This package provides the following environments. For information on any of these environments, use the help function, e.g., help(hgu95av2GO) or ?hgu95av2GO. We focus on the three environments related to GO: hgu95av2GO, hgu95av2GO2ALLPROBES, and hgu95av2GO2PROBE.
The HG-U95Av2 chip contains 12,625 probes (corresponding to the keys in the hgu95av2GO environment), with the first 10 Affymetrix probe IDs listed below.

Probes-to-most specific GO terms mappings: The hgu95av2GO environment
The hgu95av2GO environment provides key-value pairs for the mappings from Affymetrix probe IDs (keys) to GO term IDs (values). Each Affymetrix probe ID is mapped to a list of one or more elements, where each element corresponds to a particular GO term and is itself a list with the following three elements.
• "GOID": A GO term ID corresponding to the Affymetrix probe ID (key).
• "Evidence": A code for the evidence supporting the association of the GO term to the Affymetrix probe. • "Ontology": An abbreviation for the name of the ontology to which the GO term belongs: BP (Biological Process), CC (Cellular Component), or MF (Molecular Function).
Note that only the directly associated terms or most specific terms (i.e., not their less specific ancestor terms) a probe is annotated with are returned as values in hgu95av2GO. The GO package may be used to obtain more information on the GO term IDs, e.g., GO term, (all) ancestors, (immediate) parents, (immediate) children, and (all) offspring (Section 4.3). Example 5. GO terms directly associated with Affymetrix probe ID 1635 at. Let us obtain GO annotation information for the probe with Affymetrix ID 1635 at, corresponding to the ABL1 gene. The code below shows that probe 1635 at is directly annotated with 14 distinct GO terms (the same GO term ID may be returned multiple times with a different evidence code). As already noted in Example 1, one of these terms, with GO term ID GO:0004713, is protein-tyrosine kinase activity, in the Molecular Function ontology.
> probe <-"1635 at" > probe2GO <-get(probe, env = hgu95av2GO) > length(probe2GO) [1] 14 > unique(names(probe2GO)) [1] "GO:0000074" "GO:0000115" "GO:0000166" "GO:0003677" [5] "GO:0004713" "GO:0005515" "GO:0005524" "GO:0005634" [9] "GO:0006298" "GO:0006355" "GO:0006468" "GO:0007242" [13]  The hgu95av2GO environment (and analogous environments for other chip series) may be used to assemble an Affymetrix probe ID-by-GO term ID gene-annotation matrix, row by row. This may entail, however, a number of data processing steps. Firstly, only the most specific terms a probe is annotated with are returned as values in hgu95av2GO. One therefore needs to add all (less specific) ancestor terms in order to comply with the true path rule. Secondly, several probes may correspond to the same gene, i.e., several Affymetrix probe IDs may map to the same Entrez Gene ID according to the environment hgu95av2LOCUSID. Thirdly, the hgu95av2GO environment returns GO terms for all three gene ontologies at once. One may need to separate terms according to membership in the BP, CC, and MF ontologies (e.g., using the GOTERM environment from the GO package).
Alternately, one may assemble an Affymetrix probe ID-by-GO term ID geneannotation matrix, column by column, using the hgu95av2GO2ALLPROBES environment described below.

GO terms-to-directly annotated probes mappings:
The hgu95av2GO2PROBE environment The hgu95av2GO2PROBE environment provides key-value pairs for the mappings from GO term IDs (keys) to Affymetrix probe IDs (values). Values are vectors of length one or greater depending on whether a given GO term ID is mapped to one or more Affymetrix probe IDs. The value names are evidence codes for the GO term IDs. Note that the probes a particular GO term is mapped to are only those associated directly with the GO term (vs. indirectly via its immediate children or more distant offspring). For a list of all probes associated directly or indirectly with a particular GO term, one may use the hgu95av2GO2ALLPROBES environment. Example 6. Affymetrix probes directly associated with GO term ID GO:0004713. In the following example, 205 distinct Affymetrix probe IDs are associated directly with the GO term protein-tyrosine kinase activity (GO:0004713). The Affymetrix probe IDs include 1635 at, corresponding to the ABL1 gene.

GO terms-to-all annotated probes mappings:
The hgu95av2GO2ALLPROBES environment The hgu95av2GO2ALLPROBES environment provides key-value pairs for the mappings from GO term IDs (keys) to Affymetrix probe IDs (values). Values are vectors of length one or greater depending on whether a given GO term ID is mapped to one or more Affymetrix probe IDs. The value names are evidence codes for the GO term IDs. Note that, in accordance with the true path rule, the probes a particular GO term is mapped to are associated either directly with the GO term or indirectly via any of its immediate children or more distant offspring. The main difference between the hgu95av2GO2PROBE and hgu95av2GO2ALLPROBES environments is that the former considers only the GO term itself, whereas the latter considers the GO term and any of its offspring. Thus, the Affymetrix probe IDs returned by hgu95av2GO2PROBE are a subset of the probe IDs returned by hgu95av2GO2ALLPROBES. Example 7. Affymetrix probes directly or indirectly associated with GO term ID GO:0004713. In the following example, 319 distinct Affymetrix probe IDs (some with multiple evidence codes) are associated either directly or indirectly with the GO term protein-tyrosine kinase activity (GO:0004713). This list of 319 Affymetrix probe IDs indeed includes the list of 205 probe IDs associated directly with the GO term ID GO:0004713.

Assembling a GO gene-annotation matrix
This section provides R code for assembling an Entrez Gene ID-by-GO term ID gene-annotation matrix A, column by column. Specifically, rows correspond to (unique) Entrez Gene IDs mapping to probes on the HG-U95Av2 chip and columns to terms in the Molecular Function ontology mapping directly or indirectly to at least 10 Entrez Gene IDs for the HG-U95Av2 chip.
In practice, it may not be desirable to build the full G × M gene-annotation matrix, as this matrix could potentially be very large and sparse (padded with zeros). Rather, we assemble a (smaller) gene-annotation list, that provides, for each GO term ID, a list of Entrez Gene IDs annotated with the GO term. Example 8. Entrez Gene ID-by-GO term ID gene-annotation matrix for MF ontology. [1] 12625 > > ## Get all unique Entrez Gene IDs for HG-U95Av2 chip > LLID <-as.character(unique(unlist(mget(AffyID, + env=hgu95av2LOCUSID)))) > length(LLID) [1]

Acute lymphoblastic leukemia study of Chiaretti et al. [13]
Our proposed approach to tests of association with biological annotation metadata is illustrated using the acute lymphoblastic leukemia (ALL) microarray dataset of [13] and Gene Ontology (GO) annotation metadata.

Bioconductor experimental data R package ALL
The ALL dataset is available in the Bioconductor experimental data R package The expression measures have been obtained using the three-step robust multichip average (RMA) pre-processing method, implemented in the Bioconductor R package affy [11], and have been subjected to a base 2 logarithmic transformation.
For greater detail on the ALL dataset, please consult the ALL package documentation.

The BCR/ABL fusion
A number of recent articles have investigated the prognostic relevance of the BCR/ ABL fusion in adult ALL of the B-cell lineage [21].
The BCR/ABL fusion is the molecular analogue of the Philadelphia chromosome, one of the most frequent cytogenetic abnormalities in human leukemias. This t(9;22) translocation leads to a head-to-tail fusion of the v-abl Abelson murine leukemia viral oncogene homolog 1 (ABL1) from chromosome 9 with the 5' half of the breakpoint cluster region (BCR) on chromosome 22 (Figure 3). The ABL1 proto-oncogene encodes a cytoplasmic and nuclear protein tyrosine kinase that has been implicated in processes of cell differentiation, cell division, cell adhesion, and stress response. Although the BCR/ABL fusion protein, encoded by sequences from both the ABL1 and BCR genes, has been extensively studied, the function of the normal product of the BCR gene is not clear. The BCR/ABL proto-oncogene has been found to be highly expressed in chronic myeloid leukemia (CML) and acute myeloid leukemia (AML) cells [30].
An interesting question is therefore the identification of genes that are differentially expressed between B-cell ALL with the BCR/ABL fusion and cytogenetically normal NEG B-cell ALL.
In order to address this question, we consider gene expression measures for the n = 79 B-cell ALL cell samples (ALL$BT equal to B, B1, B2, B3, or B4), of the BCR/ABL or NEG molecular types (ALL$mol.biol equal to BCR/ABL or NEG).

Gene filtering
Many of the genes represented by the 12,625 probes are not expressed in B-cell lymphocytes. Accordingly, as in [44], we only retain the 2,391 probes that meet the following two criteria: (i) fluorescence intensities greater than 100 (absolute scale) for at least 25% of the 79 cell samples; (ii) interquartile range (IQR) of the fluorescence intensities for the 79 cell samples greater than 0.5 (log base 2 scale).
Furthermore, different probes may correspond to the same gene, i.e., map to the same Entrez Gene ID, according to the environment hgu95av2LOCUSID from the hgu95av2 package. In order to obtain one expression measure per gene, we choose to average the expression measures of multiple probes mapping to the same gene.
These various pre-processing steps lead to G = 2, 071 genes with unique Entrez Gene IDs.

Reduced ALL dataset: Genotypes and phenotypes of interest
The combined genotypic and phenotypic data for the n = 79 B-cell ALL cell samples of the BCR/ABL and NEG molecular types may be summarized by the set X Y n ≡ {(X i , Y i ) : i = 1, . . . , n}, of n pairs of G-dimensional gene expression profiles X i = (X i (g) : g = 1, . . . , G), G = 2, 071, and cancer class labels

Multiple hypothesis testing framework
Our primary question of interest is the identification of genes that are differentially expressed (DE) between BCR/ABL and NEG B-cell ALL. A subsequent question involves relating differential gene expression to GO annotation.
As detailed below, GO annotation metadata for the filtered list of G = 2, 071 unique genes from the HG-U95Av2 chip may be summarized by binary geneannotation profiles.
The gene-parameter profiles of interest concern differential gene expression between BCR/ABL and NEG B-cell ALL, i.e., the association between microarray gene expression measures and cancer class. Continuous gene-parameter profiles of unstandardized and standardized measures of differential expression are estimated, respectively, by (unstandardized) differences of empirical means and (standardized) two-sample t-statistics. Binary gene-parameter profiles, indicating whether genes are differentially expressed, are estimated by imposing cut-off rules on two-sample t-statistics or adjusted p-values.
The following association measures between GO gene-annotation profiles and DE gene-parameter profiles are considered: two-sample t-statistics for tests of association between binary GO gene-annotation profiles and continuous DE geneparameter profiles; χ 2 -statistics for tests of association between binary GO geneannotation profiles and binary DE gene-parameter profiles.
Significant associations between differential gene expression and GO annotation are identified by applying FWER-controlling single-step maxT Procedure 1, based on the non-parametric bootstrap null shift and scale-transformed test statistics null distribution of Procedure 2.

Gene-annotation profiles
Gene Ontology annotation metadata for the HG-U95Av2 chip series are obtained as described in Sections 4.2-4.5, from the following Bioconductor R packages: the GOspecific metadata package GO (Version 1.10.0, Bioconductor Release 1.7) and the Affymetrix chip-specific metadata package hgu95av2 (Version 1.10.0, Bioconductor Release 1.7).
For each of the three gene ontologies, binary gene-annotation matrices A BP , A CC , and A MF , are assembled for the GO terms annotating at least 10 of the G = 2, 071 filtered genes (sample R code provided in Section 4.

Gene-parameter profiles
Definition of gene-parameter profiles Consider a data structure (X, Y ) ∼ P , where X = (X(g) : g = 1, . . . , G) is a G = 2, 071-dimensional vector of microarray gene expression measures and Y ∈ {N EG, BCR/ABL} is a binary cancer class label. Let η k ≡ Pr(Y = k) denote the proportion of cancers of class k ∈ {N EG, BCR/ABL}. Define conditional G-dimensional mean vectors and G×G covariance matrices for the expression measures of class k ∈ {N EG, BCR/ABL} cancers by (43) µ

respectively.
Gene-parameter profiles, concerning differential gene expression between BCR/ ABL and NEG B-cell ALL, may be specified in various ways. Continuous DE geneparameter profiles may be defined in terms of the following unstandardized and standardized measures of differential gene expression between BCR/ABL and NEG B-cell ALL, and Absolute values of λ d (g) and λ t (g) may be used for measuring two-sided differential expression, i.e., either over-or under-expression in BCR/ABL compared to NEG B-cell ALL.
Binary DE gene-parameter profiles may be defined in terms of indicators for two-sided and one-sided differential expression.
Estimation of gene-parameter profiles The above DE gene-parameter profiles may be estimated as follows, based on the sample X Y n of gene expression measures for the n = 79 B-cell ALL cell samples of the BCR/ABL and NEG molecular types.
The resubstitution estimators of the continuous gene-parameter profiles of Equation (44) are based, respectively, on differences of empirical means and two-sample Welch t-statistics (up to the multiplier 1/ √ n). That is, and λ t n (g) ≡ 1 √ n µ BCR/ABL,n (g) − µ N EG,n (g) σ BCR/ABL,n (g,g) n BCR/ABL + σNEG,n(g,g) nNEG , where µ k,n (g) ≡ i I (Y i = k) X i (g)/n k and σ k,n (g, g) ≡ i I (Y i = k) (X i (g) − µ k,n (g)) 2 /(n k − 1) denote, respectively, the empirical means and variances of the gene expression measures for cancers of class k ∈ {N EG, BCR/ABL}.
Estimating the two-sided binary gene-parameter profile λ = of Equation (45) involves the two-sided test of the G null hypotheses H 0 (g) = I (µ BCR/ABL (g) = µ N EG (g)), of no differences in mean gene expression measures between BCR/ABL and NEG B-cell ALL. Likewise, estimating the one-sided binary gene-parameter profiles λ + and λ − involves, respectively, the one-sided test of the G null hypotheses of no over-expression (H 0 (g) = I (µ BCR/ABL (g) ≤ µ N EG (g))) and no underexpression (H 0 (g) = I (µ BCR/ABL (g) ≥ µ N EG (g))) in BCR/ABL compared to NEG B-cell ALL. For single-step common-cut-off maxT Procedure 1, adjusted p-values produce the same gene ranking as the test statistics defined in Equation (46). Simple and naive estimators of the three sets of differentially expressed genes (i.e., false null hypotheses), represented by the gene-parameter profiles λ = , λ + , and λ − , are therefore given, respectively, by the sets of genes with the largest γG values of |λ t n (g)|, λ t n (g), and −λ t n (g). That is, Analogous estimators may also be based on other test statistics, such as unstandardized difference statistics λ d n . More sophisticated estimators, that translate the proportion γ of rejected null hypotheses into a Type I error rate such as the gFWER, TPPFP, or FDR, could be based on adjusted p-values for the multiple test of the G null hypotheses H 0 (g). For example, one could estimate λ = by where P = 0n (g) are adjusted p-values for a suitably chosen multiple testing procedure, such as, FWER-controlling single-step maxT Procedure 1 or a TPPFP-controlling augmentation multiple testing procedure (Chapter 6 in [14], [40]). One-sided binary gene-parameter profiles λ + and λ − could be estimated likewise.

Association measures for gene-annotation and gene-parameter profiles
The association between continuous DE gene-parameter profiles, as in Equation (44), and binary GO gene-annotation profiles may be measured by two-sample Welch t-statistics (or corresponding p-values). Specifically, given a continuous Gvector x and a binary G-vector y, define the following association measure, where y k ≡ g I (y(g) = k),x k ≡ g I (y(g) = k) x(g)/y k , and v[x] k ≡ g I (y(g) = k) (x(g) −x k ) 2 /(y k − 1), k ∈ {0, 1}.
The association between binary DE gene-parameter profiles, as in Equation (45), and binary GO gene-annotation profiles may be measured by χ 2 -statistics (or corresponding p-values) for the test of independence of rows and columns in a 2 × 2  contingency table, such as Table 2. Specifically, given binary G-vectors x and y, define the following association measure, (50) ρ χ (x, y) ≡ G(g 00 g 11 − g 01 g 10 ) 2 (g 00 + g 01 )(g 00 + g 10 )(g 11 + g 01 )(g 11 + g 10 ) , where g kk ′ ≡ g I (x(g) = k) I (y(g) = k ′ ), k, k ′ ∈ {0, 1}. Given an association measure 3 ρ : R G×M × R G → R M , a G × M GO geneannotation matrix A, and a G-dimensional DE gene-parameter profile λ = Λ(P ), the M -dimensional association parameter vector ψ = Ψ(P ) of primary interest is defined as The corresponding resubstitution estimator ψ n =Ψ(P n ) is simply obtained by replacing the gene-parameter profile λ by an estimator thereof λ n =Λ(P n ), that is, (52) ψ n ≡ ρ(A, λ n ).

Test statistics
One-sided and two-sided tests of null hypotheses concerning any of the association parameters defined above may be based on (unstandardized) difference statistics T n (m), defined as in Equation (40).
For one-sided tests, large values of the test statistics T n (m) provide evidence against the corresponding null hypotheses H 0 (m), that is, rejection regions are of the form C n (m) = (c n (m), +∞). For two-sided tests, large values of the absolute test statistics |T n (m)| provide evidence against the corresponding null hypotheses H 0 (m).

Multiple testing procedures
For the purpose of illustration, we focus on control of the family-wise error rate, using single-step maxT Procedure 2.9, based on the non-parametric bootstrap null shift-transformed test statistics null distribution of Procedure 2.9 (null shift values λ 0 (m) = 0 and no scaling).
Let O n (m) denote indices for the ordered adjusted p-values, so that P 0n (O n (1)) ≤ · · · ≤ P 0n (O n (M )). GO terms with adjusted p-values less than or equal to α are declared significantly associated with differential gene expression at nominal FWER level α. That is, the list of GO terms found to be associated with differential gene expression is where R n (α) ≡ |R n (α)| denotes the number of identified GO terms.

Summary of testing scenarios
This section summarizes our approach for identifying GO terms associated with BCR/ABL vs. NEG differential gene expression. For each of the three gene ontologies (i.e., BP, CC, and MF), we consider the following three types of testing scenarios, each corresponding to a different association parameter ψ = ρ(A, λ) for GO annotation and BCR/ABL vs. NEG differential gene expression. Scenarios MT[t, t] and MT[d, t] are very similar and correspond, respectively, to continuous geneparameter profiles of standardized and unstandardized measures of differential gene expression. In contrast, Scenario MT[ =, χ] corresponds to a binary gene-parameter profile of differential gene expression indicators.
Scenario MT[t, t]: Association parameter ψ t,t = ρ t (A, |λ t |), for standardized continuous DE gene-parameter profile λ t . Consider the two-sided test of vs.
where the association parameter vector of interest is defined as based on Equations (44) and (49), and the null values are ψ t,t 0 (m) ≡ 0. The continuous DE gene-parameter profile λ t is estimated by λ t n , as in Equation (46), and the association parameter ψ t,t is estimated by the resubstitution estimator ψ t,t n ≡ ρ t (A, |λ t n |), as in Equation (52). The test statistics are defined as (unstandardized) difference statistics, T t,t n (m) ≡ √ n(ψ t,t n (m) − ψ t,t 0 (m)), and the null hypotheses H t,t 0 (m) are rejected for large absolute values of T t,t n (m). vs.
where the association parameter vector of interest is defined as based on Equations (44) and (49), and the null values are ψ d,t 0 (m) ≡ 0. The continuous DE gene-parameter profile λ d is estimated by λ d n , as in Equation (46), and the association parameter ψ d,t is estimated by the resubstitution estimator ψ d,t n ≡ ρ t (A, |λ d n |), as in Equation (52). The test statistics are defined as (unstandardized) difference statistics, In what follows, the G-dimensional gene-parameter profiles λ correspond to the G = 2, 071 genes with unique Entrez Gene IDs, obtained as described in Section 5.1. For each of the three gene ontologies, binary gene-annotation matrices are assembled for the GO terms annotating at least 10 of the G = 2, 071 genes of interest: G = 2, 071×M BP = 367 gene-annotation matrix A BP for the BP ontology, G = 2, 071 × M CC = 81 gene-annotation matrix A CC for the CC ontology, and G = 2, 071 × M MF = 185 gene-annotation matrix A MF for the MF ontology.

Differentially expressed genes between BCR/ABL and NEG B-cell ALL
In order to identify differentially expressed genes between BCR/ABL and NEG B-cell ALL, two-sided tests of the G null hypotheses H 0 (g) = I (µ BCR/ABL (g) = µ N EG (g)) are performed using the two-sample t-statistics λ t n (g) of Equation (46) and FWERcontrolling bootstrap-based single-step maxT Procedure 2. Adjusted p-values P = 0n (g) are obtained using the MTP function from the multtest package (Version 1.8.0, Bioconductor Release 1.7), with B = 5, 000 non-parametric bootstrap samples and other arguments set to their default values. Figure 4 displays a normal quantile-quantile plot of the test statistics λ t n (g) (Panel (a)) and a plot of the sorted bootstrap-based single-step maxT adjusted p-values P = 0n (g) (Panel (b)). A handful of genes stand out in terms of their large absolute test statistics and small adjusted p-values.
For control of the FWER at nominal level α = 0.05, Procedure 2.9 identifies 16 differentially expressed genes, i.e., 16 genes with P = 0n (g) ≤ α. Table 3 provides   (46)), and adjusted p-values P = 0n (g), for the 16 genes found to be significantly differentially expressed between BCR/ABL and NEG B-cell ALL, at nominal FWER level α = 0.05, according to the bootstrap-based single-step maxT procedure, with two-sample t-statistics λ t n (g) and B = 5, 000 bootstrap samples. A more detailed hyperlinked table, including information on gene function, chromosomal location, links to GenBank, Entrez Gene, NCBI Map Viewer, UniGene, PubMed, AmiGO, and KEGG, is provided on the website companion (Supplementary Table 10 Unsurprisingly, the gene showing the most over-expression in BCR/ABL cell samples, as measured by the t-statistics λ t n , is the ABL1 gene (v-abl Abelson murine leukemia viral oncogene homolog 1), located on the long arm of chro-mosome 9 (9q34.1). As mentioned in Section 5.1, the BCR/ABL phenotype is indeed defined in terms of the ABL1 gene.
Furthermore, many of the DE genes seem to be related to apoptosis or oncogenesis. For example, the Kruppel-like factor 9 (KLF9) gene encodes a transcription factor that binds GC-box elements in gene promoter regions. The Krüppel-like factor (KLF) family is comprised of highly related zinc-finger proteins, that are important components of the eukaryotic cellular transcriptional machinery and that take part in a wide range of cellular functions (e.g., cell proliferation, apoptosis, differentiation, and neoplastic transformation). In particular, KLFs have been linked to various cancers [25]. The intron-less gene AHNAK nucleoprotein (desmoyokin) (AHNAK), located on the long arm of chromosome 11 (11q12.2), encodes an unusually large protein (≅ 700 kiloDalton (kDa)) that is typically repressed in cell lines derived from human neuroblastomas and several other types of tumors [35]. Yet another example, the caspase 8, apoptosis-related cysteine peptidase (CASP8) gene, encodes a key enzyme at the top of the apoptotic cascade and has been linked to neuroblastoma [3]. Likewise, other genes listed in Table 3, including MX1, FYN, ACTN1, FHL1, and TRAM2, appear to be related to the molecular biology of cancer.
Our results are in general agreement with those of [44], slight differences being due, most likely, to our preliminary gene filtering, which involves averaging the expression measures of multiple probes mapping to the same Entrez Gene ID.
For greater detail, the interested reader is invited to consult Supplementary Table  10.1 on the website companion and follow links to PubMed and other databases. Further exploration of the DE genes may be performed in R using the Bioconductor packages annotate and annaffy.

GO terms associated with differential gene expression between BCR/ABL
and NEG B-cell ALL Figure 5 displays, for each of the three gene ontologies and each of the three testing scenarios, plots of the sorted adjusted p-values, P 0n (O n (1)) ≤ · · · ≤ P 0n (O n (M )), for FWER-controlling bootstrap-based single-step maxT Procedure 2 (B = 5, 000 bootstrap samples). The smaller the adjusted p-values, the less conservative the procedure, and the longer the list R n (α) = {m : P 0n (m) ≤ α} of identified GO terms at any given nominal Type I error level α. Table 4 summarizes the results in terms of the numbers R n (α) = |R n (α)| of GO terms found to be significantly associated with BCR/ABL vs. NEG differential gene expression for different nominal FWER levels α.
In general, adjusted p-values tend to be quite large, with only a handful of GO terms identified as being significantly associated with BCR/ABL vs. NEG differential gene expression for nominal FWER level α ∈ {0.05, 0. 10 Table 4. GO terms associated with differential gene expression between BCR/ABL and NEG B-cell ALL. This table reports, for each of the three gene ontologies and each of the three testing scenarios, the numbers Rn(α) = |Rn(α)| of GO terms found to be significantly associated with BCR/ABL vs. NEG differential gene expression for different nominal FWER levels α. and MT[ =, χ : γG = 20] are due to the fact that while the estimators λ = n,0.05 and λ = n,20 coincide on the original sample, they may differ on bootstrap samples of these data.
Next, the three testing scenarios are compared in terms of the contents of the lists R n (α) of identified GO terms. Specifically, let O n (r) ≡ {O n (1), . . . , O n (r)} denote the set of indices corresponding to the r smallest adjusted p-values for a given gene ontology and testing scenario. Measures of agreement between testing scenarios are provided by the numbers of common GO terms among sets of ordered GO terms O n (r) of various cardinality r, i.e., by the cardinality of intersections of sets O n (r) for different testing scenarios. Figure 6 displays plots of numbers of common GO terms for pairs of testing scenarios. As expected, there is substantial overlap between the GO terms identified by Scenarios MT[t, t] and MT[d, t] for continuous DE gene-parameter profiles (blue plotting symbols in top panels). This suggests that, for the ALL dataset, standardized (λ t ) and unstandardized (λ d ) continuous measures of differential gene expression have similar properties. In contrast, there is much less overlap between the GO terms identified by Scenario MT[ =, χ], for binary Fig 7. GO terms associated with differential gene expression between BCR/ABL and NEG B-cell ALL, conditional distribution of λ t n given A. Conditional boxplots of the estimated continuous DE gene-parameter profile λ t n given the gene-annotation profiles A(·, m) for the top two GO terms m ∈ {On(1), On(2)} identified according to each of the three testing scenarios. Rows correspond to gene ontologies and columns to testing scenarios. In each panel, the white and gray boxplots correspond, respectively, to the GO terms with the smallest and second smallest adjusted p-values; boxplots for unannotated and annotated estimated gene-parameter profiles, (λ t n (g) : A(g, m) = 0) and (λ t n (g) : A(g, m) = 1), are labeled as 0 and 1, respectively. Non-overlapping notches (informally) represent large differences in medians.  Figure 7 displays conditional boxplots of λ t n given A(·, m), that is, boxplots of the unannotated and annotated estimated gene-parameter profiles, (λ t n (g) : A(g, m) = 0) and (λ t n (g) : A(g, m) = 1), respectively. Although the boxplots reveal clear differences (non-overlapping notches) between unannotated and annotated profiles for some of the GO terms (e.g., MF term GO:0003735), the differences can be subtle for other terms (e.g., MF term GO:0003924 0n (m) for t-tests of association between estimated continuous gene-parameter profiles λ t n and geneannotation profiles A(·, m). In contrast, and also as expected, Scenario MT[ =, χ] tends to identify GO terms with small p-values P =,χ 0n (m) for χ 2 -tests of association between estimated binary gene-parameter profiles λ = n,α and gene-annotation profiles A(·, m). Furthermore, the tables corroborate our earlier observation that Table 5. GO terms associated with differential gene expression between BCR/ABL and NEG Bcell ALL, top two BP GO terms. This table provides association measures between the estimated DE gene-parameter profiles λ t n and λ = n,α and the gene-annotation profiles A(·, m) for the top two BP GO terms m ∈ {On(1), On(2)} identified according to each of the three testing scenarios. A 1 (m) = g A(g, m): Number of genes directly or indirectly annotated with GO term m (out of G = 2, 071 genes, GOALLLOCUSID environment in GO package). P t,t 0n (m): Nominal unadjusted p-value for the two-sample t-test comparing the unannotated and annotated estimated continuous DE gene-parameter profiles, (λ t n (g) : A(g, m) = 0) and (λ t n (g) : A(g, m) = 1), respectively (t.test function from the R package stats, with default argument values). P =,χ 0n (m): Unadjusted p-value for the χ 2 -test of independence between the estimated binary DE gene-parameter profile λ = n,α , α = 0.05, and the gene-annotation profile A(·, m) (chisq.test function from the R package stats, with arguments simulate.p.value = TRUE, correct=FALSE). P 0n (m): Bootstrap-based single-step maxT adjusted p-value, according to which the top two GO terms are identified for each testing scenario.

Scenario
GO term    The plots indicate that more terms tend to be identified in the BP ontology compared to either the CC or MF ontologies, and fewer terms tend to be identified in the MF ontology compared to either the BP or CC ontologies. Note that comparisons based on adjusted p-values  Tables 8, 9, and 10 list the 20 GO terms with the smallest adjusted p-values for Scenario MT[t, t], applied to the BP, CC, and MF gene ontologies, respectively. Figures 9, 10, and 11 display portions of the directed acyclic graphs for the top 20 GO terms in each ontology. The figures suggest that GO terms associated with BCR/ABL vs. NEG differential gene expression tend to concentrate in certain branches of the DAGs, i.e., differential expression is associated with related properties of gene products. Although it is known that many of the effects of the BCR/ABL fusion are mediated by tyrosine kinase activity, the MF GO term protein-tyrosine kinase activity (GO:0004713) does not appear to be significantly associated with differential gene expression between BCR/ABL and NEG B-cell ALL (adjusted p-value of 0.8890 for Scenario MT[t, t]).
For illustration purposes, we further investigate two of the GO terms from Tables 8 and 10: GO term anti-apoptosis (GO:0006916), with ninth smallest adjusted p-value for Scenario MT[t, t] applied to the BP gene ontology, and GO term structural constituent of ribosome (GO:0003735), with the smallest adjusted p-value for Scenario MT[t, t] applied to the MF gene ontology. Tables 11 and 12 list genes directly or indirectly annotated with GO terms GO:0006916 and GO:0003735, respectively. Figure 12 displays mean-difference plots of the average expression measures in BCR/ABL and NEG cell samples for genes annotated with GO terms GO:0006916 and GO:0003735. Panel (a) in Figure 12 indicates that genes annotated with BP GO term antiapoptosis (GO:0006916) tend to be over-expressed in BCR/ABL compared to NEG cell samples. Among these 21 genes, only SOCS2 is significantly differentially expressed between BCR/ABL and NEG B-cell ALL (nominal FWER level α = 0.05, Table 3). However, a brief survey of the literature reveals that several of the genes in Table 11 interact with the BCR/ABL proto-oncogene. For instance, [27] investigate mechanisms for the BCR/ABL-mediated activation of the transcription factor NF-κB/Rel encoded by the NFKB1 gene. Their findings suggest that NF-κB/Rel may be a potential target for molecular therapies of leukemia. [30] demonstrate that ectopic expression of BCR/ABL interferes with the tumor necrosis factor (TNF) signaling pathway through the down-regulation of TNF receptors. The TNF gene encodes a multifunctional proinflammatory cytokine involved in the regulation of a wide spectrum of biological processes, including cell proliferation, differentiation, apoptosis, lipid metabolism, and coagulation. The TNF gene has been implicated in a variety of diseases, including autoimmune diseases, insulin resistance, and cancer.
As seen in Table 12, 22 of the 24 genes annotated with MF GO term structural constituent of ribosome (GO:0003735) code for ribosomal proteins. Although none of the 24 annotated genes is identified as being significantly differentially expressed between BCR/ABL and NEG B-cell ALL (nominal FWER level α = 0.05, Table  3), Panel (b) in Figure 12 suggests that these genes tend to be under-expressed in BCR/ABL cell samples. GO terms associated with differential gene expression between BCR/ABL and NEG B-cell ALL, DAG for top 20 BP GO terms. Portion of the directed acyclic graph for the 20 GO terms with the smallest adjusted p-values for Scenario MT[t, t] applied to the BP gene ontology. Nodes for the top 20 terms are shaded in lavender; black and red edges indicate, respectively, "is a" and "part of a" relationships among terms. The figure was produced using the QuickGO browser. According to QuickGO, the GO term IDs GO:019882 and GO:0030333 listed in Table 8 correspond to the same term, antigen processing and presentation. (Higher-resolution color version on website companion.)

Discussion
We have proposed a general and formal statistical framework for multiple tests of association with biological annotation metadata. A key component of our approach is the systematic and precise translation of a generic biological question into a corresponding multiple hypothesis testing problem, concerning association measures between known gene-annotation profiles and unknown gene-parameter profiles. This general and rigorous formulation of the statistical inference question allows us to apply the multiple hypothesis testing methodology developed in [14] and related articles, to control a broad class of Type I error rates, in testing problems involving general data generating distributions (with arbitrary dependence structures among variables), null hypotheses, and test statistics.
The flexibility of our approach was illustrated using the ALL microarray dataset of [13], with the aim of relating GO annotation to differential gene expression between BCR/ABL and NEG B-cell ALL. This analysis demonstrates the importance of selecting a suitable DE gene-parameter profile λ and measure ρ for the association between this gene-parameter profile and GO gene-annotation profiles A. Indeed, for the ALL dataset, the choice of gene-parameter profile for measur- Nodes for the top 20 terms are shaded in lavender; black and red edges indicate, respectively, "is a" and "part of a" relationships among terms. The figure was produced using the QuickGO browser. (Higher-resolution color version on website companion.) ing differential expression between BCR/ABL and NEG B-cell ALL has a large impact on the list of identified GO terms. Testing scenarios based on binary DE GO terms associated with differential gene expression between BCR/ABL and NEG B-cell ALL, BP GO term GO:0006916 and MF GO term GO:0003735. This figure displays meandifference plots of average expression measures in BCR/ABL and NEG cell samples, i.e., plots of µBCR/ABL,n(g) − µNEG,n(g) vs. (µBCR/ABL,n(g) + µNEG,n(g))/2, for genes directly or indirectly annotated with GO terms GO:0006916 (Panel (a)) and GO:0003735 (Panel (b)). The term antiapoptosis (GO:0006916) has the ninth smallest adjusted p-value for Scenario MT[t, t] applied to the BP gene ontology (Tables 8 and 11) and the term structural constituent of ribosome (GO:0003735) has the smallest adjusted p-value for Scenario MT[t, t] applied to the MF gene ontology (Tables  10 and 12).  Our analysis of the ALL microarray dataset clearly shows the limitations of binary gene-parameter profiles of differential expression indicators, which are still the norm for combined GO annotation and microarray data analyses. Our proposed statistical framework, with general definitions for the gene-annotation and gene-parameter profiles, allows consideration of a much broader class of inference problems, that extend beyond GO annotation and microarray data analysis. Geneannotation profiles may be continuous or polychotomous and may correspond, for example, to exon/intron counts/lengths/nucleotide distributions, gene pathway membership, or gene regulation by particular transcription factors. Likewise, geneparameter profiles may be continuous or polychotomous and may correspond, for example, to regression coefficients relating possibly censored biological and clinical outcomes to genome-wide transcript levels, DNA copy numbers, and other covariates.
This first application of our proposed methodology only considered control of the family-wise error rate using single-step common-cut-off maxT Procedure 1, based on the non-parametric bootstrap null shift and scale-transformed test statistics null distribution of Procedure 2. Adjusted p-values tended to be quite large, with only a handful of GO terms identified as being significantly associated with smaller than the dimension J of the data structure X, sample sizes have considerably increased in recent genomic applications. In addition, simulation studies have indicated that our proposed MTPs have good finite sample properties in terms of both Type I error control and power.
Ongoing efforts include consideration of more general and biologically pertinent multivariate association measures ρ. For instance, for GO annotation metadata, the association parameter for a given GO term could take into account the structure of the DAG by considering the gene-annotation profiles of offspring or ancestor terms. We are also interested in developing better numerical and graphical approaches for representing and interpreting the multiple testing results, e.g., the lists of GO terms and associated adjusted p-values. Finally, we are planning on implementing the proposed methods in an R package to be released as part of the Bioconductor Project.
The experimental data (ALL) and annotation metadata (annaffy, annotate, GO, hgu95av2) packages used in the analysis of Section 5 may also be obtained from the Bioconductor Project website.