seek assembles its human gene expression compendium by obtaining data sets from ncbi s gene expression omnibus ( geo ) database and the cancer genome atlas ( tcga ) .
the compendium consists of data sets from 41 platforms including 32 platforms from affymetrix , agilent , and illumina , and 9 rna sequencing platforms ( supplementary data 1 ) .
these platforms were chosen based on the number of available data sets and the availability of raw data to perform consistent processing for each platform .
the data sets were processed consistently by applying platform - specific procedures on their raw measurements ( supplementary note 5 and supplementary data 4 ) to remove the systematic differences among data sets .
the normalized data sets containing gene - level expression values can be accessed through the seek website .
the next step of data processing is calculating the pearson correlations rd(x , y ) between all pairs of genes x and y in each data set d. as correlation values arising from different genome - wide distributions are not directly comparable across data sets , a fisher transform procedure is applied to convert the distribution of correlations to a normal - like distribution : 
 fd(x , y)=12ln1+rd(x , y)1rd(x , y ) 
 where fd(x , y ) is the fisher - transformed score .
then , the data are translated to z - scores for standardization : 
 zd(x , y)=1std(fd)[fd(x , y)avg(fd ) ] 
 where avg(fd ) is the average of fd for all ( x , y ) pairs , and std(fd ) is the standard deviation of fd . the normalization procedure has been used in previous studies , and has been found successful in transforming most correlation distributions that are generated from different platforms and technologies into a comparable normal distribution with mean 0 and variance 1 ( supplementary fig .
note that seek also works well with other correlation measures , such as spearman and bicor ( supplementary fig . 9 ) .
we found that the normalized pearson correlation provides performance better or comparable to that of spearman and bicor in the search setting , likely because the normalization procedure and the seek algorithm itself reduce the effects of outliers in search performance ( supplementary fig . 9 ) .
the search algorithm takes two inputs 1 ) a set of query genes q = { q1 ,  , qx } , and 2 ) the set of correlation z - scores containing the query : zd(g , q ) , for each data set d in the data compendium d , for all genes q in q and for all genes g in the genome g. the outputs of the algorithm are a prioritized list of data sets and co - expressed genes relevant to q. the search algorithm consists of four steps .
the first step is to load pre - computed z - scores of pearson correlations ( in the normalization step above ) containing the query across d. the second step is to perform hubbiness correction on each data set d. the correction procedure is motivated by the observation that hubby or well - connected genes in the co - expression network represent global , well - co - expressed processes , and can contaminate the search results regardless of query composition due to the effect of unbalanced gene connectivity in a scale - free co - expression network , and can lead to non - specific results in search or clustering approaches . to avoid the bias created by hubby genes that are not related to the user s query or pathway of interest , our method corrects each gene g s correlation to q in each data set d : 
 ( eq . 1)zd(g , q)=zd(g , q)1|g|xgzd(g , x ) 
 where z is the hubbiness - corrected z - score . by subtracting g s average correlation from the correlation of ( g , q ) , we expect the resulting score to emphasize g s co - expression specifically with the query rather than its general connectivity .
the control of co - expression hubbiness enables the detection of specific biological signals in the data that would otherwise be swamped by broad co - expression patterns of the most well - connected genes .
the goal is to rank data sets based on each data set s relevance to the query .
the result will be the first output of the search system and will also be used to compute the final gene - score vector for the last step .
the main idea is to upweight data sets where a subset of the query genes can retrieve the remaining query genes well based on normalized , hubbiness - corrected co - expression in that data set .
thus , it is analogous in spirit to the cross - validation procedures commonly used in machine learning , where a subset of the standard ( in this case query )  hides  from the system to assess how well the method can predict these hidden genes . to describe the weighting method
the data set d is implicit in each formula below and omitted for brevity , thus z(g , q ) is the corrected z - score for g to a query gene q in q in data set d. let rq = ( g , g , g ,  , g ) be the sequence of genes at rank 1 , 2 , 3 ,  ,
that is , rq satisfies : z(g , q )  z(g , q )  z(g , q )  let r(t , rq ) be the rank of gene t in the ranking rq minus 1 ( for example , r(g , rq ) = 0 ) , and let p < 1be a rate parameter , which we set at 0.99 based on empirical analysis ( supplementary fig .
2)w=1|q|qq[(1p)tqqpr(t , rq ) ] 
  the weighting formula performs cross - validations on q in the set q. the goal is to detect which query genes q can best retrieve the remainder query q  q ; such q s have a high contribution to w. we shorten r(t , rq ) in eq . 2 as r(t ) .
the exact form of this expression for weight ( i.e. sum of p ) is inspired by rank - biased precision , and is adapted to our setting to robustly measure the effectiveness of the data set in retrieving q 
q. here , p < 1 is the rate parameter in rank - biased precision , and is the parameter of geometric distribution , since r(t ) assumes discrete values . when it is employed , p upweights ranks for genes t in the set q
 q that are high in the rank list ( i.e. , r(t ) is small ) , which intuitively emphasizes those genes in the query that are highly co - expressed to each other .
the measure has the desired property of upweighting pairs of query genes that are well correlated while not allowing the correlations between the uninformative , non - coherent part of the query to affect the weight of the data set because the query genes may only be partially co - expressed in a given data set .
compared to previous methods , our method gains robustness to heterogeneous query signals , because the reward on the highly coherent query genes far outweighs the damaging effect of a few non - coherent query genes , which are poorly ranked to other query genes , have high r(t ) , and have scores p tending to zero .
the last step of the algorithm calculates the final integrated gene scores to generate a master ranking of co - expressed genes that is the second output of the system ( in addition to data set relevance weighting ) .
we obtain the gene - to - query score matrix mg , d , where the entry mg , d is the average corrected z - score of gene g to the query in data set d : 
 mg , d=1|q|qqzd(g , q ) 
 with the data set weight vector from the previous step w = [ w1 , w2 ,  ] , a simple formulation of the final gene - score vector f is given by : 
 f = mg , dwt,=1/ddwd 
  although previous research had some success with this formulation , our findings show that it works well only in the presence of complete gene information with no missing genes in mg , d . when there are heterogeneous sources of data in the compendium ( e.g. different microarray and rnaseq platforms ) , the confounding factor of missing genes and partial gene rankings must be accounted for .
our approach is to modify the procedure above by employing threshold parameters to exclude a data set from weighting if it does not contain enough query genes , and exclude a gene from the final ranking if it is not assayed by a sufficient number of data sets in the compendium ( supplementary note 6 ) .
the pseudocode for the entire seek search algorithm can be found in supplementary note 6 .
1 , 2 ) and data set quality , including automatically down weighting data sets with substantial batch effects ( supplementary note 4 and supplementary fig .
, the search algorithm performs the same steps above except that in the data set weighting step the algorithm assigns equal weight to all data sets .
thus , for single - gene queries , the search system will treat each data set equally and retrieve genes that are generally correlated with the query in the hubbiness - corrected space .
if users wish to perform their single - gene searches in a tissue - specific or disease - specific manner , they can manually define a category of data sets using the extensive  refine search 
interface on the seek website , which will restrict d in the search system input .
we estimate a p value for each retrieved gene by comparing the integrated score of each gene with scores from a pool of 10,000 randomly generated queries with diverse query sizes varying from 1 to 100 genes .
the random pool allows seek to estimate the significance of gene - score as well as evaluate the specificity of that gene to the query genes ( as opposed to random queries ) . for a given gene g and its final co - expression score sq(g ) generated from the user s query q
, the p value of g is estimated as the number of random queries r in which sr(g ) >
the seek algorithm is implemented in c++ and has been integrated into the open - source c++ sleipnir library , enabling other computational users to use and expand seek without website tie - in .
the backend employs the efficient data structures from the sleipnir library to facilitate the process of handling large query sets of over 100 genes without memory overflow .
seek s jobs are parallelized to make full use of the multiprocessor resources and their processing power .
the seek web server is constructed with some of the latest web technologies including jquery and qtip2 libraries .
dynamic pages are generated with java servlets running behind the apache tomcat server on a red hat centos linux operating system .
in addition , ajax technology is deployed to send and retrieve data from the server asynchronously such that users can receive instant feedback on their gene enrichment analysis , expression zoom - in function , and data set selection module without having to leave or refresh the page .
seek categorizes data sets into tissue and disease groups by mining the description , title , and sample level characteristic fields in data sets metadata .
the text mining procedure utilizes the umls metathesaurus and brenda controlled vocabularies to extract predefined concept names that are present in the individual fields . to ensure that tissue groups are accurate , we manually reviewed annotations to the frequently appearing terms generated by text mining .
meta data set groups , such as cancer and non - cancer groups and the multi - tissue profiling group ( supplementary data 5 ) , to provide users with the ability to limit their search to such groups under the  refine search  feature of the website .
we conducted a comprehensive evaluation of seek in comparison with existing algorithms gene recommender , mem , and meta - data set correlation search ( supplementary note 1 ) .
we tested each system s ability to retrieve genes from the same biological process given some chosen genes from the process as queries .
for the evaluation , we partitioned the genes in each of the 995 go biological process terms ( supplementary data 2 ) into a query building set and a testing set .
the query building set consists of a random sample of 25 genes from each term if the term has more than 40 genes , or else it is made of half of the number of genes in the term .
queries were formed by repeatedly sampling genes from the set , so that each query size has 10 different queries of that size represented , and we iteratively generated queries for sizes 2 , 3 , 4 ,  up to q genes , where q=0.8|query building set| .
the testing set consists of the remaining genes in the term ( after subtracting the query building set ) , and is used for evaluating the queries retrieval results .
a precision recall ( pr ) curve is computed on a per - query basis , averaged over all queries of a term , and finally averaged over all evaluated terms to derive an overall system performance plot for each method .
fold precision over random is calculated at an indicated recall ( 10% ) , and uses a random ranking of genes where genes rank positions are shuffled . by selecting genes randomly from each process in building the queries , we mimic the situation in which the query genes are functionally related but not well co - expressed . by keeping the two sets ( query building and testing ) separate in the evaluation
, we can reduce the performance variation between the queries of the same size within a process . for building gold standard go gene - set for evaluation , we used gene annotations with experimental evidence codes ( imp , igi , ipi , ida , iep , exp ) as well as tas ( traceable author statements ) and nas ( non - traceable author statement ) . to select the
go slim set ( supplementary data 3 ) used for studying the effect of compendium size , we carefully examined the title and description of the go terms in the context of the go hierarchy and arrived at a non - redundant subset of go terms that are both specific enough to be informative , and diverse enough to represent the hierarchy , similar to our approach in . to evaluate seek s performance as a function of the query size ,
we pooled together previously built biological process queries from 995 processes , then binned them by query size ( 220 genes ) .
we examined 3 categories of biological processes based on process size : 2040 genes , 40100 genes , and 100300 genes .
performance refers to the fold precision over random at 10% recall ( fold pr10% ) in using each query to retrieve remaining genes from its corresponding process . to evaluate the search system s robustness to noisy query genes , we selected over 1,800 5-gene and 10-gene queries from 90 kegg pathways with 50100 genes per pathway .
we established a  no noise  case , where each query was purely made of genes belonging to the same kegg pathway , and a noisy case , where 1- , 2- , 4- random genes were added to each query .
the fraction ( fold pr@10% of each noisy query ) / ( fold pr@10% of the corresponding  no noise  query ) was calculated , where fold pr@10% refers to the performance of retrieving kegg pathway genes using the queries .
seek assembles its human gene expression compendium by obtaining data sets from ncbi s gene expression omnibus ( geo ) database and the cancer genome atlas ( tcga ) .
the compendium consists of data sets from 41 platforms including 32 platforms from affymetrix , agilent , and illumina , and 9 rna sequencing platforms ( supplementary data 1 ) .
these platforms were chosen based on the number of available data sets and the availability of raw data to perform consistent processing for each platform .
the data sets were processed consistently by applying platform - specific procedures on their raw measurements ( supplementary note 5 and supplementary data 4 ) to remove the systematic differences among data sets .
the normalized data sets containing gene - level expression values can be accessed through the seek website .
the next step of data processing is calculating the pearson correlations rd(x , y ) between all pairs of genes x and y in each data set d. as correlation values arising from different genome - wide distributions are not directly comparable across data sets , a fisher transform procedure is applied to convert the distribution of correlations to a normal - like distribution : 
 fd(x , y)=12ln1+rd(x , y)1rd(x , y ) 
 where fd(x , y ) is the fisher - transformed score .
then , the data are translated to z - scores for standardization : 
 zd(x , y)=1std(fd)[fd(x , y)avg(fd ) ] 
 where avg(fd ) is the average of fd for all ( x , y ) pairs , and std(fd ) is the standard deviation of fd . the normalization procedure has been used in previous studies , and has been found successful in transforming most correlation distributions that are generated from different platforms and technologies into a comparable normal distribution with mean 0 and variance 1 ( supplementary fig .
note that seek also works well with other correlation measures , such as spearman and bicor ( supplementary fig . 9 ) .
we found that the normalized pearson correlation provides performance better or comparable to that of spearman and bicor in the search setting , likely because the normalization procedure and the seek algorithm itself reduce the effects of outliers in search performance ( supplementary fig . 9 ) .
the search algorithm takes two inputs 1 ) a set of query genes q = { q1 ,  , qx } , and 2 ) the set of correlation z - scores containing the query : zd(g , q ) , for each data set d in the data compendium d , for all genes q in q and for all genes g in the genome g. the outputs of the algorithm are a prioritized list of data sets and co - expressed genes relevant to q. the search algorithm consists of four steps .
the first step is to load pre - computed z - scores of pearson correlations ( in the normalization step above ) containing the query across d. the second step is to perform hubbiness correction on each data set d. the correction procedure is motivated by the observation that hubby or well - connected genes in the co - expression network represent global , well - co - expressed processes , and can contaminate the search results regardless of query composition due to the effect of unbalanced gene connectivity in a scale - free co - expression network , and can lead to non - specific results in search or clustering approaches . to avoid the bias created by hubby genes that are not related to the user s query or pathway of interest , our method corrects each gene g s correlation to q in each data set d : 
 ( eq .
1)zd(g , q)=zd(g , q)1|g|xgzd(g , x ) 
 where z is the hubbiness - corrected z - score . by subtracting g s average correlation from the correlation of ( g , q ) , we expect the resulting score to emphasize g s co - expression specifically with the query rather than its general connectivity .
the control of co - expression hubbiness enables the detection of specific biological signals in the data that would otherwise be swamped by broad co - expression patterns of the most well - connected genes .
the goal is to rank data sets based on each data set s relevance to the query
. the result will be the first output of the search system and will also be used to compute the final gene - score vector for the last step .
the main idea is to upweight data sets where a subset of the query genes can retrieve the remaining query genes well based on normalized , hubbiness - corrected co - expression in that data set .
thus , it is analogous in spirit to the cross - validation procedures commonly used in machine learning , where a subset of the standard ( in this case query )  hides  from the system to assess how well the method can predict these hidden genes . to describe the weighting method
the data set d is implicit in each formula below and omitted for brevity , thus z(g , q ) is the corrected z - score for g to a query gene q in q in data set d. let rq = ( g , g , g ,  , g ) be the sequence of genes at rank 1 , 2 , 3 ,  ,
that is , rq satisfies : z(g , q )  z(g , q )  z(g , q )  let r(t , rq ) be the rank of gene t in the ranking rq minus 1 ( for example , r(g , rq ) = 0 ) , and let p < 1be a rate parameter , which we set at 0.99 based on empirical analysis ( supplementary fig . 10 ) .
then the weight of the data set is 
 ( eq . 2)w=1|q|qq[(1p)tqqpr(t , rq ) ] 
  the weighting formula performs cross - validations on q in the set q. the goal is to detect which query genes q can best retrieve the remainder query q  q ; such q s have a high contribution to w. we shorten r(t , rq ) in eq . 2 as r(t ) .
the exact form of this expression for weight ( i.e. sum of p ) is inspired by rank - biased precision , and is adapted to our setting to robustly measure the effectiveness of the data set in retrieving q 
q. here , p < 1 is the rate parameter in rank - biased precision , and is the parameter of geometric distribution , since r(t ) assumes discrete values . when it is employed , p upweights ranks for genes t in the set q 
q that are high in the rank list ( i.e. , r(t ) is small ) , which intuitively emphasizes those genes in the query that are highly co - expressed to each other .
the measure has the desired property of upweighting pairs of query genes that are well correlated while not allowing the correlations between the uninformative , non - coherent part of the query to affect the weight of the data set because the query genes may only be partially co - expressed in a given data set .
compared to previous methods , our method gains robustness to heterogeneous query signals , because the reward on the highly coherent query genes far outweighs the damaging effect of a few non - coherent query genes , which are poorly ranked to other query genes , have high r(t ) , and have scores p tending to zero .
the last step of the algorithm calculates the final integrated gene scores to generate a master ranking of co - expressed genes that is the second output of the system ( in addition to data set relevance weighting ) .
we obtain the gene - to - query score matrix mg , d , where the entry mg , d is the average corrected z - score of gene g to the query in data set d : 
 mg , d=1|q|qqzd(g , q ) 
 with the data set weight vector from the previous step w = [ w1 , w2 ,  ] , a simple formulation of the final gene - score vector f is given by : 
 f = mg , dwt,=1/ddwd 
  although previous research had some success with this formulation , our findings show that it works well only in the presence of complete gene information with no missing genes in mg , d .
when there are heterogeneous sources of data in the compendium ( e.g. different microarray and rnaseq platforms ) , the confounding factor of missing genes and partial gene rankings must be accounted for .
our approach is to modify the procedure above by employing threshold parameters to exclude a data set from weighting if it does not contain enough query genes , and exclude a gene from the final ranking if it is not assayed by a sufficient number of data sets in the compendium ( supplementary note 6 ) .
the pseudocode for the entire seek search algorithm can be found in supplementary note 6 .
1 , 2 ) and data set quality , including automatically down weighting data sets with substantial batch effects ( supplementary note 4 and supplementary fig .
the search algorithm performs the same steps above except that in the data set weighting step the algorithm assigns equal weight to all data sets .
thus , for single - gene queries , the search system will treat each data set equally and retrieve genes that are generally correlated with the query in the hubbiness - corrected space .
if users wish to perform their single - gene searches in a tissue - specific or disease - specific manner , they can manually define a category of data sets using the extensive  refine search 
interface on the seek website , which will restrict d in the search system input .
we estimate a p value for each retrieved gene by comparing the integrated score of each gene with scores from a pool of 10,000 randomly generated queries with diverse query sizes varying from 1 to 100 genes .
the random pool allows seek to estimate the significance of gene - score as well as evaluate the specificity of that gene to the query genes ( as opposed to random queries ) . for a given gene g and its final co - expression score sq(g ) generated from the user s query q
, the p value of g is estimated as the number of random queries r in which sr(g ) > sq(g ) divided by the random pool size .
the seek algorithm is implemented in c++ and has been integrated into the open - source c++ sleipnir library , enabling other computational users to use and expand seek without website tie - in .
the backend employs the efficient data structures from the sleipnir library to facilitate the process of handling large query sets of over 100 genes without memory overflow .
seek s jobs are parallelized to make full use of the multiprocessor resources and their processing power .
the seek web server is constructed with some of the latest web technologies including jquery and qtip2 libraries .
dynamic pages are generated with java servlets running behind the apache tomcat server on a red hat centos linux operating system .
in addition , ajax technology is deployed to send and retrieve data from the server asynchronously such that users can receive instant feedback on their gene enrichment analysis , expression zoom - in function , and data set selection module without having to leave or refresh the page .
seek categorizes data sets into tissue and disease groups by mining the description , title , and sample level characteristic fields in data sets metadata .
the text mining procedure utilizes the umls metathesaurus and brenda controlled vocabularies to extract predefined concept names that are present in the individual fields . to ensure that tissue groups are accurate , we manually reviewed annotations to the frequently appearing terms generated by text mining .
meta data set groups , such as cancer and non - cancer groups and the multi - tissue profiling group ( supplementary data 5 ) , to provide users with the ability to limit their search to such groups under the  refine search  feature of the website .
we conducted a comprehensive evaluation of seek in comparison with existing algorithms gene recommender , mem , and meta - data set correlation search ( supplementary note 1 ) .
we tested each system s ability to retrieve genes from the same biological process given some chosen genes from the process as queries .
for the evaluation , we partitioned the genes in each of the 995 go biological process terms ( supplementary data 2 ) into a query building set and a testing set .
the query building set consists of a random sample of 25 genes from each term if the term has more than 40 genes , or else it is made of half of the number of genes in the term .
queries were formed by repeatedly sampling genes from the set , so that each query size has 10 different queries of that size represented , and we iteratively generated queries for sizes 2 , 3 , 4 ,  up to q genes , where q=0.8|query building set| .
the testing set consists of the remaining genes in the term ( after subtracting the query building set ) , and is used for evaluating the queries retrieval results . a precision recall
( pr ) curve is computed on a per - query basis , averaged over all queries of a term , and finally averaged over all evaluated terms to derive an overall system performance plot for each method .
fold precision over random is calculated at an indicated recall ( 10% ) , and uses a random ranking of genes where genes rank positions are shuffled . by selecting genes randomly from each process in building the queries , we mimic the situation in which the query genes are functionally related but not well co - expressed . by keeping the two sets ( query building and testing ) separate in the evaluation
, we can reduce the performance variation between the queries of the same size within a process . for building gold standard go gene - set for evaluation
, we used gene annotations with experimental evidence codes ( imp , igi , ipi , ida , iep , exp ) as well as tas ( traceable author statements ) and nas ( non - traceable author statement ) . to select the
go slim set ( supplementary data 3 ) used for studying the effect of compendium size , we carefully examined the title and description of the go terms in the context of the go hierarchy and arrived at a non - redundant subset of go terms that are both specific enough to be informative , and diverse enough to represent the hierarchy , similar to our approach in . to evaluate seek s performance as a function of the query size , we pooled together previously built biological process queries from 995 processes , then binned them by query size ( 220 genes ) .
we examined 3 categories of biological processes based on process size : 2040 genes , 40100 genes , and 100300 genes .
performance refers to the fold precision over random at 10% recall ( fold pr10% ) in using each query to retrieve remaining genes from its corresponding process . to evaluate the search system s robustness to noisy query genes , we selected over 1,800 5-gene and 10-gene queries from 90 kegg pathways with 50100 genes per pathway .
we established a  no noise  case , where each query was purely made of genes belonging to the same kegg pathway , and a noisy case , where 1- , 2- , 4- random genes were added to each query .
the fraction ( fold pr@10% of each noisy query ) / ( fold pr@10% of the corresponding  no noise  query ) was calculated , where fold pr@10% refers to the performance of retrieving kegg pathway genes using the queries .
 we present seek ( http://seek.princeton.edu ) , a query - based search engine across very large transcriptomic data collections , including thousands of human data sets from almost 50 microarray and next - generation sequencing platforms . 
 seek uses a novel query - level cross - validation - based algorithm to automatically prioritize data sets relevant to the query and a robust search approach to identify query - coregulated genes , pathways , and processes . 
 seek provides cross - platform handling , multi - gene query search , iterative metadata - based search refinement , and extensive visualization - based analysis options . 
Online methods
Data preparation and correlation normalization
Search algorithm
Estimating the significance of gene scores
Algorithm and interface implementations
Metadata processing
Large-scale functional evaluation setup
Supplementary Material