- Open Access
MIDAS: software for analysis and visualisation of interallelic disequilibrium between multiallelic markers
BMC Bioinformaticsvolume 7, Article number: 227 (2006)
Various software tools are available for the display of pairwise linkage disequilibrium across multiple single nucleotide polymorphisms. The HapMap project also presents these graphics within their website. However, these approaches are limited in their use of data from multiallelic markers and provide limited information in a graphical form.
We have developed a software package (MIDAS – Multiallelic Interallelic Disequilibrium Analysis Software) for the estimation and graphical display of interallelic linkage disequilibrium. Linkage disequilibrium is analysed for each allelic combination (of one allele from each of two loci), between all pairwise combinations of any type of multiallelic loci in a contig (or any set) of many loci (including single nucleotide polymorphisms, microsatellites, minisatellites and haplotypes). Data are presented graphically in a novel and informative way, and can also be exported in tabular form for other analyses. This approach facilitates visualisation of patterns of linkage disequilibrium across genomic regions, analysis of the relationships between different alleles of multiallelic markers and inferences about patterns of evolution and selection.
MIDAS is a linkage disequilibrium analysis program with a comprehensive graphical user interface providing novel views of patterns of linkage disequilibrium between all types of multiallelic and biallelic markers.
Gametic disequilibrium (widely known as linkage disequilibrium or LD) is a genetic phenomenon which occurs when alleles at different loci are non-randomly associated in a given population. This correlation between polymorphisms is caused and/or influenced by their shared history of mutation and recombination, and by many other factors including genetic drift, population growth, admixture or migration, population structure, the ages of the polymorphisms, the physical distance separating them and the effects of selective pressure . The characterization of LD is an important issue in both evolutionary and medical genetics, since it is informative in association mapping of trait or disease loci, and an indicator of the interaction between genes, the relative influence of different evolutionary forces in the generation/disruption of genetic variability, and the genetic history of populations .
The theory of estimation of LD has been substantially developed in recent years. Relevant advances have been made in the knowledge of the properties of LD coefficients and LD statistical tests, which are used respectively to measure the magnitude and to estimate the significance of LD. LD is said to exist when the frequency of a haplotype observed in a population sample is significantly greater or lesser than the frequency expected from the product of the allele frequencies, the magnitude of LD correlating with such difference. There are a variety of measures and statistical tests available for the estimation of LD (D', ρ, r, r2, d, d2, and chi-square and Fisher exact tests, being the most used LD coefficients and statistical tests), and many programs exist for that purpose (including Haploview, 2LD, Arlequin, GDA, DNAsp, ALLASS, DISEQ, DMAP, etc., reviewed in  and ). Some software, such as GOLD , GOLDsurfer  and Haploview , also include graphical displays enabling quick overviews of large regions. However, most packages are intended for use with single-nucleotide polymorphism (SNP) data in a pairwise fashion. This focus on biallelic markers makes both LD estimation and graphical representation straightforward compared with multiallelic markers such as microsatellites.
The analysis of LD between a pair of multiallelic loci represents a conceptual difference in relation to the analysis of LD between a pair of biallelic loci. In both instances, LD can be analysed at two different levels. One is the overall LD between the pair of loci, and the other is the interallelic LD between each of the alleles at the first locus and each of the alleles at the second one.
The magnitude and the significance of both overall and interallelic LD are the same for pairwise analyses involving two biallelic loci. This does not apply, however, for LD between multiallelic loci. Given a pair of multiallelic loci with k and l alleles respectively, there are k × l possible interallelic associations. In theory, pairwise combinations of alleles at different loci can differ in parameters such as magnitude, significance and patterns of LD. This has been confirmed empirically in the characterization of interallelic LD between pairs of dinucleotide repeat loci spanning human chromosome 11p15 , and in the analysis of LD between the TH01 microsatellite and IGF2 SNP haplotypes in the context of the identification of microsatellite loci tagging haplotypes relevant to association mapping of complex disease traits [8, 9]. The analysis of interallelic associations is therefore necessary for a complete description of LD between multiallelic loci .
Despite the existence of alternative estimation theory [2, 10, 11], LD between multiallelic loci is often estimated by pooling alleles into two groups in order to reduce the system to a two-allele two-locus model. This approach does not allow the analysis of all possible interallelic associations. In contrast, it reduces the LD between multiallelic loci to a single estimate of overall LD. It has been shown that the overall measure obtained by pooling alleles of multiallelic loci tends to underestimate LD, may complicate discrimination among the evolutionary forces generating LD in populations, and may decrease the success of association mapping of trait or disease loci ( and references therein). In addition to the number of alleles, the magnitude and the power to detect LD depend on other factors, including the sample size, the statistical tests and coefficients used, the allele frequencies and the sign of the association . This latter issue has been shown to be of special importance. A sign-based LD estimation method recently developed for multiallelic systems , has been shown to considerably increase both the statistical power and the accuracy of estimation of the intensity of LD [2, 13]. On the other hand, the task of presenting a graphical overview of interallelic disequilibrium between alleles of multiallelic markers is rather more challenging than for biallelic markers (with colour intensity indicating the magnitude of linkage disequilibrium between a pair of markers) and has not been previously attempted.
In this work, we have developed an integrated LD analysis software (MIDAS: Multiallelic Interallelic Disequilibrium Analysis Software) that computes interallelic LD from genotypic data incorporating the latest advances in the theory of estimation of LD, and represents graphically the intensity and significance of pairwise non-random associations between any combination of microsatellites, SNPs, haplotypes or other multi-allelic markers.
MIDAS was written in the Python programming language v2.4 , using the Tkinter module for generating a graphical user interface (GUI). The Tkinter "Canvas" widget was used for plotting of graphical data, whilst other Tkinter widgets were used for creation of menus, buttons and other aspects of the interface. All modules used were part of the standard Python distribution  and include: Tkinter, tkFileDialog, math, copy, cPickle, os and webbrowser. The program reads and writes standard tab-delimited text files, and has an additional option to save a binary analysis file (using the cPickle module) which stores all variables and allows the user to reload a previous analysis.
Figure 1 shows a flow-chart representing the program structure. Data (raw genotypes, with marker IDs and positions) are imported, and then the user selects analysis (this is a separate step to enable incorporation of different analyses in future versions). Analysis begins with an assessment of Hardy-Weinberg equilibrium (HWE) as previously described [15, 16]. Markers out of HWE are flagged for highlighting in the final outputs. The next step is estimation of LD. Finally, the results of the analysis are plotted (figure 2).
The program has been designed for simple installation and use by any computer user, and requires only the prior installation of the standard Python distribution  to function on a Microsoft® Windows® 2000/XP computer. Operation is mouse and menu-driven with optional hotkeys for scroll and zoom. Input files can be prepared in most spreadsheet programs and exported as tab-text. Results output is tab-text format and can be imported into most spreadsheet programs.
All parts of the program were scripted de novo, but the algorithm for LD calculation was based on previous programs developed by two of the authors (CZ and SR).
Estimation of LD
Given two multiallelic loci, A and B, we estimated the LD for each pair of alleles defining a two-locus haplotype. The accurate computation of all possible interallelic associations requires that each of the two-locus haplotypes defining an interallelic combination represents only the observed and expected counts for the pair of alleles under consideration. This is not attained when alleles (and therefore haplotype counts) are pooled arbitrarily. MIDAS computes interallelic disequilibrium between multiallelic loci following an approach previously described  which avoids both losing and pooling interallelic information. Both this approach and its underlying theory have been applied and discussed in detail previously ( and references therein). In brief, if locus A has k alleles A i (i = 1,......, k) and locus B has l alleles B j (j = 1,......, l), then the complete array of possible two-locus haplotypes was partitioned into k × l separate 2 × 2 contingency tables by collapsing the data into A i vs. not-A i () at the A locus, and B j vs. not-B j () at the B locus. Estimates of two-locus haplotype frequencies were obtained from genotype data by the Hill method, an expectation-maximisation (EM) algorithm . The magnitude of disequilibrium between pairs of alleles at different loci was measured by = D ij /D max , where D ij = X ij - p i q j , p i and q j are the frequencies of alleles i and j, respectively, X ij is the observed frequency of the haplotype A i B j and D max = min [p i (1 - q j ),(1 - p i )q j ] when D ij > 0 or D max = min [p i q j ,(1 - p i )(1 - q j )] when D ij < 0 [10, 18, 19]. Significance test of the null hypothesis of random association between pairs of alleles at the two loci (D ij = 0) was tested by , which approximates a χ2 distribution with one degree of freedom, where n is the number of individuals sampled [20, 21]. Yates's correction was also computed.
Estimation of the magnitude and significance of pairwise LD involving biallelic loci was performed in the same way, but considering that p i and q j are the frequencies of the commonest alleles for each biallelic locus. This establishes a homogeneous criterion for the construction of 2 × 2 contingency tables, (i.e., consideration of haplotype A i B j as the one constituted by the two more frequent alleles). This criterion was uniformly followed for the estimation of the observed haplotype frequency and for computation of pairwise LD magnitude and significance in all SNP/SNP analyses. This criterion is consistent with respect to the sign-based LD estimation method recently developed for multiallelic systems [2, 13] and establishes consistency of some biological basis. By placing the most frequent allele of each of a pair of SNPs in the top left of the 2 × 2 table then if the minor alleles coincide on some haplotype, the display shows a 'main diagonal' excess (D' positive) whereas if minor allele at locus A predominate with major allele of locus B (and vice versa) the display shows a minor diagonal pattern (D' negative). When |D'| = 1, the haplotype patterns depicted (either three or two of the possible four) give information which is relevant to their possible history not fully evident from D' nor from r2 nor any other coefficient (see figure 5).
For pairwise analyses involving two multiallelic loci, haplotype A i B j was considered to comprise the two alleles of interest. This is also consistent with the sign-based LD estimation method [2, 13] in most situations, except in rare circumstances when haplotype A i B j is constituted by one allele with frequency higher than 0.5 and another allele with frequency lower than 0.5. For pairwise analyses involving one multiallelic locus and one biallelic locus, both LD estimation and representation were performed twice for each microsatellite allele of interest: A i B j was considered to comprise the microsatellite allele of interest and the commonest allele at the biallelic locus in one analysis, and the microsatellite allele of interest and the rarest allele at the biallelic locus in the other.
For users that wish to perform the analysis of multiallelic markers by dichotomising the marker into most common allele versus all other alleles combined we provide data in the output file to indicate that analysis for each combination of markers. This is provided in rows where there is a "Y" for "MostFreq1" (first marker) and "Y" for "MostFreq2" (second marker). For users who wish to collapse multiallelic markers to biallelic markers in other ways the software will accept that data in the form of input files with multiallelic markers recoded as if they were SNPs. However, it should be noted that no dichotomization represents the actual overall LD between two multiallelic loci, but only one of the possible interallelic associations.
The MIDAS output file provides D', r2, expected and estimated haplotype frequencies, allele frequencies, χ2 and distance between markers.
Results and discussion
An example dataset is shown in figure 2, comprising a set of microsatellites and SNPs from the 11p chromosome region  and Zapata et al (in prep) (subset of 50 samples). Figure 2a shows the typical unzoomed view, with pairs of markers in a grid of 1 to n columns and 1 to n rows. The plot begins at top left with 1 versus 2 and continues to n-1 versus n at bottom right. The distance between markers is represented by the intensity of background colour (closer = darker yellow). The image can be zoomed and scrolled, and positioning the mouse over a plotted result gives the statistics, a magnified plot and a plot of |D'| (as shown in the right-hand panel of figure 2b–d).
For pairwise SNP analyses the LD is represented as in figure 2d. The vertical and horizontal lines split the black square into four rectangles, the areas of which represent the expected haplotype frequencies for each allele combination (upper left is A i B j frequency = A i × B j ) (figure 3). Each quadrant then has a coloured rectangle to represent the "observed" (i.e. estimated using EM) haplotype frequency. The dimensions of the rectangle are in proportion to the two allele frequencies it represents, and its colour intensity represents the significance (by χ2) of LD or the magnitude of D' (user option on the View menu, figures in this paper show use of the significance option). Blue rectangles represent a less frequent haplotype than expected (D' < 0), and red a more frequent haplotype than expected (D' > 0). Alleles are re-ordered to ensure that the most common alleles for each marker are represented by the top-left quadrant.
For multiallelic versus biallelic or multiallelic the plot is slightly different (figure 2c). For each marker combination there are multiple pairs of vertical and horizontal lines (matching the upper left quadrant of the biallelic display). Each pair represents one allele combination, with the black rectangle indicating the expected haplotype frequency and the coloured rectangle the "observed" (i.e. estimated using EM) haplotype frequency (figure 3). The colour scheme is the same as for biallelic markers.
The markers are arrayed with marker 1 versus marker 2 at top left and marker n-1 versus marker n at bottom right, forming a right-angle triangle of plots (figure 2a). To the bottom left of the display is a line parallel to the long side of this triangle representing a map of the genomic region in which the markers are situated. Each marker is represented by a green line from their relative position on this map to the row and column in which they are plotted (figure 2a). Placing the mouse over this line (or the circle at its end) gives marker name and position.
A typical session involves preparation of an input file of genotypes (figure 4) using any mainstream spreadsheet program and exporting as tab-text format. MIDAS is then run by double-clicking the script file. The window shows basic instructions – briefly: (1) "Open genotype file" from "Open" on the "File" menu, then (2) Select "Analysis" – "LD and haplotypes". Zoom can be operated by mouse-click, key-stroke ("i" and "o") or menu, while scrolling can be operated by cursor key or scroll-bar. At minimum zoom ("View" – "fit to screen" eg figure 2a) the user can rapidly spot statistically significant results and patterns. Placing the mouse cursor over a feature displays its statistics and detail. The high levels of zoom shown in figures 2b,2c and 2d enable the user to analyse the LD in more detail, and read the statistics by using mouse-over. Export functions include a tab-text file of all results (which can be opened by any mainstream spreadsheet program) and a binary file format that stores all analysis variables and can be used to store the whole analysis for future use. The latter format speeds up viewing of previous analyses with regions containing many markers. Finally a postscript export option is available to save graphical view, although standard screen captures (using Alt-PrintScreen in Microsoft® Windows®) are adequate for most uses.
Menu options include file operations, analysis, options to customise the view and a help option, which provides simple information and instructions for usage.
Application of MIDAS
Evolutionary relationships between SNPs
The SNP/SNP plots (which can be seen at low levels of zoom) provide a quick way of inferring evolutionary relationships between markers. Figure 5 shows how three different type of plot provide this type of information.
LD between a multiallelic microsatellite and several SNPs
We have previously described the association between allele groups of a highly polymorphic microsatellite in the Growth Hormone/chorionic somatomammotrophin (GH/CSH) gene region on chromosome 17 (CSH1.01) and phenotypes of the metabolic syndrome . For these analyses we dichotomised the microsatellite on the basis of size and distribution . Figure 6 shows the interallelic LD between SNPs in the GH/CSH gene region and the CSH1.01 microsatellite. These analyses confirm the validity of our dichotomisation, indicating two major clades of alleles within the microsatellite. However, in most cases the analysis of all alleles of a multiallelic marker rather than a dichotomisation of alleles provides the maximum information, with no necessity to make biological assumptions.
LD between input haplotypes and other markers
Figure 6 also demonstrates the potential to input haplotype data and analyse it as a multi-allelic marker. In this case a 4-SNP haplotype is analysed for LD with a multi-allelic microsatellite (figure 6f). This approach can provide an overview of how biallelic or multiallelic markers interact with haplotypes, and also how haplotypes in two different regions interact with each other.
LD between two multiallelic tandem repeat loci
Our work in the IGF2-INS-TH region of chromosome 11 has identified associations between SNPs, haplotypes and other markers and obesity [23–28]. Two multi-allelic markers within this region (included in our haplotype analyses) are the insulin gene VNTR (INS VNTR) and the tyrosine hydroxylase tetranucleotide microsatellite (TH 01). Figure 7 shows the interallelic linkage disequilibrium between these two markers, with the patterns indicating which alleles are associated and also suggesting that the VNTR mutates more rapidly than the microsatellite (one TH 01 allele associates with multiple INS VNTR alleles).
Regions of "perfect" LD
Recent work using data from HapMap [29, 30] and Celera has indicated the presence of extended regions of perfect LD (where only two major haplotypes exist) . Figure 8a shows the characteristic MIDAS pattern for this type of region, with only two haplotypes existing for each pair of SNPs (data from HapMap [29, 30]). Figure 8b and 8c demonstrates an alternative approach which we have developed (SNPFrequencyViewer) to rapidly scan for these regions, in which extended regions of isofrequent SNPs correlate with regions of perfect LD (data from HapMap [29, 30]). These can then be confirmed and examined in more detail using MIDAS. Admixture of two populations highly differentiated in these genomic regions is one possible explanation, selection another.
MIDAS is a new program that presents the novel approach of analysing and graphically representing the interallelic linkage disequilibrium (LD) between multiple pairs of bi- and multiallelic markers. The graphical representation of LD incorporates information on expected haplotype frequency (under no LD), estimated haplotype frequency and D' or significance. Distance information and statistics are also presented in the interface. This enables rapid visual interpretation and inference of evolutionary and functional relationships between SNPs and microsatellites across large genomic regions. Applications to data-sets we have analysed previously demonstrate the effectiveness of viewing patterns in the data graphically rather than numerically.
Availability and requirements
Project name: MIDAS: Multiallelic Interallelic Disequilibrium Analysis Software.
Project home page: http://www.genes.org.uk/software/midas
Operating system(s): Microsoft® Windows® 2000/XP
Programming language: Python 2.4/Tkinter
Other requirements: Python 2.4 or later  must be installed before MIDAS
License: MIDAS licence supplied with program
Any restrictions to use by non-academics: royalty-free use allowed within terms of licence
Ardlie KG, Kruglyak L, Seielstad M: Patterns of linkage disequilibrium in the human genome. Nat Rev Genet 2002, 3: 299–309. 10.1038/nrg777
Zapata C, Rodríguez S, Visedo G, Sacristán F: Spectrum of nonrandom associations between microsatellite loci on human chromosome 11p15. Genetics 2001, 158: 1235–1251.
Jorde LB: Linkage disequilibrium and the search for complex disease genes. Genome Res 2000, 10: 1435–1444. 10.1101/gr.144500
Mueller JC: Linkage disequilibrium for different scales and applications. Brief Bioinform 2004, 5: 355–364. 10.1093/bib/5.4.355
Abecasis GR, Cookson WO: GOLD – graphical overview of linkage disequilibrium. Bioinformatics 2000, 16: 182–183. 10.1093/bioinformatics/16.2.182
Pettersson F, Jonsson O, Cardon LR: GOLDsurfer: three dimensional display of linkage disequilibrium. Bioinformatics 2004, 20: 3241–3243. 10.1093/bioinformatics/bth341
Barrett JC, Fry B, Maller J, Daly MJ: Haploview: analysis and visualization of LD and haplotype maps. Bioinformatics 2005, 21: 263–265. 10.1093/bioinformatics/bth457
Rodríguez S, Gaunt TR, O'Dell SD, Chen XH, Gu D, Hawe E, Miller GJ, Humphries SE, Day IN: Haplotypic analyses of the IGF2-INS-TH gene cluster in relation to cardiovascular risk traits. Hum Mol Genet 2004, 13: 715–725. 10.1093/hmg/ddh070
Rodríguez S, Gaunt TR, Dennison E, Chen XH, Syddall HE, Phillips DI, Cooper C, Day IN: Replication of IGF2-INS-TH*5 haplotype effect on obesity in older men and study of related phenotypes. Eur J Hum Genet 2006, 14: 109–116.
Weir BS, Cockerham CC: Testing hypothesis about linkage disequilibrium with multiple alleles [abstract]. Genetics 1978, 88: 633.
Karlin S, Piazza A: Statistical methods for assessing linkage disequilibrium at the HLA-A, B, C loci. Ann Hum Genet 1981, 45: 79–94.
Zapata C, Carollo C, Rodriguez S: Sampling variance and distribution of the D' measure of overall gametic disequilibrium between multiallelic loci. Ann Hum Genet 2001, 65: 395–406. 10.1046/j.1469-1809.2001.6540395.x
Zapata C, Nunez C, Velasco T: Distribution of nonrandom associations between pairs of protein loci along the third chromosome of Drosophila melanogaster. Genetics 2002, 161: 1539–1550.
The Python Programming Language[http://www.python.org]
Curie-Cohen M: Estimates of inbreeding in a natural population: a comparison of sampling properties. Genetics 1982, 100: 339–358.
Robertson A, Hill WG: Deviations from Hardy-Weinberg proportions: sampling variances and use in estimation of inbreeding coefficients. Genetics 1984, 107: 703–718.
Hill WG: Estimation of linkage disequilibrium in randomly mating populations. Heredity 1974, 33: 229–239.
Lewontin RC: The interaction of selection and linkage. I. General considerations; heterotic models. Genetics 1964, 49: 49–67.
Hedrick PW: Gametic disequilibrium measures: proceed with caution. Genetics 1987, 117: 331–341.
Weir BS: Inferences about linkage disequilibrium. Biometrics 1979, 35: 235–254.
Day IN, Chen XH, Gaunt TR, King TH, Voropanov A, Ye S, Rodriguez S, Syddall HE, Sayer AA, Dennison EM, Tabassum F, Barker DJ, Cooper C, Phillips DI: Late life metabolic syndrome, early growth, and common polymorphism in the growth hormone and placental lactogen gene cluster. J Clin Endocrinol Metab 2004, 89: 5569–5576. 10.1210/jc.2004-0152
O'Dell SD, Miller GJ, Cooper JA, Hindmarsh PC, Pringle PJ, Ford H, Humphries SE, Day IN: Apal polymorphism in insulin-like growth factor II (IGF2) gene and weight in middle-aged males. Int J Obes Relat Metab Disord 1997, 21: 822–825. 10.1038/sj.ijo.0800483
O'Dell SD, Bujac SR, Miller GJ, Day IN: Associations of IGF2 ApaI RFLP and INS VNTR class I allele size with obesity. Eur J Hum Genet 1999, 7: 821–827. 10.1038/sj.ejhg.5200381
Gu D, O'Dell SD, Chen XH, Miller GJ, Day IN: Evidence of multiple causal sites affecting weight in the IGF2-INS-TH region of human chromosome 11. Hum Genet 2002, 110: 173–181. 10.1007/s00439-001-0663-5
Gaunt TR, Cooper JA, Miller GJ, Day IN, O'Dell SD: Positive associations between single nucleotide polymorphisms in the IGF2 gene region and body mass index in adult males. Hum Mol Genet 2001, 10: 1491–1501. 10.1093/hmg/10.14.1491
Rodriguez S, Gaunt TR, Dennison E, Chen XH, Syddall HE, Phillips DI, Cooper C, Day IN: Replication of IGF2-INS-TH(*)5 haplotype effect on obesity in older men and study of related phenotypes. Eur J Hum Genet 2006, 14: 109–116.
The International HapMap Project Nature 2003, 426: 789–796. 10.1038/nature02168
Altshuler D, Brooks LD, Chakravarti A, Collins FS, Daly MJ, Donnelly P: A haplotype map of the human genome. Nature 2005, 437: 1299–1320. 10.1038/nature04226
Lawrence R, Evans DM, Morris AP, Ke X, Hunt S, Paolucci M, Ragoussis J, Deloukas P, Bentley D, Cardon LR: Genetically indistinguishable SNPs and their influence on inferring the location of disease-associated variants. Genome Res 2005, 15: 1503–1510. 10.1101/gr.4217605
Costas J, Salas A, Phillips C, Carracedo A: Human genome-wide screen of haplotype-like blocks of reduced diversity. Gene 2005, 349: 219–225. 10.1016/j.gene.2004.12.042
TRG is funded by a BHF (British Heart Foundation) Intermediate Fellowship (FS/05/065/19497), SR by a HOPE (Wessex Medical Trust) fellowship and work in our laboratory by the Medical Research Council (UK) (Programme Grant G9800748).
MIDAS was written in Python 2.4 with a Tkinter graphical interface by TRG with design suggestions and testing by all authors. Algorithms for estimating haplotype frequencies from Hill method and for LD were adapted from a BASIC code program developed by CZ (who also suggested some statistical improvements). The manuscript was drafted by TRG and SR with inputs from all authors.
Tom R Gaunt, Santiago Rodriguez contributed equally to this work.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.