Cáceres and Gonzalez BMC Genomics (2019) 20:26 https://doi.org/10.1186/s12864-018-5340-3 METHODOLOGY ARTICLE Open Access A measure of agreement across numerous conditions: assessing when changes in network structures are tissue-specific Alejandro Cáceres1,2* and Juan R. Gonzalez1,2,3 Abstract Background: There is great interest to study how gene pathways change their structure across different tissues. The assessment of inter-study reliability of pathway changes across tissues can inform on the fraction of tissues with specific functional changes in network structure. However, there is a lack of agreement measures among studies that independently observe how a group of observations change across conditions. We, therefore, propose λ, a new inter-study reliability measure that determines the consistency to distinguish observations by condition. Results: We derived λ’s distributional characteristics, determine its reliability properties and compared it with Cohen’s κ. We studied the co-expression structure of 287 gene pathways across four brain regions in two transcriptomic studies and applied λ to assess the inter-study reliability of the pathways’ brain-regional changes. Brain-related pathways showed highest λ; the top value was for the nicotine addiction pathway whose structure was reliably distinguishable among regions with dopaminergic projections. Conclusion: Our results offer novel substantial evidence that changes in network structure across tissues can be inferred independently of samples, algorithms and experiments (RNA-sequencing or microarrays). Reliability measures, such as λ, can inform on the tissues where changes in a network’s structure are likely functional. An R package is available at https://github.com/isglobal-brge/lambda. Keywords: Brain, Transcriptome, Co-expression networks, Reliability, Coehn’s kappa, GTEx, RNA-sequencing, Addiction, Nicotine Background [4]. Inferred networks or established gene pathways may Reproducibility is a pressing issue in biomedical research be additionally assessed under different and numerous that particularly worries a large number of researchers conditions to understand the physiological changes of in the field [1]. Research guidelines from leading jour- their associated biological functions. The interest here nals and the American Statistician Association urge for is, for instance, to determine for which tissues a path- the need of confirmation studies and accurate statisti- way has a specific function [5]. Therefore, if the pathway cal reporting [2, 3]. In systems biology, gene interaction is physiological only in specific tissues, the inter-study networks are often derived from the integration of high- reproducibility of its structure across other tissues, where throughput data with the aim to determine gene struc- it is not functional, is expected to be low. As such, neu- tures with probable biological functionality. Networks, robiological networks would be expected to meaning- therefore, also need to be reproduced between indepen- fully change in structure between nervous tissues but not dent studies in order to contain valid scientific content between connective tissues like blood [6]. While pair-wise comparisons of a network reproducibility between two *Correspondence: alejandro.caceres@isglobal.org studies can be applied on a single tissue [7], there is a need 1 ISGlobal, Doctor Aiguader, 08003 Barcelona, Spain to measure the reproducible changes in network structure 2 Centro de Investigación Biomédica en Red en Epidemiología y Salud Pública (CIBERESP), Barcelona, Spain across multiple tissues which, in particular, informs on Full list of author information is available at the end of the article the fraction of tissues for which the changes may underlie © The Author(s). 2019 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Cáceres and Gonzalez BMC Genomics (2019) 20:26 Page 2 of 13 tissue-specificity. We, therefore, aimed to propose a mea- A first question that can be asked is the degree to which sure that allows us to determine the degree to which the the most similar observations between studies correspond changes in network structure among a range of tissues are to those under the same condition, and hence have a reproducible across studies; Fig. 1 shows a schematic rep- measure of how reliably the observations can be distin- resentation of the data and the type of network reliability guished by condition. Similarity can be assessed by any we intent to assess. statistic, namely a preservation statistic, that compares Classical statistics include numerous ways to measure two sets of observations at any particular property. A cor- the reliability of an observation [8]. Reliable observa- relation of the observations between two conditions is an tions are reproducible and accurate. Agreement measures example of one preservation statistic. If a preservation between two independent experiments are used to assess statistic for observations between studies is defined [4], the consistency of the observations being made. If obser- we can ask whether the preservation is maximum when vations are categorical, Cohen’s κ and its generalizations the conditions between studies match. The question can are typically used [9, 10]; if observations are continuous be addressed mathematically by computing a preserva- then a number of correlation measures can be used, such tion matrix, whose elements are the preservation statistics as intra-class correlations [11]. These and other agree- between conditions across studies, and assess the extent ment measures are suitable when experiments are per- to which the diagonal terms are maxima across rows and formed under comparable, unchanged conditions. When columns. In particular, a preservation matrix across con- studies are designed to test how a group of individuals ditions can be defined for cross-tabulated tables used in change under a range of varying conditions, it is of inter- classical inter-study reliability studies and, therefore, the est to assess how reliable the changes across conditions ability to correctly pair observations by condition can are. However, there is a lack of recommended inter-rater be studied in those cases. We thus aimed to construct measures that assess reliable changes among conditions a reliability measure for the correct condition pairing of [8, 12]. observations between studies, compare it with Cohen’s κ Fig. 1 Definition of λ for a co-expression network inferred in 6 tissues in two different studies. Schematic representation of a m-gene pathway (A) whose co-expression structure is inferred from the transcriptomic data of different tissues/conditions (1 to 6) across two studies (E1, E2). The inference for tissue 1 from data is illustrated. The figure shows the preservation matrix Z (blue) for the network A that is constructed from the pair-wise preservation statistics between tissues across different experiments. The arrows show the pairing between tissues across experiments given by the maximum Z found in the rows and columns that cross each diagonal term of Z. Network structures in tissues 1, 2, 4 and 5 in E1 are correctly matched by tissue in E2, as their diagonal elements are row and column maxima. Tissues 3 and 6 are not matched correctly. The fraction of tissues for which network structures in E1 can be correctly paired to those in E2 is 4/6, illustrated by the red stars that fall in the diagonal. λ is a reliability measure defined as the probability that the diagonal terms in Z are their row and column maxima. λ informs on how diagonal Z is Cáceres and Gonzalez BMC Genomics (2019) 20:26 Page 3 of 13 and apply it to assess the inter-study reproducibility of and 5 are correctly paired between experiments, as the changes in co-expression network structure across tissues. corresponding diagonal terms are their row and column The GTEx project is an unprecedented effort to study maxima and, therefore, their differences are reproduced the gene expression in tens of tissues in hundreds of sub- between experiments. jects using RNA-sequencing [13]. It is, therefore, a strong We define a measure λ on the preservation matrix that candidate for becoming a preferred benchmark for inter- tells us how diagonal Z is. In general, we assume that the action networks across different tissues. As the number of units’ observations corresponding to row i can be paired studies with expression data in multiple tissues is expected by similarity with those in column j if the preservation to increase, agreement measures against GTEx may serve statistic zij is maximum across row i and column i . Note to assess the reproducibility of network changes across that, as a consequence, Z is non-symmetric. We then tissues [3]. Some studies that measure gene expression define λ as the expected fraction of correct condition in different tissues are publicly available, one of which is pairings between experiments (R), given by the BRAINEAC project [14]. Here, the gene expression using microarray data was measured in hundreds on indi- 1 λ = E[R] = pii (1) viduals, free of neurodegenerative disorders at the time k i of death, in nine different brain tissues, four of which 1  overlap with those in GTEx. Using our agreement mea- σ 2 = Var[R] = 2 pii (1 − pii ) (2) sure between the studies, we assessed the fraction of k i brain regions for which 287 KEGG pathways [15] reliably changed across the cerebellum, frontal cortex, hippocam- where pii is the probability that the diagonal term zii pus and putamen. is maximum within its row and column elements. λ measures the fraction of conditions that can be reli- Results ably distinguished by structural changes, taking into For two studies that measure the same group of units account a probability for random pairing. We mod- (subjects, co-expression pairs in a network, etc.) across k eled pii assuming independence among the similarities conditions, we propose to construct the k × k preserva- of observations between conditions, and normality in tion matrix Z between conditions. The matrix element zij the preservation statistic (see “Materials and methods” (i, j = 1, ...k) is the correlation (or any preservation statis- section). As an agreement measure, we derived its tic) of the units’ observations between condition i in the reliability properties and compare it with a classical first study and condition j in the second study. When the measure of inter-study agreement, Cohen’s κ. We then observations are the connections between all node-pairs applied it to determine which of 287 KEGG pathways pre- in a network then zij may be any normally distributed sented higher rates of correct condition pairing across statistic where pairwise topologies are compared. Figure 1 four brain regions between two independent transcrip- shows the case in which a co-expression network A of tomic studies. m gene transctipts is inferred for n subjects in 6 tissues (1,...6) in 2 different experiments (E1 and E2). The matrix Comparison between λ and Cohen’s κ in simulation studies Z is given by all pair-wise comparisons between tissues We studied the properties of λ with extensive simulations. for the network connections across both experiments. The While λ is applicable to more general situations than those arrows show the values of Z that are maxima across the covered by Cohen’s κ, we compared the performance of rows and columns that intersect at each diagonal term. λ with κ in simulated cases where both measures can be We are interested in measuring how reproducible the used. We simulated numerous agreement tables between changes of units observations, such as network edges, are. studies, in which 500 units were classified by both stud- In the case of perfect inter-study reproducibility across ies into numerous categories/conditions, under different all conditions/tissues then the matrix Z looks diagonal. scenarios. From the cross-tabulation table of agreement, That is, zii is maximum across row i and column i , for see expression 14 in “Materials and methods” section, all i. In general, this is not the case and the preservation we computed κ: the inter-study reliability of classification corresponding to row i and column j  = i may be maxi- that accounts for random agreement. From the corre- mum across the elements in row i and column i . Figure 1 sponding preservation matrix obtained from Eqs. 17 and shows, for instance, that the element of Z at tissue 3 in 18 in “Materials and methods” section, we computed λ: E1 and tissue 2 in E2 is the maximum in row 3 and col- the fraction of conditions that can be correctly distin- umn 3 of Z. As the network structure of tissue 3 in E1 is guished between experiments. We also computed P0 : the paired with that of tissue 2 in E2, network connections in proportion of units correctly classified between experi- tissue 3 have not been reliably reproduced between exper- ments; and r: the observed fraction of conditions that are iments. By contrast, all the structures at tissues 1, 2, 4 correctly paired between experiments. Cáceres and Gonzalez BMC Genomics (2019) 20:26 Page 4 of 13 We recreated the reproducibility assessment between Table S1 and “Materials and methods” section for further two raters/experiments across varying number of condi- details. tions/tissues k = (5, 7, 15), where three types of cross- Simulations showed that when null agreement is tabulation tables between raters across conditions were expected (P0 = 0) λ is 0, while λ is 1 for full agree- allowed (expression 14 in “Materials and methods” section). ment (P0 = 1). Consistently with this, we observed that Those were tables with: 1) marginal equiprobabilities 1/k λ increased monotonically with κ for all the simulation of both experiments across all tissues (scenario 1), 2) scenarios, see Additional file 2: Figure S1. The functional marginal probabilities ∼ 1/j for tissue j (scenario 2) and 3) dependence was highly stable under different scenarios, marginal probabilities ∼ 1/j2 (scenario 3), where j = 1, ...k. revealing, as expected, high λ agreement for fair values Simulation scenarios 2 and 3 were designed to test situ- of κ (0.2, 0.4), given that the latter is a measure of exact ations where agreement between studies tend to concen- agreement rather than discriminative agreement. There- trate around one tissue j = 1. We performed 10,000 sim fore, the agreement measured by λ is consistently higher ulations per number of conditions and scenario. Simulations than the agreement measured by κ. We also observed that were extracted from the permutation of 50 diagonal for low values, λ tends to zero when κ takes small negative cross-tabulation tables (perfect agreement) and 50 cross- values, a situation already described in Cohen’s work [9]. tabulation tables with null diagonals (perfect disagree- For a given κ, we observed a sizable range of λ values, in ment). Conditions across subjects were permuted a total particular, as tissues become less equiprobable (scenario 3). of 30,000 for each table, preserving their marginal fre- We noted that if the number of tissues is small (k = 5) quencies, given by each simulated scenario. Every 30 per- and the marginal distribution greatly concentrates around mutations, agreement measures were computed obtaining one single tissue (j = 1 for scenario 3), then λ tends to 1/k a total of 5000 values for perfect agreement (Fig. 2) and (0.2) because the experiments can clearly distinguish that 5000 for perfect disagreement. Simulations thus covered tissue from the rest. In this case, κ tends to zero. the entire interval between perfect agreement (P0 = 1) We then studied the relationship between the agree- and null agreement (P0 = 0). See Additional file 1: ment measures that account for random agreement with Fig. 2 Simulation scheme for reliability measures under one of the proposed scenarios (scenario 2). In the figure, four conditions/tissues are considered (A, B, C, D) in which two raters (1, 2) measure 500 subjects. The cross-tabulation tables between raters are shown, all with marginal frequencies across raters given by ∼ 1/j (j = 1, 2, 3, 4 and A = 1, B = 2, C = 3, D = 4), corresponding to the simulation scenario number 2 described in “Materials and methods” section. Initially, all 500 individuals are correctly classified between raters. 3000 permutations are performed, such that the marginals frequencies are conserved. Every 30 permutations reliability measures (P0 , κ, λ and r) are computed. The process is repeated 50 times. A second simulation is performed from an initial cross-tabulation of perfect disagreement, where diagonal terms are all 0. A total of 10,000 values of each measure are then obtained for simulations that cover the whole reliability interval (0,1) under each scenario Cáceres and Gonzalez BMC Genomics (2019) 20:26 Page 5 of 13 those that do not. In our simulations, we confirmed that determined with high precision, low σ and therefore small κ is lower than the proportion of agreement P0 (Fig. 3); CIs for λ. From a practical point of view, low or moderate which illustrates the initial motivation for κ’s definition as values of λ can be used to select the tissues for which the a measure that corrects for random agreement [9]. Simi- network changes can be reliably measured, and are likely larly, r, the observed fraction of times the diagonal terms tissue-specific. Consider, for instance, a situation where in the preservation matrix are row and column maxima the structure of a gene pathway is inferred in 50 differ- (Eq. 20 “Materials and methods” section), is higher than λ, ent tissues. If we expect that the network is biologically a distributional estimate of such fraction. Note also that r, functional in only 5 of the tissues, then λ will be at most as a realization of the random variable R, is discontinuous 0.1. In addition, low σ will also indicate that the network with k + 1 possible values, while λ is a continuous variable structure in 90% of tissues could be discarded on grounds ranging from 0 to 1. of measurement reliability. By contrast, κ is maximally We additionally studied the variance σ of the fraction R, informative only as it tends to 1, that is when all the con- for which r is a realization and λ its expected value. We ditions are distinguishable and the preservation matrix is observed that σ decreased with the number of tissues and fully diagonal. with departure from marginal equiprobability (Additional We finally performed a power test for λ, where we simu- file 3: Figure S2). We observed that for a given λ, multi- lated a true agreement scenario in which 3/5 tissues could ple values of σ are allowed (Additional file 4: Figure S3). be reliably paired across two studies. We simulated a true However, the mean of κ, for a single draw of the binomial scenario using a multivariate distribution for 10 random process in Eq. (16), is a function of its variance. In partic- variables corresponding to 5 networks per experiment; 3 ular, κ’s variance is minimum at extreme values (0,1) and networks between experiments correlated with strength is inversely proportional to the number of subjects [16]. of 0.5 while all other possible network correlations, across In contrast, we observed that σ decreases towards zero tissues and studies, varied from 0 to 0.5. We also varied the for a range of values of λ. This occurs when λ tends to r, number of genes N in the network (10, 20, 40) correspond- that is, when the probabilities that the diagonal terms of ing to N(N − 1)/2 non-redundant connections (45, 190, the preservation matrix are maxima tend either to zero or 780). We observed that λ had sufficient power (> 80%) to one, see Eq. (13) in “Materials and methods” section. to detect true agreement when the difference between We also computed the variability of λ from the 5% con- matching tissue correlations and background correlations fidence intervals (CI) obtained from bootstrapping the is 0.05 for 40 gene networks or 0.2 for networks consisting units observations in the cross-tabulation tables. For sim- of 10 genes (Additional file 6: Figure S5). ulations corresponding to 5 equiprobable conditions and varying number of subjects, we observed that the length KEGG pathways of λ’s CI was proportional to σ and independent of the The Kyoto encyclopedia of genes and genomes (KEGG) number of subjects (Additional file 5: Figure S4). offers a list of experimentally characterized biochemichal Overall, these observations show that the fraction of tis- pathways. We computed λ between two independent sues for which network changes are reproducible can be transcriptomic studies, GTEx (RNA-sequencing) and a b Fig. 3 Comparison between agreement measures that correct for random agreement with those that do not. a Cohen’s κ is compared with P0 (the total fraction of agreement). b λ is compared with r (the observed fraction of diagonal elements that in the preservation matrix are their row and column maxima) Cáceres and Gonzalez BMC Genomics (2019) 20:26 Page 6 of 13 BRAINEAC (microarray), that measured gene expres- Figure 4 illustrates the structure of the strongest edges, sion levels with different technologies in the cerebellum, across tissues and studies, of the nicotine addiction path- frontal cortex, hippocampus and putamen. way, which was the top hit with λ = 0.67 and σ 2 = 0.02. For each of the 287 KEGG pathways, we selected the We observed that the structure of the network was pre- commonly annotated genes in each study. We computed served between studies in which, for instance, GABRB3 the co-expression networks, where genes were hubs and was an important hub across tissues and studies while the correlations between the expression levels were edges some genes like GABRR2 and CHRNA6 were consistently between hubs. We then derived the pathways’ adjacency unlinked from the network. While some reliable differ- matrices from the absolute value of the gene-pair corre- ences between tissues can be noticed, the preservation lations in each tissue. As a pairwise preservation statistic, matrix in Fig. 5 clearly shows that the frontal cortex, the we used the adjacency correlation (cor.cor) that quanti- hippocampus and the putamen can be correctly paired fies how similar the connectivity is between two network between studies, while the cerebellum cannot. The lower structures [4], see “Materials and methods” section. We value of λ = 0.67 from the observed R = 0.75 is due then computed the preservation matrices of network con- to the variability of the preservation values and, perhaps, nectivities for each pathway. As a result, we obtained a to the closeness in structure between the hippocampus and λ per pathway. Table 1 shows the pathways with the top the frontal cortex. λ values (> 0.5). While no pathway correctly matched We also computed benchmark-λ as a measure of its structure between studies across all four tissues, we benchmarking the pathway’s structures obtained from observed seven pathways (2%) with agreement between BRAINEAC against GTEx, (Table 2 and Additional file 8: (0.5, 0.75); those are pathways that correctly paired their Table S3). Benchmark-λ corresponded to the expected structures across two to three tissues between studies. value of tissues from BRAINEAC that were correctly Remarkably, five of these pathways are directly linked paired to those in GTEx, given by the maximum over with signaling processes specific to brain, suggesting that rows of the preservation matrix at a given column (see the differences between network structures across brain “Materials and methods” section). We observed higher regions may be tissue-specific. We further observed that values of benchmark-λ, and in particular, 20 pathways in the complete list of all KEGG pathways (Additional with benchmark-λ > 0.7, 15 of which were brain-related. file 8: Table S2), brain pathways were highly ranked by λ. Seven of those pathways are involved in synaptic signaling In the top-ranked pathways, we also observed MicroRNAs (glutamatergic, cholinergic, GABAergic and dopaminer- in cancer and adrenergic signaling in cardiomyocytes. gic synapse, and long-term potentiation, retrograde endo- While these are not brain-specific pathways, numerous cannabinoid and calcium signaling pathways), and four MicroRNAs have been shown to express in brain [17] and in addiction processes (nicotine, morphine, amphetamine adrenergic signaling plays an important role in long-term and alcohol). potentiation in hippocampus [18]. We also computed the bootstrap 5% CI of λ and estimated σ for each pathway. Transcriptome-wide gene network We observed a linear relationship between the length of We inferred the transcriptome-wide co-expression net- the interval and σ , confirming our previous observation work for 9071 genes across the GTEx and BRAINEAC (Additional file 5: Figure S4) that the variability of λ is a studies in the four brain tissues. The network was fully function of σ and independent of the number of genes in characterized by 8.2 × 107 gene-pair correlations. We the pathways; see Additional file 7: Figure S6). observed λ = 0.5 for the transcriptome-wide network across all four tissues. Figure 5 illustrates the preservation matrix between studies across tissues. We observed that Table 1 Agreement measure λ > 0.5 between BRAINEAC and all the preservation statistics in the matrix were similar in GTEx across 4 brain regions in 287 KEGG pathways size, between 0.38 and 0.49, and that the cerebellum and λ CI: 2.5% CI: 97.5% σ Ref Description frontal cortex diagonals were correctly paired between 0.68 0.72 0.75 0.02 hsa05033 Nicotine addiction studies. We noted that the diagonal terms at the hip- pocampus and putamen were the second maxima after 0.67 0.61 0.88 0.04 hsa04720 Long-term potentiation their correlations with the frontal cortex in GTEx. There- 0.58 0.61 0.77 0.04 hsa05206 MicroRNAs in cancer fore, the transcriptome-wide networks cannot clearly dis- 0.55 0.50 0.59 0.01 hsa04080 Neuroactive ligand-receptor entangle the frontal cortex from the hippocampus and interaction putamen, suggesting that a large amount of the 8.2 × 107 0.53 0.42 0.69 0.03 hsa04020 Calcium signaling pathway gene-pair correlations may not be brain-region specific. 0.52 0.45 0.65 0.03 hsa04261 Adrenergic signaling in We finally benchmarked the transcriptome-wide net- cardiomyocytes work obtained in BRAINEAC against GTEx across the 0.51 0.49 0.62 0.02 hsa04912 GnRH signaling pathway four brain regions. In this case, we confirmed that all Cáceres and Gonzalez BMC Genomics (2019) 20:26 Page 7 of 13 Fig. 4 Co-expression structure of the nicotine addiction pathway in the cerebellum (CRB), frontal cortex (FCTX), hippocampus (HIPP) and putamen (PUTM) across the GTEx and BRAINEAC studies. The figure illustrates the strongest correlations between genes only (greater than the last decile) diagonal terms were their column maxima (see Fig. 5) of agreement between studies that measure the effects and therefore benchmark-λ = 1. These results show of varying conditions on a set of units (subjects, co- that the transcriptome-wide network can correctly pair expression pairs, etc). We showed that the measure con- BRAINEAC tissues with those of GTEx. As the bench- formed to the properties of agreement measures and used marking only takes into account the variability between simulations to compare it with Cohen’s κ. Our results tissues of GTEx and not BRAINEAC, the observation: illustrate the large potential of λ, in particular, for studying benchmark-λ = 1, is consistent with the consideration the changes in gene-network structure across numerous that RNA-sequencing (GTEx) may provide better charac- tissues. The measure is particularly useful to determine terization of the transcriptome than microarray technol- the fraction of tissues for which network structures can ogy (BRAINEAC) [19]. be reliably distinguished, and thus informs on the ratio of structural changes that may be tissue-specific. Discussion In our application to co-expression networks in brain, We proposed a new measure, λ, of agreement between we observed that the structure of a number of brain- studies. The motivation of the measure is the assessment related pathways could be reliably matched by tissue Fig. 5 Preservation matrices for the nicotine addiction pathway and the transcriptome-wide network. The elements in the preservation matrices are Fisher’s Z-transformed correlations of the co-expression relationships among genes between studies Cáceres and Gonzalez BMC Genomics (2019) 20:26 Page 8 of 13 Table 2 Agreement measure λ > 0.7 for the benchmarking of [20]. The reliability of other addiction pathways (mor- BRAINEAC against GTEx across 4 brain regions in 287 KEGG phine, amphetamine and alcohol) supports these findings pathways and indicates that structural changes of those pathways λ CI: 2.5% CI: 97.5% σ Ref Description in regions with dopamine projections, which are essential 0.97 0.99 1.00 0.01 hsa04724 Glutamatergic in addiction processes (frontal cortex, hippocampus and synapse putamen), are likely tissue-specific [21]. 0.95 0.97 1.00 0.01 hsa05033 Nicotine addiction The agreement measure λ is based on the preservation 0.94 0.93 1.00 0.01 hsa04360 Axon guidance matrix between tissues. Langfelder and colleagues thor- 0.91 0.93 1.00 0.02 hsa05414 Dilated oughly studied a number of pair-wise preservation statis- cardiomyopathy tics of Cholesterol Biosynthetic Process network across 0.89 0.97 1.00 0.02 hsa04931 Insulin resistance sex and tissue differences in mice [4]. Their approach 0.84 0.88 0.99 0.03 hsa04725 Cholinergic to assessing the reproducibility of the data was to test synapse the preservation of network structure against a reference 0.84 0.82 0.97 0.02 hsa05032 Morphine set of conditions (female/liver), allowing them to identify addiction the condition (male/liver) for which structures are most 0.82 0.78 0.93 0.02 hsa04720 Long-term similar. However, their data show, more generally, that potentiation sex changes are much smaller than tissue changes and, 0.82 0.76 0.89 0.02 hsa05031 Amphetamine therefore, that differences in network structures can reli- addiction ably pair tissues between sex changes. In this case, the 0.81 0.76 0.95 0.02 hsa04912 GnRH signaling inter-study reliability for changes in the network structure pathway across a range of tissues is clearly high, as it would be sum- 0.79 0.70 0.91 0.03 hsa04261 Adrenergic marized by a by high λ. While a pair-wise preservation signaling in cardiomyocytes statistic assesses the reliability of a network’s structure between two conditions [4], the preservation matrix reg- 0.76 0.75 0.78 0.01 hsa04727 GABAergic synapse isters all pair-wise changes of the structure among numer- ous conditions. λ is a measure on the preservation matrix 0.75 0.75 0.75 0.00 hsa04080 Neuroactive ligand-receptor that allows the assessment of agreement between stud- interaction ies to reliably discriminate between conditions, measuring 0.75 0.75 0.75 0.00 hsa04723 Retrograde the probability that the diagonal terms of the matrix are endocannabinoid their row and column maxima. signaling We are unaware of similar measures of agreement 0.74 0.75 0.75 0.00 hsa04020 Calcium signaling on akin preservation matrices that test the changes in pathway the structure of a network across tissues. Other func- 0.74 0.73 0.76 0.01 hsa04728 Dopaminergic tions of the preservation matrix could also be investi- synapse gated. We observed, in particular, that λ is a measure of 0.74 0.76 0.90 0.03 hsa05206 MicroRNAs in inter-observer agreement that is consistently higher than cancer Cohen’s κ. Perfect agreement for κ is exclusively given by 0.73 0.75 0.76 0.01 hsa04015 Rap1 signaling pathway diagonal tables, while perfect agreement for λ is given by maxima diagonal terms in tables where their elements, 0.71 0.67 0.95 0.04 hsa04971 Gastric acid secretion irrespective of their magnitude, are estimated with suffi- cient low variability. λ does not test the level of pair-wise 0.71 0.68 0.80 0.02 hsa05034 Alcoholism preservation but the extent to which structural patterns are correctly classified by conditions. This is an impor- across GTEx and BRAINEAC. We observed that the tant difference between the measures, which allows λ to nicotine addiction pathway presented the highest λ in be used in more general situations where the elements which the frontal cortex, the hippocampus and the of cross-tabulated tables between raters/studies are infer- putamen were reliably distinguished. Interestingly the ences. In particular, we observed that λ is informative for cerebellum, with a limited role in addiction, could not intermediate values, as it can be used to select the frac- be reliably distinguished. Benchmarking of BRAINEAC tion of tissues where network changes are reproducible against GTEx, confirmed that the differences in the nico- and likely to be tissue-specific. tine pathway structure across relevant brain regions were Currently, the validity of a gene or protein network highly reliable, as were those in glutamatergic, choliner- derived from high-throughput data is often benchmarked gic, GABAergic, dopaminergic synapse pathways, whose against networks derived from current knowledge of complex interactions are critical to nicotine addiction specific interactions, given by curated pathways, specific Cáceres and Gonzalez BMC Genomics (2019) 20:26 Page 9 of 13 experiments or even text mining of published articles The problem [22]. This type of confirmatory analysis extracts networks Let us assume that we have two studies that measure the that are a mixture of interactions that have been indi- expression of a set of genes in two different population vidually reported on different tissues. Therefore, while samples, in the same range of tissues and using different validity is investigated, in terms of consistency with pre- experimental setups. For instance, one experiment may vious knowledge, reproducibility of network structure in use RNA-seq and the other a microarray technology. We an independent confirmatory study is not being assessed. may be interested in inferring the co-expression relation- We observed important preservation of network patterns ships between the genes across tissues and determine between two equivalent, independent studies, based on whether the changes in the relationships are consistent different technologies (RNA-seq and microarray), analysis between experiments. The co-expression between two methods and subjects. As such, our results fully test inter- genes in the network is given by their correlation over the observer reproducibility and provide support that inferred subjects’ gene expression levels. Figure 1 illustrates the sit- network structures are biological entities independent of uation where the co-expression between the genes of a numerous sources of variability and heterogeneity in the network is inferred for 6 tissues in two different studies. studies. Previous work has tested the agreement of tran- Our aim is then to propose a measure that tells us the scription measurements between RNA-sequencing and m fraction of tissues that can be reliably distinguished by the icroarray experiments on the same subject samples ([23, 24]). network structure across studies. The measure will then These are important studies to validate experimental tech- inform on the fraction of tissues for which the changes in niques. However, fully testing inter-observer reproducibil- gene relationships may be tissue-specific. ity of network inferences, as we have done, must be based A network A of m genes (nodes) can be represented on independent confirmatory experiments on different by a m × m adjacency matrix between genes, where the population samples, experiments and analyses. Finally, entries are the edges linking a pair of nodes, see Fig. 1. We recent efforts have been made to assess Bayesian network are interested in gene-network modules, such as pathways reproducibility, where intra-study preservation statistics or gene-sets associated with a biological function. There have been proposed [25]. While there is still a need to are numerous structural properties of networks that can apply those metrics to assess inter-study reproducibility, n be used to test the preservation of a network between ote that λ will also be applicable in such a case. two studies. Preservation statistics quantify given aspects of within-network topology that are preserved between Conclusions studies, see Langfelder et al. for a review [4]. Given a We propose a novel inter-study measure of agreement λ to preservation statistics z, we aim to assess the inter-study determine the fraction of conditions for which structural reliability of a network’s structure between studies E1 and changes in a set of observations may be deemed reliable. E2 across numerous conditions. The measure revealed that changes in the co-expression For studies that compute a network across different structure of addiction-related pathways can correctly dis- conditions like tissues (i = 1, ...k), we want to assess tinguish the frontal cortex, hippocampus and putamen the preservation between studies of the network changes but not cerebellum between studies. Reliable tissue dif- across tissues. We thus form the preservation matrix Z ferences in network structure can help to identify tissue- across conditions  where the entries are pair-wise preser- specific pathway-biology and increase the reproducibility vation statistics z Ai , Aj = zij of network A between of network inferences. More generally, our agreement tissues i and j in E1 and E2, respectively. For three tissues measure can be used in a set of independent studies that we thus have measure how the same group of units changes across numerous conditions. ⎛ ⎞ z11 z12 z13 Materials and methods Z = ⎝ z21 z22 z23 ⎠ (3) We propose an inter-study measure of agreement to dis- z31 z32 z33 criminate conditions. While the measure can be applied to different research studies, such as those with factorial We would then like to have a measure of inter-study agree- designs, we illustrate how its need arises from an exam- ment from Z, which can tell how diagonal Z is or, more ple in current functional genomic research. We consider explicitly, the extent to which i can be correctly paired two studies that measure the same units (subjects, gene- with i , for all i. pair correlations, etc.) across multiple conditions. We aim to construct a measure that informs on the fraction of A solution conditions that can be reliably distinguished by structural In Z, we can pair by similarity two network structures cor- changes in the observations. responding to condition i in E1 and condition j in E2 if Cáceres and Gonzalez BMC Genomics (2019) 20:26 Page 10 of 13 the element zij is maximum within row i and column j. As The probability that the diagonal term zii is the maximum such, we are interested in measuring the extent to which in the column i follows a similar form. two structures can be paired when the conditions match Our agreement measure then follows from the fraction between studies, that is when j = i . We then propose R of diagonal terms in Z that are their row and column to measure the probability that the diagonal terms in the maxima. The fraction R distributes as the fraction of suc- preservation matrix are their row and column maxima. cesses in k Bernoulli trials each of which has its own For Z in equation 3, we thus compute probability pii ; that is, as a Poisson binomial distribution for the success fraction with mean and variance • p11 = Pr(z11 > z12 , z13 , z21 , z31 ), 1 • p22 = Pr(z22 > z21 , z23 , z12 , z32 ) and λ = E[R] = pii (9) • p22 = Pr(z33 > z31 , z32 , z13 , z23 ), k i 1  where pii is the probability that the network structures in σ 2 = Var[R] = 2 pii (1 − pii ). (10) k tissue i between studies are correctly paired by maximum i similarity. Taken the preservation statistic as a similar- The agreement measure λ is then the expected frac- ity measure, we can assume that the similarity of one tion of conditions with the correct matching of network structure (i) to another (i ) is independent of its simi- structures between studies. In the case that E1 is the larity to a third (j ). Therefore, the probabilities pii can benchmark for experiment E2, then one is interested in be computed as the product of the individual pair-wise testing whether the diagonal terms are the maxima of their probabilities columns only, generalizing the concepts of sensitivity and specificity. In this case λ can be computed by simply set- k ting Pr(zii > zji ) = 1 for j = 1, ...k. Note also that it is pii = Pr(zii > zij ) ∗ Pr(zii > zji ), (4) straightforward to generalize the measure for more than j=1 two studies by expanding the products in Eq. (4). where the first factor is the probability that zii is the max- imum over row i of Z, and the the second factor is the Reliability properties of λ probability that zii is the maximum over column i . We Suitable reliability measures satisfy three basic properties: can write down an explicit form for the pair-wise prob- i) their values range from 0 to 1; ii) they tend to 0 when no abilities, further assuming that the preservation statistic agreement is expected and tend to 1 when a full agreement distributes normally is expected and iii) they account for random agreement.   Regarding property i), λ clearly ranges from 0 to 1, as zij ∼ N μij , σij2 . (5) it is given by the fraction between the sum of k proba- bilities pii and k. To show properties ii) and iii), we look Following the independence of structural similarities at the explicit model for pii given by the substitution of between experiments in different tissues, we then have Eq. (7) in (4) that the pair-wise probability Pr(zii > zij ) can be derived ⎡ ⎛ ⎞⎤ ⎡ ⎛ ⎞⎤ from the suitable integration of the joint distribution of zii k 1⎢ ⎜ 1 μij − μii ⎟⎥ 1 ⎢ ⎜ 1 μji − μii ⎟⎥ and zij pii = ⎣1 − erf ⎝ √  ⎠⎦ ∗ ⎣1 − erf ⎝ √  ⎠⎦ . 2 2 σ2 + σ2 2 2 σ2 + σ2 j=1 ij ii ji ii ∞ ∞   Pr(zii > zij ) = N μii , σii2 N μij , σij2 dzii dzij . (11) −∞ zij erf(x) is the error function which, importantly, defines (6) the Heaviside step function in the limit: H(x) = 1/2 ∗ ⎡ ⎛ ⎞⎤ limt→0 erf(−x/t) [26]. Therefore, when {σii , σji , σij } → 0, 1⎢ ⎜ 1 μij − μii ⎟⎥ pii tends to = ⎣1 − erf ⎝ √  ⎠⎦ . (7) k 2 2 σ 2 + σ 2 ij ii pii = H(μii > μij ) ∗ H(μii > μji ) (12) j=1 where erf is the error function [26]. Therefore, we have  that the probability that the diagonal term zii is the maxi- 1, if μii > {μij , μji }∀j = 1...k = (13) mum in the row i is 0, otherwise. ⎛ ⎡ ⎤⎞ That is, when there is no variability in zij , pii is 1 if zii k k 1 ⎜ ⎢ 1 μij − μii ⎥⎟ is maximum across rows and columns and 0 otherwise. Pr(zii > zij ) = ⎝1 − erf ⎣ √  ⎦⎠ . Therefore, regarding property ii), null and perfect relia- 2 2 σ 2 + σ 2 j=1 j=1 ij ii bilities are clearly expected with null variability. For such (8) cases all pii are either 0 or 1 according to Eq. (13). In the Cáceres and Gonzalez BMC Genomics (2019) 20:26 Page 11 of 13 case of null reliability none of the μii is maximum, all pii For for large ni the binomial distribution tends to a nor- are 0 and λ = 0. Whereas, in the case of perfect reliabil- mal distribution. Therefore, the values in Eqs. 17 and 18 ity all μii are maxima, all pii are 1 and λ = 1. Regarding can be used in Eq. 7. With a similar computation for property iii), we can see that under null variability, λ tends the column elements, the measure λ can be obtained to the observed fraction of conditions r that are correctly for a table in the form (14) and compared with the paired between experiments. The fraction r, which is a value of κ for varying values of the total fraction of realization of R, can be considered as a reliability measure agreement P0 . but does not account for random agreement. By contrast, As κ is an agreement measure that corrects P0 for ran- λ, as an expected value of R under a probability distribu- dom agreement, we compared λ with the observed frac- tion, incorporates a variability that accounts for correct tion of diagonal elements that are their row and column random pairing. maxima r, explicitly given by Comparison between measures of agreement    Let us take two raters who observe l units in k categories. 1, if nii = max nij , nji j Then, for instance, if k = 3 and categories are labeled a, ri = (19) b and c, the cross-tabulation matrix of observaions takes a 0, otherwise similar form of expression (3), 1 r= ri . (20) k a b c i a naa nab nac (14) Simulations b nba nbb nbc We performed a series of simulations to study the proper- c nca ncb ncc , ties of λ with respect to κ and r. Simulations were obtained for three total number of conditions/tissues k = (5, 7, 15), where ni,j is the number of units observed in cate- and three fifferent forms for the marginal frequencies gory i and j by the first and second rater respectively,  P1 (i) and P2 (i) (i = 1, ...k) across studies 1 and 2 and ki,j=1 nij = l. Agreement is typically measured by Cohen’s κ • Scenario 1 (equiprobable): P1 (i) = P2 (i) = k1 , ∀i k   P(i, i) − ki=1 P1 (i)P2 (i) • Scenario 2: P1 (i) = P2 (i) = 1i / kj=1 1j κ = i=1 k , (15) • Scenario 3 (the least equiprobable): 1 − i=1 P1 (i)P2 (i)  P1 (i) = P2 (i) = i12 / kj=1 j12 where P(i, i) = nii /l is the observed frequency of units that were measured in category i by both raters and Pd (i) We set the number of observations to l = 500. For is the frequencyof units in i observed by rater d (d = 1, 2). each scenario, we simulated 50 cases of perfect agree- The sum P0 = i P(i, i) is the total fraction of agreement: ment tables (P0 = 1), i.e. diagonal matrices, and 50 cases the proportion of observations that falls in the diagonal, of perfect disagreement (P0 = 0); those are tables with which does not account for random agreement. Cohen’s κ zeros on the diagonal terms except for the cell of maxi- measures the fraction of discordant observations expected mum joint probability. For each case, we permuted a pair by chance that are actually observed in agreement. of observations 3000 times, such that the original marginal While λ can be applied to a wider range of reliability frequencies were conserved. After each set of 30 permu- studies than κ, we compared the two measures in cases tations, we computed the four agreement measures. This where both of them can be computed. Note that a matrix procedure allowed the assessment of 10,000 simulations, Z in Eq. (3) can be computed from the cross-tabulation in each scenario and tissue level, covering the whole agree- table in expression (14). Given that in row i, in expres- ment interval of P0 . The simulation scheme is shown in sion (14), the number of observed units is ni = P1 (i) ∗ l, Fig. 2 and further details are given in Additional file 1: we can then assume that nij is one draw of a binomial Table S1. process We used R.3.30 and the package psych to perform calculations and compute Cohen’s κ. nij ∼ Binomial(ni , θij ) (16) Gene expression data where θ̂ij = nij /ni , and the mean and variance of the mean We downloaded expression data from the GTEx project are given by obtained from RNA-seq [27]. Reads per kilobase per mil- lion mapped reads (RPKM) of version 6 were obtained μij = nij (17) for all brain tissues. Covariates for each tissue were also σij2 = nij (1 − nij /ni ), (18) downloaded. Cáceres and Gonzalez BMC Genomics (2019) 20:26 Page 12 of 13 We also downloaded the brain expression data of the Pearson’s correlation between the vectorized matrices of BRAINEAC project [28] obtained from winsorized val- pathway A ues of exon array data (Affymetrix Human Exon 1.0 ST   array). Downloaded data had been previously normalized cor.cor Ai , Aj and corrected for batch effects.    We identified four brain tissues common in both data- = cor vectorizeMatrix(Ai ), vectorizeMatrix Aj sets and for which GTEx had covariates’ information. (22) Those were cerebellum (CRBL) with 125 individuals in GTEx and 130 in BRAINEAC, frontal cortex (FCTX) with between tissues i and j, across experiments E1 and E2, 108 and 135 individuals, (HIPP) hippocampus with 94 and respectively. The elements zij of the preservation matrix 130 individuals, and putamen (PUTM) with 82 and 135 Z for each pathway A were obtained  by the Fisher’s Z- individuals, respectively. transformation of cor.cor Ai , Aj . Between the two studies, we mapped 9071 common genes to compute all gene-pair correlations of expression Additional files levels. Additional file 1: Table S1. Simulation details. (52.6 KB) Co-expression networks Additional file 2: Figure S1. Relationship between λ and Cohen’s κ For each tissue in GTEx, we downloaded its correspond- across all simulation scenarios. (406 KB) ing table of covariates. These included 20 variables, 3 Additional file 3: Figure S2. Relationship between λ and σ across all of which are genome-wide principal components that simulated scenarios. (401 KB) inform on ancestry, 15 surrogate transcriptomic variables Additional file 4: Figure S3. Left) Relationship between κ and its variance for the simulated scenario with 5 conditions and a single subject, that inform on batch effects, one variable containing gen- corresponding to a single draw of the binomial process in Eq. (16). The der and another platform. As expression levels for GTEx functional dependence of κ variance with κ is unchanged but reduced in are derived from count data, we performed Spearman’s magnitude with increasing N. Right) Relationship between λ and σ in the same scenario. (120 KB) partial correlations of the expression levels between each Additional file 5: Figure S4. Bootstrap 5% confidence intervals for λ for a gene-pair, among 9071 genes. We used the par.cor simulation with 5 conditions and a varying number of individuals. We see function of the ppcor R package to test the rank cor- that the length of the CI is proportional to σ and independent of the relations adjusting for covariates. A transcriptome-wide number of subjects. (7.61 KB) adjacency matrix was constructed for each tissue, where Additional file 6: Figure S5. Power estimations for true agreement r = 3/5 when the true correlations between matching tissues (diagonal the entries, corresponding to the edges of the genome- terms in the preservation matrix) are fixed at cor = 0.5 and all other wide network, were the absolute value of all gene-pair background correlations c0, given by all other elements in the preservation correlations. matrix, varied from 0 to 0.5. The figure shows the differences between fixed and varying correlations (corDiff = cor − cor0) by a different number of For the BRAINEAC data, similar adjacency matrices per genes N in the network. Networks were simulated as random variables tissue were obtained. In this case, however, the down- from multivariate distributions. (4.90 KB) loaded gene expression values had already been normal- Additional file 7: Figure S6. Bootstrap 5% confidence intervals for λ, ized and corrected for batch effects. As gene expression between GTEx and BRAINEAC, for 287 KEGG pathways. We see that the length of the CI is proportional to σ , regardless of the different number of levels are close to normality, we performed a Pearson’s genes of each pathway. (13.4 KB) correlation between all gene-pairs and used these values Additional file 8: Tables S2-S3. Agreement measure λ between GTEx to compute the adjacency matrices. and BRAINEAC for 287 KEGG pathways across four brain regions. (56.0 KB) We computed λ to test the inter-study agreement of Additional file 9: Computer Code. Entire computer code in R to the co-expression networks of 287 pathways, downloaded reproduce the results discussed in the manuscript. (12.1 KB) from https://www.genome.jp/kegg/. For each pathway, we Abbreviations extracted a total of eight adjacency matrices (4 tissues Binomial(n, p): Binomial distribution of n trials with success probability p; × 2 studies) whose elements were the absolute value of BRAINEAC study: The brain eQTL almanac study; CRBL: Cerebellum; erf(x): Error the expression correlations between all gene pairs within function of x; FCTX: Frontal cortex; GTEx study: Genotype-tissue expression study; HIPP: Hippocampus; KEGG: The Kyoto encyclopedia of genes and the pathway. We computed the pathways’ vectorized adja- genomes; netA: Network A; PUTM: Putamen; RPKM: Reads per kilobase per cency matrices, which are the vectors with non-redundant million mapped reads; RNA-seq: RNA sequencing; N(μ, σ 2 ): Normal components distribution with mean μ and variance σ 2 vectorizeMatrix(A) Acknowledgments (21) We thank the reviewers of the manuscript for their contribution. = (a2,1 , a3,1 , a3,2 , a4,1 , a4,2 , a4,3 , ..., an,1 , am,m−1 ), where m in the number of genes in the pathway A. Funding This work was partially supported by Red Española de Supercomputación As a preservation statistic between studies, we used the (BCV-2016-3-0002) and by Ministerio de Economía e Innovación (Spain) adjacency correlation (cor.cor in [4]), based on the (MTM2015-68140-R). Cáceres and Gonzalez BMC Genomics (2019) 20:26 Page 13 of 13 Availability of data and materials 15. Kyoto encyclopedia of genes and genomes. http://www.genome.jp/ The entire code to fully reproduce this work is given in Additional file 9: kegg. Accessed 12 Apr 2018. Computer Code and all required pre and post-processed data has been 16. Fleiss JL, Cohen J, Everitt B. Large sample standard errors of kappa and deposited in Figshare and assigned the identifier https://figshare.com/ weighted kappa. Psychol Bull. 1969;72(5):323. account/home#/projects/20065, see README file within computer code. The 17. Krichevsky AM, King KS, Donahue CP, Khrapko K, Kosik KS. A microrna function lambda.agreement has been implemented in an R package array reveals extensive regulation of micrornas during brain available in https://github.com/isglobal-brge/lambda. development. Rna. 2003;9(10):1274–81. 18. O’Dell TJ, Connor SA, Guglietta R, Nguyen PV. β-adrenergic receptor Authors’ contributions signaling and modulation of long-term potentiation in the mammalian AC designed the study, analyzed the data and wrote the manuscript, JRG hippocampus. Learn Mem. 2015;22(9):461–71. supervised the work and revised the manuscript. All authors have read and 19. Zhao S, Fung-Leung W-P, Bittner A, Ngo K, Liu X. Comparison of rna-seq approved the manuscript. and microarray in transcriptome profiling of activated t cells. PLoS ONE. 2014;9(1):78644. Ethics approval and consent to participate 20. Markou A. Phil Trans R Soc B Biol Sci. 2008;363(1507):3159–68. Not Applicable. 21. Volkow ND, Wang G-J, Fowler JS, Tomasi D, Telang F, Baler R. Addiction: decreased reward sensitivity and increased expectation sensitivity Consent for publication conspire to overwhelm the brain’s control circuit. Bioessays. 2010;32(9): Not Applicable. 748–55. 22. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Competing interests Huerta-Cepas J, Simonovic M, Roth A, Santos A, Tsafou KP, et al. String The authors declare that they have no competing interests. v10: protein–protein interaction networks, integrated over the tree of life. Nucl Acids Res. 2014;43:D447–52. Publisher’s Note 23. Trost B, Moir CA, Gillespie ZE, Kusalik A, Mitchell JA, Eskiw CH. Springer Nature remains neutral with regard to jurisdictional claims in Concordance between rna-sequencing data and dna microarray data in published maps and institutional affiliations. transcriptome analysis of proliferative and quiescent fibroblasts. R Soc Open Sci. 2015;2(9):150402. Author details 24. Guo Y, Sheng Q, Li J, Ye F, Samuels DC, Shyr Y. Large scale comparison 1 ISGlobal, Doctor Aiguader, 08003 Barcelona, Spain. 2 Centro de Investigación of gene expression levels by microarrays and rnaseq using tcga data. PloS Biomédica en Red en Epidemiología y Salud Pública (CIBERESP), Barcelona, ONE. 2013;8(8):71462. Spain. 3 Department of Mathematics, Universitat Autònoma de Barcelona, 25. Cohain A, Divaraniya AA, Zhu K, Scarpa JR, Kasarskis A, Zhu J, Chang R, Bellaterra, Barcelona 08193 Spain. Dudley JT, Schadt EE. Exploring the reproducibility of probabilistic causal molecular network models. In: PACIFIC SYMPOSIUM ON BIOCOMPUTING Received: 27 April 2018 Accepted: 29 November 2018 2017. Singapore: World Scientific; 2017. p. 120–31. 26. Weisstein ew. heaviside step function. in mathworld–a wolfram web resource. 2018. http://mathworld.wolfram.com/heavisidestepfunction. html. Accessed 12 Apr 2018. References 27. Genotype tissue-expression project. http://www.gtexportal.org. Accessed 1. Baker M. 1500 scientists lift the lid on reproducibility. Nature. 12 Apr 2018. 2016;533(7604):452–4. 28. Braineac project. http://www.braineac.org/. Accessed 12 Apr 2018. 2. Wasserstein RL, Lazar NA. The asa’s statement on p-values: context, process, and purpose. Am Stat. 2016;70(2):129–33. 3. Mogil JS, Macleod MR. No publication without confirmation. Nature. 2017;542(7642):409–11. 4. Langfelder P, Luo R, Oldham MC, Horvath S. Is my network module preserved and reproducible?. PLoS Comput Biol. 2011;7(1):1001057. 5. Oldham MC, Konopka G, Iwamoto K, Langfelder P, Kato T, Horvath S, Geschwind DH. Functional organization of the transcriptome in human brain. Nat Neurosci. 2008;11(11):1271. 6. Cai C, Langfelder P, Fuller TF, Oldham MC, Luo R, van den Berg LH, Ophoff RA, Horvath S. Is human blood a good surrogate for brain tissue in transcriptional studies?. BMC Genomics. 2010;11(1):589. 7. Chen C, Cheng L, Grennan K, Pibiri F, Zhang C, Badner JA, Gershon ES, Liu C. Two gene co-expression modules differentiate psychotics and controls. Mol Psychiatry. 2013;18(12):1308. 8. Hallgren KA. Computing inter-rater reliability for observational data: an overview and tutorial. Tutor Quant Methods Psychol. 2012;8(1):23. 9. Cohen J. A coefficient of agreement for nominal scales. Educ Psychol Meas. 1960;20(1):37–46. 10. Banerjee M, Capozzoli M, McSweeney L, Sinha D. Beyond kappa: A review of interrater agreement measures. Can J Stat. 1999;27(1):3–23. 11. Shrout PE, Fleiss JL. Intraclass correlations: uses in assessing rater reliability. Psychol Bull. 1979;86(2):420. 12. Gwet KL. Handbook of Inter-rater Reliability: The Definitive Guide to Measuring the Extent of Agreement Among Raters, 4th edn. Gaithersburg, USA: Advanced Analytics, LLC; 2014. 13. Lonsdale J, Thomas J, Salvatore M, Phillips R, Lo E, Shad S, Hasz R, Walters G, Garcia F, Young N, et al. The genotype-tissue expression (gtex) project. Nat Genet. 2013;45(6):580–5. 14. Trabzuni D, Ryten M, Walker R, Smith C, Imran S, Ramasamy A, Weale ME, Hardy J. Quality control parameters on a large dataset of regionally dissected human control brains for whole genome expression studies. J Neurochem. 2011;119(2):275–82.