extfieldsofstudy:
Biology
Medicine
Computer Science

text:
Modelling and visualizing fine-scale linkage disequilibrium structure

Background Detailed study of genetic variation at the population level in humans and other species is now possible due to the availability of large sets of single nucleotide polymorphism data. Alleles at two or more loci are said to be in linkage disequilibrium (LD) when they are correlated or statistically dependent. Current efforts to understand the genetic basis of complex phenotypes are based on the existence of such associations, making study of the extent and distribution of linkage disequilibrium central to this endeavour. The objective of this paper is to develop methods to study fine-scale patterns of allelic association using probabilistic graphical models. Results An efficient, linear-time forward-backward algorithm is developed to estimate chromosome-wide LD models by optimizing a penalized likelihood criterion, and a convenient way to display these models is described. To illustrate the methods they are applied to data obtained by genotyping 8341 pigs. It is found that roughly 20% of the porcine genome exhibits complex LD patterns, forming islands of relatively high genetic diversity. Conclusions The proposed algorithm is efficient and makes it feasible to estimate and visualize chromosome-wide LD models on a routine basis.

Background
Alleles at two loci are said to be in linkage disequilibrium (LD) when they are correlated or statistically dependent. The term refers to the idea that in a large homogeneous population subject to random mating, recombination between two loci will cause any initial association between them to vanish over time. In observed data, however, non-zero allelic associations are pervasive, particularly at short distances, but also at long distances and even between chromosomes. These associations arise in a complex interplay between processes such as mutation, selection, genetic drift and population admixture, and are broken down by recombination. The patterns of association are of interest, partly because they underpin the relation of genotype to phenotype at the population level, and partly because they reflect population history.
Patterns of LD may be represented in different ways. A common method is to display pairwise measures of LD as triangular heatmaps [1,2]: in these displays, LD blocks (genomic intervals within which all loci are in high LD) stand out clearly. Early work in the HapMap project led http://www.biomedcentral.com/1471-2105/ 14/179 (1, 1, 1), (1, 2, 1), (2,1,1) and (1,2,2). Observe that in general s 1 and s 3 are marginally associated (are in LD), but in the subpopulations corresponding to s 2 = 1 and s 2 = 2 they are unassociated: in other words, they are conditionally independent given s 2 . More complex mutation histories give rise to more complex patterns of conditional independences that can be represented as graphical models [8].
Other authors have used graphical models for the joint distributions of allele frequencies. Usually, in highdimensional applications, attention is restricted to a tractable subclass, the decomposable graphical models [10]. In the first use of decomposable models in this context [11], models were selected using a greedy algorithm based on significance tests. In [8,12] methods and programs for selecting decomposable graphical models using Monte Carlo Markov Chain (MCMC) sampling were described. These methods are computationally feasible for modest numbers of markers (say, several hundreds), but not for modern SNP arrays with hundreds of thousands of SNPs per chromosome. To improve efficiency, the search space may be restricted to graphical models whose dependence graphs are interval graphs [13,14]. These are graphs for which each vertex may be associated with an interval of the real line such that two vertices are connected by an edge if and only if their intervals overlap. In this context the ordering of SNPs along the real line is their physical ordering along the chromosome. MCMC sampling from this model class may be performed more efficiently [13,14]. This work was extended in [15] to a more general subclass of decomposable models, namely those in which distant marker pairs (i.e., with more than a given number of intervening markers) are conditionally independent given the intervening markers.
In an alternative approach [16][17][18] latent mixtures of forests have been applied, in order to accommodate short-, medium-and long-range LD patterns. Also directed graphs (Bayesian networks) have been applied, selecting edges and their directions using causal discovery algorithms [19]. There are close links between decomposable models and Bayesian networks ( [10], Sect. 4.5.1).
A rather different approach to modelling the joint distribution of allele frequencies [20,21] is implemented in the software package BEAGLE [22], which is widely used to process data from SNP arrays. The approach is based on a class of models arising in the machine learning literature called acyclic probabilistic finite automata (APFA) [23]. These are related to time-variant variable length Markov chains. For phase estimation and imputation BEAGLE uses an iterative scheme analogous to the EM algorithm, alternating between sampling from a haplotype-level model given the observed genotype data (the E-step) and selecting a haplotype-level model given the samples (the M-step). A similar computational scheme for decomposable graphical models has been described and implemented in the FitGMLD program [15].
Characterization of genetic variation at the population level is of fundamental importance to understanding how phenotypes relate to genotypes. Some specific uses to which joint models for allele frequencies have been put include 1. Insight into the population history of different genomic intervals. Under simplifying assumptions, the ancestral history of a short genomic interval can be reconstructed from a decomposable graphical model for the SNPs in the interval [8]. 2. Quality control of genome assembly, in that some motifs may suggest errors in SNP positioning. 3. Phase estimation and imputation as described above [15,21]. 4. Derivation of more informative covariates involving multiple loci to plug into genomic prediction models [20]. 5. Use in pedigree simulation to model LD in founders, for example in connection with gene drop simulation [14] and assessment of SNP streak statistics [24].
In this paper decomposable graphical models are used to model fine-scale, local LD patterns. By local is meant that the proposed methods are designed to capture shortrange associations between loci, but not long range ones. The same is true of other approaches [13][14][15] mentioned above. Here, an efficient, linear-time algorithm is developed to select a model using a penalized likelihood criterion, and it is shown how such a model may conveniently be displayed, allowing fine-scale LD structure to be visualized.

Graphs and graphical models
The following notation and terminology is mainly based on [9]. A graph is defined as a pair G = (V , E), where V is a set of vertices or nodes and E is a set of edges. Each   duroc  1  4249  2863  4762  0  0  4340  0  0  35084  1153  37  37  15   duroc  2  4249  2113  3676  0  3  3339  2  2  18450  920  26  27  10   duroc  3  4249  1535  2873  1  13  2643  4  15  8695  679  20  20  8   duroc  4  4249  2211  4392  0  0  4012  0  0  24005  1007  35  36  11   duroc  5  4249  1255  2357  7  2  2166  8  2  6030  559  16  16  10   duroc  6  4249  1937  3358  0  0  3069  0  0  15396  778  24   The table shows the performance of the algorithms applied to the pig data. N and p denote the numbers of observations in the data and number of SNPs available after filtering. The edge sets found in steps 1 and 2 of the standard algorithm are denoted E s and E sp , respectively: for the fast algorithm they are denoted E f and E fp . The numbers of undershoots and overshoots, i.e. |E s \ E f | and |E f \ E s | for step 1, and |E sp \ E fp | and |E fp \ E sp | for step 2, are also shown. The corresponding running times in seconds are denoted t s , t f , t sp and t fp , respectively. The backward step is much faster than the forward step in both cases (t sp and t fp ). Overall, the fast algorithm is 23.4 times faster than the standard algorithm, and its inaccuracy, assessed as The last column shows the maximum width of G = (V, E s ). The computations were run under Redhat Fedora 10 Linux on a Intel i7 four-core 2.93GHz machine with 48 GB RAM.
edge is associated with a pair of nodes, its endpoints. Here only undirected graphs are considered, that is, with graphs undirected edges only. Two vertices α and β are said to be adjacent, written α ∼ β, if there is an edge between them. The neighbours of a vertex is the set of nodes that are adjacent to it. A subset A ⊆ V is complete if all vertex pairs in A are connected by an edge. A clique is a maximal complete subset, that is to say, a complete subset that is not contained in a larger complete subset.
A path (of length n) between vertices α and β in an undirected graph is a set of vertices α = α 0 , α 1 , . . . , α n = β where α i−1 ∼ α i for i = 1, . . . , n. If a path α = α 0 , α 1 , . . . , α n = β has α = β then the path is said to be a cycle of length n. If a cycle α = α 0 , α 1 , . . . , α n = α has adjacent elements α i ∼ α j with j ∈ {i − 1, i + 1} then it is said to have a chord. If it has no chords it is said to be chordless. A graph with no chordless cycles of length ≥ 4 is called triangulated or chordal.
A subset D ⊂ V in an undirected graph is said to sepa- The boundary bd(A) of a vertex set A ⊆ V is the set of vertices adjacent to a vertex in A but not in A, that is Let G = (V , E) be an undirected graph with cliques C 1 , . . . C k . Consider a joint density f () of the variables in V. If this admits a factorization of the form If all the densities in a model factorize according to G, then the model is said to be G-Markov. When this is true G encodes the conditional independence structure of the model, through the following result (the global Markov property): whenever sets A and B are separated by a set C in G, A and B are conditionally independent given C under the model. This is written as A ⊥ ⊥ B | C. A decomposable graphical model is one whose dependence graph is triangulated.

Selecting graphical models for chromosome-wide LD
Suppose that N observations of p SNPs from the same chromosome are sampled from some population. The variable set is written V = (v 1 , . . . v p ), and it is assumed that these are ordered by physical position on the chromosome. The variables may be either observed genotypes or inferred haplotypes, if these have been imputed: the former are trichotomous and the latter binary. Here an algorithm to use these data to select a graphical model for the distribution of V is described. It is based on a penalized likelihood criterion is the number of free parameters under G, and α is a penalization constant. For the AIC, α = 2 and for the BIC, α = log(N). The latter penalizes complex models more heavily and so selects simpler models. Under suitable regularity conditions, the BIC is consistent in the sense that for large N it will select the simplest model consistent with the data ( [25], Sect. 2.6). http://www.biomedcentral.com/1471-2105/14/179 A technical but nonetheless important point is that calculation of the correct model dimension is not straightforward for high-dimensional models, since not all parameters may be estimable. A simple expression exists however for the difference in model dimension between two decomposable models that differ by a single edge ( [10], pp. [37][38][39][40]. This is useful when the search space is restricted to decomposable models and the search algorithm only involves comparison of models differing by single edges (as here).
A forward-backward approach to estimate the graphical model for V is used. The first (forward) step is based on a greedy forward search algorithm. To take account of the physical ordering of the SNPs, the algorithm starts from a model G 0 = (V , E 0 ) where E 0 is the set of edges between physically adjacent SNPs: this model is called the skeleton. To seek the minimum BIC (or AIC) model, the algorithm repeatedly adds the edge associated with the greatest reduction in BIC (AIC): only edges whose inclusion results in a decomposable model are considered. The process continues until no more edges improve the criterion. The search space in this step is the class of decomposable models that include the skeleton. Note that this algorithm -as with almost all greedy algorithms -is not guaranteed to find the global optimum.
There are several advantages to initially constraining the search space to include the skeleton. An unconstrained search starting off from the null model (with no edges) would not take the physical ordering into account. Since the graphs considered are consistent with this ordering,  In (b), a corresponding plot of median entropy versus height is shown. The entropy was calculated for each interval of 7 adjacent SNPs using the phased data.
they are conveniently displayed as LD maps, as illustrated below. Because decomposable models contain no chordless cycles of length ≥ 4, two distal loci cannot be linked by an edge unless there are sufficiently many intermediate (chordal) edges to prevent the occurrence of such cycles. Thus the algorithm proceeds by filling out around the skeleton by incrementally adding edges. The effect of this is that only local LD is modelled. Since both association and (implicitly) proximity are required, the edges included are more likely to be for real (that is, less likely to be due to chance), so in this sense the model selection process is more reliable. In addition, restricting the search space in this way improves computational efficiency. The linear-time algorithm described in the following section also builds on the physical vertex ordering.
In the second (backward) step the graph found in the first step is pruned, again using a greedy algorithm that seeks to minimize the BIC (AIC) criterion by removing edges, without requiring that the skeletal edges are retained. Keeping these in the model would be appropriate if adjacent SNPs are in high LD, but this is not always the low complexity interval  It is seen that the interval shown in (a) has low haplotype diversity, dominated by a single common haplotype, whereas that shown in (b) has high haplotype diversity, with one or two relatively common haplotypes.
case. For example, there may be recombination hotspots, or genome assembly errors may have led to errors in positioning the SNPs. The graphs may be used to examine whether this has occurred. To display the resulting graph, a suitable graph layout is needed. After some experimentation, one such was obtained by modifying the horizontal dot layout in Rgraphviz [26], by exponentially smoothing the ycoordinates and replacing the x-coordinates with the SNP-number. In this way the y-coordinates are chosen so as to clarify the graphical structure: a consequence of this is that nodes with large absolute y-coordinates tend to signal genomic regions of high structural complexity. For some purposes it is helpful to use more objective measures of complexity, and two such measures are used here. Consider a chromosome with p SNPs. For each i = 1, . . . p−1, the height h i of the interval between SNPs i and i+1 is defined to be the number of edges of the form (j, k), where j ≤ i and i + 1 ≤ k, and the width w i is defined to be the maximum value of |k − j| of such edges. Note that by the global Markov property, when h i = 0 or equivalently As a measure of haplotype diversity, the entropy [27] is used here. Suppose that in some genomic interval there are k distinct haplotypes with relative frequencies f 1 , . . . f k . Then the entropy is defined as H = − i=1...k f i log f i . It is large when there are many distinct haplotypes that are equally frequent.
A useful way to display marker data over short genomic intervals is the sample tree [23]. This summarizes a set of discrete longitudinal data of the form A rooted tree is constructed in which the non-root nodes represent partial outcomes of the form (x 1 , . . . x k ) for k ≤ q present in the sample data. Edges may be coloured with the marker value and labelled with the corresponding sample counts, or drawn with width proportional to the sample counts.

A fast selection algorithm
For high-dimensional models, the first step of the algorithm described above can be computationally demanding. To address this a much faster algorithm for the same task is now described. This involves combining a series of overlapping marginal models of lower dimension. Call the block length (model dimension) L and the overlap K.
Thus the first block is . . v 2L−K } and so on. In the applications described here L = 100 and K = 20 are used. Suppose that the true model for V = (v 1 , . . . v p ) is G, and thatĜ is the estimate of G obtained by performing the first step of the algorithm described in the previous section. The goal is to construct an approximationG toĜ by combining modelsĜ i = (V i ,Ê i ) obtained by applying that same algorithm to blocks V i for i = 1, 2 . . . .
A way to construct a modelG 12 for V 1 ∪V 2 by combininĝ G 1 andĜ 2 is now described. Let m * = max{w : The rationale for this is that marginal models may include spurious associations on the boundaries. For example, let G 1 and G 2 be the subgraphs of G induced by V 1 and V 2 , respectively. Then the marginal distribution of V 2 will not in general be G 2 -Markov, but it will be G * 2 -Markov for a graph G * 2 derived from G 2 by completing the boundary of each connected component of G V \V 2 in G [28]. SoĜ 2 will tend to contain edges not in G 2 . To estimate the boundary of V 1 \V 2 in G,Ĝ 1 is used: the boundary is estimated to be contained in the set L − K + 1, . . . , m * . Hence by only including edges (v, w) ∈Ê 2 with w > m * edges in this boundary are omitted. Similarly, the boundary of V 2 \ V 1 in G 1 may contain spurious associations, soÊ 0 1 may contain unnecessary edges. To avoid this the overlap length K is chosen to be sufficiently large so that m * < min{v : ∃ (v, w) ∈ E 2 with w ≥ L + 1}. If the maximum width was known in advance this inequality could be ensured by setting K to be twice the maximum width, but so large a value appears to be unnecessary in practice.
The algorithm proceeds in the obvious fashion com-biningG 12 withĜ 3 and so on. Since the construction may result in chordless cycles it is necessary to triangulate the resulting graphG. This can be done using the maximum cardinality search algorithm [29] which can be implemented in linear time.
Assuming that the time for deriving the estimatesĜ i for fixed L is bounded, the algorithm described here is by construction linear in p. But the proportionality constant can http://www.biomedcentral.com/1471-2105/14/179 be expected to be depend on various factors, such as the density of the minimum AIC/BIC model.

Implementation
The methods are implemented in a set of R functions which are provided as Additional files 1 and 2 to the paper. The selection algorithms build on the forward search algorithm implemented in the stepw function in the gRapHD package [30]. Currently this function requires that there are no missing data. The backward search algorithm was developed by the author in the course of preparing this paper, and has subsequently been implemented in the gRapHD package. The functions to perform the selection algorithms, and others to work with and display the graphs, build on the following packages: igraph [31], graph [32], Rgraphviz [26] and gRapHD. The triangular heatmaps were produced using the LDheatmap package [33]. Simulation from a fitted decomposable model was carried out using the gRain package [34]. The package rJPSGCS provides an interface between R and the Fit-GMLD program which was used to perform the algorithm of [15].

Results
To illustrate the methods described above, they were applied to SNP data obtained from three commercial http://www.biomedcentral.com/1471-2105/14/179 pig breeds. In all 4239, 1979 and 2123 pigs of the Duroc, Landrace and Yorkshire breeds were genotyped using the Illumina Porcine SNP60 BeadChip [35]. After preprocessing on the basis of call rate, minimal allele frequency and other quality criteria, missing values were imputed using BEAGLE [22]. Using the methods described above, decomposable graphical models were selected for each chromosome and breed, using data at the genotype level. The BIC penalizing constant α = log(N) was used throughout. Figure 1 compares the running times of the two algorithms, and confirms that the fast algorithm is approximately linear in the number of SNPs. Table 1 gives more detailed information. It is seen that the estimateG is indeed a good approximation toĜ, and that the algorithm is considerably faster than the standard greedy algorithm. Figure 2 shows the LD graph for one chromosome (Duroc chromosome 1): the same graph is obtained with both algorithms. It is a long thin graph with 2863 nodes and 4340 edges. It has 32 connected components, of which 20 are isolated vertices -suggesting SNP positioning errors -and four are long intervals with 997, 569, 786 and 452 SNPs. Curiously, the subgraph induced by the last 35 SNPs contains 7 connected components, which suggests positioning errors in this region. The most striking feature of the graph is that for about 80% of the chromosome, a simple serial or near-serial association structure is found, but for the remaining 20% more complex patterns of LD are observed. Similar results are found for all chromosomes and breeds.
Genomic intervals with simple association structure tend to be associated with low haplotype diversity, and The table summarizes the inaccuracy of the fast selection algorithm applied to data sets simulated under two models for Duroc chromosome 1: the skeleton, and the more complex LD graph shown in Figure 2. The columns labelled undershoot and overshoot show |E \ E j | and |E j \ E|, where E and E j are the edge sets of the true model and the model found by the algorithm, respectively. For the data sets simulated under the skeleton, the edgewise false negative rate 1 10|E| j |E \ E j | is 0.029 and the false positive rate 1 10|E| j |E j \ E| is 0.00024. Under the more complex LD graph, the corresponding rates are 0.031 and 0.019. intervals with complex structure with high diversity (Figure 3). To further compare the low and high complexity regions, two representative intervals of the same length were selected: SNPs numbered 1800-1825 (low complexity), and SNPs numbered 2470-2495 (high complexity). These have sample entropies of 0.79 and 4.05, respectively. Their subgraphs are shown in Figure 4, and their sample trees in Figure 5. These show the low diversity of the low complexity interval, and the relatively high diversity of the high complexity interval. Corresponding triangular heatmaps are shown in Figures 6 and 7. The low complexity interval has high LD whereas the high complexity interval shows a more mosaic structure.
Thus low complexity regions tend to consist of series of haplotype blocks with high LD and low haplotype complexity, and may be dominated by a few common haplotypes. In contrast, regions of high complexity tend to have high haplotype diversity and little or no haplotype block structure. It must be stressed that the SNPs in such regions are not generally in linkage equilibrium: on the contrary, complex patterns of association, not marginal independences, are observed. Figure 2 and Table 1 show the results of applying the algorithms to unphased, genotype data, but as mentioned above they may also be applied to phased haplotype-level data. This implicitly regards the inferred haplotypes as a random sample of size 2N from an underlying population of haplotypes. Analyses based on inferred haplotypes may be subject to a loss of efficiency [36]. For the current data phase imputation was carried out using the BEAGLE software [22]: applying the fast algorithm to the phased data for Duroc chromosome 1 resulted in the graph shown in Figure 8. This is slightly denser than the graph in Figure 2, with 375 more edges. This may be ascribed to the reduced model dimension due to the use of binary rather than trichotomous variables, which leads to a weaker penalization of complex models in (1). Haplotype-and genotype-level graphical models are related but in general distinct: it has been shown that they are identical only when the haplotype-level graph is acyclic [37].
To assess the accuracy of the algorithm, twenty simulated data sets each with the same number of observations as the data sets analyzed above (N = 4239) were constructed. The first ten were generated by taking N random samples from the skeleton, that is, the graphical model with edges between physically adjacent SNPs only, fitted to the Duroc chromosome 1 genotype data. The second ten were generated by taking N random samples from the model shown in Figure 2, fitted to the same data. The results of applying the fast algorithm to each data set are summarized in Table 2. For the skeleton, the edgewise false negative rate was 0.029 and the edgewise false inclusion rate was 0.00024. For the model in Figure 2 the corresponding rates were 0.031 and 0.019. http://www.biomedcentral.com/1471-2105/ 14/179 To visually represent the latter results, Figure 9 shows a graph G = (V , ∪ j E j ), in which the edge colours represent the frequency that the edge was found in the 10 models. The graph suggests that model uncertainty is primarily restricted to the genomic intervals with complex dependence structure.
Finally, for comparison purposes the algorithm of [15] was applied to the Duroc chromosome 1 genotype data. This algorithm automatically imputes phase and missing marker data, cycling between imputation and model selection as sketched in the background section. The selected graph is shown in Figure 10. Like the graphs shown in Figures 2 and 8, it is long and stringy, but rather more dense. This greater density may be explained by the use in [15] of a Metropolis acceptance rule that is based on the penalized likelihood (1), but with a smaller penalizing constant α = (1/8) log(N). A detailed comparison of the two algorithms would be valuable but is not attempted here.

Discussion
This paper has introduced an efficient, linear-time forward-backward algorithm to estimate chromosomewide probabilistic graphical models of fine-scale linkage  Figure 9 The union of LD graphs for Duroc chromosome 1 derived from ten simulated data sets. The figure represents the results of applying the fast algorithm to ten data sets simulated under the LD graph shown in Figure 1. The union of the ten selected graphs is shown. Edges present in all ten graphs are black, and the colour of the others is interpolated between light blue (those found once) and blue (those found nine times).

Figure 10
Haplotype-level graph found using the algorithm of [15]. The figure shows the result of applying the algorithm of [15] to the same data set as used with Figure 2. The algorithm estimates a haplotype-level graph for the data. disequilibrium, and has described a convenient way to display these models. In illustration, the methods have been applied to data obtained from three commercial breeds of pigs using the Illumina Porcine SNP60 BeadChip. The forward part of the algorithm proceeds by combining a series of overlapping marginal decomposable models of dimension L and overlap length K. This implicitly assumes that the maximum width of the true graph is at most K. The resulting model is then triangulated and backward selection performed. The search space closely resembles that of [15], which samples from decomposable models with a given maximum width, using sliding windows. The difference in approaches lies primarily in the search method: here greedy search to optimize a penalized likelihood criterion is used, whereas in [15] MCMC sampling methods are applied.
The R function used for greedy forward search currently requires that the input data contain no missing values, so it was necessary to impute missing values prior to performing the algorithm, and the BEAGLE software [22] was used for this. This raises the possibility of circularity, or more precisely, that the model selection is influenced by constraints or assumptions implicit in the models used by BEAGLE. But given BEAGLE's high imputation accuracy with such data it seems unlikely that this plays an important role here. http://www.biomedcentral.com/1471-2105/14/179 The algorithm was found to have high accuracy when applied to simulated data sets of the dimensions considered here (N ∼ 4000; p ∼ 3000), with edgewise false positive and negative rates of around 3%. That it is good at reconstructing the model generating the data is reassuring. Needless to say, the algorithm does not necessarily identify the "true" model, which may not be in the search space. As noted previously, the approach captures short range associations, but not long range ones. Moreover, higher edgewise false negative rates may occur when observed data are used if there are infinitesimal departures from a model that are not detectable at the given sample size, as these are filtered out when data are simulated under a selected model.
A comparison of triangular heatmaps with LD graphs is instructive. The former are compact graphical representations of all pairwise marginal associations for a set of SNPs. They are particularly well-suited to identify LD blocks, which stand out as highlighted triangles. LD graphs supplement heatmaps by showing patterns of association but not their strength. Since they are parametric models they can be put to a number of quantitative uses as described in the background section. At one level, heatmaps and LD graphs convey similar information, since intervals with simple dependence structures tend to appear as series of LD blocks in the heatmap, whereas those with complex structures tend to occur in the interblock regions and have a more mosaic appearance. But the graphs provide a more incisive characterisation of genetic variation, building on the concept of conditional rather than marginal dependence. In this regard, it may be helpful to regard conditional independence statements as expressing the notion of irrelevance, in the sense that A ⊥ ⊥ B | C implies that if we know C, information about A is irrelevant for knowledge of B. Thus the graphs say something about connections between specific SNPs, for example about which SNPs are required to predict a specific SNP.
The graphs may also be useful in fine mapping. Genome-wide association studies seek to find the genetic basis of complex traits, typically using single locus methods -that is, by identifying SNPs with strong marginal association with the complex trait. Due to LD, many SNPs in a genomic region may exhibit strong association with the trait, making it hard to identify the causal loci. A way to address this is to assess the effect of a putative causal SNP on the complex trait in a linear model that also includes terms for the two flanking SNPs, in order to adjust for confounding with nearby effects due to LD. This implicitly assumes a simple serial dependence structure between the SNPs, and when the dependence structure is more complex such adjustment might be insufficient, leading to false positives. This may be prevented by including terms for the neighbours of the putative causal SNP, not just for the flanking SNPs. Any SNP is separated from the remaining SNPs by its neighbours, so by the global Markov property it is independent of the remaining SNPs given its neighbours. A similar method has been proposed based on Bayesian networks [38].
Like graphical models, the APFA models that underlie BEAGLE may be represented as graphs that encode a set of conditional independence relations, but in a very different way. In the present context, for example, nodes in APFA graphs represent haplotype clusters rather than SNPs. Where the model classes intersect, APFA graphs are much more complex than the corresponding dependence graphs, and so less amenable to visualization.
A striking feature of the LD graphs for the pig data was that for roughly 80% of the genome, simple serial or nearserial LD patterns were found, but for the remaining 20%, more complex patterns were observed. The regions with the simple serial structure tend to have low haplotype diversity, which is to be expected in livestock breeds with small effective population sizes [39]. Perhaps more unexpected is that roughly 20% of the porcine genome exhibits complex LD patterns, forming islands of relatively high genetic diversity. This information may be useful in an animal breeding context, to identify regions with high genetic variation. It will also be interesting to compare graphs obtained using different SNP densities in a given breed or species to examine whether and how their topologies change with varying marker densities.

Conclusions
The proposed algorithm is efficient and makes it feasible to estimate and visualize chromosome-wide LD models on a routine basis.


extfieldsofstudy: Medicine Biology

text:
Microarray and Pathway Analysis Reveal Distinct Mechanisms Underlying Cannabinoid-Mediated Modulation of LPS-Induced Activation of BV-2 Microglial Cells

Cannabinoids are known to exert immunosuppressive activities. However, the mechanisms which contribute to these effects are unknown. Using lipopolysaccharide (LPS) to activate BV-2 microglial cells, we examined how Δ9-tetrahydrocannabinol (THC), the major psychoactive component of marijuana, and cannabidiol (CBD) the non-psychoactive component, modulate the inflammatory response. Microarray analysis of genome-wide mRNA levels was performed using Illumina platform and the resulting expression patterns analyzed using the Ingenuity Pathway Analysis to identify functional subsets of genes, and the Ingenuity System Database to denote the gene networks regulated by CBD and THC. From the 5338 transcripts that were differentially expressed across treatments, 400 transcripts were found to be upregulated by LPS, 502 by CBD+LPS and 424 by THC+LPS, while 145 were downregulated by LPS, 297 by CBD+LPS and 149 by THC+LPS, by 2-fold or more (p≤0.005). Results clearly link the effects of CBD and THC to inflammatory signaling pathways and identify new cannabinoid targets in the MAPK pathway (Dusp1, Dusp8, Dusp2), cell cycle related (Cdkn2b, Gadd45a) as well as JAK/STAT regulatory molecules (Socs3, Cish, Stat1). The impact of CBD on LPS-stimulated gene expression was greater than that of THC. We attribute this difference to the fact that CBD highly upregulated several genes encoding negative regulators of both NFκB and AP-1 transcriptional activities, such as Trib3 and Dusp1 known to be modulated through Nrf2 activation. The CBD-specific expression profile reflected changes associated with oxidative stress and glutathione depletion via Trib3 and expression of ATF4 target genes. Furthermore, the CBD affected genes were shown to be controlled by nuclear factors usually involved in regulation of stress response and inflammation, mainly via Nrf2/Hmox1 axis and the Nrf2/ATF4-Trib3 pathway. These observations indicate that CBD, and less so THC, induce a cellular stress response and that this response underlies their high immunosuppressant activities.

Introduction
Cannabis sativa extracts (marijuana and hashish) have been used for centuries both as therapeutic agents and recreational drugs, primarily due to their ability to regulate neurobehavioral processes including memory, mood and appetite [1,2]. The spectra of possible therapeutic uses of marijuana and its active constituents, the cannabinoids, range from handling nausea, vomiting and cachexia (in cancer chemotherapy and AIDS patients) to treatment of chronic pain, glaucoma, epileptic seizures, Parkinsonian tremor as well as multiple sclerosis [1,[2][3][4][5][6]. Amongst the 60 different cannabinoids identified from Cannabis preparations, D 9tetrahydrocannabinol (THC), the major psychoactive Cannabis constituent, and the non-psychoactive cannabidiol (CBD) are the most abundant and important.
Many of the beneficial effects of cannabinoids have been attributed to their potent immunosuppressive and anti-inflammatory properties [7][8][9][10][11][12][13][14]. Additionally, cannabinoids are known to possess pro-apoptotic, neuroprotective and anti-tumor properties [15][16][17]. To date, two cannabinoid receptors have been characterized: the CB 1 and the CB 2 receptors. The CB 1 receptor is highly expressed in neural cells and mediates the psychoactive and addictive activities of cannabinoids, while the CB 2 receptor is abundantly present in the periphery including the immune system and is involved in cannabinoid immunomodulation [18][19][20]. THC binds to both of these receptors with similar efficiency [21] and has been reported to have effects on both the nervous and the immune systems [7,11,18,22].
Phytocannabinoids were reported to affect various populations of immune cells [7,[11][12][13][14]31]. Both THC and CBD have been shown to decrease cytokine production in human immune cell lines [32] and to suppress T cell proliferation and their effector functions [12,33,34]. In response to stimulation with the bacterial endotoxin lipopolysaccharide (LPS), both monocytes and microglial cultures treated with either THC or CBD produce lower levels of cytokines such as tumor necrosis factor a (TNFa), interleukin-1a (IL-1a), IL-1b and IL-6 [11,35]. However, the molecular mechanisms involved in these cannabinoid-mediated effects are not yet fully characterized. Eljaschewitsch et al., [36] showed that activation of CB 1 and/or CB 2 receptors in the murine microglial BV-2 cell line leads to rapid induction of mitogen-activated protein kinase phosphatase-1 (MKP-1) and that this event switches off MAPK signal transduction which was activated by LPS stimulation. Altered adenosine signaling (through inhibition of its uptake) has been reported as a potential non-cannabinoid receptor mechanism by which CBD, but less so THC, can decrease inflammation [37]. Other studies identified the nuclear receptor peroxisome proliferator-activated receptor c (PPARc) as a novel intracellular target, which mediates the cannabinoid-associated immunosuppression in a manner that is independent of the known cannabinoid receptors CB 1 and CB 2 [38,39]. Other targets including the G-protein-coupled receptors GPR55 and GPR18 as well as the transient receptor potential (TRP) channels were also suggested [19,20,30,40].
Microglial cells are the resident macrophage-like cells of the CNS. They are highly ramified cells and their processes are very dynamic under non-pathological conditions, actively scanning their environment. These cells have important roles in brain's innate immunity and neuronal homeostasis as well as in neuroinflammatory pathologies [41]. Microglia can be activated by infection, injury or by endogenously released neurotoxic factors and their activation is associated with the production and secretion of a variety of compounds such as cytokines, reactive oxygen species (ROS), reactive nitrogen species, matrix metalloproteinases and prostaglandins. Although microglial activation is considered a protective mechanism involved in the clearance of pathogen infection and in regulating tissue repair and recovery, excessive or chronic activation can lead to harmful effects [42]. Interestingly, the mechanisms that give rise to either the protective or the damaging microglial phenotypes are not fully elucidated. Enhancing the microglial-mediated innate immunity in the CNS and/or preventing the harmful effects associated with their chronic activation may offer new therapeutic approaches for the treatment of brain injury and neurodegenerative diseases [43].
One of the most potent stimuli for microglia activation is the bacterial endotoxin LPS, that mimics infection by Gram-negative bacteria. LPS activates intracellular signaling pathways in a complex way leading to secretion of cytokines and to overexpression of several markers of the immune response. Previous studies reported that LPS stimulation induces gene expression of TNFa, IL-1b, IL-6, iNOS and COX-2 as well as the production of NO and PGE 2 in primary and BV-2 microglial cell cultures [11,44,45]. We reported that CBD reduces the activity of the NF-kB pathway and upregulates the activation of STAT3 transcription factor in LPSstimulated BV-2 cells, and that both CBD and THC decrease the activation of the LPS-induced STAT1 transcription factor, a key player in IFNb-dependent pro-inflammatory processes [11]. Moreover, performing comparative microarray analysis of genome-wide mRNA levels in the BV-2 cells, we reported that CBD, but less so THC, shows a specific gene expression profile associated with oxidative stress and glutathione depletion involving the GCN2/eIF2a/p8/ATF4/Chop-Trib3 pathway [13]. Furthermore, the CBD-stimulated genes were shown to be controlled by nuclear factors known to be involved in regulation of stress response and inflammation, mainly via the (EpRE/ARE)-Nrf2/ Atf4 system and the Nrf2/Hmox axis. We reported that CBD, but less so THC, affects the expression of genes involved in zinc homeostasis, suggesting that the regulation of zinc levels could have an important role through which CBD may exert its antioxidant and anti-inflammatory effects [14].
Although the inhibitory functions of cannabinoids on LPSactivated NF-kB and IFN-b/STAT proinflammatory pathways and on the secretion of selected cytokines in BV-2 microglial cells has been studied [11], a genome-wide search for cannabinoid molecular targets in LPS-activated BV-2 cells has not yet been performed. We have, therefore, performed gene array studies and comparative gene profiling analysis of BV-2 cells treated with LPS, CBD+LPS or THC+LPS. This approach allowed us to analyze the changes induced by CBD and THC on gene expression patterns in LPS-treated BV-2 cells and to explore the genomewide interaction network affected by these treatments. In this regard, a structured network knowledge-based approach to analyze genome-wide transcriptional responses in the context of known functional interactions among proteins, small molecules and phenotypes has been established [46]. We applied this analysis to show the interactions and signaling networks elicited by the cannabinoids in LPS-stimulated BV-2 cells. Our results show that CBD is a potent modulator of microglial activation. Identification of the CBD-and THC-regulated genes and related networks provides a molecular basis for understanding the effects of these cannabinoids on LPS-activated microglia.

Materials
LPS from Escherichia coli (serotype 055:B5) was purchased from Sigma (St. Louis, MO, USA). THC and CBD were obtained from the National Institute on Drug Abuse (NIDA; Baltimore, MD, USA). Stock solutions of cannabinoids were prepared in ethanol and diluted into culture medium before experiments. Final ethanol concentration in the medium did not exceed 0.1%.

Total RNA extraction
Extraction, quantification and quality of extracted RNA were performed and analyzed as described previously [13].

Microarray transcript analysis
Comparative microarray analysis was performed using the Illumina MouseRef-8 BeadChip platform, as described previously [13]. Four replicates of each experiment were carried out, resulting in four independent microarrays for each individual treatment and controls (24 total arrays). Gene products that were affected by 2-fold or more were further analyzed as described in the Result section. Differentially expressed genes were classified according to their gene ontology (GO; http://www.geneontology. org/; [47]), using DAVID Bioinformatics online tools (Database for Annotation, Visualization and Integrated Discovery; http:// david.abcc.ncifcrf.gov/; [48]). Cellular pathway association was analyzed according to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/) and pathway maps according to BioCarta (http://www.biocarta. com/genes/index.asp).

Ingenuity pathway analysis
Pathway and global functional analyses were performed using Ingenuity Pathway Analysis 6.0 (IPA; IngenuityH Systems, www. ingenuity.com). A data set containing gene identifiers and corresponding expression values was uploaded into the application, and each gene identifier was mapped using the Ingenuity Pathways Knowledge Base (IPKB). The IPKB analyses identify the biological functions as well as the pathways from the IPA library that are most significant to the data set. Genes from the data sets associated with biological functions or with a canonical pathway in the IPKB, that met the p-value cutoff of 0.005 were used to build the interactome as described below. Fisher's exact test was used to calculate a p-value determining the probability that each biological function and/or canonical pathway assigned to this data set was not due to chance alone.

Interactome/network analysis of differentially expressed genes
The accession number and the ratio (CBD+LPS versus LPS and THC+LPS versus LPS) for each gene whose expression was changed by at least 2-fold by one of the treatments (focused genes), were uploaded into the IPA system. This IPA application was used to query databases for interactions between sets of differentially expressed genes and all other genes stored in the knowledge base to generate a set of interactive networks taking into consideration canonical pathways, relevant biological interactions as well as cellular and disease processes. The IPA system computes a score for each network according to the fit of the set of the supplied ''focus genes''. The p value scores for such a set indicate the likelihood of focused genes to belong to a network versus those obtained by chance.

Microarray data validation by quantitative real time PCR
Many of the gene products that were found to be affected using the microarray gene analysis were validated by quantitative real time reverse transcription polymerase chain reaction (qPCR). Primer sets for qPCR were designed using Primer Quest, an online tool provided by Integrated DNA Technologies (http://eu.idtdna. com/Scitools/Applications/Primerquest/) ( Table S1). Wherever possible, designs with at least one of the primer sequences located on an intron-exon boundary were chosen, thus avoiding coamplification of minor contaminating amounts of genomic DNA that could be present in the RNA samples. All primers were analyzed using the nucleotide program BLAST to ensure primer specificity for the gene of interest (http://blast.ncbi.nlm.nih.gov/ Blast.cgi). cDNA was generated by using the QuantiTect Reverse Transcription kit containg ''genomic DNA wipe out'' (to eliminate contamination with genomic DNA), according to the manufacturer's instructions (Qiagen, AG, Basel, Switzerland). qPCR was carried out as detailed by Juknat et al., [13] using the Rotor-Gene 3000 instrument (Corbett Research, Sydney, Australia). Expression levels of genes of interest were normalized to the reference gene, b 2 -microglobulin (B2m), whose expression was found not to be affected by the various treatments, and are expressed as fold change using the calculation method described previously [13]. The qPCR experiments were repeated 3 to 4 times using different mRNA batches from independent experiments and reactions were performed in duplicates for each cDNA sample.
Statistics qPCR data were plotted as the mean6SEM of 3-4 independent experiments. Statistical significance was assessed using a one-way or two-way ANOVA, followed by Bonferroni post hoc multiple comparison test as implemented with the Statistics Toolbox Software MATLAB, Version 6.1 (R2007b), MathWorks (http:// www.mathworks.com/help/toolbox/stats/rn/brasjn_.html). A p value ,0.05 was defined as statistically significant.

Results
Effect of cannabinoids on LPS-stimulated gene expression in BV-2 microglial cells BV-2 microglial cells were pretreated for 2 h with either THC or CBD (both at 10 mM) followed by the addition of LPS (100 ng/ ml) to the incubation medium for another 4 h. The control treatment with cannabinoids alone lasted for 6 h and with LPS alone for 4 h. The choice of these time points for transcriptional profiling was guided by our previous studies [11,13] as well as by other reports [44,45] which investigated the general temporal pattern of microglial activation by LPS. Moreover, according to our previous results, neither THC nor CBD treatments (both at 10 mM) significantly affected the viability of the BV-2 microglial cells during this 6 h period [11].
The RNA prepared from these samples (in sextuplicates) was analyzed for changes in transcriptional levels using the MouseRef-8 v1.1 Expression BeadChip Illumina Arrays. Each of these arrays has .24,000 mouse targets based on the NCBI mouse Reference Sequence Database, including 16,287 constitutive exons/islands based on the splice variants in the mouse transcriptome (Molecular Signature Database; MouSDB3) and NCBI LocusLink databases. The results of the analyses of the arrays showed that 32% of the transcripts were consistently ''present'' in the BV-2 RNA samples across all arrays. Moreover, clustering based on inter-array Pearson correlation coefficient indicated no batch effects. Microarray analysis based on a threshold of p# 0.005, revealed that a total of 22% (5338 out of 24,000 transcripts) of the Illumina gene set was differentially regulated across treatments. Of these, 1319 gene probe sets were upregulated and 1829 transcripts were downregulated by the LPS treatment ( Figure 1A); and from these numbers of genes, 400 transcripts were found to be upregulated and 145 downregulated by LPS by 2-fold or more. When the fold change was set on $3-fold (p# 0.005), we found that 226 gene products were upregulated and 33 were downregulated by LPS ( Figure 1A). The vast majority of the LPS-affected transcripts (89% or 1181 gene probe sets of the upregulated genes, p# 0.005 and 90% or 1641 gene probes of the downregulated transcripts, p# 0.005) represented genes that were exclusively responsive to LPS stimulation, and not to treatment with CBD alone or THC alone ( Figure 1B).
the CBD+LPS treatment was analyzed, 1381 gene probe sets (25.8%) were upregulated and 1666 transcripts (31.2%) were downregulated by CBD+LPS ( Figure 1A). From these numbers of genes, 379 gene products were found to be upregulated and 489 downregulated exclusively as a response to the combination of CBD+LPS and not affected by LPS or CBD alone ( Figure 1B). When the fold change was set on $2-fold, 502 transcripts were upregulated and 297 gene probe sets were downregulated by the CBD+LPS treatment ( Figure 1A).
When LPS was applied in the presence of THC, 1216 gene probe sets (22.8%) were upregulated and 1638 transcripts (30.7%) were downregulated; and from these, 424 transcripts were found to be upregulated and 149 downregulated by 2-fold or more ( Figure 1A). From the 1216 upregulated genes, 157 transcripts were exclusively responsive to THC+LPS treatment alone and from the 1638 downregulated transcripts, 285 gene probes were affected only when THC and LPS were together ( Figure 1B). When the fold change was set on $3-fold, we found that 230 gene products were upregulated by CBD+LPS and 236 transcripts by THC+LPS whereas 51 gene products were downregulated by CBD+LPS and 31 transcripts by THC+LPS ( Figure 1A).
Our results also reveal that from the 5338 transcripts found to be differentially regulated by the various treatments (p # 0.005), 680 gene probe sets were found to be upregulated by CBD alone and 58 gene products by THC alone. CBD had also a much larger effect compared to THC on the number of downregulated genes, 524 gene products were downregulated by CBD and 36 by THC ( Figure 1A).
The groups of LPS-upregulated and downregulated genes that showed a change in expression of $2-fold (p#0.005) in either direction were subjected to gene ontology (GO) analysis allowing functional annotation using the DAVID Bioinformatics Resources and the KEGG pathways ( Figure S1). All major Biological Processes, Cellular Components and Molecular Functions within GO for the LPS-upregulated transcripts, were for the most part, genes associated with the immune and defense responses as well as apoptosis and cell death. KEGG database analysis included genes related to Toll-like receptor pathways and antigen processing as well as to the MAPK pathways. IPA global functional analysis of the LPS-upregulated and downregulated genes ($2-fold; p# 0.005) confirmed the DAVID representation of genes within IPAspecific GO categories ( Figure S2). IPA global functional analysis of the LPS-upregulated transcripts show gene products mainly implicated in cell-to-cell signaling, cellular movements, growth and proliferation, as well as cell death, immune response and signaling, including genes involved in the NF-kB pathway and in the activation of the liver X receptor/retinoic X receptor (LXR/ RXR). The largest subsets of downregulated transcripts included genes known to be involved in the regulation of macromolecules and nucleotide metabolism, differentiation and development as well as regulation of gene expression and transcription (Figures S1 and S2).
A list of 183 LPS-upregulated genes (stimulated by 2-fold or more; p# 0.005) are presented in Table S2. This list includes LPS-upregulated genes that were (in most cases) significantly affected by CBD+LPS or by THC+LPS. They appear in the Table according to their categories and specific annotations based on their GO and IPA. The highly upregulated transcripts were related to inflammation, host defense and adaptive response. Inflammatory cytokines and chemokines together with their cognate receptors comprised a large group of more than 20 genes. Interferon related transcripts formed a significant population comprised of 12 genes. Based on GO and relevant literature references, genes were also classified into other functional categories including metabolic enzymes, membrane transport and secretion, kinases, phosphatases as well as transmembrane Gprotein coupled receptors. The remainder of the LPS-upregulated genes was grouped according to their involvement in several other cellular processes such as apoptosis, proliferation and cell cycle progression, transcriptional and translational control as well as stress response. Additionally, we have identified a cluster of genes associated with regulation of extracellular matrix that are known to be involved in molecular recognition between cells, cell adhesion and migration. Genes implicated in cytoskeleton remodeling and control of cell motility and morphogenesis form a separate category.
When CBD was applied without LPS stimulation, it did not have a significant effect on these genes, except for Ccl7, Ccl2 and Ccl9 which were downregulated by CBD by 71.4%, 71.4% and 50%, respectively. Careful evaluation of the list of genes reveals that the inhibitory effect of CBD on LPS-upregulated genes was concentrated mainly within gene categories related to inflammation. These included cytokines and chemokines, interferon regulated components of the immune response as well as molecules involved in host defense and adaptive immune response. The inhibitory effect of CBD on LPS-upregulated genes was less strong on other gene products such as those for regulatory molecules, kinases, phosphatases, metabolic enzymes, G proteincoupled receptors and molecules related to regulation of transcription.
From the 82 LPS-downregulated genes presented in Table S3, CBD pretreatment counteracted the reduction in mRNA levels of 16 genes (by 30% or more) and 2 were affected by more than 60%. These genes almost exclusively belong to the functional categories of cell adhesion and migration, motility and morphogenesis, regulation of transcription as well as of membrane transport and secretion. These genes ( Table 3) include neurophilin 1(Nrp1; whose reduction was completely abolished), neural precursor cell expressed developmentally down-regulated gene 9 (Nedd9; reversed by 69%), basic helix-loop-helix domain containing class B2 (Bhlhb2; by 56%), unc-5 homolog B (Unc5b; by 53%), Cd28 (by 49%), SRY-box containing gene 4 (Sox4; by 42%) and potassium channel, subfamily K, member 13 (Kcnk13; by 42%). Treatment with THC resulted in a much weaker attenuation of the LPS-induced downregulation of gene expression and only in the case of Nrp1 exceeded 30%. DAVID Bioinformatics Resources and KEGG pathway analyses identify Nrp1 and Unc5b as genes related to axon guidance and Nedd9 as a gene related to adherens junction dynamics.

Validation of microarray results
A large number of genes that were identified by microarray analysis as differentially regulated were subjected to validation by qPCR using B2m as a reference gene. Many of them were selected to be validated according to the distinct effects of cannabinoids on the LPS affected genes. The qPCR assays were repeated several times using at least 3 mRNA preparations from independent experiments. The results are expressed as fold change relative to control levels. We found that in almost all cases, there has been a very good agreement between the microarray and the qPCR data in terms of direction of change as well as of its magnitude ( Table 5). For example, the qPCR data show that LPS upregulates the expression of the inflammatory chemokines Ccl2, Ccl12 and Cxcl10, the transporters and channels Aqp9, Mcoln2 (TRPML2) and Slc7a11, the kinase inhibitor Cdkn1a, the phosphatases Dusp1 and Dusp2, the regulatory molecule Sqstm1/p62, the G protein-coupled receptor Ptgir and the receptor Tlr2. The qPCR data also show that LPS downregulates the expression of the G protein-coupled receptors Ebi2/Gpr183, Cnr2 and Gpr55, and the transcription factors Pparg1 and Pparg2. In addition, in agreement with the microarray data, the stress response genes Trib3 and Hmox1 are upregulated by CBD+LPS but not by LPS alone or by THC+LPS.

Network analysis and signaling pathways -Functional associations between the effects of cannabinoids and LPS
The differential mRNA expression patterns as revealed by the arrays of BV-2 cells treated with LPS alone or in combination with either CBD or THC, were analyzed using the Ingenuity System Database, software that includes the Ingenuity Knowledge Base (IKB) and the Global Molecular Network (GMN). These databases integrate published findings on biologically meaningful genetic or molecular gene/gene product interactions and identify functionally related gene networks. The differential gene expression values (e.g.; CBD+LPS versus LPS treatment or THC+LPS versus LPS) were entered into IPA to determine the most highly regulated networks of gene interactions and to highlight the biological processes that are relevant to each of the treatments. The IPA computation gives a score for each network according to the fit of the set of supplied affected genes  that are also contained in the GMN (focus genes: network eligible molecules). The scores indicate the likelihood that ''focus genes'' belong to a network versus those that are obtained by chance. A score of .2 indicates a $99% confidence that a ''focus gene'' network was not generated by chance. According to the degree of interconnectedness among the molecules, a higher or lower IPA network score is assigned. Only networks with a score of 15 or higher, have been selected for further analysis. Networks describing the relationships between a subset of genes and their neighboring genes are presented in Figures 2, 3, 4 for CBD and  (Figure 2) received the highest score (48) and contains 34 differentially regulated genes. This network has in its central position two transcription factors, ATF4 and NFkB, interconnecting several important canonical pathways. The later include Nrf2/ATF4-mediated oxidative stress response (including the highly upregulated Trib3, Didit3/Chop/Gadd153, Asns, Gfd15), aminoacyl-tRNA synthetases (including Sars, Nars, Gars), and ERK/MAPK and related phosphatases (such as Dusp2 and EIF2AK1). In addition, 3 members of the Bcl-2 family, Bik, Bcl2l2 and Bag4, related to apoptotic processes were found to be present in the NFkB interactome. The Nrf2/ATF4-regulated transcripts were mainly upregulated (red symbols) whereas the genes downstream of NFkB were mostly downregulated (green symbols). The genes Trib3 and Gdf15 seem to serve as a point of interaction between these two families of gene transcripts. The ATF4 and NFkB networks are involved in cellular and molecular functions of cell death, protein synthesis and cancer [49,50].
Network 2 ( Figure 3) with a score of 32 contains 27 genes. This network shows Fos, a member of the AP-1 transcription factor family, as a hub of connectivity, as it appears to play a central role possessing 15 edges (via direct connections). Two other transcripts belonging to the AP-1 superfamily, FosL2/Fra-2 and FosL1/Fra-1 (a gene whose product is known to decrease AP-1-dependent transcription) are also present in the Fos interactome. FosL2/Fra-2 was found to be downregulated while FosL1/Fra-1 is upregulated. A series of canonical pathways including those involved in glutathione and xenobiotic metabolism (Mgst2, Gstm5, Gstm6, Gsto2, Gstm4, Gsta5, zinc metallothionein 2a (Mt2a), Znf467, folate transporter Slc19a1/Rfc1, Pdk1, Ube1c/Uba3) as well as those involved in LPS/IL-1 mediated inhibition of either RXR or the peroxisome proliferator activated receptor (PPAR) signaling (Rar, Rdh5, Ncoa5, Ap1, Fos, Ets) were found to participate in this network. Fos directly interacts with Gsta5 (from the glutathione-Stransferase family), known to belong together with cytochrome P450 enzymes, NAD(P) H:quinone oxido-reductase 1 (Nqo1) and aldehyde dehydrogenase 3 to a group of genes encoding phase 1 and phase 2 xenobiotic metabolizing enzymes [51]. Altogether, Network 2 is in correlation with a set of cellular and molecular functions which include amino acid-, vitamin-and mineralmetabolism and transport as well as cell signaling. Network 3 ( Figure 4) with a score of 28, includes 25 genes. This network presents p38MAPK in a central position and overlaps with canonical pathways for calcium signaling (Nfat5, Itpr1, Pka, ATP2A2, Ank2, Ogt, Mef2b), FGF pathway (Fgf21, ERK1/2), B cell signaling (Mef2c, Itpr1, ERK1/2, Sh3bp2, Creb) and p38MAPK signaling (Dusp8, Mef2c, Creb, Nfat5, Esrra). The IPA interactome analysis has linked the CBD-stimulated expression of Dusp8 (via the negative modulation of p38MAPK) to the attenuation of LPSmediated induction of Nfat5, a member of the inflammatory transcription family. In addition, 3 members of the complement component 1q (C1qa, C1qb and C1qc) that recognize antibodybound pathogens and trigger the complement classical pathway, were found to be present in the p38MAPK interactome. This network seems to be related to functions of cell signaling, tissue morphology and neurological diseases.
The levels of changes in gene expression induced by combined treatment of LPS and THC (versus that induced by LPS alone) are smaller than those observed for the combined treatment of LPS with CBD (versus LPS alone). These LPS+THC affected genes encompassed three networks of specific molecular interactions. Network 4 ( Figure 5) with a score of 32 contains 16 differentially regulated focus molecules and presents the transcription factors AP-1, Akt and NFkB in its central positions. The gene products within the displayed network were found to be part of a number of canonical pathways. For example, the genes X-box binding protein 1 (Xbp1) and Tgfbr2/3 are related to the glucocorticoid receptor signaling. The progressive multifocal leukoencephalopathy (Pml) is related to retinoic acid receptor (RAR) activation and PPAR signaling, and the neuraminidase 3 (Neu3) to ceramide pathway and to p38MAPK. The gene Trem1 is involved in immune response and Vegfa and Thbs1 are related to the angiogenesis pathway. Altogether, this network shows that THC affects several biological processes including immune and lymphatic system development and function as well as tissue morphology.
THC also led to changes in the LPS effects on several other groups of genes. Network 5 received a score of 23 and is populated with 12 genes, with retinoic acid as a transcription regulator positioned at its center ( Figure 6). A number of canonical pathways related to axonal guidance (Sema3a, Rock2, Dpysl2/Crmp2, Dpysl5/Cram, Spg3a/Atl1, Foxd3), synaptic plasticity (Raver1), cell adhesion (Adam2), actin cytoskeleton (Rock2, Foxd3) as well as gene products involved in the activation of the nuclear receptor family of transcription factors LXR/RXR (Es1, Rnf113a1), were found to participate in this network. Cell to cell interaction as well as tissue development are the major biological processes related to this network.
Finally, Network 6 with a score of 16 and 9 genes shows Tp53, Cdkn1A and beta-estradiol at the central nodes of the interactome (Figure 7). This network overlaps with the aryl hydrocarbon and Figure 2. Ingenuity Interactome analysis of the effect of CBD on LPS-affected gene expression. Network analysis of the genes whose expression was affected by LPS+CBD vs LPS was performed using the Ingenuity software. Each network displays the genes/gene products as nodes (different shapes representing the functional classes of gene products) and the biological relationships between the nodes as lines. The length of each line reflects the amount of literature evidence supporting this node-to-node relationship. The color intensity of each node indicates the degree of upregulation (red) or downregulation (green) of the respective gene transcript. Genes in gray have been shown (in the literature) to interact with the colored gene products that appear in this scheme. Network 1: The NFkB-ATF4 interactome -Interaction between the CBD-upregulated Trib3 and the attenuation of the NF-kB transcription factor (responsible for the transcription of many proinflammatory genes) shows the connection of this network to anti-inflammatory responses. doi:10.1371/journal.pone.0061462.g002

Effect of cannabinoids on LPS-stimulated gene expression in BV-2 cells
This study addresses the effects of the two major cannabinoids present in cannabis, CBD and THC, on mRNA expression in LPS-stimulated BV-2 microglial cells. Our results show that pretreatments with CBD or THC differentially affect LPSregulated gene transcription. The modulation of LPS-regulated gene expression by CBD was more pronounced as compared with THC. This result is in line with the more profound effects of CBD on mRNA regulation in surveillant (resting) microglial cells, as previously described by our group [11,13].
Here we show that CBD, more than THC, suppresses the LPSinduction of many proinflammatory genes in selected functional categories ( Table 1 and Table S2) including cytokines (e.g., Il-1b, Il-1a, Il-27, Il-6, Il-18, Tnfrsf5/Cd40) and chemokines (e.g.; Ccl12, Ccl7, Ccl9 and Ccl2). Moreover, CBD also repressed the basal expression of Ccl2, Ccl7 and Ccl9. Interestingly, Ccl2, Ccl7 and Ccl12 bind to the same receptor (CCR2), thus, showing binding promiscuity while they activate different signal transduction pathways in different cell populations. From these chemokines, Ccl2 plays a critical role in multiple sclerosis and in its murine model, experimental autoimmune encephalomyelitis (EAE). Several reports show that increased expression of Ccl2 in immune cells is closely associated with the clinical activity of EAE ( [52] and references therein). Thus, the CBD downregulation of Ccl2 mRNA is in line with the finding that CBD is ameliorating the EAE disease symptoms [12].
Another hallmark of inflammation is the increased expression of proinflammatory mediators like matrix metalloproteinases (Mmps). Here, we show that CBD downregulates the expression of the LPS-upregulated matrix metalloproteinase13 (Mmp13). MMP13 level and activity are enhanced in correlation with the degenerative changes in osteoarthritis cartilage and this molecule co-localizes with its specific type II collagen cleavage products. In agreement with this result, Malfait et al., [53] reported that in the murine collagen-induced arthritis (CIA) CBD treatment blocks the progression of the disease.
CD69 antigen is another mediator found to be involved in the CIA model of rheumatoid arthritis. This antigen is highly expressed in the leukocyte infiltrates of various chronic inflammatory diseases. In this regard, CIA disease severity is increased by antibody-induced blockade of the transforming growth factor b (TGF-b) in wild type but not in CD69 -/mice, suggesting that CD69 is a negative modulator of autoimmune response and inflammatory reaction [54]. Our current results show that CBD downregulates the mRNA expressions of the LPS-upregulated Cd69 antigen and of Mmp13 suggesting a possible mechanism for the therapeutic activity of CBD in the murine CIA model.
The expression of several genes in the inflammatory cytokine functional category was found to be enhanced by CBD (and to a lower extent by THC). These include the LPS-upregulated expression of Gdf15/MIC-1 and Il-15. From these genes, Gdf15 mRNA expression is upregulated (8.2-fold) in surveillant/ resting BV-2 cells by CBD (and 2.4-fold by THC) but synergistically upregulated by CBD+LPS (55.7-fold). Gdf15/ MIC-1 is a member of the TGF-b superfamily that plays key roles in the regulation of cellular responses to stress signals, inflammation and tissue repair. Gdf15 mRNA expression was found to be upregulated in activated macrophages by secreted proinflammatory cytokines (including TNFa, IL-1b and IL-6, as well as IL-2 and the macrophage colony-stimulating factor) suggesting that Gdf15 may act through an autocrine loop as an inhibitory factor in the late phases of inflammation by suppressing inflammation through the inhibition of macrophage activation [55].
Other genes highly upregulated by CBD are Lcn2 and Aqp9. Lcn2 has been implicated in many cellular processes, such as cell death/survival, cell migration/invasion, cell differentiation, inflammation and iron sequestration, including a role in the acute phase response ( [56] and references therein). With regard to Aqp9, we have previously shown that CBD upregulates Aqp9 mRNA in BV-2 cells (14-fold; [13]). Moreover, CBD synergistically increases (by 877%) the LPS-upregulated Aqp9 (3.1-fold). Aqp9 participates in the transport of small solutes (such as glycerol, lactate and urea) and takes part in osmotic swelling induced by apoptotic stimuli [57].
Two major kinase-mediated signaling pathways are activated following TLR4 stimulation by LPS. These include the MAPKs and IkB kinase complexes, which lead to activation of AP-1 and NF-kB transcription factors, respectively. Here we show that LPSupregulated expression of Dusp1 (MKP-1) was significantly upregulated by CBD (5.1-fold) but not by THC. An increase in the expression of Dusp1 was also observed following incubation of resting BV-2 cells with CBD, but not with THC. A similar CBDmediated regulation of mRNA level was observed for Dusp8, another negative regulator of MAPKs. The IPA interactome analysis has linked this CBD-stimulated expression of Dusp8, via negative modulation of p38 MAPK, to the attenuation of LPSmediated induction of NFAT5, a member of the inflammatory transcription complex. On the other hand, the LPS-stimulated expression of Dusp2 (Pac-1), a key positive regulator of the inflammatory MAPKs signaling [58], was significantly repressed by CBD as well as by THC. Moreover, CBD but not THC decreased basal levels of Dusp2/Pac-1 mRNA in resting BV-2 cells.
IPA and interactome analysis indicate that pretreatment with CBD, but less so with THC, result in attenuation of LPSstimulated activation of NF-kB and its dependent gene transcription pathways. Indeed, our previous data show that CBD decreases the activity of the NFkB signaling pathway in BV-2 microglial cells via the partial reversal of the LPS-induced degradation of IRAK-1 intermediate kinase, thus reversing IkB degradation and reducing NF-kB p65 subunit phosphorylation. This effect of CBD in decreasing NF-kB signaling activity is in line with the CBD-induced diminished transcription of NF-kBdependent genes, e.g., IL-1b and IL-6 [11].
We and others previously reported that BV-2 cells express CB 1 , CB 2 , GPR55, GPR18 and TRPV2, G protein-coupled receptors and a channel that are known to interact with cannabinoids [20,30,31,[59][60][61]. Using gene array analysis [13,31], we have observed that the relative levels of CB 1 , CB 2 , GPR18 and TRPV2 as well as of the fatty acid amide hydrolase (FAAH; known to cleave the endocannabinoid anandamide) gene transcripts were not significantly affected by the cannabinoid treatments and their levels did not exceed the 2-fold induction or 50% reduction by either CBD or THC treatment. On the other hand, we show here that LPS markedly downregulates CB 2 and GPR55 and that this downregulation is not affected by either CBD or THC pretreatment. This result is in agreement with our previous report showing that LPS markedly downregulates CB 2 and GPR55 mRNAs in BV-2 microglial cells and in microglial primary cultures [59].

Network analysis and signaling pathways
A relationship between CBD-mediated oxidative stress response and glutathione depletion was previously reported ( [13] and references therein). More recently, we showed that CBD-specific gene expression profile in BV-2 cells displays changes normally occurring under either nutrient limiting conditions or proteasome inhibition, and that are attributed to activation of GCN2/eIF2a/ p8/ATF4/CHOP-Trib3 pathway leading to autophagy as well as to apoptotic cell death [13]. The Trib3 gene product seem to be of high importance to the CBD effect due to its ability to serve as a master regulator of an array of pathways including AP-1 [62], ER stress [63], Akt/PKB [64] and NF-kB [65]. Trib3 expression is significantly upregulated by CBD (16-fold) as well as by THC (3.4fold) [13] and as observed here, remains upregulated after LPS treatment (CBD + LPS, 14-fold; THC + LPS, 2.3-fold). According to these gene array studies and the qPCR results, LPS by itself does not significantly affect the expression of Trib3 mRNA. IPA interactome analysis of the microarrays data reveals an interaction between the CBD-upregulated Trib3 and the NF-kB transcription factor pathway (as shown in Figure 2). This interaction seems to be responsible for the attenuation by CBD of the transcription of many proinflammatory genes. There are several indications suggesting interaction between these two pathways. First, a direct interaction between p65/RelA and Trib3 protein which induces inhibition of PKA dependent p65 phosphorylation, was described [65]. Second, Trib3 protein can negatively regulate the serinethreonine kinase Akt/PKB [64], a downstream effector of PI3K that has been implicated in the potentiation of NF-kB-induced transcription of proinflammatory mediators [66]. This negative regulation of Akt activity by the highly induced Trib3 gene product could point to the mechanism for the CBD-mediated regulation of LPS-stimulated gene expression. Indeed, the effect of CBD treatment on a number of LPS-stimulated genes as reported here is reminiscent of the effects described for the PI3K inhibitor (LY294002) and for the NFkB inhibitor (BAY 11-7082) in the murine macrophage cell line RAW264.7 activated with LPS [67]. Both PI3K and NFkB signaling pathways exert important roles in gene expression in response to LPS, but they are not overlapping. Specifically, treatment with CBD repressed a number of typical proinflammatory genes (e.g.; CD40/Tnfrsf5, Tlr2 and Ccl2/Mcp-1) stimulated by LPS, which are known to be NFkB dependent and of other genes including Csf3, Il-1b, Il-1a and Cox2/Ptgs2, which are under the control of both PI3K and NFkB pathways.
Finally, Trib3 was documented to interfere with the inflammatory MAPK signaling via direct interaction with MEK-1 and MKK7 leading to attenuation of AP-1 mediated transcriptional activity in cancer HeLa cells [62]. AP-1 is a transcription factor involved in the regulation of inflammation-mediated cellular functions and has been shown to be inhibited by Nrf2-activating agents [68]. Indeed, our IPA network analysis indicates that the observed decrease in mRNA levels for a number of genes is probably related to a reduction in AP-1 dependent transcription. Additionally, according to these IPA results, this repression is reinforced by combined treatments of CBD and LPS as observed by the induction of FosL1 gene product, another negative regulator of AP-1 (see Figure 3).
Trib3 has been shown to downregulate PPARc transcription and serve as a potent negative regulator of adipocyte differentiation [69] and PPARc is a molecular target for CBD that could be involved in mediating transcriptional effects in BV-2 microglial cells. Indeed, CBD has been shown to bind to PPARc in vitro as well as to activate its transcriptional activity in 3T3L1 fibroblast and in HEK293 transfected cells [70]. In addition, Necela et al., [71] described a regulatory feedback loop in which PPARc represses NF-kB-mediated inflammatory signaling in unstimulated macrophages. Moreover, they show that upon activation of TLR4 in LPS-stimulated macrophages, NF-kB drives down PPARc expression. These results are in agreement with our results showing that LPS highly downregulates the expression of Pparg1 and Pparg2 (Table 5) in BV-2 cells.
The profiles of CBD-induced gene expression with either resting or LPS-activated BV-2 cells, show that CBD stimulates the transcription of several anabolic genes encoding amino acid biosynthetic enzymes, amino acid transporters and aminoacyl-tRNA synthetases known to be activated by ATF4, a basic leucine zipper transcription factor, that is increased when cultured cells are deprived of amino acids or subjected to endoplasmic reticulum (ER) stress ( [13] and references therein). The divergent types of stress (nutrient deficiency, ER stress, hypoxia and oxidative stress) converge on a single event-phosphorylation of the translation initiation factor eIF2a, resulting in a general translational pause followed by selective increase in ATF4 mRNA translation and subsequent stimulation of expression of ATF4 target genes [72]. Many of the CBD-affected transcripts (e.g., Herpud1, Hmox1, Nqo1, Gclm, Gstmu1, Gstmu6 upregulated by CBD and by CBD+LPS treatments) are indeed classified as Nrf2-mediated oxidative stress response genes, including enzymes involved in the biosynthesis of glutathione [13]. Thus, the observed CBD-mediated induction of ATF4-dependent anabolic genes may serve to replenish the amino acids reduced during the elevated turnover of GSH (cysteine, glutamate and glycine). The mechanism underlying CBD action presumably engages generation of ROS which in turn depletes intracellular GSH. Perturbations in redox tone and GSH levels activate the ''phase 2 response'', a mechanism used by cells to mitigate oxidative stress [73]. As we have previously shown, many of the ''phase 2'' gene products are significantly upregulated by CBD [13].
Our present results show that CBD, and less so THC, have immunosuppressive and protective activities that are reminiscent of other clinically applied drugs such as glucocorticoids (GCs; [74]), rexinoids and synthetic triterpenoids [75]. GCs are immunomodulatory agents known to act as suppressive and protective mediators against inflammation. GCs are known to clear antigens by stimulating cell trafficking as well as scavenger systems and matrix metalloproteinases while they stop cellular immune responses by inhibiting antigen presentation and T cell activation [74]. Synthetic oleanane triterpenoids (SO) were shown to be highly effective (much more than the natural triterpenoids) in many in vivo models in the prevention and treatment of cancer and other diseases with an inflammatory component. Molecular targets of SO include KEAP1 (the inhibitor of the transcription factor Nrf2), PPARc, IkB kinase, TGF-b signaling and STAT signaling [75]. SO are among the most potent known inducers of the phase 2 response both in vivo and in vitro and affect the expression of several key cell cycle proteins (e.g., cyclin D1, Cdkn1a/p21 and caveolin 1). In some cancer cells, SO signal through PPARc to inhibit proliferation. The rexinoids bind almost exclusively to the RXRs and are involved in regulation of development, cell proliferation, differentiation and apoptosis [76]. Because RXRs heterodimerize with other receptors (such as RARs, vitamin D receptor, thyroid hormone receptor, PPARa, PPARc and LXRs), rexinoids modulate the actions of many steroid-like molecules that control metabolism and cellular energetics. In view of these results, triterpenoids and rexinoids are defined as multifunctional drugs. Their targets are either regulatory proteins that control the activity of transcription factors (triterpenoids) or transcription factors themselves (rexinoids). These complex modulatory activities exerted by GCs, rexinoids and SO display a panorama of effects that closely resembles the complex actions of CBD.

Conclusions
We found that CBD induces robust responses related to oxidative stress and inflammation and that CBD, and to a lesser extent THC, have immunosuppressive and protective activities. The anti-inflammatory effects are mediated by mainly downregulating the expression of proinflammatory genes but also by upregulating anti-inflammatory mediators. The impact of pretreatment with CBD on LPS-stimulated gene expression was much greater than that of THC.
We show here, that CBD highly upregulates genes encoding negative regulators of NF-kB and AP-1 transcriptional activities (e.g.; Trib3, Dusp1). CBD response to oxidative stress and inflammation is controlled by Nrf2 and ATF4 transcription factors and seem to involve the Nrf2-Hmox1 and the Nrf2/ATF4 systems. The central role of phase 2-response genes and the ATF4dependent Trib3 induction offer an explanatory mechanism by which CBD (but less so THC) lowers the LPS-induced activation of microglial cells and leads to attenuation of transcription of many proinflammatory genes.
All together, our results show that CBD and to a lesser extent THC, are potent modulators of microglial activation and affect several signaling pathways and related networks including induction of cellular stress responses that underlies their high immunosuppressant activity. Figure S1 Functional and pathway annotations according to DAVID Bioinformatics Resources of affected gene products following LPS treatment. Gene Ontology (GO) analysis was performed separately for upregulated (red) and downregulated (green) genes, to identify functional categories and cluster of genes significantly affected by LPS treatment. The GO terms were arbitrarily chosen at various levels within the ontology to avoid redundant allocation. Categories shown are significantly represented at p#0.05. The bottom pie charts represent the cellular pathways according to KEGG database.


Advancements in Medical Statistics Promote Economies of Scale

Modern medical studies, especially those intended for high impact journals, can easily cost millions of dollars (USD). For example, cost of drug development, including phase I to III trials, averages about USD 800 million [1,2]. Studies searching for single nucleotide polymorphisms (SNPs) and copy number variants (CNVs) that contribute to disease are also costly. The full processing cost of an Affymetrix SNP 6.0 microarray, for example, can be $600 [3] multiplied by the 6000 to 9000 patients typical for recent high impact publications of GenomeWide Association Studies (GWAS) [4-9]. Still, researchers aim for even larger sample sizes because clear answers of variants’ relationships with disease are not forthcoming given current sample sizes. For example, exome sequencing studies of 5300 and 6700 patients only yielded either very rare variants or variants with small effect size [10]. More traditional clinical studies, which measure clinical parameters to differentiate between groups of patients, may also have high costs due to recruitment, treatment, testing, and follow-up. They may also have high cost if a large sample size is required to achieve sufficient power for detection of significance in the questions of interest [11,12].

Introduction
Modern medical studies, especially those intended for high impact journals, can easily cost millions of dollars (USD). For example, cost of drug development, including phase I to III trials, averages about USD 800 million [1,2]. Studies searching for single nucleotide polymorphisms (SNPs) and copy number variants (CNVs) that contribute to disease are also costly. The full processing cost of an Affymetrix SNP 6.0 microarray, for example, can be $600 [3] multiplied by the 6000 to 9000 patients typical for recent high impact publications of Genome-Wide Association Studies (GWAS) [4][5][6][7][8][9]. Still, researchers aim for even larger sample sizes because clear answers of variants' relationships with disease are not forthcoming given current sample sizes. For example, exome sequencing studies of 5300 and 6700 patients only yielded either very rare variants or variants with small effect size [10]. More traditional clinical studies, which measure clinical parameters to differentiate between groups of patients, may also have high costs due to recruitment, treatment, testing, and follow-up. They may also have high cost if a large sample size is required to achieve sufficient power for detection of significance in the questions of interest [11,12].
The burden of funding required for these studies strains the resources of even large pharmaceutical companies [4], wealthy corporations [13], and government funded groups, The Cancer Genome Atlas (TCGA) [14], the National Heart, Lung and Blood Institute and Exome Sequencing Project (ESP) [15], who host some of the largest current medical studies. To fund studies of increasing size, it is necessary to implement economies of scale. Economies of scale maintain or lower study cost while increasing the study's number of patients, number of potential analyses, and power. The majority of articles in this special issue increase economies of scale in one or more of these ways.

Recent Advances
Four studies in this issue increase their economy of scale by lowering the cost per patient. Two studies accomplished this by inexpensive data collection. Rogus et al. [16] retrieved data that was obtained from a health system in the normal course of business. In this way, over 600,000 measurements from over 39,000 patients were obtained at low additional cost. Such a large number of patients gave the study large power to detect the difference between case and control groups, which was found to be highly significant. Maetani and Gamel [12] used existing data from the Cancer Institute Hospital in Tokyo, which included surgery, follow-up, and survival of about 3600 patients for up to 50 years. The relatively low cost of obtaining existing data enabled them to obtain sufficient sample size to run analyses on three subgroups by age, 50 subgroups by demographic and prognostic factors, and on patients treated in different years with different regimens. Lifelong follow-up data was also cheaply and readily available, as required for validation of their Boag model [17]. Another two studies lowered their cost per patient by creating new study designs. Matched case-control studies can have high cost due to issues like follow-up costs until the end of recruitment. Sugihara et al. [18] created a study design that avoids this cost. Their dynamic registration method draws from a pool of potential patients while balancing the two groups of patients for all prognostic factors. Erbas et al. [19] avoided the cost of control recruitment and follow-up by implementing a case-crossover design for their asthma study. Given typical recruitment difficulties, the larger the study, the more resources this design will save.
Advancements in medical statistics also increase economies of scale by increasing power and/or the number of feasible analyses without increasing patient size. Meta-analysis increases power with very little increase in cost by combining existing studies into a single dataset of larger patient size. Meta-analysis also makes additional analyses feasible, when in a single study these analyses either cannot be run or cannot reach statistical significance. Sandoval and Zarate [11] list these and six more ways that economy of scale is increased by meta-analysis.

Rogus et al. [16]
, Ketchum et al. [20], and Maetani and Gamel [12] increased the ability of their studies to achieve significant results without increasing cost. The advancement of Rogus et al. [16] was to use Monte Carlo resampling to correct for intra-subject correlation. This correlation kept their highly negative z-statistic from being significant despite the large number of subjects in the study. Maetani and Gamel [12] developed the Boag model in such a way that the expense of longterm follow-up is not necessary to predict the overall survival curve and the long-term effects of the measured variables in a study. Ketchum et al. [20] developed a within-subject normal-mixture model, minimizing costs by making only a few subjects necessary to obtain significant results. This was achieved by incorporating data from multiple tests per subject into the model. The model has accurate estimation with as few as 20 patients. It uses all the data and can include all covariates.
Several studies in this issue increase the number of potential analyses without increasing cost. The study design of Erbas et al. [19] allows more analyses than a case-control design study of the same size. Particularly, they had more power to detect interactions of time-related risk factors with asthma. Halabi [21] developed an adjustment for type I error rate that preserves the rate at the selected alpha level through multiple intermediate analyses. The adjustment also allows flexibility in the allocation of the type I error rate. In this way, the cost of the study is maintained while maximizing the study's ability to detect meaningful changes. Also, the study becomes adaptable to intermediate discoveries without the need to make costly changes to the study or to design an entirely new study. Bergemann et al. [22] combined three test scores over four time points to make global normalized z-scores. These scores have more power, detecting significant differences at all time points when single test scores only detected significance in some time points. The global scores also yield lower p-values than in previous analyses.