BACKGROUND
increasingly it is a commonplace for researchers to conduct large scale genomic studies involving multiple experimental groups along with a control group, also called the normal or the reference group.
the goal is to determine features that are differentially expressed in a given experimental group  and to determine if a differentially expressed feature is up or down regulated.
for example, a toxicologist may be interested in identifying differences in the gene expression profile of spontaneous tumors and chemically induced tumors, relative to normal tissues .
there is considerable interest among cancer researchers to understand the gene expression profile of tumors according to tumor size .
accordingly, this has been an active area of research for various cancers over the past decade  <cit> .
for example, gieseg et al.
<cit>  investigated gene expression of  <dig> different sizes of colon cancer tumors and found gene expressions to be constant with tumor size.
however, diaz et al.
<cit>  demonstrated that the gene expression of β <dig> integrin is correlated with both size and nuclear grade of the breast cancer tumors, with increased expression in large tumors and in higher tumor grades.
more recently riis et al.
<cit>  conducted an extensive gene expression analysis of breast cancer tumors according to size categories and identified genes and activation pathways that were specific to tumor size.
they concluded that the molecular signatures derived from their analysis will help clinicians determine whether to treat tumors of a given size category aggressively or not.
earlier ciarmela et al.
<cit>  conducted an exhaustive study of the role of various growth factors on the development of uterine fibroids .
thus it is clear there is considerable interest among clinicians and biologists to investigate the expression of genes according to the tumor size or category.
in all such investigations, one is typically interested in performing several pairwise comparisons, of thousands of features, relative to a reference group .
often researchers are not only interested in determining if a feature is differentially expressed but are also interested in determining whether it is up or down-regulated in the experimental group .
for simplicity of exposition, throughout this paper we shall replace the term “feature” by “gene”.
multiple testing problems involving multiple pairwise comparisons of high dimensional data along with directional decisions has not received much attention in the literature, yet such testing problems are commonly encountered in practice.
when the number of genes is very small  several methods have been proposed that control the directional errors as well as the family wise error rate .
however, such methods are very conservative when the number of genes is very large as in a microarray data or cpg methylation data .
several ad-hoc methods and strategies are used in the literature when the number of genes is large.
for example, some researchers apply multiple testing procedures  or bonferroni procedure) within each pairwise comparison and ignore the fact that they are conducting several pairwise comparisons.
once a differentially expressed gene  for a pairwise comparison is identified then they are declared to be up or down-regulated by looking at the direction of the fold change  without accounting for the statistical error associated with such a directional decision.
such strategies result in an inflated overall error rate due to multiple pairwise comparisons and directional decisions.
this overall error rate is called the mixed directional false discover rate .
another procedure that is commonly used by researchers is to first perform anova across all experimental groups for each feature.
using the resulting p-values, a multiple testing procedure is applied to obtain an initial list of degs.
these degs are then subject to a tukey’s test or a bonferroni test to identify the pairwise comparisons where a deg is significant.
the directional decision regarding whether a gene is up or down regulated, is made on the basis of the direction of the fold change.
while this strategy is more intuitive than the earlier one, it is potentially a conservative procedure resulting in a substantial loss of power, firstly because intrinsically anova compares all treatment groups.
secondly, the thresholds or the critical values at the two steps are not derived with the goal to control the overall error rate.
the only formal methodology available in the literature that controls the mdfdr for the above directional multiple testing problems is the method by guo et al.
<cit> , which is designed to make decisions on thousands of features when making multiple pairwise comparisons and deciding on the direction of comparison.
thus the method controls the false discovery rate when making multiple pairwise comparisons on thousands of genes while also controlling the directional errors committed when falsely declaring a deg to be up-regulated  when it is not.
guo et al.
<cit>  procedure generalizes the procedure of benjamini and yekutieli  <cit>  which was designed for the case when there were only two groups to compare.
the procedure of guo et al.
<cit>  is available in the software oriogen  that can be accessed from http://www.niehs.nih.gov/research/atniehs/labs/bb/staff/peddada/.
while guo et al.
<cit>  methodology is useful for making several multiple pairwise comparisons; it is relatively conservative since it relies on the bonferroni procedure to deal with multiple pairwise comparisons within each gene.
in this paper we develop a general mdfdr controlling testing procedure that allows us to use any mixed directional familywise error rate  controlling procedure in place of the bonferroni procedure, for conducting pairwise comparisons in high dimensional data that is broadly applicable to a wide range of genomic data including gene expression microarray data, cpg methylation data, rna-seq data and others.
based on this general procedure, we develop a specific methodology using the dunnett’s test  which is designed for making comparisons of several experimental groups with the control or the reference group.
not only that the resulting methodology is practically relevant but as demonstrated in the numerical simulations, the resulting methodology not only controls the mdfdr but is more powerful relative to some potential alternative methods.
we illustrate the proposed methodology by analyzing a gene expression data obtained from the niehs’ fibroid growth study   <cit> .
using our methodology we gain deeper insights into molecular characteristics of uterine fibroids according to the tumor size, which was not understood until now.
we have identified several differentially expressed genes and pathways that are specifically enriched according to the tumor size.
while researchers and clinicians who study fibroids are well aware of many of the genes and pathways described in this paper, we have provided a characterization of these genes and pathways according to the tumor size.
our data can be further mined to gain deeper insights regarding fibroids.
methods
in this section we present the general testing procedure that controls mdfdr while making multidimensional directional decisions, and develop the statistical methodology based on dunnett’s test that is used to analyze fgs gene expression data.
the relevant notations, definitions, hypotheses and other statistical concepts are described in the online additional file  <dig>  and additional file 1: figure s <dig>
the general mdfdr controlling procedure
the general testing algorithm:use any global testing method to obtain the screening p-values.
apply bh procedure  <cit>  to identify genes that are differentially expressed in at least one pairwise comparison.
let r denote the number of genes so discovered.
for each gene discovered in step  <dig>  use any two-sample testing method to obtain the p-values for each pairwise comparison and apply any mixed directional family wise error  controlling procedure, such as holm, hochberg etc.
<cit> , to the pairwise p-values at level rα/m.
for a given gene discovered in step  <dig>  if a pairwise hypothesis is rejected in step  <dig>  then we declare the gene to be up or down regulated in the pairwise comparison according to the sign of the corresponding test statistic.
this procedure is general in the sense that any global testing method can be used to obtain the screening p-values in step  <dig> of the procedure, any pairwise comparison testing method can be used to obtain the pairwise p-values in step  <dig> and any mdfwer controlling procedure can be used in steps 2- <dig> of the procedure.
it is important to note here that the method of guo et al.
<cit>  is a special case of the proposed general procedure in which bonferroni global test is used for testing the screening hypotheses and bonferroni method along with additional directional decisions works as the mdfwer controlling procedure.
this general algorithm is proved to control mdfdr through the following theorem with the mathematical proof provided in the online additional file 6: section s <dig>
theorem  <dig>
under assumption of independence of p-value vectors pj, j= <dig> …,m, the mdfdr of the general testing algorithm is strongly controlled at level α.
statistical methodology for fgs gene expression data
suppose we are interested in comparing “q” experimental groups with a reference group  on the basis of the mean expression of “m” genes.
for example, suppose we are interested in comparing “small”, “medium” and “large” fibroids with a “normal” tissue  from uterus on the basis of “m” genes.
our goal is not only to identify differentially expressed genes in any given pairwise comparison but also to determine if the mean expression is up or down-regulated in the tumor tissue compared to the normal myometrium.
our statistical methodology for the fgs gene expression data analysis proceeds in three steps as follows.
specific details of implementation of each step are described in the online additional file 6: section s <dig>
for each gene we obtain a dunnett  based screening p-value from all “q” pairwise comparisons with the reference group.
this results in one p-value for each gene.
applying the bh procedure  <cit>  on these screening p-values, we obtain genes that are differentially expressed in at least one pairwise comparison.
suppose we discover r genes in this step.
thus there are r genes which are differentially expressed in at least one pairwise comparison with the reference group.
for the jth gene discovered in step  <dig>  we compute dunnett’s p-value pij for each pairwise comparison i, i= <dig> ,…,q, with the reference group.
if pij≤rαm then we declare that the ith pairwise comparison with the reference group is significant.
if a gene j is found to be significant in the ith pairwise comparison with the reference group then we declare it to be up-regulated in the ith group relative to the reference group if tij> <dig>  otherwise it is declared to be down-regulated.
here tij denotes the test statistic associated with the jth gene in the ith pairwise comparison.
remarks.
it is important to note in steps  <dig> and  <dig> that it is possible for a gene to be differentially expressed on the basis of the screening p-value, yet it may not be significant in any of the pairwise comparisons.
this happens when the bh-adjusted p-values are only marginally significant.
in such cases we declare that gene to be not differentially expressed.
we expect this phenomenon to occur with very small frequency.
RESULTS
in this section we present two sets of numerical results.
first, the operating characteristics of the proposed procedure based on a simulation study are presented.
next, we apply the proposed methodology to a gene expression data obtained from the niehs’ fgs to identify tumor size specific degs.
numerical simulation and discussion
in this section we present the operating characteristics of the proposed procedure based on a simulation study.
this study evaluates the power of the proposed dunnett based methodology and compares its performance with some competing procedures, the details of which are given in the online additional file 6: section s <dig>
in our simulation we considered p= <dig> groups, the first three were taken to be experimental groups and the last group was taken to be the reference group.
thus all pairwise comparisons are made with the last group.
our simulated microarray chip consisted of m= <dig> genes per chip with n= <dig> chips per group.
as often done in simulation studies for microarray data , we generated the expression of each gene in each chip using a normal distribution.
more precisely, for the jth gene, j= <dig> ,…, <dig>  in the sth sample, s= <dig> ,…, <dig> in the gth group, g= <dig> , <dig> , we generated its expression zgjs from a normal distribution with mean value ezgjs=μgj and variance vzgjs= <dig>  for the reference group  we set μ4j= <dig> for all j. to create the null data we set μgj= <dig> for j= <dig> ,…,m <dig>  g= <dig> , <dig> and non-null data were created by generating μgj∼independentu, for g= <dig> , <dig> and j=m0+ <dig> m0+ <dig> …,m0+m <dig>  where m1=m−m <dig> and u represents the uniform distribution.
it is important to note that, for the non-null means considered here, the standard deviation used in this simulation study is large.
consequently, all methods considered in this simulation study are expected to have small power.
we considered three patterns of correlation structure as follows:  independent gene expressions: the correlation coefficient ρ between any pair of genes and any pair of sample is  <dig>  thus, the data are completely independent.
gene expressions within sample are dependent: for a given sample s, the correlation coefficient between any pair of genes is ρ but the correlation coefficient between any pair of samples is  <dig>   gene expressions within genes are dependent: for a given gene j, the correlation coefficient between any pair of samples is ρ but the correlation coefficient between any pair of genes is  <dig>  although we considered a variety of patterns of ρ, due to space limitations we only present the results corresponding to independence.
similar results are obtained for other patterns of ρ .
we simulated the mdfdr and average power using  <dig> simulation runs.
our nominal mdfdr level was taken to be α= <dig> .
details regarding the calculations of the mdfdr and power are provided in the accompanying online additional file  <dig>  for comparison purposes, we compared the proposed procedure with the guo et al.
<cit>  procedure and variants of the proposed procedure in which the mdfwer controlling procedures used in steps  <dig> and  <dig> of the algorithm are respectively, the holm’s procedure  <cit> , the hochberg’s procedure  <cit>  and the bonferroni procedure.
these procedures are described in detail in the online additional file  <dig>
we summarize the results of the simulation studies in fig.
<dig> for the independence case.
the horizontal axis denotes the proportion of truly differentially expressed genes on the array and the vertical axis denotes the average mdfdr  and average power .
as desired, all five procedures control the mdfdr on average at α= <dig> .
however, the proposed dunnett based procedure has highest power compared to all other methods considered
fig.
<dig> mdfdr , average power .
mdfdr , average power  with the proposed methodology and three variants using holm, hochberg and bonferroni procedures, respectively, in steps  <dig> and  <dig> along with guo et al.
<cit>  procedure, under independence among genes and within genes
we summarize the results for the case of dependence within genes across groups in fig.
<dig>  additional file 2: figure s <dig> and additional file 3: figure s <dig> and the results for the case of dependence among genes in fig.
<dig>  additional file 4: figure s <dig> and additional file 5: figure s <dig>  the details of these dependence structures are given in the online additional file 6: section s <dig>  although all five procedures control the mdfdr on average at  <dig> , the power of our proposed methodology is much higher compared to all four competing procedures under the different dependence structures considered.
fig.
<dig> mdfdr , average power .
mdfdr , average power  with the proposed methodology and three variants using holm, hochberg and bonferroni procedures, respectively, in steps  <dig> and  <dig> along with guo et al.
<cit>  procedure, under dependence  within genes
fig.
<dig> mdfdr , average power .
mdfdr , average power  with the proposed methodology and three variants using holm, hochberg and bonferroni procedures, respectively, in steps  <dig> and  <dig> along with guo et al.
<cit>  procedure, under dependence  among genes
analysis of gene expressions in uterine fibroids
uterine fibroids are benign smooth muscle tumors occurring in the uterus of women with a very high prevalence rates in many populations.
these tumors are considered to be hormonally mediated and therefore occur during the pre-menopausal years.
according to some estimates, the cumulative incidence by age  <dig> of these tumors among caucasian women exceeds  <dig> % and it is much higher among women of african american descent.
the direct and indirect annual cost of fibroids in the us is as high as  <dig> billion dollars, yet these tumors are not well-studied perhaps because they are considered to be benign.
the niehs conducted a large prospective study tracking the growth of fibroids in  <dig> pre-menopausal women over  <dig> months and discovered difference in growth patterns of fibroids according to age for white and black women  <cit> .
a total of  <dig> tumor samples and  <dig> normal myometrium samples were collected from  <dig> women who opted for either hysterectomy or myomectomy.
the rna from these tissue samples were processed to obtain gene expression data using affymetrix chips .
a total of  <dig> probe sets were obtained.
using these gene expression data davis et al.
<cit>  studied the molecular characteristics of fibroids.
several important degs and canonical pathways were identified in davis et al.
<cit> , which were also experimentally validated using rtpcr.
in this paper we take the next step towards the understanding of molecular characteristics of fibroids by investigating the changes in gene expression according to tumor size.
as done in the literature for other tumors such as breast cancer tumors and others , we are interested in identifying genes that are specific to size of uterine fibroids.
we classified fibroids into three size groups according to their volume, small , medium  and large  and compared each size category with the normal myometrium, the reference group, using the proposed methodology.
the cut points between size categories were chosen so that there is a natural separation between the size categories.
the sorted volumes of the  <dig> tumors are provided in fig.
<dig>  due to large differences in the volumes between large and small tumors, for clarity we provided volumes in three panels according to their size.
fig.
<dig> volume of fibroids according to size category
we identified a total of  <dig> probe sets to be differentially expressed in at least one pairwise comparison  at mdfdr of  <dig> .
these  <dig> probe sets map to  <dig> genes.
the venn diagram of the degs by tumor size is in fig.
<dig>  based on the  <dig> genes, using the ingenuity pathway analysis , we discovered a total  <dig> distinct enriched canonical pathways.
the venn diagram of the number of enriched canonical pathways by tumor size is provided in fig.
<dig>
fig.
<dig> venn diagram of differentially expressed genes by tumor size
fig.
<dig> venn diagram of enriched pathways by tumor size
it is not surprising that a third of the degs and nearly  <dig> % of all enriched pathways are common to tumors of all sizes since tumor tissues are fundamentally different from the normal myometrium.
however, we discover several degs and pathways that are uniquely enhanced according to the tumor size, suggesting that changes in molecular characteristics might be taking place as tumors grow.
we highlight some important groups of genes and pathways discovered by our methodology and provide the complete list of all genes and pathways in the additional files  <dig> and  <dig>  respectively, for researchers to mine our data.
the “gene symbol” and “gene name” columns in these tables are obtained from the website david .
growth factors, vascularization and related pathways: the netrin signaling pathway is well-known as a versatile pathway with multiple functions.
one of its functions is to promote vascular networks and branching of blood vessels  <cit>  and angiogenesis  <cit> .
according to our analysis this pathway is uniquely enriched in small tumors only.
the interleukin- <dig>  pathway is known to induce inflammatory response and the production of prostaglandins and expression of il- <dig> which may play a critical role in the fibroid initiation and early development.
for example, prostaglandins play a critical role in the promotion of growth factors involved in angiogenesis, such as vegf, basic fibroblast growth factor  and platelet-derived growth factor   <cit>  and development of fibroid requires vascularization and blood supply.
thus interleukin- <dig>  pathway likely plays an important role during fibroid development.
interestingly, according to our analysis, this pathway is uniquely enriched in small size tumors but not in the medium or large tumors .
furthermore, the fibroblast growth factors  <dig> and  <dig>  which belong to the regulation of the epithelial-mesenchymal transition pathway and are well-known to be involved in vascularization and angiogenesis, are both uniquely down-regulated in small tumors and not differentially expressed in medium or large tumors.
our analysis further implies that the epithelial-mesenchymal transition pathway was enriched only in small and medium size tumors.
according to ciarmela et al.
<cit> , estrogen may promote fibroid growth through up-regulation of epidermal growth factor receptor .
however, we found egfr to be down-regulated in fibroids and that too only in the medium size tumors.
similarly, the fibroblast growth factors   were differentially expressed only in the medium size tumors.
the acidic fgf was up-regulated whereas the basic fgf was down regulated.
these findings are consistent with ciarmela et al.
<cit>   in that they are expressed during tumor progression.
similarly, insulin like growth factor  was only differentially expressed  in medium size tumors.
additionally, growth factor signaling pathways such as vegf, pdgf, tgf β and egf are uniquely enriched in medium size tumors and not in small or large tumors.
while the above results suggest that some growth factors and related pathways are very specific to small and medium tumors, we discovered several growth factors to be differentially expressed in tumors of all sizes.
these included, insulin like growth factor  <dig> , insulin like growth factor binding protein  <dig> , and platelet derived growth factor c  which were up-regulated whereas insulin-like growth factor binding protein  <dig> , connective tissue growth factor , heparin-binding egf-like growth factor , transforming growth factor beta receptor ii , fibroblast growth factors  <dig> and  <dig>  were down-regulated.
similarly, growth factor signaling pathways such as the human growth factor  and igf- <dig> were enriched in tumors of all sizes.
thus it appears that the differential expression of these genes and the enrichment of the above pathways is necessary for tumor onset and progression.
estrogen and related genes: fibroids are hormonally mediated and it is also well documented in the literature that accordingly estrogen and progesterone receptors and prostaglandins promote proliferation of fibroids .
not surprisingly, we found the estrogen receptor esr- <dig> to be up-regulated in tumors of all sizes.
this result was confirmed using the rtpcr data .
our rtpcr data further indicates that the expression of esr- <dig> in small tumors is almost  <dig> % higher than in large tumors, suggesting that esr- <dig> has a large effect on small tumors but is also needed for the continued growth of the tumor.
interestingly, the progesterone receptor  was up-regulated in only medium size tumors and not differentially expressed in small or large tumors, a result that is also confirmed by the rtpcr .
this suggests that perhaps pgr may not be involved in tumor initiation  but is only involved in growth of the tumor.
however, its function ends once the tumor becomes large enough.
most prostaglandins were generally down-regulated in tumors of all sizes.
for example, prostaglandin e receptor  <dig> , prostaglandin f receptor  and prostaglandin-endoperoxide synthase  <dig>   are down-regulated in tumor of all sizes.
the result of ptgs <dig> was also confirmed using rtpcr .
however, some prostaglandins were differentially expressed according to the size of the tumor.
for example, prostaglandin-endoperoxide synthase  <dig>   was differentially expressed only in the medium sized tumors where it was down-regulated and prostaglandin e synthase  <dig>  was down-regulated in large tumors only but prostaglandin e receptor  <dig>  was down-regulated in both medium and large tumors.
similar to estrogen and progesterone receptors, the α and γ isoforms of peroxisome proliferator-activated  receptors have been associated with the regulation of proliferation of uterine fibroids  <cit> .
in our data both these isoforms are down-regulated in all tumor sizes compared to the normal myometrium.
we also discovered the related retinoid x receptor gamma to be down-regulated in the medium size tumors but was not significant in small or large tumors.
collagens: there is a vast amount of literature implicating collagens to smooth muscle tumors such as the fibroids .
consequently, it is not surprising that several collagens  and extra cellular matrix proteins are differentially expressed in tumors of all sizes.
apart from col4a <dig>  col4a <dig> and col21a <dig>  which were down-regulated in all tumor size groups, the remaining  <dig> collagens were up-regulated in tumor samples.
we investigated a subset of these collagens using rtpcr which confirmed the above findings for those genes .
there were other collagens that were differentially expressed depending upon the size of the tumor.
for example, col2a <dig>  col4a <dig> and col5a <dig> were not differentially expressed in small tumors but were significantly up-regulated in medium and large tumors.
collagen col23a <dig> was significantly up-regulated in small and medium tumors but not in the large fibroids.
interestingly, some of the collagens were differentially expressed in only one of the size categories.
for example, col4a3bp, col6a <dig>  col9a <dig>  col10a <dig> and col11a <dig> were differentially expressed only in the large tumors.
again, except for col4a3bp, the remaining  <dig> collagens were up-regulated in the tumor samples.
interestingly, col5a <dig> was only differentially expressed  in small tumors.
other genes: leptin receptor is well-known to be negatively associated with the obesity and obesity is a potential risk factor for fibroids.
interestingly, we discover leptin receptor  to be significantly down-reglated in all tumor sizes.
this finding is also confirmed by the rtpcr data reported in additional file  <dig>  as noted earlier, diaz at al.
<cit>  demonstrated that β <dig> integrin had an increased expression in larger breast tumors and in higher tumor grades.
in our fibroid data, however, we notice a down-regulation in β <dig> integrin in medium and large tumors and was not differentially expressed in small tumors.
comparison of various procedures with dunnett’s test based methodology
we analyzed the niehs fgs data using various methods considered in this paper at mdfdr level of  <dig> .
we compared the methods in terms of number of probe sets that were differentially expressed  in;  all tumors,  small tumors,  medium tumors, and  large tumors.
results are summarized using venn diagrams.
in fig.
<dig> we compare our dunnett test based methodology with guo et al.
<cit>  procedure.
not only did the dunnett based methodology identified many of the probe sets identified by the guo et al.
<cit>  procedure, i.e.
shared several probe sets , it identified several more that were not identified by guo et al.
<cit>  procedure .
on the contrary, there were very few probe sets that were uniquely identified by guo et al.
<cit>  procedure .
we see similar results when we compare dunnett screening and holm procedure , dunnett screening and hochberg procedure  and dunnett screening and bonferroni procedure .
the venn diagram in fig.
<dig> compares the number of probe sets identified by dunnett based method to each of the other methods.
in each case the dunnett based method tends to uniquely identify many probe sets than each of the other methods.
the total number of probe sets identified by the various methods are summarized in table  <dig>
fig.
<dig> venn diagram of identified probe sets by the dunnett based methodology and guo et al.
<cit>  method.
the venn diagrams compare the two procedures in terms of the probe sets identified in at least one tumor size, large tumor size , medium tumor size  and small tumor size
fig.
<dig> venn diagram of identified probe sets by the dunnett based methodology and other methods discussed in numerical studies.
the other methods considered in the numerical simulation are dunnett screening and holm procedure , dunnett screening and hochberg procedure  and dunnett screening and bonferroni procedure
CONCLUSIONS
this paper offers a general statistical methodology for conducting pairwise comparisons in high dimensional data that is broadly applicable to a wide range of genomic data including gene expression microarray data, cpg methylation data, rna-seq data and others and proposes a specific testing methodology for the analysis of fgs gene expression data.
the proposed methodology not only controls the false discoveries made when making several pairwise comparisons on the basis of high dimensional data, but it controls for directional errors, i.e.
falsely declaring a gene to be up-regulated  when it is not.
thus it controls the mixed directional false discovery rate .
the framework we proposed here is general enough that, depending upon the application, a researcher can use different testing procedures within our framework.
in all situations the procedure strongly controls mdfdr.
for the application under consideration in this paper, the experimental design calls for the comparison of several experimental groups against a particular group.
for such designs, in the classical multiple comparison problem, the dunnett’ test is known to be powerful.
so, it is not surprising that our dunnett’s test based procedure is most powerful among all the procedures we considered in this paper.
our methodology based on dunnett procedure works best here in terms of power because  the dunnett’s p-values are based on the joint distribution of the test statistics and  our procedure uses the maximum order statistic to find the screening p-value.
using our methodology we gain deeper insights into molecular characteristics of uterine fibroids according to the tumor size.
we have identified several differentially expressed genes and pathways that are specifically enriched according to the tumor size.
while researchers and clinicians who study fibroids are well aware of many of the genes and pathways described in this paper, we have provided a characterization of these genes and pathways according to the tumor size.
our data can be further mined to gain deeper insights regarding fibroids.
in view of the recent publication  <cit>  we acknowledge that our results are potentially confounded by race.
since the sample size by race is very small, we are unable to separate out the race effect.
despite that potential weakness, we think that our findings provoke researchers to explore further research along these lines.
additional files
additional file  <dig>
figure s <dig>  a graphical display of various hypotheses of interest.
additional file  <dig>
figure s <dig>  mdfdr , average power  with the proposed methodology and three variants using holm, hochberg and bonferroni procedures, respectively, in steps  <dig> and  <dig> along with guo et al.
<cit>  procedure, under dependence  within genes.
additional file  <dig>
figure s <dig>  mdfdr , average power  with the proposed methodology and three variants using holm, hochberg and bonferroni procedures, respectively, in steps  <dig> and  <dig> along with guo et al.
<cit>  procedure, under dependence  within genes.
additional file  <dig>
figure s <dig>  mdfdr , average power  with the proposed methodology and three variants using holm, hochberg and bonferroni procedures, respectively, in steps  <dig> and  <dig> along with guo et al.
<cit>  procedure, under dependence  among genes.
additional file  <dig>
figure s <dig>  mdfdr , average power  with the proposed methodology and three variants using holm, hochberg and bonferroni procedures, respectively, in steps  <dig> and  <dig> along with guo et al.
<cit>  procedure, under dependence  among genes.
additional file  <dig>
supplementary materials.
section s <dig>  the mathematical notations and definitions used in the paper are described in this section.
the hypotheses tested in different steps of the testing procedure are formulated and interpreted.
section s <dig>  the mathematical proof of the theorem that states that the general testing algorithm controls mdfdr.
section s <dig>  the details of the statistical methodology developed for the testing and analysis of the fgs microarray data.
section s <dig>  the supplementary results for the simulation study performed to evaluate the power of the developed statistical methodology compared to several other competing methods.
additional file  <dig>
differentially expressed probe sets in different sizes of tumors compared to normal myometrium.
significant probe sets are grouped and color coded according to different sizes of fibroids.
the last column shows the bh-adjusted p-value for step  <dig> of the methodology.
the columns under "dunnett adjusted p-values" show the p-values for step  <dig> of the methodology.
additional file  <dig>
pathways identified in different sizes of tumors compared to normal myometrium.
in the table, the enriched pathways are grouped and color coded according to different sizes of fibroids.
the last column shows the negative of the logarithm of the bh adjusted p-value, computed by ipa.
additional file  <dig>
rtpcr data.
this table shows the rtpcr data for a subset of the genes that are identified in our analysis.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
ag, wg and sdp conceived the study.
ag performed all calculations and all numerical work.
ag, wg and sdp wrote the paper.
all authors read and approved the final manuscript.

 increasingly it is a commonplace for researchers to conduct large scale genomic studies involving multiple experimental groups along with a control group, also called the normal or the reference group. the goal is to determine features that are differentially expressed in a given experimental group and to determine if a differentially expressed feature is up or down regulated. for example, a toxicologist may be interested in identifying differences in the gene expression profile of spontaneous tumors and chemically induced tumors, relative to normal tissues. there is considerable interest among cancer researchers to understand the gene expression profile of tumors according to tumor size. accordingly, this has been an active area of research for various cancers over the past decade. for example, gieseg et al. investigated gene expression of different sizes of colon cancer tumors and found gene expressions to be constant with tumor size. however, diaz et al. demonstrated that the gene expression of β integrin is correlated with both size and nuclear grade of the breast cancer tumors, with increased expression in large tumors and in higher tumor grades. more recently riis et al. conducted an extensive gene expression analysis of breast cancer tumors according to size categories and identified genes and activation pathways that were specific to tumor size. they concluded that the molecular signatures derived from their analysis will help clinicians determine whether to treat tumors of a given size category aggressively or not. earlier ciarmela et al. conducted an exhaustive study of the role of various growth factors on the development of uterine fibroids. thus it is clear there is considerable interest among clinicians and biologists to investigate the expression of genes according to the tumor size or category. in all such investigations, one is typically interested in performing several pairwise comparisons, of thousands of features, relative to a reference group. often researchers are not only interested in determining if a feature is differentially expressed but are also interested in determining whether it is up or down-regulated in the experimental group. for simplicity of exposition, throughout this paper we shall replace the term “feature” by “gene”. multiple testing problems involving multiple pairwise comparisons of high dimensional data along with directional decisions has not received much attention in the literature, yet such testing problems are commonly encountered in practice. when the number of genes is very small several methods have been proposed that control the directional errors as well as the family wise error rate. however, such methods are very conservative when the number of genes is very large as in a microarray data or cpg methylation data. several ad-hoc methods and strategies are used in the literature when the number of genes is large. for example, some researchers apply multiple testing procedures or bonferroni procedure) within each pairwise comparison and ignore the fact that they are conducting several pairwise comparisons. once a differentially expressed gene for a pairwise comparison is identified then they are declared to be up or down-regulated by looking at the direction of the fold change without accounting for the statistical error associated with such a directional decision. such strategies result in an inflated overall error rate due to multiple pairwise comparisons and directional decisions. this overall error rate is called the mixed directional false discover rate. another procedure that is commonly used by researchers is to first perform anova across all experimental groups for each feature. using the resulting p-values, a multiple testing procedure is applied to obtain an initial list of degs. these degs are then subject to a tukey’s test or a bonferroni test to identify the pairwise comparisons where a deg is significant. the directional decision regarding whether a gene is up or down regulated, is made on the basis of the direction of the fold change. while this strategy is more intuitive than the earlier one, it is potentially a conservative procedure resulting in a substantial loss of power, firstly because intrinsically anova compares all treatment groups. secondly, the thresholds or the critical values at the two steps are not derived with the goal to control the overall error rate. the only formal methodology available in the literature that controls the mdfdr for the above directional multiple testing problems is the method by guo et al., which is designed to make decisions on thousands of features when making multiple pairwise comparisons and deciding on the direction of comparison. thus the method controls the false discovery rate when making multiple pairwise comparisons on thousands of genes while also controlling the directional errors committed when falsely declaring a deg to be up-regulated when it is not. guo et al. procedure generalizes the procedure of benjamini and yekutieli which was designed for the case when there were only two groups to compare. the procedure of guo et al. is available in the software oriogen that can be accessed from http://www.niehs.nih.gov/research/atniehs/labs/bb/staff/peddada/. while guo et al. methodology is useful for making several multiple pairwise comparisons; it is relatively conservative since it relies on the bonferroni procedure to deal with multiple pairwise comparisons within each gene. in this paper we develop a general mdfdr controlling testing procedure that allows us to use any mixed directional familywise error rate controlling procedure in place of the bonferroni procedure, for conducting pairwise comparisons in high dimensional data that is broadly applicable to a wide range of genomic data including gene expression microarray data, cpg methylation data, rna-seq data and others. based on this general procedure, we develop a specific methodology using the dunnett’s test which is designed for making comparisons of several experimental groups with the control or the reference group. not only that the resulting methodology is practically relevant but as demonstrated in the numerical simulations, the resulting methodology not only controls the mdfdr but is more powerful relative to some potential alternative methods. we illustrate the proposed methodology by analyzing a gene expression data obtained from the niehs’ fibroid growth study. using our methodology we gain deeper insights into molecular characteristics of uterine fibroids according to the tumor size, which was not understood until now. we have identified several differentially expressed genes and pathways that are specifically enriched according to the tumor size. while researchers and clinicians who study fibroids are well aware of many of the genes and pathways described in this paper, we have provided a characterization of these genes and pathways according to the tumor size. our data can be further mined to gain deeper insights regarding fibroids. methods in this section we present the general testing procedure that controls mdfdr while making multidimensional directional decisions, and develop the statistical methodology based on dunnett’s test that is used to analyze fgs gene expression data. the relevant notations, definitions, hypotheses and other statistical concepts are described in the online additional file and additional file 1: figure s the general mdfdr controlling procedure the general testing algorithm:use any global testing method to obtain the screening p-values. apply bh procedure to identify genes that are differentially expressed in at least one pairwise comparison. let r denote the number of genes so discovered. for each gene discovered in step use any two-sample testing method to obtain the p-values for each pairwise comparison and apply any mixed directional family wise error controlling procedure, such as holm, hochberg etc., to the pairwise p-values at level rα/m. for a given gene discovered in step if a pairwise hypothesis is rejected in step then we declare the gene to be up or down regulated in the pairwise comparison according to the sign of the corresponding test statistic. this procedure is general in the sense that any global testing method can be used to obtain the screening p-values in step of the procedure, any pairwise comparison testing method can be used to obtain the pairwise p-values in step and any mdfwer controlling procedure can be used in steps 2- of the procedure. it is important to note here that the method of guo et al. is a special case of the proposed general procedure in which bonferroni global test is used for testing the screening hypotheses and bonferroni method along with additional directional decisions works as the mdfwer controlling procedure. this general algorithm is proved to control mdfdr through the following theorem with the mathematical proof provided in the online additional file 6: section s theorem under assumption of independence of p-value vectors pj, j= …,m, the mdfdr of the general testing algorithm is strongly controlled at level α. statistical methodology for fgs gene expression data suppose we are interested in comparing “q” experimental groups with a reference group on the basis of the mean expression of “m” genes. for example, suppose we are interested in comparing “small”, “medium” and “large” fibroids with a “normal” tissue from uterus on the basis of “m” genes. our goal is not only to identify differentially expressed genes in any given pairwise comparison but also to determine if the mean expression is up or down-regulated in the tumor tissue compared to the normal myometrium. our statistical methodology for the fgs gene expression data analysis proceeds in three steps as follows. specific details of implementation of each step are described in the online additional file 6: section s for each gene we obtain a dunnett based screening p-value from all “q” pairwise comparisons with the reference group. this results in one p-value for each gene. applying the bh procedure on these screening p-values, we obtain genes that are differentially expressed in at least one pairwise comparison. suppose we discover r genes in this step. thus there are r genes which are differentially expressed in at least one pairwise comparison with the reference group. for the jth gene discovered in step we compute dunnett’s p-value pij for each pairwise comparison i, i=,…,q, with the reference group. if pij≤rαm then we declare that the ith pairwise comparison with the reference group is significant. if a gene j is found to be significant in the ith pairwise comparison with the reference group then we declare it to be up-regulated in the ith group relative to the reference group if tij> otherwise it is declared to be down-regulated. here tij denotes the test statistic associated with the jth gene in the ith pairwise comparison. remarks. it is important to note in steps and that it is possible for a gene to be differentially expressed on the basis of the screening p-value, yet it may not be significant in any of the pairwise comparisons. this happens when the bh-adjusted p-values are only marginally significant. in such cases we declare that gene to be not differentially expressed. we expect this phenomenon to occur with very small frequency.  in this section we present two sets of numerical results. first, the operating characteristics of the proposed procedure based on a simulation study are presented. next, we apply the proposed methodology to a gene expression data obtained from the niehs’ fgs to identify tumor size specific degs. numerical simulation and discussion in this section we present the operating characteristics of the proposed procedure based on a simulation study. this study evaluates the power of the proposed dunnett based methodology and compares its performance with some competing procedures, the details of which are given in the online additional file 6: section s in our simulation we considered p= groups, the first three were taken to be experimental groups and the last group was taken to be the reference group. thus all pairwise comparisons are made with the last group. our simulated microarray chip consisted of m= genes per chip with n= chips per group. as often done in simulation studies for microarray data, we generated the expression of each gene in each chip using a normal distribution. more precisely, for the jth gene, j=,…, in the sth sample, s=,…, in the gth group, g=,, we generated its expression zgjs from a normal distribution with mean value ezgjs=μgj and variance vzgjs= for the reference group we set μ4j= for all j. to create the null data we set μgj= for j=,…,m g=, and non-null data were created by generating μgj∼independentu, for g=, and j=m0+ m0+ …,m0+m where m1=m−m and u represents the uniform distribution. it is important to note that, for the non-null means considered here, the standard deviation used in this simulation study is large. consequently, all methods considered in this simulation study are expected to have small power. we considered three patterns of correlation structure as follows: independent gene expressions: the correlation coefficient ρ between any pair of genes and any pair of sample is thus, the data are completely independent. gene expressions within sample are dependent: for a given sample s, the correlation coefficient between any pair of genes is ρ but the correlation coefficient between any pair of samples is gene expressions within genes are dependent: for a given gene j, the correlation coefficient between any pair of samples is ρ but the correlation coefficient between any pair of genes is although we considered a variety of patterns of ρ, due to space limitations we only present the results corresponding to independence. similar results are obtained for other patterns of ρ. we simulated the mdfdr and average power using simulation runs. our nominal mdfdr level was taken to be α=. details regarding the calculations of the mdfdr and power are provided in the accompanying online additional file for comparison purposes, we compared the proposed procedure with the guo et al. procedure and variants of the proposed procedure in which the mdfwer controlling procedures used in steps and of the algorithm are respectively, the holm’s procedure, the hochberg’s procedure and the bonferroni procedure. these procedures are described in detail in the online additional file we summarize the results of the simulation studies in fig. for the independence case. the horizontal axis denotes the proportion of truly differentially expressed genes on the array and the vertical axis denotes the average mdfdr and average power. as desired, all five procedures control the mdfdr on average at α=. however, the proposed dunnett based procedure has highest power compared to all other methods considered fig. mdfdr, average power. mdfdr, average power with the proposed methodology and three variants using holm, hochberg and bonferroni procedures, respectively, in steps and along with guo et al. procedure, under independence among genes and within genes we summarize the results for the case of dependence within genes across groups in fig. additional file 2: figure s and additional file 3: figure s and the results for the case of dependence among genes in fig. additional file 4: figure s and additional file 5: figure s the details of these dependence structures are given in the online additional file 6: section s although all five procedures control the mdfdr on average at, the power of our proposed methodology is much higher compared to all four competing procedures under the different dependence structures considered. fig. mdfdr, average power. mdfdr, average power with the proposed methodology and three variants using holm, hochberg and bonferroni procedures, respectively, in steps and along with guo et al. procedure, under dependence within genes fig. mdfdr, average power. mdfdr, average power with the proposed methodology and three variants using holm, hochberg and bonferroni procedures, respectively, in steps and along with guo et al. procedure, under dependence among genes analysis of gene expressions in uterine fibroids uterine fibroids are benign smooth muscle tumors occurring in the uterus of women with a very high prevalence rates in many populations. these tumors are considered to be hormonally mediated and therefore occur during the pre-menopausal years. according to some estimates, the cumulative incidence by age of these tumors among caucasian women exceeds % and it is much higher among women of african american descent. the direct and indirect annual cost of fibroids in the us is as high as billion dollars, yet these tumors are not well-studied perhaps because they are considered to be benign. the niehs conducted a large prospective study tracking the growth of fibroids in pre-menopausal women over months and discovered difference in growth patterns of fibroids according to age for white and black women. a total of tumor samples and normal myometrium samples were collected from women who opted for either hysterectomy or myomectomy. the rna from these tissue samples were processed to obtain gene expression data using affymetrix chips. a total of probe sets were obtained. using these gene expression data davis et al. studied the molecular characteristics of fibroids. several important degs and canonical pathways were identified in davis et al., which were also experimentally validated using rtpcr. in this paper we take the next step towards the understanding of molecular characteristics of fibroids by investigating the changes in gene expression according to tumor size. as done in the literature for other tumors such as breast cancer tumors and others, we are interested in identifying genes that are specific to size of uterine fibroids. we classified fibroids into three size groups according to their volume, small, medium and large and compared each size category with the normal myometrium, the reference group, using the proposed methodology. the cut points between size categories were chosen so that there is a natural separation between the size categories. the sorted volumes of the tumors are provided in fig. due to large differences in the volumes between large and small tumors, for clarity we provided volumes in three panels according to their size. fig. volume of fibroids according to size category we identified a total of probe sets to be differentially expressed in at least one pairwise comparison at mdfdr of. these probe sets map to genes. the venn diagram of the degs by tumor size is in fig. based on the genes, using the ingenuity pathway analysis, we discovered a total distinct enriched canonical pathways. the venn diagram of the number of enriched canonical pathways by tumor size is provided in fig. fig. venn diagram of differentially expressed genes by tumor size fig. venn diagram of enriched pathways by tumor size it is not surprising that a third of the degs and nearly % of all enriched pathways are common to tumors of all sizes since tumor tissues are fundamentally different from the normal myometrium. however, we discover several degs and pathways that are uniquely enhanced according to the tumor size, suggesting that changes in molecular characteristics might be taking place as tumors grow. we highlight some important groups of genes and pathways discovered by our methodology and provide the complete list of all genes and pathways in the additional files and respectively, for researchers to mine our data. the “gene symbol” and “gene name” columns in these tables are obtained from the website david. growth factors, vascularization and related pathways: the netrin signaling pathway is well-known as a versatile pathway with multiple functions. one of its functions is to promote vascular networks and branching of blood vessels and angiogenesis. according to our analysis this pathway is uniquely enriched in small tumors only. the interleukin- pathway is known to induce inflammatory response and the production of prostaglandins and expression of il- which may play a critical role in the fibroid initiation and early development. for example, prostaglandins play a critical role in the promotion of growth factors involved in angiogenesis, such as vegf, basic fibroblast growth factor and platelet-derived growth factor and development of fibroid requires vascularization and blood supply. thus interleukin- pathway likely plays an important role during fibroid development. interestingly, according to our analysis, this pathway is uniquely enriched in small size tumors but not in the medium or large tumors. furthermore, the fibroblast growth factors and which belong to the regulation of the epithelial-mesenchymal transition pathway and are well-known to be involved in vascularization and angiogenesis, are both uniquely down-regulated in small tumors and not differentially expressed in medium or large tumors. our analysis further implies that the epithelial-mesenchymal transition pathway was enriched only in small and medium size tumors. according to ciarmela et al., estrogen may promote fibroid growth through up-regulation of epidermal growth factor receptor. however, we found egfr to be down-regulated in fibroids and that too only in the medium size tumors. similarly, the fibroblast growth factors were differentially expressed only in the medium size tumors. the acidic fgf was up-regulated whereas the basic fgf was down regulated. these findings are consistent with ciarmela et al. in that they are expressed during tumor progression. similarly, insulin like growth factor was only differentially expressed in medium size tumors. additionally, growth factor signaling pathways such as vegf, pdgf, tgf β and egf are uniquely enriched in medium size tumors and not in small or large tumors. while the above results suggest that some growth factors and related pathways are very specific to small and medium tumors, we discovered several growth factors to be differentially expressed in tumors of all sizes. these included, insulin like growth factor, insulin like growth factor binding protein, and platelet derived growth factor c which were up-regulated whereas insulin-like growth factor binding protein, connective tissue growth factor, heparin-binding egf-like growth factor, transforming growth factor beta receptor ii, fibroblast growth factors and were down-regulated. similarly, growth factor signaling pathways such as the human growth factor and igf- were enriched in tumors of all sizes. thus it appears that the differential expression of these genes and the enrichment of the above pathways is necessary for tumor onset and progression. estrogen and related genes: fibroids are hormonally mediated and it is also well documented in the literature that accordingly estrogen and progesterone receptors and prostaglandins promote proliferation of fibroids. not surprisingly, we found the estrogen receptor esr- to be up-regulated in tumors of all sizes. this result was confirmed using the rtpcr data. our rtpcr data further indicates that the expression of esr- in small tumors is almost % higher than in large tumors, suggesting that esr- has a large effect on small tumors but is also needed for the continued growth of the tumor. interestingly, the progesterone receptor was up-regulated in only medium size tumors and not differentially expressed in small or large tumors, a result that is also confirmed by the rtpcr. this suggests that perhaps pgr may not be involved in tumor initiation but is only involved in growth of the tumor. however, its function ends once the tumor becomes large enough. most prostaglandins were generally down-regulated in tumors of all sizes. for example, prostaglandin e receptor, prostaglandin f receptor and prostaglandin-endoperoxide synthase are down-regulated in tumor of all sizes. the result of ptgs was also confirmed using rtpcr. however, some prostaglandins were differentially expressed according to the size of the tumor. for example, prostaglandin-endoperoxide synthase was differentially expressed only in the medium sized tumors where it was down-regulated and prostaglandin e synthase was down-regulated in large tumors only but prostaglandin e receptor was down-regulated in both medium and large tumors. similar to estrogen and progesterone receptors, the α and γ isoforms of peroxisome proliferator-activated receptors have been associated with the regulation of proliferation of uterine fibroids. in our data both these isoforms are down-regulated in all tumor sizes compared to the normal myometrium. we also discovered the related retinoid x receptor gamma to be down-regulated in the medium size tumors but was not significant in small or large tumors. collagens: there is a vast amount of literature implicating collagens to smooth muscle tumors such as the fibroids. consequently, it is not surprising that several collagens and extra cellular matrix proteins are differentially expressed in tumors of all sizes. apart from col4a col4a and col21a which were down-regulated in all tumor size groups, the remaining collagens were up-regulated in tumor samples. we investigated a subset of these collagens using rtpcr which confirmed the above findings for those genes. there were other collagens that were differentially expressed depending upon the size of the tumor. for example, col2a col4a and col5a were not differentially expressed in small tumors but were significantly up-regulated in medium and large tumors. collagen col23a was significantly up-regulated in small and medium tumors but not in the large fibroids. interestingly, some of the collagens were differentially expressed in only one of the size categories. for example, col4a3bp, col6a col9a col10a and col11a were differentially expressed only in the large tumors. again, except for col4a3bp, the remaining collagens were up-regulated in the tumor samples. interestingly, col5a was only differentially expressed in small tumors. other genes: leptin receptor is well-known to be negatively associated with the obesity and obesity is a potential risk factor for fibroids. interestingly, we discover leptin receptor to be significantly down-reglated in all tumor sizes. this finding is also confirmed by the rtpcr data reported in additional file as noted earlier, diaz at al. demonstrated that β integrin had an increased expression in larger breast tumors and in higher tumor grades. in our fibroid data, however, we notice a down-regulation in β integrin in medium and large tumors and was not differentially expressed in small tumors. comparison of various procedures with dunnett’s test based methodology we analyzed the niehs fgs data using various methods considered in this paper at mdfdr level of. we compared the methods in terms of number of probe sets that were differentially expressed in; all tumors, small tumors, medium tumors, and large tumors. results are summarized using venn diagrams. in fig. we compare our dunnett test based methodology with guo et al. procedure. not only did the dunnett based methodology identified many of the probe sets identified by the guo et al. procedure, i.e. shared several probe sets, it identified several more that were not identified by guo et al. procedure. on the contrary, there were very few probe sets that were uniquely identified by guo et al. procedure. we see similar results when we compare dunnett screening and holm procedure, dunnett screening and hochberg procedure and dunnett screening and bonferroni procedure. the venn diagram in fig. compares the number of probe sets identified by dunnett based method to each of the other methods. in each case the dunnett based method tends to uniquely identify many probe sets than each of the other methods. the total number of probe sets identified by the various methods are summarized in table fig. venn diagram of identified probe sets by the dunnett based methodology and guo et al. method. the venn diagrams compare the two procedures in terms of the probe sets identified in at least one tumor size, large tumor size, medium tumor size and small tumor size fig. venn diagram of identified probe sets by the dunnett based methodology and other methods discussed in numerical studies. the other methods considered in the numerical simulation are dunnett screening and holm procedure, dunnett screening and hochberg procedure and dunnett screening and bonferroni procedure  this paper offers a general statistical methodology for conducting pairwise comparisons in high dimensional data that is broadly applicable to a wide range of genomic data including gene expression microarray data, cpg methylation data, rna-seq data and others and proposes a specific testing methodology for the analysis of fgs gene expression data. the proposed methodology not only controls the false discoveries made when making several pairwise comparisons on the basis of high dimensional data, but it controls for directional errors, i.e. falsely declaring a gene to be up-regulated when it is not. thus it controls the mixed directional false discovery rate. the framework we proposed here is general enough that, depending upon the application, a researcher can use different testing procedures within our framework. in all situations the procedure strongly controls mdfdr. for the application under consideration in this paper, the experimental design calls for the comparison of several experimental groups against a particular group. for such designs, in the classical multiple comparison problem, the dunnett’ test is known to be powerful. so, it is not surprising that our dunnett’s test based procedure is most powerful among all the procedures we considered in this paper. our methodology based on dunnett procedure works best here in terms of power because the dunnett’s p-values are based on the joint distribution of the test statistics and our procedure uses the maximum order statistic to find the screening p-value. using our methodology we gain deeper insights into molecular characteristics of uterine fibroids according to the tumor size. we have identified several differentially expressed genes and pathways that are specifically enriched according to the tumor size. while researchers and clinicians who study fibroids are well aware of many of the genes and pathways described in this paper, we have provided a characterization of these genes and pathways according to the tumor size. our data can be further mined to gain deeper insights regarding fibroids. in view of the recent publication we acknowledge that our results are potentially confounded by race. since the sample size by race is very small, we are unable to separate out the race effect. despite that potential weakness, we think that our findings provoke researchers to explore further research along these lines. additional files additional file figure s a graphical display of various hypotheses of interest. additional file figure s mdfdr, average power with the proposed methodology and three variants using holm, hochberg and bonferroni procedures, respectively, in steps and along with guo et al. procedure, under dependence within genes. additional file figure s mdfdr, average power with the proposed methodology and three variants using holm, hochberg and bonferroni procedures, respectively, in steps and along with guo et al. procedure, under dependence within genes. additional file figure s mdfdr, average power with the proposed methodology and three variants using holm, hochberg and bonferroni procedures, respectively, in steps and along with guo et al. procedure, under dependence among genes. additional file figure s mdfdr, average power with the proposed methodology and three variants using holm, hochberg and bonferroni procedures, respectively, in steps and along with guo et al. procedure, under dependence among genes. additional file supplementary materials. section s the mathematical notations and definitions used in the paper are described in this section. the hypotheses tested in different steps of the testing procedure are formulated and interpreted. section s the mathematical proof of the theorem that states that the general testing algorithm controls mdfdr. section s the details of the statistical methodology developed for the testing and analysis of the fgs microarray data. section s the supplementary results for the simulation study performed to evaluate the power of the developed statistical methodology compared to several other competing methods. additional file differentially expressed probe sets in different sizes of tumors compared to normal myometrium. significant probe sets are grouped and color coded according to different sizes of fibroids. the last column shows the bh-adjusted p-value for step of the methodology. the columns under "dunnett adjusted p-values" show the p-values for step of the methodology. additional file pathways identified in different sizes of tumors compared to normal myometrium. in the table, the enriched pathways are grouped and color coded according to different sizes of fibroids. the last column shows the negative of the logarithm of the bh adjusted p-value, computed by ipa. additional file rtpcr data. this table shows the rtpcr data for a subset of the genes that are identified in our analysis. competing interests the authors declare that they have no competing interests. authors’ contributions ag, wg and sdp conceived the study. ag performed all calculations and all numerical work. ag, wg and sdp wrote the paper. all authors read and approved the final manuscript. 
our powerful methodology allowed us to draw several interesting conclusions regarding the molecular characteristics of uterine fibroids. often researchers are interested in comparing multiple experimental groups with a reference group on the basis of thousands of features and determine if a differentially expressed feature is up or down regulated in a pairwise comparison. a complete list of all degs and the corresponding enriched pathways according to tumor size is provided for researchers to mine these data. there are two sources of false discoveries, one due to multiple testing involving several pairwise comparisons and the second due to falsely declaring a feature to be up regulated when it is not. thus we are providing a general framework for making multiple pairwise comparisons. together, the total error rate is called the mixed directional false discovery rate. our analysis of the uterine fibroid growth gene expression data suggests that molecular characteristics of a fibroid changes with size. it appears that the molecular apparatus necessary for fibroid growth and development is established during tumor development. for example, il- signaling pathway, associated with inflammation and known to activate prostaglandins that are implicated in the progression of fibroids, is significantly enriched only in small tumors. the proposed general procedure strongly controls mdfdr. several specific methodologies can be derived from this general methodology by using appropriate testing procedures at different steps of the general procedure. mdfdrbh proceduremultiple testingdifferential gene expressiontumor sizeissue-copyright-statement© the author 2016
 often researchers are interested in comparing multiple experimental groups with a reference group on the basis of thousands of features and determine if a differentially expressed feature is up or down regulated in a pairwise comparison. there are two sources of false discoveries, one due to multiple testing involving several pairwise comparisons and the second due to falsely declaring a feature to be up regulated when it is not. together, the total error rate is called the mixed directional false discovery rate.  we develop a general powerful mdfdr controlling testing procedure and illustrate the methodology by analyzing uterine fibroid gene expression data. we identify several differentially expressed genes and pathways that are specifically enriched according to the size of a uterine fibroid.  the proposed general procedure strongly controls mdfdr. several specific methodologies can be derived from this general methodology by using appropriate testing procedures at different steps of the general procedure. thus we are providing a general framework for making multiple pairwise comparisons. our analysis of the uterine fibroid growth gene expression data suggests that molecular characteristics of a fibroid changes with size. our powerful methodology allowed us to draw several interesting conclusions regarding the molecular characteristics of uterine fibroids. for example, il- signaling pathway, associated with inflammation and known to activate prostaglandins that are implicated in the progression of fibroids, is significantly enriched only in small tumors. it appears that the molecular apparatus necessary for fibroid growth and development is established during tumor development. a complete list of all degs and the corresponding enriched pathways according to tumor size is provided for researchers to mine these data. identification of these degs and the pathways may potentially have clinical implications. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords mdfdrbh proceduremultiple testingdifferential gene expressiontumor sizeissue-copyright-statement© the author 2016 