BACKGROUND
one of the key aspects in the study of cancer is to understand the principles and mechanisms of gene expression variation contributing to cancer genesis and progression.
the identification of genes differentially expressed between normal and cancer cells/tissues is not only helpful for designing diagnostic and therapeutic procedures, but also for understanding cancer biology as a whole.
in this regard, dna microarrays have been the dominating platform in the high-throughput study of cancer transcriptomes since their emergence in the mid-1990s  <cit> .
however, there are several drawbacks, which include: high background level signals resulting from cross-hybridization  <cit> ; difference in hybridization properties due to different probe sequences; limited dynamic range due to background level and saturation, and difficulty in detecting splicing isoforms and unknown genes.
for these reasons, with the advancement of the next generation sequencers, we are seeing high-throughput transcriptome mapping and quantifying method, also known as rna-seq, to begin to supersede microarray in expression profiling.
however, rna-seq experiments are relatively demanding in terms of time, cost, and computation equipment.
experimental differences between different sequencing platforms may complicate transcriptome analysis with multiple tissue sources.
since exploring meta-analysis from traditional digital expression data such as est derived from cdnas  <cit>  is more feasible, this study may serve as a precursor to more complicated experiments.
originally primarily aimed for cataloging of transcript repertoire, ests from large-scale cdna sequencing projects such as cancer genome anatomy project , human cancer genome project , and cancer genome project  also allow searching for differentially expressed genes  in specific tissue types or in whole genomes  <cit> .
several in silico analysis tools such as ncbi unigene cdna xprofiler  <cit> , cgap digital differential display   <cit> , and cgap digital gene expression displayer   <cit>  are available online allowing the analysis of publicly available data.
while standard statistical methods such as fisher's exact test for finding degs in two-class problems  or pearson's correlation are commonly used  <cit> , there are also specially developed methods for finding degs in the landscape of digital signals for two-library problems  <cit>  or for multiple libraries  <cit> .
the online tools as well as the statistical methods remain useful to this day in est or even rna-seq projects  <cit> .
aside from searching for degs, the searches for gene transcript isoforms specific to particular libraries were also demonstrated and many of these attribute differentially expressed isoforms to human cancers  <cit> .
in spite of the successful applications, these tools or methods are not without limitations.
xprofiler reports differential expression in an all-or-none manner where only a list, but not statistical quantification, of candidates is reported.
ddd allows quantification using fisher's exact test.
however, the nature of the test dictates that comparisons of three or more libraries involve multiple pair-wise comparisons, and thus there are no easy comparisons of library specific genes.
dged uses a bayesian approach to find degs, but it is also pair-wise.
the reported "odds ratio" is perhaps better described as "relative risk" and may be biased with unequal sampling.
another popular and useful bayesian-based method originally developed for est analysis by audic and cleverie  <cit>  is also popular for rna-seq data.
it is less conservative than fisher's exact test, but it also does not apply to multi-class problems.
the multi-class comparison method established by stekel et al.
<cit>  finds specificity in one condition out of all and is useful in application such as finding degs in multi-tissue libraries.
however, in the search for cancer degs, a subsequent analysis of differential expression between cancer and normal libraries of the tissue of interest may not yield fruitful results due to the possible scarcity of est sampling in the particular tissue type.
on the other hand, the naïve method of pooling all data into the two-class problem of normal versus cancer when searching for differentially expressed genes or differentially splice variants  <cit>  risks introducing bias.
in extreme cases, one may encounter the fallacy of simpson's paradox  <cit>  where genes in reality more active in the normal condition appear to be more so in the cancer condition .
we now report on the application of a computational and integrative approach to analyze cancer differentially expressed genes .
the statistical method we employed is cochran-mantel-haenszel statistics   <cit>  and to the best of our knowledge has not been applied in this context.
instead of pooling all normal and all cancer ests from different tissue types to fit into a two-class problem as by using the  <dig> by  <dig> contingency χ <dig> test or the fisher's exact test, cmh allows original stratification of libraries in their respective tissue types, yet exhaustively analyzes expression between cancer and normal conditions across all tissue types.
the method is an extension to χ <dig> test which, in our application, measures the association between cancer and gene expressions, adjusting for the tissue confounding factor.
this approach allows one to find genes that are overall differentially expressed in cancer, or multiple-cancer genes, irrespective to a specific tissue type.
the method is demonstrated in this paper to exhaustively analyze ests from the dbest database  <cit> .
to the best knowledge of the authors, such an all-inclusive, whole-transcriptome analysis has not been redone in recent years now that more ests than ever are available.
our filtering of est libraries was also more rigorous than many previous studies.
notably, we excluded the orestes  libraries  <cit>  on which a normalization procedure had been applied.
libraries from cell line were also excluded owning to their unrepresentativeness of primary cancer cell transcriptomes.
our analysis pipeline further focused on enrichment of the degs by cross examination with expression data of a different platform, i.e.
the microarray data, and selecting for membrane and secretory associated protein genes since we intend to find therapeutic targets or biomarkers, and conducting string  network analysis to show the cancer enriched clusters  <cit> .
with real-time qpcr validation, we have identified three candidates that are inclined to express in cancer across more than one tissue types.
we hope such a meta-analytical and multiple-tissue comparison can serve as an exploratory experiment for future multi-library or multi-tissue study of other digital sources such as rna-seq.
methods
overview
our approach was to exploit the entire collection of human est sequences from dbest  <cit>  to obtain transcripts from different type of cells/tissues/organs.
the assumption was that the activities of the genes can be represented by their transcripts, and also reflected by the number of representing ests in the ncbi dbest database, given that a large number of mrnas  were sequenced.
pertinent sequences from different sources were matched to genes and tallied together.
through the annotation of each est record, we obtained the tissue type and condition type  from which it was derived.
with the information, we then had the entire gene transcription profile for all the tissues and conditions.
next, cross examining data of other sources including microarray data, secretory and membrane associations as well as analyzing protein associations with string  <cit>  allowed us to narrow down the list of candidate genes.
the process is illustrated in figure  <dig>
human gene reference sequence preparation
the ncbi reference sequences   <cit>  were downloaded from its ftp site  <cit> .
homo sapiens refseq records were selected and subjected to repeat masking via repeatmasker  <cit> .
human est sequence preparation and library filtering
human est data  and their cdna library information were downloaded from ncbi dbest database  <cit>  and cgap  <cit> .
program in python language was written to mark for discard the unsuitable libraries when the keywords such as "enrichment", "subtract", "pcr", and "normalized" were found in the descr, unique_protocol, or keywords fields of the library information.
an arbitrary cutoff of >  <dig> was chosen to the highly unrepresentative libraries .
to curb from incorrect inclusions or exclusions, we finalized the process with manual curation.
libraries made from mixed tissues or cell lines were also discarded.
the final libraries from cgap were manually classified into  <dig> different tissue types and two different conditions, normal and cancer.
est to gene assignment
the blat alignment tool was used to align ests to refseqs as a mean to assign ests to genes  <cit> .
the criteria of having an identity of 95% or above and the minimum length of  <dig> nucleotides were set for a match.
the refseq match with the highest identity was assigned for the est.
if two refseq matches shared exactly the same identity, the program chose the first encountered.
est count and summarization
the procedure attributes each transcript represented by refseq its expression profile across different tissue and condition types based on est assignment counts.
each est has its corresponding tissue type and condition type classification, based on its source clone library.
for example, a transcript with an aligned est from a lung cancer clone library is one expression count each in tissue type lung and condition type cancer.
this way, after all ests were counted, each transcript has a profile of expression across various libraries and conditions.
expressions from different transcript variants of the same gene were pooled to obtain a single gene expression.
the raw counts were thus made into transcription profile for each gene for further statistical analysis.
statistical evaluation of cancer candidates
cochran-mantel-haenszel statistics  was applied to evaluate cancer differential expression of each gene.
to evaluate each gene, other genes were pooled as "other genes" to create a  <dig> ×  <dig> × k table consisting of data from tissue-condition cross, where k was the number of tissues ×  <dig> .
a contrived example of  <dig> ×  <dig> × k table where k is  <dig> is shown in table  <dig>  gene a is the gene under study while other genes are pooled together as "other genes".
only tissue i and tissue ii columns are calculated in cmh.
the pooled ones are not part of the analysis.
akin to fisher's exact test, the test assumes that "other genes" should consist mostly of genes not differentially expressed between normal or cancer conditions.
or, some of them are degs for one condition, but they are at least partly canceled out by degs for the other.
in any case, the imbalances of cancer counts to normal counts in the second row is regarded as owning to sample bias and it serves as a metric against which gene a is measured.
by continuously isolate values for gene currently under study while pooling all other genes to the second row, an odds ratio and a confidence interval is calculated for each gene.
genes with a p-value <  <dig>  and an lower bound of confidence interval of odds ratio ≧  <dig>  are selected for further analyses.
this hypothetical case serves both as an example of how cochran-mantel-haenszel  is applied as well as the occurrence of simpson's paradox.
gene a is the gene under investigation.
expressions from all other genes are pooled into the "other genes" row.
bold typeface indicates columns showing higher cancer vs. normal propensities.
cmh is applied on the stratified tissue columns .
a casual observation involving only the pooled data would suggest gene a as having higher expression in cancer .
however, a closer inspection on each of the tissue columns reveals otherwise.
the observed difference between cancer and normal of the "other genes" is theoretically mostly due to sampling bias.
microarray cross reference
human u <dig> plus  <dig>  genechip array cel data were downloaded from gene expression omnibus   <cit> .
when computing power allows, the data were processed with affyplm  <cit>  using its three-step procedure of processing background signals with gcrma, normalizing signals with quantile normalization, and summarize probe signals with medium polish.
for large experimental datasets that were computationally infeasible for us, we used justrma from the affy package  <cit> .
for experimental dataset without raw cel data, we obtained the pre-processed matrix files via geoquery  <cit> .
regardless of the source of array signal processing, we analyzed the genes for differential expression with limma  <cit> .
differentially expressed gene candidates with p-value <  <dig>  and logfc >  <dig>  were selected and crossed with genes from est profiling with statistical evaluation.
for each array, the significant genes were crossed with our est profiling results.
the union of these intersecting genes was selected for further evaluation.
annotation of secretory proteins
to identify our differentially expressed genes with secretory annotation, a list of  <dig>  proteins with secretory annotation originated from the conglomeration of data from uniprot   <cit> , human plasma proteome organization  , and secreted protein database    <cit> .
this list was matched against degs to give them secretory annotation.
annotation of membrane proteins
membrane protein annotations were gathered from five sources - topdb   <cit> , locate   <cit> , pdb_tm   <cit> , opm   <cit> , and mpdb   <cit>  - to generate a unique list of  <dig>  membrane proteins.
any degs on this list would confer it a membrane annotation.
validation of tissue expression profiles of candidate genes
tissuescan™ cancer survey panel 96-i qpcr array panel  containing the cdnas of  <dig> normal and  <dig> cancer tissues each from  <dig> organs  was used to examine the expression profiles of selected cancer differentially expressed gene candidates.
real-time qpcr analyses with the taqman® gene expression assay kits  and fam- and vic-labeled target genes and hprt <dig> internal control primers, respectively, were performed according to the manufacturer's suggested procedure on an applied biosystems prism  <dig> system.
relative specific gene expression was quantified by normalization against the hprt <dig> with the Δct method.
gene expression changes were quantified as  <dig> - .
RESULTS
human ests selection and tissue distribution
the basic steps of our analysis are illustrated in figure  <dig>  a total of  <dig> , <dig> human est sequences  were downloaded from the ncbi.
despite the size of the data, not all ests are relevant for our gene expression analysis.
after screening the  <dig>  est libraries as described in the methods section above,  <dig>  unsuitable libraries, the preparation of which involved pcr amplification, normalization, subtraction, etc.
or originated from cell lines, were discarded.
the remaining  <dig> libraries consisted of  <dig> , <dig> est sequences representing approximately a third of all the downloaded human ests.
after blat alignment of the  <dig> , <dig> ests to  <dig>  gene transcripts from refseqs, approximately  <dig> , <dig>  ests with at least  <dig> nucleotides matched to refseqs were detected.
an examination of the sources of the matched ests indicated that the representativeness of each tissue is skewed and that the brain is the most represented out of all tissues.
among the  <dig> different tissues, brain ests constituted 26% of all matched ests, uterus  ranked second, followed by testis , placenta , pancreas , muscle , liver , kidney  and others each below 3% .
similarly, condition type  representation was also skewed.
normal tissue type had  <dig> , <dig> ests combined, and cancer tissue had  <dig>  ests in the ratio of roughly  <dig> to  <dig>  originally before filtering out those from the cell lines, there were more cancer ests and the ratio of normal ests to cancer ests was roughly  <dig> to  <dig>  this showed how much more rigorous our filtering was.
unfortunately, this also meant we had a much smaller dataset to work with.
the unequal distribution of the  <dig> , <dig> matched ests in different tissue types caused some tissue types to be ill-represented.
for example, the number of brain est hits dominated over other tissue types.
on the other hand, spinal cord had the least count with  <dig> est hits.
the latter had little value for our application.
therefore, we only took a tissue type into consideration when its total est hit count was above the cut-off of  <dig> .
considering that the human genome has approximately  <dig>  genes, the cut-off still did not allow "deep" probe into gene expression.
nevertheless, the method we employed did not attempt to identify specific gene expression in one particular tissue; therefore, the problem was mitigated.
we also categorized ests according to their clone library classification, to either be from normal or from cancer.
sometimes a certain tissue-condition type was so under-represented that the information was not trustworthy.
for example, adipose had  <dig>  normal hits but only  <dig> cancer hits, and heart tissue had  <dig>  normal hits but no cancer hits.
for these cases, data was kept throughout the analysis.
but these data did not make contribution to our analysis.
since our est assignments were made to transcripts represented by refseq sequences, when the entire assignment procedure was done, each transcript variant had its expression profile across all tissue-condition types.
due to the lack of enough ests data, differentiating between different splicing variants of the same gene was not feasible.
we had to pool expression from different splicing variants into a single expression profile representing the gene.
analysis of differentially expressed genes
due to the small sample size , it was only realistic to evaluate gene expression based on all ests of all tissues.
however, tissue type was a confounder.
if all counts for each gene were pooled as "normal" or "cancer" regardless of the tissue of origin, the count would be incorrect.
to solve both the sample size and the tissue confounder problems, cochran-mantel-haenszel statistical method was employed to identify genes with differential expression as described in the method.
we used the arbitrary cut-offs of p-value <  <dig>  and odds ratio ≧  <dig>  to obtain a primary set of candidates.
as a result, a total of  <dig> cancer differentially expressed gene candidates were selected.
the  <dig>  cut-off is chosen based on a good coverage to a list of well-known biomarkers or genes known to associate with cancer .
to show that this list of  <dig> genes was enriched for cancer and thus obtains credibility for our methodology, we looked for cancer related pathways associated with them in genego  <cit>  pathways, which covered  <dig> signaling and metabolic networks .
among the  <dig> most significantly matched pathways, several are cancer related - pathways number  <dig>   <dig>  and  <dig> involve immune response; number  <dig> and  <dig> involve cytoskeleton remodeling; number  <dig> is transition and termination of dna replication; and number  <dig> and number  <dig> are adhesion related.
in addition, the result of genego disease enrichment analysis  indicates our set of genes as neoplasm enriched: seven out of the  <dig> most associated diseases are related to cancer.
the disease ranks the highest is neoplasms, followed by neoplasm by site, and digestive systems neoplasm.
this list reveals that our  <dig> degs covers general neoplasm related functions, and not specific to any particular neoplasm, as digestive, urogenital and breast are all covered.
to narrow down this list of biomarkers, we crossed examined the expression profiles of the candidates with the differentially expressed genes in  <dig> microarray experiments, i.e.
two each of ovary and uterus, and one each of pancreas and colon .
these tissue types were selected based on the following reasons.
we noticed that many of our candidate genes had the most expression in ovary tissue .
the other concern was the number of ests.
since our candidate genes were derived from est sampling of various tissue types, they were influenced more heavily by tissue types with more est representation due to deeper sampling from them.
therefore, the rest of the tissue types were selected based on their representativeness.
of the  <dig> degs,  <dig> candidates were also found to be differentially expressed genes in our microarray analysis.
n: normal, c: cancer
gse <dig> has two entries since we compared pair-wise between normal vs. benign and normal vs. malignant.
since membrane and secretory proteins could be potential therapeutic target or serum biomarkers, the subcellular location of the  <dig> degs were examined against the secretory and membrane protein lists consolidated from public databases.
among these,  <dig> degs were putative membrane or secretory proteins -  <dig> had only secretory annotation,  <dig> had only membrane annotation and  <dig> had both.
literature search and string analysis of the  <dig> degs
to further examine whether the  <dig> membrane/secretory degs identified in our est database mining had enriched cancer-related genes, we searched the literatures for known associations with cancers.
in additions, they were also analyzed with string for interactions, which are based on experimental evidence or prediction, such as conserved genomic neighborhood, gene fusion, co-occurrence across genomes, pathways, protein complex, co-regulation, or other literature sources such as co-mentioning.
the network of the string interactions of the  <dig> degs together with the literature search results were plotted based on the combined string score with cytoscape  <cit>  .
approximately  <dig> proteins formed a big cluster of interacting proteins and a large proportion of the degs  had published cancer association with clinical or non-clinical experimental supports.
this demonstrates the value of our integration strategy since we had an ample of literature supports.
the  <dig> degs were selected out of their general cancer propensity without necessarily referring to any particular tissue type.
however, we can still assess the general tissue distributions shown in figure  <dig>  a gene has a tissue representation if any est from a clone library of the tissue type is matched to it.
we can see that some genes are observed across many tissue types.
a gene could be observed across a variety of tissue types if it is pan-tissue, and its expression measure is relatively abundant.
separately, woolf's test for heterogeneity can also give hints to whether a gene is pan-cancer.
those that were found as significant in this test were considered having unequal representation in different genes; although whether they are pan-cancer require further evaluation.
three candidates had higher expression in several cancer tissues
three cancer differentially expressed secreted protein gene candidates, col3a <dig>  chain), dlg <dig> , and rnf <dig> , which had an odds ratio of  <dig> ,  <dig> , and  <dig> , respectively, and with limited or no clinical support were selected for real-time qpcr analysis using the taqman® gene expression assay kits  .
with the hprt <dig> as the reference, higher expressions of these genes were noticed in at least some of the cancer tissues.
apparently, the average relative expression levels of col3a <dig> in breast, liver, thyroid cancer samples were higher than their normal counterparts.
the average expression levels of dlg <dig> in breast, kidney, liver, lung, and ovarian cancers, and rnf <dig> in colon, liver, lung, ovarian, and prostate cancers were also found to be higher than their normal tissues.
the expression of col3a <dig> in approximately  <dig> of the liver cancers, dlg <dig> in  <dig> of the liver,  <dig> lung and  <dig> ovarian cancers, and rn <dig> in  <dig> of the colon,  <dig> ovarian and  <dig> prostate cancers seemed to have higher expressions than the normal tissues.
in light of the limited sample size, the three candidates appear to have an overall higher expression in cancer tissues.
discussion
reported here is an integrative, meta-analytic approach for the discovery of pan-cancer differentially expressed gene candidates.
our primary enrichment included a set of  <dig> degs with cancer associations supported by genego disease and pathway analysis.
further integrative evaluations with cancer differentially expressed genes suggested by microarray data narrowed the list down to  <dig> genes, and among these there were  <dig> degs likely belonged to either secretory and membrane protein genes.
further string protein network analysis and literature reviewing indicated 71% of the  <dig> degs were highly connected and many of them were associated with cancers in previous publications.
simpson's paradox
the meta-analytic nature of our study brought us the opportunities as well as challenges to study the digital signatures of various transcriptomes in a new perspective.
comparing to experimental methods that focus on a single tissue type or limited tissue types, our approach allows us to find genes inclined to express in cancer in a pan-tissue manner.
an important challenge of our approach is to avoid simpson's paradox which can occur in a meta-analysis study  <cit> .
simpson's paradox is where the association between two variables may show a correlation that is reversed in direction from what is observed from stratified subgroups.
a contrived example is shown in table  <dig> in which gene a appears to have a higher cancer expression when pooled, but it is in fact not so under individual stratified sub-tables.
this may be somewhat of an extreme case, where directionality of the ratios actually differs between the sub-tables and the pooled table.
however, the tissue confounder still introduces bias, large or small, that may throw our judgment off.
in this study, we used cmh to analyze the data based on stratified sub-tables to avoid running into this paradox.
one could also analyze only one tissue at a time for differential expression, but this means one has a smaller dataset to work with.
cmh could avoid this problem since it uses est counts from all tissues instead of analyzing just the normal and cancer propensity under each individual tissue type.
in our actual data, the odds ratio of the pooled table is also different from that of the stratified table.
for example, the gene ptrf, a polymerase i and transcript release factor, has a pooled odds ratio of  <dig>  and a cmh odds ratio of  <dig>  calculated from stratified sub-tables.
in this particular case, both odds ratios indicated an inclination toward a higher normal expression and are both statistically significant although at different degree .
for the gene vcan , the pooled odds ratio is  <dig>  and χ <dig> test yields a significant p-value of  <dig> e- <dig>  however, cmh gives an insignificant result for this gene with p-value of  <dig> .
as an extreme case, gbp <dig>  has a pooled odds ratio of  <dig>  and χ <dig> test gives a p-value smaller than  <dig> e- <dig> , whereas with cmh the odds ratio is  <dig> , actually indicating a higher normal counts, although cmh p-value of  <dig>  is insignificant.
this indicates simpson's paradox in action.
careful inspection showed that all cancer counts and most normal counts of gbp <dig> were contributed by the tongue tissue source.
out of a total of  <dig> cancer counts and  <dig> est normal counts, tongue accounts for cancer and normal counts of  <dig> and  <dig>  respectively.
for this gene, the tongue cancer count  <dig> is not influential under a total of  <dig>  cancer counts and  <dig>  counts for the tongue.
thus pooling loses information in this respect and gives a false impression that its cancer expression is much higher when summing all cancer counts from all tissues.
stratifying by tissue type guards against this bias.
heterogeneity of odds ratios
in the strictest application, the use of cochran-mantel-haenszel method requires the odds ratios of the sub-tables be homogeneous.
in our context, it means the ratio of gene expressions between cancer and normal tissues are probably the same among all tissue types under study and any observed variability is most likely due to sampling bias.
also, the calculated odds ratio would be the estimated common odds ratio across the tissue strata.
in our case, however, not all genes had similar ratios under each tissue , and this was of course expected.
in spite of this, we were interested in the overall expression patterns of the genes in cancer conditions.
we were not interested in an estimate of common odds ratio across the strata, which often does not exist.
we were interested in hypothesis testing - to give us leads to the genes that had higher cancer expression in general.
in this regard, the test could be applied  <cit> .
the cmh odds ratio is a weighted average of the odds ratio in each tissue classification and can give us a summary measure  <cit> , which we used to prioritize and followed up with subsequent biological analyses.
in other word, an odds ratio in our data was merely a value that "average up" across all tissue types.
from these ratios we were able to reveal the preferential cancer expressions, since the list covered a number of important known biomarkers, and enrichment of cancer-related genes were supported by knowledge-based genego analyses and previous publications.
lower bound of confidence interval
another distinctive tactic we used is the selection of degs among the statistically significant genes  base on lower bounds of the confidence interval of the odds ratio estimates.
the popular approach to search for degs is to select genes base on p-value first, and then select the subset base on parameter estimators such as odds ratio or fold change values.
the p-value criterion selects the statistically significant ones .
the subsequent criterion is based on prior domain knowledge.
however, among those with statistically significant p-values and similar parameter estimators, the ranges of the estimations can vary widely.
using our dataset as an example, the two genes tuba1b and fam60a both have odd values of  <dig>  .
however, for tuba1b, it is within the 95% confidence that its true odds ratio is between  <dig>  and  <dig> .
yet for fam60a it is between  <dig>  and  <dig> .
based on our background knowledge and for future application, if we must select genes having odds ratios greater than  <dig> , then using odds ratio as cutoff would not serve this purpose since it is quite possible that the real odds ratio  is below  <dig> .
choosing genes based on their confidence intervals would be more precise, but this has not been much appreciated.
multi-cancer biomarkers
the multi-cancer approach compares genes that are overall differentially expressed among multiple cancer types comparing to their respective normal tissue types.
although many biomarker studies focus on gene differentially expressed in a particular tissue type, wu et al.
found  <dig> proteins in the conditioned media of  <dig> cell lines showing negative or weak tissue staining in the human protein atlas, suggesting them to be potential pan-cancer markers  <cit> .
sahin et al., found that claudin- <dig> splice variant  <dig> had the ectopic activations in pancreatic, esophageal, ovarian, and lung tumors while its expression in normal tissue only occurred in differentiated epithelial cells of the gastric mucosa, confirmed by rt-pcr  <cit> .
these studies suggested that relatively multi-cancer genes or multi-cancer splice variants exist.
the three candidates col3a <dig>  chain), dlg <dig>  , and rnf <dig>  are putative secreted or plasma membrane proteins with the potential of developing serum diagnostic reagents.
in reviewing the involvement of these genes with cancers in previous studies, hint for pan-cancer marker was surfaced as the expression of the extracellular matrix protein col3a <dig> gene in brain cancer  <cit>  and angiofibroma  <cit>  was elevated.
while secreted membrane bound rnf <dig> protein gene was known to be up-regulated in colorectal cancer  <cit> .
interestingly, upon the real-time qpcr analysis of three cancer differentially expressed secreted protein gene candidates, col3a <dig>  dlg <dig>  and rnf <dig> identified in this study, higher cancer expression levels of these genes in multiple cancer types were verified.
this does not only indicate the usefulness of our computational approach and filtering procedure but also encourages us to devote further resources for assessing the clinical usages of these three candidates.
pooling of gene expression
earlier in this discussion, we mentioned that naïve pooling of data may introduce bias and at worst may produce simpson's paradox.
we also mentioned that we have tackled this problem with cmh.
nevertheless, two other occasions of pooling actually took place.
we pooled expression from different splicing variants from the same gene to make one gene expression.
we also pooled different libraries of the same tissue into one tissue classification.
in both of these cases, we may encounter expression bias, since different splicing variants and different tissue libraries  might have differences in expression patterns.
this is an unfortunate limitation in this and similar studies, since dbest data consists of many different sources, and given the relative lack of data after the very stringent criteria we have used in our library selection compare to previous studies .
we opted for pooling since we had comparatively limited number of sequences to work with .
nonetheless, future digital expression profiling can be made better with the rna-seq methodology that offers a greater depth of coverage than ests obtained from traditional cdna sequencing.
it gives a much larger sampling size that makes more realistic the differentiation among isoforms and also makes pooling of different libraries of the same tissue less necessary.
as for discovery of pan-cancer genes or isoforms when studying multiple tissue types, similar idea as outlined in this study would be just as applicable.
CONCLUSIONS
we have demonstrated that the use of the cochran-mantel-haenszel statistic in the integrative approaches allowed us to identify potential biomarkers or therapeutic targets via exhaustive search of various est libraries from dbest.
as shown in previous study, splice variant could be useful target of antibody therapy  <cit> .
the method can be easily extended over to searching cancer differential splicing variants had there been enough data.
the issues involved in the analysis, such as the simpson's paradox and the pan-cancer markers, may also be encountered in other multi-class digital analysis.
the three targets confirmed by real-time qpcr, col3a <dig>  dlg <dig>  and rnf <dig>  are worthy of further evaluation for clinical applications.
competing interests
the authors declare that they have no competing interests.
authors' contributions
tw and lc selected and filtered est libraries as well as string analysis.
tw conceived the deployment of the statistical method and implemented and ran the est analysis pipeline.
microarray slides were selected by lc and analyzed by tw.
tw, lc, tc were involved with literature search.
lc selects the genes for expression validation and string analysis.
jw performed the real-time qpcr analysis.
yt is involved with the statistical interpretation of the results.
wl and wn provided direction and guidance.
all authors read and approved the final manuscript.
supplementary material
additional file 1
tissue and library distributions of  <dig> , <dig> ests.
this table shows the number of ests assigned to each tissue type prior to matching to reference sequences.
click here for file
additional file 2
est pipeline raw data.
this is the raw est count from the est pipeline imported into excel.
the columns are the condition type, tissue, and condition-tissue type stratifications.
the rows represent the est counts that are assigned to genes.
click here for file
acknowledgements
this work was supported by grants nsc 99-3112-b-010- <dig>  from the national science council and an intramural grant derived from the aim for the top university grant awarded to national yang ming university from the ministry of education, taiwan, the republic of china.
this article has been published as part of bmc genomics volume  <dig> supplement  <dig>  2012: eleventh international conference on bioinformatics : computational biology.
the full contents of the supplement are available online at http://www.biomedcentral.com/bmcgenomics/supplements/13/s <dig>

 one of the key aspects in the study of cancer is to understand the principles and mechanisms of gene expression variation contributing to cancer genesis and progression. the identification of genes differentially expressed between normal and cancer cells/tissues is not only helpful for designing diagnostic and therapeutic procedures, but also for understanding cancer biology as a whole. in this regard, dna microarrays have been the dominating platform in the high-throughput study of cancer transcriptomes since their emergence in the mid-1990s. however, there are several drawbacks, which include: high background level signals resulting from cross-hybridization; difference in hybridization properties due to different probe sequences; limited dynamic range due to background level and saturation, and difficulty in detecting splicing isoforms and unknown genes. for these reasons, with the advancement of the next generation sequencers, we are seeing high-throughput transcriptome mapping and quantifying method, also known as rna-seq, to begin to supersede microarray in expression profiling. however, rna-seq experiments are relatively demanding in terms of time, cost, and computation equipment. experimental differences between different sequencing platforms may complicate transcriptome analysis with multiple tissue sources. since exploring meta-analysis from traditional digital expression data such as est derived from cdnas is more feasible, this study may serve as a precursor to more complicated experiments. originally primarily aimed for cataloging of transcript repertoire, ests from large-scale cdna sequencing projects such as cancer genome anatomy project, human cancer genome project, and cancer genome project also allow searching for differentially expressed genes in specific tissue types or in whole genomes. several in silico analysis tools such as ncbi unigene cdna xprofiler, cgap digital differential display, and cgap digital gene expression displayer are available online allowing the analysis of publicly available data. while standard statistical methods such as fisher's exact test for finding degs in two-class problems or pearson's correlation are commonly used, there are also specially developed methods for finding degs in the landscape of digital signals for two-library problems or for multiple libraries. the online tools as well as the statistical methods remain useful to this day in est or even rna-seq projects. aside from searching for degs, the searches for gene transcript isoforms specific to particular libraries were also demonstrated and many of these attribute differentially expressed isoforms to human cancers. in spite of the successful applications, these tools or methods are not without limitations. xprofiler reports differential expression in an all-or-none manner where only a list, but not statistical quantification, of candidates is reported. ddd allows quantification using fisher's exact test. however, the nature of the test dictates that comparisons of three or more libraries involve multiple pair-wise comparisons, and thus there are no easy comparisons of library specific genes. dged uses a bayesian approach to find degs, but it is also pair-wise. the reported "odds ratio" is perhaps better described as "relative risk" and may be biased with unequal sampling. another popular and useful bayesian-based method originally developed for est analysis by audic and cleverie is also popular for rna-seq data. it is less conservative than fisher's exact test, but it also does not apply to multi-class problems. the multi-class comparison method established by stekel et al. finds specificity in one condition out of all and is useful in application such as finding degs in multi-tissue libraries. however, in the search for cancer degs, a subsequent analysis of differential expression between cancer and normal libraries of the tissue of interest may not yield fruitful results due to the possible scarcity of est sampling in the particular tissue type. on the other hand, the naïve method of pooling all data into the two-class problem of normal versus cancer when searching for differentially expressed genes or differentially splice variants risks introducing bias. in extreme cases, one may encounter the fallacy of simpson's paradox where genes in reality more active in the normal condition appear to be more so in the cancer condition. we now report on the application of a computational and integrative approach to analyze cancer differentially expressed genes. the statistical method we employed is cochran-mantel-haenszel statistics and to the best of our knowledge has not been applied in this context. instead of pooling all normal and all cancer ests from different tissue types to fit into a two-class problem as by using the by contingency χ test or the fisher's exact test, cmh allows original stratification of libraries in their respective tissue types, yet exhaustively analyzes expression between cancer and normal conditions across all tissue types. the method is an extension to χ test which, in our application, measures the association between cancer and gene expressions, adjusting for the tissue confounding factor. this approach allows one to find genes that are overall differentially expressed in cancer, or multiple-cancer genes, irrespective to a specific tissue type. the method is demonstrated in this paper to exhaustively analyze ests from the dbest database. to the best knowledge of the authors, such an all-inclusive, whole-transcriptome analysis has not been redone in recent years now that more ests than ever are available. our filtering of est libraries was also more rigorous than many previous studies. notably, we excluded the orestes libraries on which a normalization procedure had been applied. libraries from cell line were also excluded owning to their unrepresentativeness of primary cancer cell transcriptomes. our analysis pipeline further focused on enrichment of the degs by cross examination with expression data of a different platform, i.e. the microarray data, and selecting for membrane and secretory associated protein genes since we intend to find therapeutic targets or biomarkers, and conducting string network analysis to show the cancer enriched clusters. with real-time qpcr validation, we have identified three candidates that are inclined to express in cancer across more than one tissue types. we hope such a meta-analytical and multiple-tissue comparison can serve as an exploratory experiment for future multi-library or multi-tissue study of other digital sources such as rna-seq. methods overview our approach was to exploit the entire collection of human est sequences from dbest to obtain transcripts from different type of cells/tissues/organs. the assumption was that the activities of the genes can be represented by their transcripts, and also reflected by the number of representing ests in the ncbi dbest database, given that a large number of mrnas were sequenced. pertinent sequences from different sources were matched to genes and tallied together. through the annotation of each est record, we obtained the tissue type and condition type from which it was derived. with the information, we then had the entire gene transcription profile for all the tissues and conditions. next, cross examining data of other sources including microarray data, secretory and membrane associations as well as analyzing protein associations with string allowed us to narrow down the list of candidate genes. the process is illustrated in figure human gene reference sequence preparation the ncbi reference sequences were downloaded from its ftp site. homo sapiens refseq records were selected and subjected to repeat masking via repeatmasker. human est sequence preparation and library filtering human est data and their cdna library information were downloaded from ncbi dbest database and cgap. program in python language was written to mark for discard the unsuitable libraries when the keywords such as "enrichment", "subtract", "pcr", and "normalized" were found in the descr, unique_protocol, or keywords fields of the library information. an arbitrary cutoff of > was chosen to the highly unrepresentative libraries. to curb from incorrect inclusions or exclusions, we finalized the process with manual curation. libraries made from mixed tissues or cell lines were also discarded. the final libraries from cgap were manually classified into different tissue types and two different conditions, normal and cancer. est to gene assignment the blat alignment tool was used to align ests to refseqs as a mean to assign ests to genes. the criteria of having an identity of 95% or above and the minimum length of nucleotides were set for a match. the refseq match with the highest identity was assigned for the est. if two refseq matches shared exactly the same identity, the program chose the first encountered. est count and summarization the procedure attributes each transcript represented by refseq its expression profile across different tissue and condition types based on est assignment counts. each est has its corresponding tissue type and condition type classification, based on its source clone library. for example, a transcript with an aligned est from a lung cancer clone library is one expression count each in tissue type lung and condition type cancer. this way, after all ests were counted, each transcript has a profile of expression across various libraries and conditions. expressions from different transcript variants of the same gene were pooled to obtain a single gene expression. the raw counts were thus made into transcription profile for each gene for further statistical analysis. statistical evaluation of cancer candidates cochran-mantel-haenszel statistics was applied to evaluate cancer differential expression of each gene. to evaluate each gene, other genes were pooled as "other genes" to create a × × k table consisting of data from tissue-condition cross, where k was the number of tissues ×. a contrived example of × × k table where k is is shown in table gene a is the gene under study while other genes are pooled together as "other genes". only tissue i and tissue ii columns are calculated in cmh. the pooled ones are not part of the analysis. akin to fisher's exact test, the test assumes that "other genes" should consist mostly of genes not differentially expressed between normal or cancer conditions. or, some of them are degs for one condition, but they are at least partly canceled out by degs for the other. in any case, the imbalances of cancer counts to normal counts in the second row is regarded as owning to sample bias and it serves as a metric against which gene a is measured. by continuously isolate values for gene currently under study while pooling all other genes to the second row, an odds ratio and a confidence interval is calculated for each gene. genes with a p-value < and an lower bound of confidence interval of odds ratio ≧ are selected for further analyses. this hypothetical case serves both as an example of how cochran-mantel-haenszel is applied as well as the occurrence of simpson's paradox. gene a is the gene under investigation. expressions from all other genes are pooled into the "other genes" row. bold typeface indicates columns showing higher cancer vs. normal propensities. cmh is applied on the stratified tissue columns. a casual observation involving only the pooled data would suggest gene a as having higher expression in cancer. however, a closer inspection on each of the tissue columns reveals otherwise. the observed difference between cancer and normal of the "other genes" is theoretically mostly due to sampling bias. microarray cross reference human u plus genechip array cel data were downloaded from gene expression omnibus. when computing power allows, the data were processed with affyplm using its three-step procedure of processing background signals with gcrma, normalizing signals with quantile normalization, and summarize probe signals with medium polish. for large experimental datasets that were computationally infeasible for us, we used justrma from the affy package. for experimental dataset without raw cel data, we obtained the pre-processed matrix files via geoquery. regardless of the source of array signal processing, we analyzed the genes for differential expression with limma. differentially expressed gene candidates with p-value < and logfc > were selected and crossed with genes from est profiling with statistical evaluation. for each array, the significant genes were crossed with our est profiling results. the union of these intersecting genes was selected for further evaluation. annotation of secretory proteins to identify our differentially expressed genes with secretory annotation, a list of proteins with secretory annotation originated from the conglomeration of data from uniprot, human plasma proteome organization, and secreted protein database. this list was matched against degs to give them secretory annotation. annotation of membrane proteins membrane protein annotations were gathered from five sources - topdb, locate, pdb_tm, opm, and mpdb - to generate a unique list of membrane proteins. any degs on this list would confer it a membrane annotation. validation of tissue expression profiles of candidate genes tissuescan™ cancer survey panel 96-i qpcr array panel containing the cdnas of normal and cancer tissues each from organs was used to examine the expression profiles of selected cancer differentially expressed gene candidates. real-time qpcr analyses with the taqman® gene expression assay kits and fam- and vic-labeled target genes and hprt internal control primers, respectively, were performed according to the manufacturer's suggested procedure on an applied biosystems prism system. relative specific gene expression was quantified by normalization against the hprt with the Δct method. gene expression changes were quantified as -.  human ests selection and tissue distribution the basic steps of our analysis are illustrated in figure a total of, human est sequences were downloaded from the ncbi. despite the size of the data, not all ests are relevant for our gene expression analysis. after screening the est libraries as described in the methods section above, unsuitable libraries, the preparation of which involved pcr amplification, normalization, subtraction, etc. or originated from cell lines, were discarded. the remaining libraries consisted of, est sequences representing approximately a third of all the downloaded human ests. after blat alignment of the, ests to gene transcripts from refseqs, approximately, ests with at least nucleotides matched to refseqs were detected. an examination of the sources of the matched ests indicated that the representativeness of each tissue is skewed and that the brain is the most represented out of all tissues. among the different tissues, brain ests constituted 26% of all matched ests, uterus ranked second, followed by testis, placenta, pancreas, muscle, liver, kidney and others each below 3%. similarly, condition type representation was also skewed. normal tissue type had, ests combined, and cancer tissue had ests in the ratio of roughly to originally before filtering out those from the cell lines, there were more cancer ests and the ratio of normal ests to cancer ests was roughly to this showed how much more rigorous our filtering was. unfortunately, this also meant we had a much smaller dataset to work with. the unequal distribution of the, matched ests in different tissue types caused some tissue types to be ill-represented. for example, the number of brain est hits dominated over other tissue types. on the other hand, spinal cord had the least count with est hits. the latter had little value for our application. therefore, we only took a tissue type into consideration when its total est hit count was above the cut-off of. considering that the human genome has approximately genes, the cut-off still did not allow "deep" probe into gene expression. nevertheless, the method we employed did not attempt to identify specific gene expression in one particular tissue; therefore, the problem was mitigated. we also categorized ests according to their clone library classification, to either be from normal or from cancer. sometimes a certain tissue-condition type was so under-represented that the information was not trustworthy. for example, adipose had normal hits but only cancer hits, and heart tissue had normal hits but no cancer hits. for these cases, data was kept throughout the analysis. but these data did not make contribution to our analysis. since our est assignments were made to transcripts represented by refseq sequences, when the entire assignment procedure was done, each transcript variant had its expression profile across all tissue-condition types. due to the lack of enough ests data, differentiating between different splicing variants of the same gene was not feasible. we had to pool expression from different splicing variants into a single expression profile representing the gene. analysis of differentially expressed genes due to the small sample size, it was only realistic to evaluate gene expression based on all ests of all tissues. however, tissue type was a confounder. if all counts for each gene were pooled as "normal" or "cancer" regardless of the tissue of origin, the count would be incorrect. to solve both the sample size and the tissue confounder problems, cochran-mantel-haenszel statistical method was employed to identify genes with differential expression as described in the method. we used the arbitrary cut-offs of p-value < and odds ratio ≧ to obtain a primary set of candidates. as a result, a total of cancer differentially expressed gene candidates were selected. the cut-off is chosen based on a good coverage to a list of well-known biomarkers or genes known to associate with cancer. to show that this list of genes was enriched for cancer and thus obtains credibility for our methodology, we looked for cancer related pathways associated with them in genego pathways, which covered signaling and metabolic networks. among the most significantly matched pathways, several are cancer related - pathways number and involve immune response; number and involve cytoskeleton remodeling; number is transition and termination of dna replication; and number and number are adhesion related. in addition, the result of genego disease enrichment analysis indicates our set of genes as neoplasm enriched: seven out of the most associated diseases are related to cancer. the disease ranks the highest is neoplasms, followed by neoplasm by site, and digestive systems neoplasm. this list reveals that our degs covers general neoplasm related functions, and not specific to any particular neoplasm, as digestive, urogenital and breast are all covered. to narrow down this list of biomarkers, we crossed examined the expression profiles of the candidates with the differentially expressed genes in microarray experiments, i.e. two each of ovary and uterus, and one each of pancreas and colon. these tissue types were selected based on the following reasons. we noticed that many of our candidate genes had the most expression in ovary tissue. the other concern was the number of ests. since our candidate genes were derived from est sampling of various tissue types, they were influenced more heavily by tissue types with more est representation due to deeper sampling from them. therefore, the rest of the tissue types were selected based on their representativeness. of the degs, candidates were also found to be differentially expressed genes in our microarray analysis. n: normal, c: cancer gse has two entries since we compared pair-wise between normal vs. benign and normal vs. malignant. since membrane and secretory proteins could be potential therapeutic target or serum biomarkers, the subcellular location of the degs were examined against the secretory and membrane protein lists consolidated from public databases. among these, degs were putative membrane or secretory proteins - had only secretory annotation, had only membrane annotation and had both. literature search and string analysis of the degs to further examine whether the membrane/secretory degs identified in our est database mining had enriched cancer-related genes, we searched the literatures for known associations with cancers. in additions, they were also analyzed with string for interactions, which are based on experimental evidence or prediction, such as conserved genomic neighborhood, gene fusion, co-occurrence across genomes, pathways, protein complex, co-regulation, or other literature sources such as co-mentioning. the network of the string interactions of the degs together with the literature search results were plotted based on the combined string score with cytoscape. approximately proteins formed a big cluster of interacting proteins and a large proportion of the degs had published cancer association with clinical or non-clinical experimental supports. this demonstrates the value of our integration strategy since we had an ample of literature supports. the degs were selected out of their general cancer propensity without necessarily referring to any particular tissue type. however, we can still assess the general tissue distributions shown in figure a gene has a tissue representation if any est from a clone library of the tissue type is matched to it. we can see that some genes are observed across many tissue types. a gene could be observed across a variety of tissue types if it is pan-tissue, and its expression measure is relatively abundant. separately, woolf's test for heterogeneity can also give hints to whether a gene is pan-cancer. those that were found as significant in this test were considered having unequal representation in different genes; although whether they are pan-cancer require further evaluation. three candidates had higher expression in several cancer tissues three cancer differentially expressed secreted protein gene candidates, col3a chain), dlg, and rnf, which had an odds ratio of,, and, respectively, and with limited or no clinical support were selected for real-time qpcr analysis using the taqman® gene expression assay kits. with the hprt as the reference, higher expressions of these genes were noticed in at least some of the cancer tissues. apparently, the average relative expression levels of col3a in breast, liver, thyroid cancer samples were higher than their normal counterparts. the average expression levels of dlg in breast, kidney, liver, lung, and ovarian cancers, and rnf in colon, liver, lung, ovarian, and prostate cancers were also found to be higher than their normal tissues. the expression of col3a in approximately of the liver cancers, dlg in of the liver, lung and ovarian cancers, and rn in of the colon, ovarian and prostate cancers seemed to have higher expressions than the normal tissues. in light of the limited sample size, the three candidates appear to have an overall higher expression in cancer tissues. discussion reported here is an integrative, meta-analytic approach for the discovery of pan-cancer differentially expressed gene candidates. our primary enrichment included a set of degs with cancer associations supported by genego disease and pathway analysis. further integrative evaluations with cancer differentially expressed genes suggested by microarray data narrowed the list down to genes, and among these there were degs likely belonged to either secretory and membrane protein genes. further string protein network analysis and literature reviewing indicated 71% of the degs were highly connected and many of them were associated with cancers in previous publications. simpson's paradox the meta-analytic nature of our study brought us the opportunities as well as challenges to study the digital signatures of various transcriptomes in a new perspective. comparing to experimental methods that focus on a single tissue type or limited tissue types, our approach allows us to find genes inclined to express in cancer in a pan-tissue manner. an important challenge of our approach is to avoid simpson's paradox which can occur in a meta-analysis study. simpson's paradox is where the association between two variables may show a correlation that is reversed in direction from what is observed from stratified subgroups. a contrived example is shown in table in which gene a appears to have a higher cancer expression when pooled, but it is in fact not so under individual stratified sub-tables. this may be somewhat of an extreme case, where directionality of the ratios actually differs between the sub-tables and the pooled table. however, the tissue confounder still introduces bias, large or small, that may throw our judgment off. in this study, we used cmh to analyze the data based on stratified sub-tables to avoid running into this paradox. one could also analyze only one tissue at a time for differential expression, but this means one has a smaller dataset to work with. cmh could avoid this problem since it uses est counts from all tissues instead of analyzing just the normal and cancer propensity under each individual tissue type. in our actual data, the odds ratio of the pooled table is also different from that of the stratified table. for example, the gene ptrf, a polymerase i and transcript release factor, has a pooled odds ratio of and a cmh odds ratio of calculated from stratified sub-tables. in this particular case, both odds ratios indicated an inclination toward a higher normal expression and are both statistically significant although at different degree. for the gene vcan, the pooled odds ratio is and χ test yields a significant p-value of e- however, cmh gives an insignificant result for this gene with p-value of. as an extreme case, gbp has a pooled odds ratio of and χ test gives a p-value smaller than e-, whereas with cmh the odds ratio is, actually indicating a higher normal counts, although cmh p-value of is insignificant. this indicates simpson's paradox in action. careful inspection showed that all cancer counts and most normal counts of gbp were contributed by the tongue tissue source. out of a total of cancer counts and est normal counts, tongue accounts for cancer and normal counts of and respectively. for this gene, the tongue cancer count is not influential under a total of cancer counts and counts for the tongue. thus pooling loses information in this respect and gives a false impression that its cancer expression is much higher when summing all cancer counts from all tissues. stratifying by tissue type guards against this bias. heterogeneity of odds ratios in the strictest application, the use of cochran-mantel-haenszel method requires the odds ratios of the sub-tables be homogeneous. in our context, it means the ratio of gene expressions between cancer and normal tissues are probably the same among all tissue types under study and any observed variability is most likely due to sampling bias. also, the calculated odds ratio would be the estimated common odds ratio across the tissue strata. in our case, however, not all genes had similar ratios under each tissue, and this was of course expected. in spite of this, we were interested in the overall expression patterns of the genes in cancer conditions. we were not interested in an estimate of common odds ratio across the strata, which often does not exist. we were interested in hypothesis testing - to give us leads to the genes that had higher cancer expression in general. in this regard, the test could be applied. the cmh odds ratio is a weighted average of the odds ratio in each tissue classification and can give us a summary measure, which we used to prioritize and followed up with subsequent biological analyses. in other word, an odds ratio in our data was merely a value that "average up" across all tissue types. from these ratios we were able to reveal the preferential cancer expressions, since the list covered a number of important known biomarkers, and enrichment of cancer-related genes were supported by knowledge-based genego analyses and previous publications. lower bound of confidence interval another distinctive tactic we used is the selection of degs among the statistically significant genes base on lower bounds of the confidence interval of the odds ratio estimates. the popular approach to search for degs is to select genes base on p-value first, and then select the subset base on parameter estimators such as odds ratio or fold change values. the p-value criterion selects the statistically significant ones. the subsequent criterion is based on prior domain knowledge. however, among those with statistically significant p-values and similar parameter estimators, the ranges of the estimations can vary widely. using our dataset as an example, the two genes tuba1b and fam60a both have odd values of. however, for tuba1b, it is within the 95% confidence that its true odds ratio is between and. yet for fam60a it is between and. based on our background knowledge and for future application, if we must select genes having odds ratios greater than, then using odds ratio as cutoff would not serve this purpose since it is quite possible that the real odds ratio is below. choosing genes based on their confidence intervals would be more precise, but this has not been much appreciated. multi-cancer biomarkers the multi-cancer approach compares genes that are overall differentially expressed among multiple cancer types comparing to their respective normal tissue types. although many biomarker studies focus on gene differentially expressed in a particular tissue type, wu et al. found proteins in the conditioned media of cell lines showing negative or weak tissue staining in the human protein atlas, suggesting them to be potential pan-cancer markers. sahin et al., found that claudin- splice variant had the ectopic activations in pancreatic, esophageal, ovarian, and lung tumors while its expression in normal tissue only occurred in differentiated epithelial cells of the gastric mucosa, confirmed by rt-pcr. these studies suggested that relatively multi-cancer genes or multi-cancer splice variants exist. the three candidates col3a chain), dlg, and rnf are putative secreted or plasma membrane proteins with the potential of developing serum diagnostic reagents. in reviewing the involvement of these genes with cancers in previous studies, hint for pan-cancer marker was surfaced as the expression of the extracellular matrix protein col3a gene in brain cancer and angiofibroma was elevated. while secreted membrane bound rnf protein gene was known to be up-regulated in colorectal cancer. interestingly, upon the real-time qpcr analysis of three cancer differentially expressed secreted protein gene candidates, col3a dlg and rnf identified in this study, higher cancer expression levels of these genes in multiple cancer types were verified. this does not only indicate the usefulness of our computational approach and filtering procedure but also encourages us to devote further resources for assessing the clinical usages of these three candidates. pooling of gene expression earlier in this discussion, we mentioned that naïve pooling of data may introduce bias and at worst may produce simpson's paradox. we also mentioned that we have tackled this problem with cmh. nevertheless, two other occasions of pooling actually took place. we pooled expression from different splicing variants from the same gene to make one gene expression. we also pooled different libraries of the same tissue into one tissue classification. in both of these cases, we may encounter expression bias, since different splicing variants and different tissue libraries might have differences in expression patterns. this is an unfortunate limitation in this and similar studies, since dbest data consists of many different sources, and given the relative lack of data after the very stringent criteria we have used in our library selection compare to previous studies. we opted for pooling since we had comparatively limited number of sequences to work with. nonetheless, future digital expression profiling can be made better with the rna-seq methodology that offers a greater depth of coverage than ests obtained from traditional cdna sequencing. it gives a much larger sampling size that makes more realistic the differentiation among isoforms and also makes pooling of different libraries of the same tissue less necessary. as for discovery of pan-cancer genes or isoforms when studying multiple tissue types, similar idea as outlined in this study would be just as applicable.  we have demonstrated that the use of the cochran-mantel-haenszel statistic in the integrative approaches allowed us to identify potential biomarkers or therapeutic targets via exhaustive search of various est libraries from dbest. as shown in previous study, splice variant could be useful target of antibody therapy. the method can be easily extended over to searching cancer differential splicing variants had there been enough data. the issues involved in the analysis, such as the simpson's paradox and the pan-cancer markers, may also be encountered in other multi-class digital analysis. the three targets confirmed by real-time qpcr, col3a dlg and rnf are worthy of further evaluation for clinical applications. competing interests the authors declare that they have no competing interests. authors' contributions tw and lc selected and filtered est libraries as well as string analysis. tw conceived the deployment of the statistical method and implemented and ran the est analysis pipeline. microarray slides were selected by lc and analyzed by tw. tw, lc, tc were involved with literature search. lc selects the genes for expression validation and string analysis. jw performed the real-time qpcr analysis. yt is involved with the statistical interpretation of the results. wl and wn provided direction and guidance. all authors read and approved the final manuscript. supplementary material additional file 1 tissue and library distributions of, ests. this table shows the number of ests assigned to each tissue type prior to matching to reference sequences. click here for file additional file 2 est pipeline raw data. this is the raw est count from the est pipeline imported into excel. the columns are the condition type, tissue, and condition-tissue type stratifications. the rows represent the est counts that are assigned to genes. click here for file acknowledgements this work was supported by grants nsc 99-3112-b-010- from the national science council and an intramural grant derived from the aim for the top university grant awarded to national yang ming university from the ministry of education, taiwan, the republic of china. this article has been published as part of bmc genomics volume supplement 2012: eleventh international conference on bioinformatics: computational biology. the full contents of the supplement are available online at http://www.biomedcentral.com/bmcgenomics/supplements/13/s 
the three candidates identified in this study, col3a dlg and rnf are valuable targets for further evaluation with a larger sample size of normal and cancer tissue or serum samples. genego analysis results indicated this set as neoplasm enriched. finally, the differential expression patterns of three selected cancer biomarker candidates were confirmed by real-time qpcr analysis. after examined the candidates in protein interaction network, public tissue expression databases, and literatures, we selected three genes for further evaluation by real-time qpcr with eight major normal and cancer tissues. searching digitized transcriptome using cmh enabled us to identify multi-cancer differentially expressed gene candidates. although the availability of public databases make possible the comprehensive mining of degs among the ests from multiple tissue types, existing studies usually employed statistics suitable only for two categories. on the other hand, simple pooling of cancer and normal genes from multiple tissue types runs the risk of simpson's paradox. our methodology demonstrated simultaneously analysis for cancer biomarkers of multiple tissue types with the est data. the higher-than-normal tissue expression of col3a dlg and rnf in some of the cancer tissues is in agreement with our in silico predictions. seven hundred and twenty three primary deg candidates were selected from a curated est database with the application of cochran-mantel-haenszel statistic. here we presented a different approach which searched for multi-cancer deg candidates by analyzing all pertinent ests in all categories and narrowing down the cancer biomarker candidates via integrative analysis with microarray data and selection of secretory and membrane protein genes as well as incorporation of network analysis. cross-examination with microarray data further narrowed the list down to genes, among which had membrane or secretory annotations. 3- october asia pacific bioinformatics network eleventh international conference on bioinformatics bangkok, thailand multi-class test has been developed to enable the finding of tissue specific genes, but subsequent search for cancer genes involves separate two-category test only on the ests of the tissue of interest.
 researches have been conducted for the identification of differentially expressed genes by generating and mining of cdna expressed sequence tags for more than a decade. although the availability of public databases make possible the comprehensive mining of degs among the ests from multiple tissue types, existing studies usually employed statistics suitable only for two categories. multi-class test has been developed to enable the finding of tissue specific genes, but subsequent search for cancer genes involves separate two-category test only on the ests of the tissue of interest. this constricts the amount of data used. on the other hand, simple pooling of cancer and normal genes from multiple tissue types runs the risk of simpson's paradox. here we presented a different approach which searched for multi-cancer deg candidates by analyzing all pertinent ests in all categories and narrowing down the cancer biomarker candidates via integrative analysis with microarray data and selection of secretory and membrane protein genes as well as incorporation of network analysis. finally, the differential expression patterns of three selected cancer biomarker candidates were confirmed by real-time qpcr analysis.  seven hundred and twenty three primary deg candidates were selected from a curated est database with the application of cochran-mantel-haenszel statistic. genego analysis results indicated this set as neoplasm enriched. cross-examination with microarray data further narrowed the list down to genes, among which had membrane or secretory annotations. after examined the candidates in protein interaction network, public tissue expression databases, and literatures, we selected three genes for further evaluation by real-time qpcr with eight major normal and cancer tissues. the higher-than-normal tissue expression of col3a dlg and rnf in some of the cancer tissues is in agreement with our in silico predictions.  searching digitized transcriptome using cmh enabled us to identify multi-cancer differentially expressed gene candidates. our methodology demonstrated simultaneously analysis for cancer biomarkers of multiple tissue types with the est data. with the revived interest in digitizing the transcriptomes by ngs, cancer biomarkers could be more precisely detected from the ests. the three candidates identified in this study, col3a dlg and rnf are valuable targets for further evaluation with a larger sample size of normal and cancer tissue or serum samples. 3- october asia pacific bioinformatics network eleventh international conference on bioinformatics bangkok, thailand 