Introduction
============

Gene coexpression analysis, which is a guilt-by-association approach based on gene expression profiles, can uncover functionally related gene pairs.  Of particular importance, coexpression information based on a large amount of publicly available transcriptome data is not highly affected by a specific experimental condition and provides a fundamental view of gene functional networks.  A number of databases have been developed to provide condition-independent coexpression information for various applications, from gene prioritization to the delineation of global relationships among multiple network modules ([@pcx191-B1], [@pcx191-B36], [@pcx191-B32]).  In addition, comparisons of gene networks among multiple species can be used to identify conserved and specific gene modules in an evolutionary context.  Conserved gene relationships suggest a core function for a cell system, whereas species-specific relationships are linked to the differentiation of species ([@pcx191-B35], [@pcx191-B31], [@pcx191-B25], [@pcx191-B19], [@pcx191-B30], [@pcx191-B33]). 

Although the importance of coexpression data is widely recognized, related methodologies have various limitations.  For example, the determination of coexpression consists of multiple steps, from the preparation of gene expression data to the calculation of a gene-to-gene matrix ([@pcx191-B38]).  To optimize this procedure, we should find the optimal combination of algorithms for each step.  In the case of the detection of differentially expressed genes from microarray or RNA sequencing (RNAseq) data, which is related to the task of coexpression calculation, we can use a benchmark dataset, such as MicroArray Quality Control datasets ([@pcx191-B7], [@pcx191-B13]).  However, determining a gold standard for gene coexpression is quite difficult because gene coexpression is a summarization metric of gene expression data and thus it depends on the composition of individual experimental samples ([@pcx191-B36]). 

We have been continuously developing a coexpression database for plants, ATTED-II (<http://atted.jp>).  In previous reports, we have described associations between gene coexpression and *cis*-elements ([@pcx191-B26]), the development of a coexpression measure with sample weights and the mutual rank (MR) index ([@pcx191-B22]), the analysis of condition-specific coexpression ([@pcx191-B27]) and the evaluation of coexpression ([@pcx191-B28]).  In the most recent report associated with ATTED-II, we demonstrated the detection of lineage-specific gene coexpression using eight coexpression platforms for four species ([@pcx191-B3]).  However, we could not apply statistical tests to the identification of lineage-specific coexpression for two reasons.  First, the correspondence of sample conditions among the four species was not clear.  To account for potential differences in the compositions of samples, we focused on lineage-specific coexpression, rather than species-specific coexpression, in which detected coexpression was also supported by conserved coexpression relationships.  Secondly, statistical characteristics of the coexpression index were not clear.  We have adopted the MR index of gene-to-gene correlations as a coexpression measure because it has a higher predictive power for gene function than the Pearson's correlation coefficient (PCC) ([@pcx191-B23]).  However, the MR index was not statistically characterized, thereby restricting the usage of this index for meta-analyses of coexpression data.  Here, we report an update of ATTED-II that focuses on the following three aims to promote meta-coexpression analyses among species: (i) establishment of comparable assessment measures of gene coexpression for different species; (ii) investigation of the statistical properties of the MR index; and (iii) construction of high-quality coexpression data.  As a result of these revisions, ATTED-II provides multispecies coexpression data with improved accuracy and usability. 

Results and Discussion
======================

New ATTED-II coexpression data for nine species
-----------------------------------------------

We updated both the microarray-based and RNAseq-based coexpression data in ATTED-II (**[Table 1](#pcx191-T1){ref-type="table"}**), based on Affymetrix GeneChip data in ArrayExpress ([@pcx191-B15]) and Illumina RNAseq data in the DNA Data Bank of Japan (DDBJ) ([@pcx191-B29]).  We also added two sets of RNAseq-based coexpression data for barrel medick (*Medicago truncatula* RNAseq; Mtr-r) and grape (*Vitis vinifera* RNAseq; Vvi-r).  As a result, ATTED-II provides two coexpression platforms for each of nine species, except for field mustard (*Brassica rapa*) and poplar (*Populus trichocarpa*) for which only a single platform is available.  Multiple platforms for individual species are preferable for determining the reproducibility of coexpression data, as shown below.  These coexpression data are available in the search or draw tools in ATTED-II.  To assist users, examples of queries for every species are provided (<http://atted.jp/top_search.shtml> and <http://atted.jp/top_draw.shtml>).  Table 1Coexpression data in ATTED-II version 9SpeciesPlatform ID[^a^](#tblfn1){ref-type="table-fn"}VersionGenesSamplesLogit-MR errorFunction score[^b^](#tblfn2){ref-type="table-fn"}Reproducibility[^c^](#tblfn3){ref-type="table-fn"}*Arabidopsis thaliana*Ath-mc7.120,81916,0330.375.430.136Ath-rc3.022,7602,1200.715.17*Brassica rapa*Bra-rc2.128,9781881.044.77--*Glycine max*Gma-mc3.115,7461,1310.743.370.076Gma-rc3.08,3735991.027.64*Medicago truncatula*Mtr-mc3.120,3769751.044.430.021Mtr-rc1.13,753411.462.65*Oryza sativa*Osa-mc6.119,8672,2500.764.980.041Osa-rc2.124,4373361.074.06*Populus trichocarpa*Ppo-mc2.121,9107651.103.82--*Solanum lycopersicum*Sly-mc2.15,7214011.044.080.041Sly-rc2.120,5642821.013.87*Vitis vinifera*Vvi-mc3.19,4213141.144.470.028Vvi-rc1.118,5873460.903.10*Zea mays*Zma-mc3.110,7778061.114.620.055Zma-rc2.132,2741,7940.884.42[^1][^2][^3]

Reproducibility of co expression data
-------------------------------------

Determining the reproducibility of coexpression data is straightforward.  First, we checked microarray-based and RNAseq-based coexpression data for Arabidopsis, with the longest history of coexpression data in ATTED-II. **[ Table 2](#pcx191-T2){ref-type="table"}** shows overlap in the coexpression edges between each data type, i.e. microarray-based Arabidopsis coexpression data (Ath-m) and RNAseq-based data (Ath-r).  This table clearly shows the successful development of Arabidopsis coexpression data in ATTED-II, with the highest Jaccard coefficient between the current coexpression data, Ath-m.c7-1 and Ath-r.c3-0.  Reproducibility results for the other species are summarized in **[Table 1](#pcx191-T1){ref-type="table"}**.  In addition to the highest level of reproducibility with the Arabidopsis data, data for soybean (*Glycine max*) and maize (*Zea mays*) were relatively highly reproducible, whereas coexpression data for *M. truncatula* and *V. vinifera* were not consistent among data types within each species, suggesting that one or both of the coexpression datasets (i.e. microarray and RNAseq) is of low quality.  Table 2Jaccard coefficients of common edges among a series of coexpression data for Arabidopsis in ATTED-IIAth-r.c3-0Ath-r.c2-0Ath-r.c1-0Ath-m.c7-00.1340.0550.038Ath-m.c6-00.1110.0570.040Ath-m.c5-00.1060.0560.040Ath-m.c4-10.0780.0460.032Ath-m.c3-10.0610.0420.029[^4][^5]

Consistency of gene co expression data with gene function
---------------------------------------------------------

Although reproducibility can be measured without any external data, this measure mutually depends on the coexpression datasets.  In other words, this index does not indicate the quality of individual gene coexpression data.  Another concern is that it is possible to select a limited set of samples to achieve highly reproducible coexpression data.  In such cases, high reproducibility does not indicate high quality.  As another quality measure, the enrichment of functional annotations in coexpressed genes is useful because gene coexpression information is used to identify functionally related genes.  Gene Ontology (GO) annotation provides a rich resource for the functional annotation of genes ([@pcx191-B9]).  Conveniently, the directed acyclic graph structure of GO terms enables the selection of terms that have a particular range of information content ([@pcx191-B17]).  However, GO annotations are not sufficient for plant species in ATTED-II, except Arabidopsis, and thus cannot be used for the quality assessment of coexpression data among species ([@pcx191-B28]).  To resolve this problem, we adopted the KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway annotation ([@pcx191-B14]) as an alternative method.  To select highly informative annotations, we excluded pathways associated with \>100 genes in a genome, resulting in 102.2 pathways, on average, for the nine species.  To characterize the selected KEGG pathway annotations, we first checked the consistency of the KEGG evaluation results with the GO evaluation results using a series of coexpression data for Arabidopsis.  The discriminative power based on gene pairs with and without common gene function annotations showed a similar trend between the two annotation types (GO and KEGG), supporting the validity of KEGG pathway annotation for the quality assessment of coexpression ([Fig.  1](#pcx191-F1){ref-type="fig"}A).  Because KEGG pathway annotations are provided for a broad range of species with almost the same annotation density across species, this will be useful not only for comparisons among plant species, but also for comparisons of coexpression data across species in distinct kingdoms, such as microalgae and animals ([@pcx191-B30], [@pcx191-B2]). 

![ Quality assessment based on the consistency of known gene functions.  As a measure of the quality of gene coexpression data, the power to discriminate gene pairs sharing a common functional annotation from other gene pairs was used. ( A) Previous and current Arabidopsis coexpression data were assessed by this discrimination analysis.  Irrespective of the gene annotation source (GO Biological Process or KEGG pathway), the quality trend was consistent. ( B) The current 16 coexpression platforms were assessed with the discrimination analysis using the KEGG pathway annotation.  Circles and diamonds indicate microarray-based and RNAseq-based coexpression datasets, respectively.](pcx191f1){#pcx191-F1}

Using the function score based on the KEGG pathway annotation, the efficiency of coexpression for all coexpression platforms was compared (**[Table 1](#pcx191-T1){ref-type="table"}**; [Fig.  1](#pcx191-F1){ref-type="fig"}B).  The function scores were correlated with the logarithmic number of samples for each species (PCC = 0.44), as reported previously ([@pcx191-B4]).  These trends were observed for both microarray- and RNAseq-based data.  The only exception was Gma-r, which showed a high score, despite a relatively small number of samples.  In this functional scoring, we did not eliminate paralogous gene pairs, which usually have similar functional annotations and expression profiles, resulting in overestimation of scores in species having a large number of paralogous genes.  Furthermore, the function score is based on unevenly distributed annotations across genes in a species.  For example, highly expressed genes tend to be well studied and thus well annotated (**[Supplementary Fig.  S1](#sup1){ref-type="supplementary-material"}**).  Therefore, the function score for coexpression may not reflect the overall performance of coexpression data.  Nonetheless, the KEGG scores for each species were generally consistent with the reproducibility scores (**[Table 1](#pcx191-T1){ref-type="table"}**).  As described above, reproducibility is mutually dependent, and thus a lower quality of coexpression data for a species mainly limits the reproducibility score.  In fact, reproducibility scores were well correlated with the smaller KEGG scores across the two platforms in a species (PCC = 0.82). 

MR as a statistical indicator
-----------------------------

We previously reported that the MR index is a powerful indicator of the co-function of a gene pair ([@pcx191-B23]).  However, the poor characterization of this index limits downstream analyses.  To investigate the distribution and confidence interval of the MR index, we repeatedly calculated the coexpression matrix from bootstrapped samples of the experimental unit. [ Fig.  2](#pcx191-F2){ref-type="fig"}A shows the SD of the MR index for 100 sets of bootstrapped RNAseq expression data in Arabidopsis.  The horizontal axis represents the mean of 100 bootstrapped MR values expressed as a percentile from a positive correlation (small MR) to a negative correlation (large MR).  The MR ranges from '1' to 'the number of genes minus 1' (e.g. 1--22,759 for the Ath-r.c3-0 platform, which includes 22,760 genes).  The vertical axis represents the median (black curved line) and the first and third quantiles (gray curved lines) of the SD of the bootstrapped MR values.  The MR index was precise near each end of its range, but in the middle of the range, indicating no correlation, MR was easily affected by chance ([Fig.  2](#pcx191-F2){ref-type="fig"}A).  Additionally, the skewness and kurtosis of the error distribution indicated that the data were not normally distributed (**[Supplementary Fig.  S2A](#sup1){ref-type="supplementary-material"}**, **B**).  In fact, this non-normal distribution is commonly observed in the PCC (*r*) of samples, which can be converted to a normal distribution using a logit transformation known as the Fisher transformation ([@pcx191-B8]).  This logit transformation can also be applied to the range-standardized MR index (0 \< *p*~*i*~ \< 1), which is calculated by dividing by the total number of genes (*N*), $p_{i} = MR_{i}/N$. $$\text{log}\text{i}\text{t}\left( p_{i} \right) = \operatorname{}\text{log}\left( \frac{p_{i}}{1 - p_{i}} \right)$$ Note that this transformation is equivalent to the Fisher transformation with the range standardization of *r* as $p_{i} = \left( r + 1 \right)/2$. [Fig.  2](#pcx191-F2){ref-type="fig"}B shows the distribution of SD for the logit-transformed MR values (hereafter, logit-MR error).  The logit-MR error was nearly constant across the entire range of MR values.  Along with the skewness and kurtosis of the distribution (**[Supplementary Fig.  S2C](#sup1){ref-type="supplementary-material"}**, **D**), the logit-MR error can be regarded as a normal distribution with a constant SD for the entire range of MR values.  In the case of Ath-r.c3-0, shown in [Fig.  2](#pcx191-F2){ref-type="fig"}B, the average logit-MR error was 0.87.  The logit-MR errors varied across the coexpression platforms (**[Table 1](#pcx191-T1){ref-type="table"}**) and were strongly correlated with the logarithmic number of samples for each coexpression platform (PCC = −0.90) ([Fig.  2](#pcx191-F2){ref-type="fig"}C).  This characteristic is useful for estimating confidence intervals for the MR index solely based on the number of samples.  For example, the logit-MR error for coexpression data based on 500 samples will result in a value of approximately 1. **[ Table 3](#pcx191-T3){ref-type="table"}** shows the 90% confidence intervals for typical MR values with logit-MR error = 1 or 0.5.  ATTED-II shows the top 300 coexpressed genes as the default.  Taking into account the confidence intervals of the MR index, this list typically includes most of the top 100 coexpressed genes.  Note that the logit-MR error on the rightmost edge of the MR percentile is slightly larger than that on the leftmost edge ([Fig.  2](#pcx191-F2){ref-type="fig"}B), suggesting that negative coexpression relationships are more unstable than positive relationships in the MR calculation.  Table 3The 90% confidence intervals of typical MR values with different error levelsMRBootstrap SD = 0.5Bootstrap SD = 110.81.40.72.131.94.91.38.3105.917.33.530.23014.252.66.692.510056.9175.832.4308.4

![ Properties of the error distribution of MR and logit-MR values.  SDs of the bootstrapped MRs from RNAseq-based Arabidopsis coexpression data (Ath-r.c3-0) are shown. ( A, B) SDs of bootstrapped MR values (A) and of logit-MR values (B) are shown against the mean of the bootstrapped MR values.  The black lines show the median values with a sliding window corresponding to the 0.01 percentile of the MR without overlap, whereas the gray lines represent the first and third quantiles. ( C) Mean SDs for the current 16 coexpression platforms in ATTED-II are plotted against the number of samples for each platform.  Circles and diamonds indicate microarray-based and RNAseq-based coexpression datasets, respectively.](pcx191f2){#pcx191-F2}

Integration of multiple MR values
---------------------------------

Based on the normality of the logit-transformed MR values, we can integrate multiple MR values.  In general, the arithmetic average can be generalized using the transformation function $f$ and its reverse function $f^{- 1}$. $$Average_{F} = f^{- 1}\left( \sum f\left( x_{i} \right)/\sum i \right)$$ For example, using identity, logarithmic and rational functions as the transformation functions results in arithmetic, geometric and harmonic averages, respectively.  Here, to determine the average with a logit transformation, we adopted a logit function for $f$ and a sigmoid function for $f^{- 1}$. For *k* MR values (*MR*~1~, *MR*~2~, ..., *MR*~*k*~) with individual weights (*w*~1~, *w*~2~, ... *w*~*k*~; Σ*w*~*i*~ = 1), the weighted logit average of the MR values is as follows: $$MR_{average} = \frac{N\Pi MR_{i}^{w_{i}}}{\Pi\left( N - MR_{i})^{w_{i}} + \Pi MR_{i}^{w_{i}} \right.}$$ Note that this equation can be approximated as a weighted geometric mean under *MR*~*i*~ \<\< *N*.  $$MR_{geometricmean} \simeq \Pi MR_{i}^{w_{i}}$$ In ATTED-II, the CoExSearch tool provides an integrated list of coexpressed genes for multiple gene queries (<http://atted.jp/top_search.shtml#CoExSearch>).  This tool uses equal weights ($w_{i} = 1/k$) from [Eqn.  3](#E3){ref-type="disp-formula"}, simplified as follows: $$MR_{average} = \frac{N\sqrt[k]{\Pi MR_{i}}}{\sqrt[k]{\Pi\left( N - MR_{i} \right)} + \sqrt[k]{\Pi MR_{i}}}$$ Even given the high manageability of logit-transformed MR values, we have retained the original MR index in the ATTED-II database.  Because the MR index is derived from the ranking of the coexpression strengths, this index implies the maximum size of the coexpression network including the guide gene and its coexpressed genes. 

Slight, but stable improvement in coexpression quality based on the bootstrapping procedure
-------------------------------------------------------------------------------------------

Sample bootstrapping is widely used to improve the generalization ability of a model ([@pcx191-B5]).  As we proposed to use the average method of the MR index, we could use the average of the bootstrapped MR values as a coexpression index.  As a result, the bootstrapping and average approach (a simple bagging approach) returned substantially better results than the original coexpression dataset using all samples at once ('Bagging effect' in **[Table 4](#pcx191-T4){ref-type="table"}**).  Although the gains from the bagging procedure were not large for most coexpression platforms, this procedure resulted in relatively greater gains for platforms with a smaller number of samples.  The PCC between the bagging effect and the logarithmic number of samples was −0.42, indicating that the bagging procedure was not effective for platforms that include larger numbers of samples.  This may be explained by the limited number of bootstrap replicates in our procedure (i.e. 100).  Table 4Comparison of calculation methodsPlatformNo.  of genesNo.  of samplesQCQCQBCombat effectBagging effectAth-m20,81916,0335.465.425.430.991.00Gma-m15,7461,1313.293.353.371.021.01Mtr-m20,3769753.414.324.431.271.03Osa-m19,8672,2504.604.924.981.071.01Ppo-m21,9107653.493.743.821.071.02Sly-m5,7214013.313.804.081.151.07Vvi-m9,4213143.934.224.471.071.06Zma-m10,7778064.004.464.621.121.03Ath-r22,7602,1204.565.125.171.121.01Bra-r28,9781884.544.634.771.021.03Gma-r8,3735996.567.607.641.161.01Mtr-r3,753412.542.622.651.031.01Osa-r24,4373363.433.804.061.111.07Sly-r20,5642823.553.703.871.041.05Vvi-r18,5873463.213.073.100.961.01Zma-r32,2741,7943.984.354.421.091.01[^6]

Batch normalization of gene expression data
-------------------------------------------

Batch effects are a major source of technical noise in transcriptome data ([@pcx191-B16], [@pcx191-B10]).  In the previous coexpression calculation procedure in ATTED-II, we applied batch zero-centering ([@pcx191-B36]), which standardizes the mean expression level of each gene in each experiment to zero ([@pcx191-B34]).  On the other hand, in our previous investigation, direct standardization of variances of each gene in each experiment did not improve the quality of coexpression results, probably due to loss of fold-change information in gene expression.  ComBat is a method that can be used to estimate stably the mean and variance of the batch effect by the empirical Bayes method ([@pcx191-B12]) and is effective for transcriptome and coexpression analyses ([@pcx191-B20], [@pcx191-B37]).  However, batch estimation for data with an unbalanced batch-group design risks introducing another bias ([@pcx191-B21]).  Because the expression data for coexpression calculation are aggregated from different experiments, this is a case of the unbalanced batch-group design.  Therefore, we investigated the effect of variance standardization by the empirical Bayes estimation in addition to batch zero-centering for a set of the current coexpression platforms in ATTED-II (**[Table 4](#pcx191-T4){ref-type="table"}**).  As a result, batch normalization with gene- and experiment-specific noise variances had notable effects on microarray-based and RNAseq-based coexpression data ('ComBat effect' in **[Table 4](#pcx191-T4){ref-type="table"}**).  Based on this finding, we incorporated the ComBat normalization method into the latest calculation pipeline in ATTED-II. 

Materials and Methods
=====================

Construction of gene coexpression data
--------------------------------------

In ATTED-II version 9, the coexpression calculation procedure was slightly modified from the previous procedure to reduce the calculation cost sufficiently to enable bootstrapping trials for multiple species.  For microarray-based coexpression data, we downloaded the GeneChip CEL files from ArrayExpress ([@pcx191-B15]), normalized them based on the RMA method ([@pcx191-B11]) and applied ComBat normalization ([@pcx191-B12]) using an experimental unit as a batch.  Before calculating correlations among all probe pairs, we selected a single probe for each gene.  We made this selection based on the similarity of coexpression patterns for the same gene in the RNAseq coexpression data.  For the microarray platforms of *M. truncatula* (Mtr-m) and *P. trichocarpa* (Ppo-m), reliable RNAseq coexpression data were not available.  For these two platforms, the expression patterns of multiple probes were averaged to generate a single pattern of expression for that gene.  After the selection of probes, a gene-to-gene correlation matrix, which is identical to the probe-to-probe correlation matrix, was calculated by PCC, and this was then converted to the MR index, $MR_{ij} = \sqrt{R_{ij}R_{ji}}$, where $R_{ij}$ indicates that gene *j* is the $R_{ij}$-th strongest coexpressed gene for the guide gene *i* ([@pcx191-B23], [@pcx191-B24]). 

For RNAseq-based coexpression data, we downloaded the Sequence Read Archive format data from the DDBJ ([@pcx191-B29]) and mapped it to NCBI RefSeq sequences ([@pcx191-B6]).  For this mapping, we employed unique sequence signatures for each gene using Matataki software (<https://github.com/informationsea/matataki>; Okamura and Kinoshita, in preparation), which enables much faster quantification of RNAseq data than using normal mapping procedures.  After conversion to a base-2 logarithm with a pseudo-count of 0.125, ComBat normalization was applied ([@pcx191-B12]), and the average expression levels were subtracted for each gene and experiment.  Using all of the experimental data simultaneously, PCCs between all gene pairs were calculated, and these values were then converted to the MR index. 

Bootstrapping procedure for gene coexpression data
--------------------------------------------------

To characterize the MR index, we conducted a bagging procedure for gene coexpression data.  One bootstrap coexpression dataset was calculated based on randomly selected experiments with replacement to generate the same number of samples as the original dataset.  In this study, we repeated this bootstrap calculation of coexpression 100 times for every coexpression platform.  The coexpression values between any gene pair should vary with a particular mean and deviation depending on the randomly selected sample.  In general, the bootstrap distribution can be used as an estimate of the population distribution and thus can be used to estimate confidence intervals for the value.  In addition, the mean of the bootstrapping coexpression values (bagging MR value) was used as another estimate of the MR index. 

Reproducibility between microarray-based and RNAseq-based gene coexpression data
--------------------------------------------------------------------------------

To measure the reproducibility of gene coexpression data from different coexpression platforms in the same species, we used genes that were in common between the coexpression platforms.  Then, the top three coexpression relationships from every analyzed gene were used as the edge for this assessment.  This edge criterion is also used in the network drawing functions in ATTED-II.  The Jaccard coefficient between the two sets of edges was used as the reproducibility score. 

Discriminative power of gene coexpression data for the identification of gene pairs with common functional annotations
----------------------------------------------------------------------------------------------------------------------

To validate the biological significance of coexpressed gene data, we used the enrichment of gene pairs sharing common functional annotations relative to other pairs.  For coexpression data in Arabidopsis, GO Biological Process annotations with \<50 genes were used.  All of the genes within each dataset were divided into two groups, i.e. gene pairs sharing at least one GO annotation and gene pairs without any common annotation.  The difference in the distributions of the degrees of coexpression for the two groups was assessed using the partial AUROC~0.01~, which is the area under a part of the receiver operating characteristic curve with a false-positive rate of \< 0.01 ([@pcx191-B18]).  To improve visibility, the partial AUROC~0.01~ multiplied by the square value of the false-positive rate (10,000 in this study) was used as the function score, so that the function scores using randomized gene--function associations were always 0.5.  To apply this evaluation method for coexpression data to any species in ATTED-II, KEGG pathway annotations downloaded from KEGG FTP were used ([@pcx191-B14]) (June 26, 2017).  We selected KEGG pathways associated with \<100 genes with a relatively high information content.  The genes associated with at least one selected KEGG pathway were then used in this assessment. 

Supplementary Data
==================

[Supplementary data](#sup1){ref-type="supplementary-material"} are available at PCP online. 

Funding
=======

This work was supported by the Japan Society for the Promotion of Science (JSPS) \[KAKENHI grant Nos.  24114005, 15K18464 and 16HP7003 to T.O., and 15K20863 to Y.A.\]. 

Supplementary Material
======================

###### 

Click here for additional data file. 

We thank Mr. Kota Jin for the web design and species icons for ATTED-II.  Super-computing resources were provided by the Human Genome Center, Institute of Medical Science, University of Tokyo. 

Disclosures
===========

The authors have no conflicts of interest to declare. 

AUROC

:   Area under the receiver operating characteristic curve

GO

:   Gene Ontology

KEGG

:   Kyoto Encyclopedia of Genes and Genomes

MR

:   mutual rank

PCC

:   Pearson's correlation coefficien

RNAseq

:   RNA sequencing

[^1]: ^a^Xxx-m, microarray-based coexpression; Xxx-r, RNAseq-based coexpression. 

[^2]: ^b^Predictive performance of the KEGG annotation represented by partial AUROC (1E-04).  A higher score indicates a better performance. 

[^3]: ^c^Jaccard coefficient for common edges between the platforms in the same species.  The top three coexpressed genes from every gene were used as edges. 

[^4]: Xxx-m, microarray-based coexpression; Xxx-r, RNAseq-based coexpression. 

[^5]: Note that a Jaccard coefficient of 1 indicates complete overlap between the two sets of coexpression edges, whereas a Jaccard coefficient of 0 indicates no overlap. 

[^6]: Q, quantile normalization; CQ, ComBat-Quantile normalization; CQB, Bagging procedure for ComBat-Quantile normalized expression data.