Massively parallel sequencing data
R mRNA editing was measured in the specimens obtained from the human dorsolateral prefrontal cortex and rat spinal cord. The 101 human subjects comprised 45 individuals diagnosed with major depression, and 56 normal controls 
. The 19 rats comprised 7 controls and 12 rats whose spinal cord was transacted six weeks prior to the data collection 
. In these rat specimens, the mRNA levels are assumed to be unaffected by the transaction, being collected from a region above it. Overall, the analysed data included 56,690,398 human reads (an average of 561,291 per subject) and 5,659,108 rat reads (an average of 297,848 per rat) (Supplementary Table S1
Each measurement (mRNA molecule) is represented by a binary vector indicating the editing states of the five sites A, B, E, C, and D. For example, a measurement in which sites A, B, and D are edited but E and C are not is represented by the binary pattern 11001. For a collection of measurements, we denote the editing pattern
as the vector
is the number of binary vectors whose decimal representation is
First, we tested whether the editing patterns of all human normal controls were statistically indistinguishable from the editing patterns of all subjects with major depression. To this end, we conducted a conservative randomization test, whereby the
-test statistic was repeatedly computed on modified data. In each repetition, we randomly assigned subjects as normal or as depressed, keeping the total number of normal controls and the total number of subjects with major depression fixed. For each repetition, we computed the test statistic of the
are the editing patterns of normal and depressed samples, respectively, and
are the total number of measurements from normal controls and from subjects with major depression, respectively. This procedure was repeated 106
times, and the p-value of the test was computed as the number of random test statistics that were larger than the true test statistic. A similar procedure was used to compare normal rats with transacted ones. We found that the editing pattern in normal human controls was indistinguishable from the editing pattern in subjects with major depression (P
0.80), and that the editing pattern in normal rats was indistinguishable from that in transacted rats (P
0.65). This result justifies pooling together all human subjects and all rats for further analysis. Using a similar randomization procedure, we found that the editing pattern in humans is very different from that in rats (P
) (Supplementary Figure S1
Finding the probability model with the best fit to the data
Clustering, by nature, identifies groups of associated sites. However, to obtain finer resolution of the relationship between the sites, we resorted to more elaborate methods. The ultimate description of the dependency between the editing sites would be their joint probability distribution. For five editing sites, there are 8,782 possible joint distribution functions. We enumerated all the 8,782 functions, and ranked them according to how well they fit the data using both maximum-likelihood and Bayesian inference (see Methods
In both human and rat, and for both maximum-likelihood and Bayesian inference, the best model was the maximally-dependent joint probability distribution,
. We graphically represent probability models by a pDAG (partial Directed Acyclic Graph), which is a Bayesian network containing a mixture of directed and undirected edges (see Methods
). The pDAG of this maximally-dependent model is simply the fully connected undirected graph (, ). This result is consistent with our previous finding that all pairs of sites are significantly dependent.
Figure 2 BIC scores (dots) and pDAGs of human models through .
Figure 3 BIC scores (dots) and pDAGs of rat models through .
Such result is expected given the large size of the data. In order to find which edges in the graph are more strongly supported by the data, we divided all the 8,782 probability models into 11 groups according to the number of edges in the corresponding pDAG. The first group consists of all models with zero edges (which is simply the single model
), the second group consists of all (ten) models with one edge, etc. Then, we computed the best-fitting probabilistic model within each group. Hereinafter,
denotes the best-fitting model from within the group of models with
edges. The results for the maximum-likelihood Bayesian Information Criterion (BIC) score (see Methods
) in human are shown in . Adding an edge to a model always improves its score,
, but the improvement becomes smaller as
increases. The estimated parameters of each of the best-fitting models are given in Supplementary Table S3
To further explore the relative impact of the different edges, we ranked the edges by the order in which they first appear in the sequence of models
. Specifically, the rank of an edge is the smallest integer
contains this edge (). Importantly, an edge in
need not necessarily be included in
, which is the reason why two edges – (B,E) and (E,C) – have the same rank, and why no new edge appeared in
. To account for the possibility that edges may disappear or reappear as the number of edges grows, we define for each edge its support
. The support of an edge with rank
is the fraction of the models
that contain this edge (). Clearly, the higher the support, the more confident we are that the edge has a unique contribution to making the respective model better fitting the data. The edge (A,B) is the first to appear (
), and has a full support (it appears in all the models
), indicating that the dependence between A and B is obviously the strongest among all pairs of sites, in accord with the findings described above. Next appear edges (B,D) and (A,C) that both also have full support. This observation is consistent with the clustering analysis results () but provides more detail on the interdependencies among A, B, and C, D. The next edge to appear is (A,E), but it does not have full support which lowers our confidence in its unique contribution to the score of the best model. The edges that appear in models
– (C,D) and (E,D) –make (at least qualitatively) only marginal contributions to the score. Repeating the analysis with the maximum-likelihood Akaike Information Criterion (AIC) scores, or with Bayesian scores, gave the same series of best-fitting models
(Supplementary Figure S3
Ranking of the strength of dependency of each pair of editing sites in human 5-HT2CR mRNA.
We conducted the same analysis for the rat data. Here, too, adding edges kept improving the BIC score of the model (). The estimated parameters of the best-fitting models are given in Supplementary Table S4
. For rat, all edges have full support, which means that an edge with rank
appears in all the models
(). In rat, the two edges with the lowest rank – (A,B) and (B,D) – have the same rank as in human. However, the edge with
in rat is (B,C) as opposed to (A,C) in the equivalent human
model. Similarly, the edge with
is (A,E) in human, but it is (B,E) in rat. This suggests that the central role of site A in governing the editing state of sites E, C, and D in human is taken by site B in rat. Indeed, referring collectively to site A or site B as F, human and rat show very similar edge rankings (compare and ). Repeating the analysis with Bayes scores yielded identical series of best-fitting models for rat (Supplementary Figure S4b
). Using AIC scores produced only a single difference, in model
. The edge (C,D), which is present in the BIC and Bayes scores, was replaced by the edge (E,C) for the AIC score (Supplementary Figures S4a and S5). However, as we have seen, these edges anyway have marginal contribution to the best-fitting model. On the whole, the information contribution of additional edges dropped much faster for the rat data than it did for the human data (compare and ), suggestive of a more complex pattern of dependencies among editing sites and accordingly more subtle regulation of the editing process in human brain.
Ranking of the strength of dependence of each pair of editing sites in rat 5-HT2CR mRNA.
The above analysis lacks measure of score variance, thus hindering quantitative evaluation of the significance of each edge to the total score. To overcome this, we repeated the analysis for each individual separately, for both human and rat. In this way, each individual
provides its own sequence of best fitting models
, and for each number of edges
), there is now a sequence
of best models, where
is the total number of individuals. For a certain
, such that
is supported by
individuals. Let us further assume that we have sorted the sequence
according to the level of support, such that
. Next, we define a set
models that are equally supported by the different individuals (
). To this end we make a Bonferroni-corrected series of proportion tests, asking whether
is supported significantly more than the other best-fitting models
is the first model whose support is significantly lower than that of
. The results for the BIC scores in human at significance level 0.05 are given in . The results for the AIC and Bayes scores are similar, and are given in Supplementary Tables S5
. As an example, in the BIC score analysis, out of 101 individuals 78 (77.2%) support the single-edge best-fitting model (A→B) (). The second-supported model
is supported by 21 individuals (20.8%), which is significantly lower than the support for
and so in this case
and the set of best fitting models is simply (
). As another example,
is supported by 14 individuals (13.9%), but this level of support is not statistically different from the support by 3 individuals (3.0%) of the model
, and in this case
. Overall, there is a good agreement between this individual-based analysis and the pooled analysis. The pooled best-fitting model for each
is marked by asterisk in and Supplementary Tables S5
, and it is always within the group of models that are equally supported by the different individuals. Very similar results had been obtained for rat (Supplementary Tables S7
). Here too, the pooled best-fitting model for each
is always within the group of models that are equally supported by the different individuals.
Statistics on the individual best-models for BIC scores in human.
The individual-based approach can be used not only to re-evaluate the support for the different graphical models but also to perform an edge-by-edge analysis. To this end, we can look at each edge, and count how many times (in either direction) it appears in the sequence
. These counts are binomial random variables, so if an edge appears, overall, in the best-fitting model of
individuals out of a total of
individuals, its variance is
. The support of each edge for any
in human is given in . Consider, for example,
. Overall, the 101 individuals support
different best-fitting models. Yet, the edge (A,B) appears in all of them, and thus is supported by all the individuals. The edge (B,D) is supported by 98 individuals, or by 97% of the best-fitting models. For each
, we can take the first
most-supported edges as the basic set
of edges in the
model. Then, we can check how unique is this set of edges by testing (using proportion test) whether the e'th
supported edge is significantly more supported than the next edges (). From and we see that the first three edges (A,B), (B,D), and (A,C) are clearly more supported than all other edges, in this order. However, the next edges (B,E), (E,C), and (A,D) all have approximately the same support and no one is more significant than the others. The results are almost identical when using AIC or Bayes scores (Supplementary Figures S6 and S7).
Level of support of each edge in all individual best-fitting models with fixed number of edges.
Edges present in the best-fitting models in human (BIC scores).
Similar analysis for rat shows, in accord with our previous results, a more hierarchical relationship between the edges ( and ). Here, the order of importance is clear for the first six edges: (A,B), (B,D), (B,C), (B,E), (A,C), and (A,D). The edges (C,D), (E,D), and (E,C) all have approximately the same support and no one is more significant than the others. The results are almost identical when using AIC or Bayes scores (Supplementary Figures S8 and S9).
Level of support of each edge in all individual best-fitting models with fixed number of edges.
Edges present in the best-fitting models in rat (BIC scores).