Venn and Euler Data Diagrams - Semantic Scholar

0 downloads 205 Views 282KB Size Report
Most existing computer programs that implement Venn and Euler diagrams are limited to ..... 15 seconds on a 2.5 GHz MacB
Venn and Euler Data Diagrams Leland Wilkinson Abstract— Scientists conducting microarray and other experiments use Venn and Euler diagrams to analyze and illustrate their results. Most existing computer programs that implement Venn and Euler diagrams are limited to three sets. The few programs that plot more than three sets rest on ad-hoc methods, such as force-directed algorithms designed for laying out graphs. As a general solution to this problem, this article introduces a statistical model for fitting Venn and Euler diagrams to observed data. The statistical model outlined in this report includes an explicit loss function and a minimization procedure that enables formal estimation of the Venn/Euler model for the first time. A significance test of the null hypothesis is computed for the solution. Residuals from the model are available for inspection. As a result, this algorithm can be used for both exploration and inference on real datasets. A Java program implementing this algorithm is available under the Mozilla Public License. An R function venneuler() is available as a package in CRAN. Index Terms—Visualization, bioinformatics, statistical graphics.

1

I NTRODUCTION

Venn diagrams are collections of n simple closed (Jordan) curves dividing the plane into 2n regions uniquely representing all possible intersections of the interiors and exteriors of the curves [26]. The simplest form of these diagrams displayed in Venn’s original paper and in most applications today involves two or three intersecting circles of constant radius. Figure 1 shows an example. Relaxing the restriction that all possible intersections be represented results in an Euler Diagram [6]. Figure 2 shows an example. Venn and Euler diagrams have had wide use in teaching logic and probability. In almost all of these applications, their use has been confined to two or three circles of equal size. Venn diagrams based on circles do not exist for more than three circles. Higher-order Venn diagrams can be drawn on the plane with nonconvex polygons, but they are difficult to compute for more than a few sets and they are difficult to decode visually [21, 5]. Higher-order circular Euler diagrams exist for many configurations, but they are difficult to compute analytically. Recently, the microarray community has discovered a new use for these diagrams [11, 18, 16, 4]. To reveal overlaps in gene lists, researchers use Venn and Euler diagrams to locate genes induced or repressed above a user-defined threshold. Consistencies across experiments are expected to yield large overlapping areas. Published examples appear to rest on several desirables: 1) for visual simplicity, researchers prefer circles to other convex or nonconvex curves or polygons, 2) an algorithm should produce a Venn diagram for two or three sets if the data can be fit exactly by a Venn diagram, 3) an algorithm should produce an Euler diagram when data can be fit exactly by that model, 4) an algorithm should handle more than three sets, 5) an algorithm should make areas proportional to counts of elements represented by those areas, and 6) an algorithm should handle data generated by a model with a random component. 2

R ELATED W ORK

Existing software programs meet only a few of these requirements. Those generating two- and three-circle diagrams use basic trigonometry to calculate lunes of intersecting circles. [2] introduced an algorithm for producing these diagrams, when they exist, with unequalsized circles based on error-free data. Variants of this algorithm have been used in most applications [22, 19, 7]. Others have used graph layout methods to fit Euler diagrams (when they exist) to more than three sets of error-free data [20, 8]. Although • Leland Wilkinson is Executive VP of Systat Software Inc., Adjunct Professor of Statistics at Northwestern University and Adjunct Professor of Computer Science at University of Illinois at Chicago. E-mail: [email protected]

these approaches have an axiomatic graph-theoretic foundation, the layouts are not unique because they depend on the the choice of graph layout algorithm (force directed, annealing, spectral, etc.). [13, 14] developed an algorithm for the general Euler problem (more than three sets, with circles sized by set cardinality). To deal with area calculations, they use regular polygons instead of circles. Their loss function incorporates multiple parameters to incorporate aesthetic considerations and there is no explicit error term in their model. As such, the method is ad hoc. This paper features an algorithm called venneuler() that produces Venn/Euler diagrams for one or more sets based on a statistical goodness-of-fit function. The advantage of this approach is that data with error can be handled appropriately in a framework that meets all six requirements. By giving the problem a statistical foundation, we place it in a class of solutions that are well-understood. 3

T HE A LGORITHM

The venneuler() algorithm is based on a simple statistical regression model, a method for computing areas of intersections of circles, and a minimization function. We present these in sequence. 3.1

Defining the Model

We begin with a listTof finite data sets X = [X1 , X2 , ..., Xn ] varying in cardinality. Let P = n {X} be a list of all possible intersections of the sets in X, including the void set and the intersections of each Xi with itself. P has m = 2n sets as entries and is ordered as P = [0, / X1 , ..., Xn , X1 ∩ X2 , ..., X1 ∩ X2 ∩ ...Xn ]

(1)

The particular ordering we use for P induces a binary n-bit pattern on each entry of X that we use to index all of our other lists of length m. In other words, each intersection structure in X is uniquely indexed by a length-n binary string that we can use to map entries of X to entries of P [10]. Let P− = Dis joint(P), where the Dis joint function produces disjoint entries. It does so by beginning at the next-to-the-last entry of the list in P and iterating to the first entry while removing duplicate elements found in sets later in the list. The result is a list of sets in which the last entry of P− contains a set whose elements are unique to the highest-order intersection. Entries earlier in P− contain remaining elements that are unique after conditioning on entries later in the list. The entries in P− thus correspond to all the possible disjoint intersections of elements in X. Next, we construct a list of disks, D = [D1 , . . . , Dn ], each having raS dius ri = |Xi |/| {X}|. Each disk Di is centered on coordinates (xi , yi ). From D, we construct the corresponding list Q = [0, / D1 , ..., Dn , D1 ∩ D2 , ..., D1 ∩ D2 ∩ ...Dn ]

(2)

We then apply the same disjoint operation we used on P in order to produce Q− = Dis joint(Q). We now have a one-to-one correspondence between the entries of Q− (disjoint disk intersections) and the entries of P− (disjoint dataset intersections). Both are indexed by the same list of binary strings. From P− and Q− we make a vector c = (|P1− |, . . . , |Pm− |) consisting of the counts of elements in each disjoint intersection of the sets in − X and a vector a = Area(Q− 1 , . . . , Qm ) consisting of the areas of the disjoint disk intersections. Given these entities, a Venn diagram with areas proportional to counts is defined by the equation a = βc

(3)

Since we allow disk sizes to vary, there may not exist a set of coordinates (xi , yi ) for which this equation is satisfied. Moreover, we will assume the elements in the data sets Xi are generated by a process having a random component. Our model is therefore a = βc+ε

(4)

where ε is a random variable with zero expected value. Our ordinary least-squares estimate of β in this case is

3.4

Our remaining task is to find the coordinates (xi , yi ) that minimize the sum of squared residuals from data fit by Equation 4. We use the method of steepest descent with the gradient calculated from the residuals. For each disk Di , the descent step on each iteration is proportional to m

5F(x, y)i =

3.2

(5)

Computing Areas

For two or three circles, analytic computation of areas in Q− is straightforward [2]. With more than three, computations become prohibitively complex. To deal with this problem, [13, 14] worked with regular polygons instead of circles and employed standard polygon intersection algorithms. This method is not only expensive, but it also fails to deal directly with the circles that researchers want to use. A simple method for solving this problem is based on numerical quadrature and binary indexing. In order to compute areas on the entries of Q− , we “draw” circles on n binary grids (bit planes), each of resolution p × p. Each “pixel” in a bit plane has the value 1 if it is inside a circle and 0 if not. The string of 1’s and 0’s derived from passing through the corresponding pixel on each bit plane yields the same binary indexing that we use for Q− itself. We simply sum the result over all pixels to get intersection areas. The method is very fast. Since we need to run through n such grids to detect which entries of Q− are indexed by each cell in the grid, the complexity of this computation is proportional to n. In practice, we find that p = 200 is sufficient resolution to allow the iterations to converge. 3.3

Initial Circle Locations

The above algorithm will work acceptably with random starting locations for the circles. It is more efficient, however, to begin with a rational starting configuration. To accomplish this, we adopt an approach from classical multidimensional scaling [25]. We compute the similarity matrix S : si j = |Xi ∩ X j |. We convert this to a quasi-distance matrix D by negating the similarities and adding a constant to make the quasi-distance matrix positive definite [17]. We then choose an arbitrary row (col) in D and compute a matrix of scalar products on the distances conditioned on this row k. The resulting matrix is Bk : bi jk = dik d jk cos θik j ,

(6)

2 cos θik j = (dik + d 2jk − di2j )/2

(7)

where

We then compute the singular value decomposition Bk = UVU0 . The starting coordinates (xi , yi ) are found in the rows of the first two columns of U.

∑ ∑ {(xi − x j )(ak − aˆk ), (yi − y j )(ak − aˆk )},

(8)

k=1 i6= j

where Qk ∈ Qi . The venneuler() program produces a stress statistic at convergence. Stress is defined as SSE/SST (residual sum of squares divided by total sum of squares around zero). A correlation coefficient can be √ computed as r = 1 − stress. As with any regression without a constant, this correlation can differ from the Pearson. To check the extent of this difference, we added a constant term to the model and fit it to real datasets. As expected, the constant term was negligible, so interpreting the correlation as an ordinary goodness-of-fit statistic should be reasonable in practice. Theory implies null sets should have zero area, so we omit the constant. 4

βˆ = a0 c/c0 c

Minimizing Loss

T HE D ISTRIBUTION

OF THE

S TRESS S TATISTIC

We computed a Monte Carlo simulation to estimate the distribution of our stress statistic. For each number of circles (n = 3, . . . 10), we generated 100 simulations. For each simulation, we generated 2n uniform random numbers to represent the areas based on the entries in Q− . We ran venneuler() on the random data and computed order statistics on the resulting stress values. For the empirical stress fractiles s.01 and s.05 , we fit the logistic function sα = expb(n−c) /(1 + expb(n−c )

(9)

The fit for both equations was extremely close (r2 > .99). For s.01 , c = 6.105, b = 0.909 and for s.05 , c = 5.129, b = 0.900. Table 1 shows the critical values for n = 3, . . . , 10. These stress values are substantially higher than corresponding critical stress values in the MDS literature, assuming n represents the number of points [24, 15, 23, 3, 1]. The venneuler() model is much more constrained than the MDS model, however. Not only are all possible pairwise intersections included in the loss function, but also all higher-order intersections are included. Moving points (disk centers) around on the plane affects 2n−1 areas rather than n(n − 1)/2 distances as in MDS. Furthermore, the regression function on which loss is based has a zero intercept; MDS ordinarily includes an intercept parameter. These restrictions are an unusual advantage of the venneuler() model. It is so highly constrained that researchers will want to consider whether the model is appropriate for their data before they waste time trying to fit the model.

Table 1. Critical stress values for venneuler()

n

s.01

s.05

3 4 5 6 7 8 9 10

0.056 0.129 0.268 0.476 0.693 0.848 0.933 0.972

0.128 0.266 0.471 0.687 0.843 0.930 0.970 0.988

5

E XAMPLES B

Figure 1 shows a 3-ring Venn diagram produced by the input: venneuler(A = {a, ab, ac, abc}, B = {b, ab, bc, abc}, C = {c, ac, bc, abc}) The stress value for this solution is 0.103, consistent with the fact that a perfect area-based solution for 3 equal-sized circles does not exist even though the ordinary Venn diagram requirement is met [2, 9]. Nevertheless, this is as close to equal-area as we can get; Figure 1 resembles aesthetically the canonical “Ballantine” charts appearing in Venn diagram tutorials [21].

A

C

D

A

C

Fig. 2. Four-ring diagram produced by venneuler() on incomplete Venn structure.

BD AC

-10

-8

BCD ABD ABC

ACD

-6

-4

AB AD CD BC

-2 0 2 RESIDUAL

B

A C

D

4

ABCD

6

8

10

Fig. 3. Residuals for four-ring diagram in Figure 2.

B

venneuler() program uses an ampersand to represent intersection.

Fig. 1. Three-ring Venn diagram produced by venneuler() on complete Venn structure

Figure 2 shows a 4-ring diagram produced by venneuler() on data that lack some 3-way intersections: venneuler(A = {a, ab, ac, ad, abc, abd, acd, abcd}, B = {b, ab, bc, bd, abc, abd, bcd, abcd}, C = {c, ac, bc, cd, abc, acd, bcd, abcd}, D = {d, ad, bd, cd, abd, acd, bcd, abcd}) The resulting diagram is neither Venn nor Euler because two 2-way intersections (A ∩ C and B ∩ D) are missing in the plot. It nevertheless approximates the Euler diagram for this set in the way we would expect. There is a trade-off between moving the circles outward to eliminate the 3-way areas and moving them inward to represent the 4-way area. The stress for this solution is .30. Figure 3 contains a residual plot from the solution in Figure 2. This plot reveals the trade-off the venneuler() algorithm made. The two largest residuals show that AC and BD are under-predicted. This happens because the 3-way intersections (for which there are no data values) are stealing area from these 2-way intersections (for which there are data values). Residual plots are a natural complement to graphics produced by venneuler(). They are useful for diagnosing the statistical properties of the solution, a feature unavailable in ad-hoc algorithms. Figure 4 shows a 6-ring Euler diagram produced from areas calculated from a pre-existing diagram. In this version of the syntax, we directly input areas for each disjoint subset. For this type of input, the

venneuler(A = 4, B = 6,C = 3, D = 2, E = 7, F = 3, A&B = 2, A&F = 2, B&C = 2, B&D = 1, B&F = 2,C&D = 1, D&E = 1, E&F = 1, A&B&F = 1, B&C&D = 1) The stress for this solution is .005 and the reproduced diagram matches the original, confirming the ability of venneuler() to capture a moderately complex error-free structure. Figure 5 shows two diagrams for data shown in Figure 1 of [12]: venneuler(SE = 13, Treat = 28, Anti-CCP = 101, DAS28 = 91, SE&Treat = 1, SE&DAS28 = 14, Treat&Anti-CCP = 6, SE&Anti-CCP&DAS28 = 1) These diagrams depict the overlap of genes detected in four different populations. The left panel is from the original article. The original graphic shows four circles, so the sets cannot be represented by a circular Venn diagram, Nevertheless, the Euler diagram in the right panel computed by venneuler() accurately represents the data. The stress for this solution is .0002. Figure 6 shows an Euler diagram for 11 animals based on gene lists downloaded from the Agilent DNA oligo microarray database (http://www.chem.agilent.com). The analysis was based on 247,412 gene symbols and 11 animal names. The stress for this solution is .06, with corresponding correlation of .97. Many genes in these lists have yet to be classified. When this task is completed, we would expect that the Human genome would overlap more substantially with Mouse and Rat. The venneuler() solution provides an accurate portrait of this work in progress. The computation was completed in under 15 seconds on a 2.5 GHz MacBook Pro running the Java 1.5 Virtual Machine in 2GB of allocated memory. Figure 7 shows an Euler diagram for six works of English fiction. We downloaded files from the Project Gutenberg Web

site (http://www.gutenberg.org/wiki/Main Page). Stop words (a, and, the, ...) were filtered and a list of distinct words was constructed for each corpus. The combined lists (totaling 65,432 words) were submitted to venneuler(). The stress for this solution is .02, with a corresponding correlation of .99. The size of the circles is based on the number of unique words in each book. Ulysses, the King James translation of the Bible, and Moby Dick anchor the configuration; they contain the lion’s share of unique words. Ulysses is notable for its large number of unique words – a familiar aspect to anyone struggling to read that novel. The Bible has a smaller number of unique words; many of these are proper names not shared by the other literature. Not surprisingly, Shakespeare’s language in Macbeth is almost a complete subset of its contemporary, the Bible. Nevertheless, despite its relatively short length, Macbeth contains almost as many unique words as Silas Marner and The Picture of Dorian Gray. Shakespeare embedded a rich vocabulary in a small package.

E

F D

A C

B

Bible

Fig. 4. Six-ring Euler diagram almost perfectly fit to artificial dataset.

(a)

(b) Silas Marner

Ulysses

Anti-CCP

Macbeth

Dorian Gray

DAS28

Moby Dick

SE Treat

Fig. 5. Four-ring diagram of data underlying Figure 1 of [12]. Panel (a) is reproduced from the original article. Panel (b) is the venneuler() solution. The correlation between the areas of the circles and their intersections is .99.

Frog

Fig. 7. Euler diagram for six works of English literature. The sizes of the circles are proportional to the number of unique words in each of the works. The areas of the intersections are proportional to the number of shared words among the works. The correlation between the areas of the circles and their intersections and the word counts is .99.

Sheep Rabbit Horse

Cow

Pig Human Mouse

Dog Chicken Rat

Fig. 6. Euler diagram for 11 gene lists containing 247,412 unique genes from the Agilent DNA oligo microarray database. The sizes of the circles and intersections are due to the number of genes listed for each animal in the database as opposed to the count of the genome itself. The correlation between the areas of the circles and their intersections and the gene counts is .97.

6 C ONCLUSION The algorithm described in this report provides for the first time a statistical basis for estimating Venn and Euler diagrams on real data. Its distinguishing features include an explicit statistical loss function that accommodates data with error, the ability to evaluate probabilistically the goodness-of-fit of a solution, the ability to represent counts proportionally by areas, and the ability to accommodate unconnected sets. The consequences of this design are several. The venneuler() program produces a classic Venn diagram for input data sets exactly representable by a Venn diagram, as Figure 1 shows. It does the same for data sets exactly representable by an Euler diagram, as Figure 4 shows. For data sets not exactly representable by a Venn or Euler diagram, venneuler() achieves an appropriate compromise, as Figure 2 shows. Finally, for datasets possibly containing error, venneuler() minimizes a simple loss function that fits the requirements and understanding of researchers who use these diagrams, as Figure 5 shows. Given the highly-constrained nature of the Venn/Euler model, it is remarkable to find fits on real data as good as the one in Figure 6. The bioinformatics community appears to have adopted a model that reflects the underlying genetic structures they are investigating. Now there is a way to fit it. The venneuler() program is written in Java and available as a function in R under the Mozilla Public License.

ACKNOWLEDGMENTS Simon Urbanek (AT&T Labs) installed the Java code in CRAN. Funding was provided by NSF/DHS grant DMS-FODAVA-0808860 R EFERENCES [1] Ingwer Borg and Patrick J. F. Groenen. Modern Multidimensional Scaling: Theory and Applications. Springer, 2005. [2] S. Chow and F. Ruskey. Drawing area-proportional Venn and Euler diagrams. In Graph Drawing, volume 2912 of Lecture Notes in Computer Science, pages 466–477. Springer, 2004. [3] Trevor F. Cox and Michael A. A. Cox. Multidimensional Scaling, Second Edition. Chapman & Hall/CRC, 2000. [4] Jos´e R. Dinneny, Terri A. Long, Jean Y. Wang, Jee W. Jung andmDaniel Mace, Solomon Pointer, Christa Barron, Siobhan M. Brady, John Schiefelbein, and Philip N. Benfey. Cell identity mediates the response of arabidopsis roots to abiotic stress. Science, 320:942–945, 2008. [5] A.W.F. Edwards. Cogwheels of the Mind: The Story of Venn Diagrams. Johns Hopkins University Press, Baltimore, 2004. [6] L. Euler. Lettres a Une Princesse d’Allemagne, volume 2. 1761. [7] H. Fang, S.C. Harris, Z. Su, M. Chen, F. Qian, L. Shi, R. Perkins, and W. Tong. Arraytrack: An FDA and public genomic tool. Methods in Molecular Biology, 563:379–398, 2009. [8] J. Flower, A. Fish, and J. Howse. Euler diagram generation. Journal of Visual Languages & Computing, 19:675–694, 2008. [9] P. Hamburger and R. E. Pippert. Venn said it couldn’t be done. Mathematics Magazine, 73:105–110, 2000. [10] Peter Hamburger. Cogwheels of the mind: The story of Venn diagrams, a review. Mathematical Intelligencer, pages 36–38, 2005. [11] Natalia B. Ivanova, John T. Dimos, Christoph Schaniel, Jason A. Hackney, Kateri A. Moore, and Ihor R. Lemischka. A stem cell molecular signature. Science, 298:601–604, 2002. [12] Cristina Moraes M. Junta, Paula Sandrin-Garcia, Ana Lucia L. FachinSaltoratto, Stephano Spano S. Mello, Rene D. R. Oliveira, Diane Meyre M. Rassi, Silvana Giuliatti, Elza Tiemi T. Sakamoto-Hojo, Paulo Louzada-Junior, Eduardo Antonio A. Donadi, and Geraldo A. S. Passos. Differential gene expression of peripheral blood mononuclear cells from rheumatoid arthritis patients may discriminate immunogenetic, pathogenic and treatment features. Immunology, 127:365–372, 2008. [13] H. A. Kestler, A. Mller, T. M. Gress, and M. Buchholz. Generalized Venn diagrams: a new method of visualizing complex genetic set relations. Bioinformatics, 21:1592–1595, 2005. [14] Hans A. Kestler, Andr Mller, Johann M. Kraus, Malte Buchholz, Thomas M. Gress, Hongfang Liu, David W. Kane, Barry Zeeberg, and John N. Weinstein. VennMaster: Area-proportional Euler diagrams for functional GO analysis of microarrays. BMC Bioinformatics, 9, 2008. [15] D. Klahr. A Monte Carlo investigation of the statistical significance of Kruskal’s nonmetric scaling procedure. Psychometrika, 34:319–330, 1969. [16] Cheng Lu, Shivakundan Singh Tej, Shujun Luo, Christian D. Haudenschild, Blake C. Meyers, and Pamela J. Green. Elucidation of the small RNA component of the transcriptome. Science, 309:1567–1569, 2005. [17] S. M. Messick and R. P. Abelson. The additive constant problem in multidimensional scaling. Psychometrika, 21:1–15, 1956. [18] Duncan T. Odom, Nora Zizlsperger, D. Benjamin Gordon, George W. Bell, Nicola J. Rinaldi, Heather L. Murray, Tom L. Volkert, J˙org Schreiber, P. Alexander Rolfe, David K. Gifford, Ernest Fraenkel, Graeme I. Bell, and Richard A. Young. Control of pancreas and liver gene expression by HNF transcription factors. Science, 303:1378–1381, 2004. [19] Mehdi Pirooznia, Vijayaraj Nagarajan, and Youping Deng. GeneVenn - a web application for comparing gene lists using Venn diagrams. Bioinformation, 1:420–422, 2007. [20] P. Rodgers, L.S. Zhang, G. Stapleton, and A. Fish. Embedding wellformed Euler diagrams. In IEEE International Conference on Information Visualisation, pages 585–593, 2008. [21] F. Ruskey. A survey of Venn diagrams. Electronic Journal of Combinatorics, 4, 1997. [22] G. K. Smyth. Limma: linear models for microarray data. In R. Gentleman, V. Carey, S. Dudoit, R. Irizarry, and W. Huber, editors, Bioinformatics and Computational Biology Solutions using R and Bioconductor, pages 397–420. Springer, New York, 2005.

[23] I. Spence and J.C. Ogilvie. A table of expected values for random rankings in nonmetric multidimensional scaling. Multivariate Behavioral Research, 8:511–517, 1973. [24] H.H. Stenson and R.L. Knoll. Goodness of fit for random rankings in Kruskal’s nonmetric scaling procedure. Psychological Bulletin, 71:122– 126, 1969. [25] Warren Torgerson. Multidimensional scaling: I. theory and method. Psychometrika, 17(4):401–419, December 1952. [26] J. Venn. On the diagrammatic and mechanical representation of propositions and reasonings. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 9:1–18, 1880.