Uganda Genome Resource Enables Insights into Population History and Genomic Discovery in Africa.

Genomic studies in African populations provide unique opportunities to understand disease etiology, human diversity, and population history. In the largest study of its kind, comprising genome-wide data from 6,400 individuals and whole-genome sequences from 1,978 individuals from rural Uganda, we find evidence of geographically correlated fine-scale population substructure. Historically, the ancestry of modern Ugandans was best represented by a mixture of ancient East African pastoralists. We demonstrate the value of the largest sequence panel from Africa to date as an imputation resource. Examining 34 cardiometabolic traits, we show systematic differences in trait heritability between European and African populations, probably reflecting the differential impact of genes and environment. In a multi-trait pan-African GWAS of up to 14,126 individuals, we identify novel loci associated with anthropometric, hematological, lipid, and glycemic traits. We find that several functionally important signals are driven by Africa-specific variants, highlighting the value of studying diverse populations across the region.

Rare variant contribution to human disease in 281,104 UK Biobank exomes.

Genome-wide association studies have uncovered thousands of common variants associated with human disease, but the contribution of rare variants to common disease remains relatively unexplored. The UK Biobank contains detailed phenotypic data linked to medical records for approximately 500,000 participants, offering an unprecedented opportunity to evaluate the effect of rare variation on a broad collection of traits1,2. Here we study the relationships between rare protein-coding variants and 17,361 binary and 1,419 quantitative phenotypes using exome sequencing data from 269,171 UK Biobank participants of European ancestry. Gene-based collapsing analyses revealed 1,703 statistically significant gene-phenotype associations for binary traits, with a median odds ratio of 12.4. Furthermore, 83% of these associations were undetectable via single-variant association tests, emphasizing the power of gene-based collapsing analysis in the setting of high allelic heterogeneity. Gene-phenotype associations were also significantly enriched for loss-of-function-mediated traits and approved drug targets. Finally, we performed ancestry-specific and pan-ancestry collapsing analyses using exome sequencing data from 11,933 UK Biobank participants of African, East Asian or South Asian ancestry. Our results highlight a significant contribution of rare variants to common disease. Summary statistics are publicly available through an interactive portal ( http://azphewas.com/ ).

Exome sequencing and characterization of 49,960 individuals in the UK Biobank.

The UK Biobank is a prospective study of 502,543 individuals, combining extensive phenotypic and genotypic data with streamlined access for researchers around the world1. Here we describe the release of exome-sequence data for the first 49,960 study participants, revealing approximately 4 million coding variants (of which around 98.6% have a frequency of less than 1%). The data include 198,269 autosomal predicted loss-of-function (LOF) variants, a more than 14-fold increase compared to the imputed sequence. Nearly all genes (more than 97%) had at least one carrier with a LOF variant, and most genes (more than 69%) had at least ten carriers with a LOF variant. We illustrate the power of characterizing LOF variants in this population through association analyses across 1,730 phenotypes. In addition to replicating established associations, we found novel LOF variants with large effects on disease traits, including PIEZO1 on varicose veins, COL6A1 on corneal resistance, MEPE on bone density, and IQGAP2 and GMPR on blood cell traits. We further demonstrate the value of exome sequencing by surveying the prevalence of pathogenic variants of clinical importance, and show that 2% of this population has a medically actionable variant. Furthermore, we characterize the penetrance of cancer in carriers of pathogenic BRCA1 and BRCA2 variants. Exome sequences from the first 49,960 participants highlight the promise of genome sequencing in large population-based studies and are now accessible to the scientific community.

The African Genome Variation Project shapes medical genetics in Africa.

Given the importance of Africa to studies of human origins and disease susceptibility, detailed characterization of African genetic diversity is needed. The African Genome Variation Project provides a resource with which to design, implement and interpret genomic studies in sub-Saharan Africa and worldwide. The African Genome Variation Project represents dense genotypes from 1,481 individuals and whole-genome sequences from 320 individuals across sub-Saharan Africa. Using this resource, we find novel evidence of complex, regionally distinct hunter-gatherer and Eurasian admixture across sub-Saharan Africa. We identify new loci under selection, including loci related to malaria susceptibility and hypertension. We show that modern imputation panels (sets of reference genotypes from which unobserved or missing genotypes in study sets can be inferred) can identify association signals at highly differentiated loci across populations in sub-Saharan Africa. Using whole-genome sequencing, we demonstrate further improvements in imputation accuracy, strengthening the case for large-scale sequencing efforts of diverse African haplotypes. Finally, we present an efficient genotype array design capturing common genetic variation in Africa.

Whole-Genome Sequencing Coupled to Imputation Discovers Genetic Signals for Anthropometric Traits.

Deep sequence-based imputation can enhance the discovery power of genome-wide association studies by assessing previously unexplored variation across the common- and low-frequency spectra. We applied a hybrid whole-genome sequencing (WGS) and deep imputation approach to examine the broader allelic architecture of 12 anthropometric traits associated with height, body mass, and fat distribution in up to 267,616 individuals. We report 106 genome-wide significant signals that have not been previously identified, including 9 low-frequency variants pointing to functional candidates. Of the 106 signals, 6 are in genomic regions that have not been implicated with related traits before, 28 are independent signals at previously reported regions, and 72 represent previously reported signals for a different anthropometric trait. 71% of signals reside within genes and fine mapping resolves 23 signals to one or two likely causal variants. We confirm genetic overlap between human monogenic and polygenic anthropometric traits and find signal enrichment in cis expression QTLs in relevant tissues. Our results highlight the potential of WGS strategies to enhance biologically relevant discoveries across the frequency spectrum.

The 2018 Otto Aufranc Award: How Does Genome-wide Variation Affect Osteolysis Risk After THA?

BACKGROUND: Periprosthetic osteolysis resulting in aseptic loosening is a leading cause of THA revision. Individuals vary in their susceptibility to osteolysis and heritable factors may contribute to this variation. However, the overall contribution that such variation makes to osteolysis risk is unknown. QUESTIONS/PURPOSES: We conducted two genome-wide association studies to (1) identify genetic risk loci associated with susceptibility to osteolysis; and (2) identify genetic risk loci associated with time to prosthesis revision for osteolysis. METHODS: The Norway cohort comprised 2624 patients after THA recruited from the Norwegian Arthroplasty Registry, of whom 779 had undergone revision surgery for osteolysis. The UK cohort included 890 patients previously recruited from hospitals in the north of England, 317 who either had radiographic evidence of and/or had undergone revision surgery for osteolysis. All participants had received a fully cemented or hybrid THA using a small-diameter metal or ceramic-on-conventional polyethylene bearing. Osteolysis susceptibility case-control analyses and quantitative trait analyses for time to prosthesis revision (a proxy measure of the speed of osteolysis onset) in those patients with osteolysis were undertaken in each cohort separately after genome-wide genotyping. Finally, a meta-analysis of the two independent cohort association analysis results was undertaken. RESULTS: Genome-wide association analysis identified four independent suggestive genetic signals for osteolysis case-control status in the Norwegian cohort and 11 in the UK cohort (p ≤ 5 x 10). After meta-analysis, five independent genetic signals showed a suggestive association with osteolysis case-control status at p ≤ 5 x 10 with the strongest comprising 18 correlated variants on chromosome 7 (lead signal rs850092, p = 1.13 x 10). Genome-wide quantitative trait analysis in cases only showed a total of five and nine independent genetic signals for time to revision at p ≤ 5 x 10, respectively. After meta-analysis, 11 independent genetic signals showed suggestive evidence of an association with time to revision at p ≤ 5 x 10 with the largest association block comprising 174 correlated variants in chromosome 15 (lead signal rs10507055, p = 1.40 x 10). CONCLUSIONS: We explored the heritable biology of osteolysis at the whole genome level and identify several genetic loci that associate with susceptibility to osteolysis or with premature revision surgery. However, further studies are required to determine a causal association between the identified signals and osteolysis and their functional role in the disease. CLINICAL RELEVANCE: The identification of novel genetic risk loci for osteolysis enables new investigative avenues for clinical biomarker discovery and therapeutic intervention in this disease.

A novel variant in GLIS3 is associated with osteoarthritis.

OBJECTIVES: Osteoarthritis (OA) is a complex disease, but its genetic aetiology remains poorly characterised. To identify novel susceptibility loci for OA, we carried out a genome-wide association study (GWAS) in individuals from the largest UK-based OA collections to date. METHODS: We carried out a discovery GWAS in 5414 OA individuals with knee and/or hip total joint replacement (TJR) and 9939 population-based controls. We followed-up prioritised variants in OA subjects from the interim release of the UK Biobank resource (up to 12 658 cases and 50 898 controls) and our lead finding in operated OA subjects from the full release of UK Biobank (17 894 cases and 89 470 controls). We investigated its functional implications in methylation, gene expression and proteomics data in primary chondrocytes from 12 pairs of intact and degraded cartilage samples from patients undergoing TJR. RESULTS: We detect a genome-wide significant association at rs10116772 with TJR (P=3.7×10-8; for allele A: OR (95% CI) 0.97 (0.96 to 0.98)), an intronic variant in GLIS3, which is expressed in cartilage. Variants in strong correlation with rs10116772 have been associated with elevated plasma glucose levels and diabetes. CONCLUSIONS: We identify a novel susceptibility locus for OA that has been previously implicated in diabetes and glycaemic traits.

HLA-DRB1*11 and variants of the MHC class II locus are strong risk factors for systemic juvenile idiopathic arthritis.

Systemic juvenile idiopathic arthritis (sJIA) is an often severe, potentially life-threatening childhood inflammatory disease, the pathophysiology of which is poorly understood. To determine whether genetic variation within the MHC locus on chromosome 6 influences sJIA susceptibility, we performed an association study of 982 children with sJIA and 8,010 healthy control subjects from nine countries. Using meta-analysis of directly observed and imputed SNP genotypes and imputed classic HLA types, we identified the MHC locus as a bona fide susceptibility locus with effects on sJIA risk that transcended geographically defined strata. The strongest sJIA-associated SNP, rs151043342 [P = 2.8 × 10(-17), odds ratio (OR) 2.6 (2.1, 3.3)], was part of a cluster of 482 sJIA-associated SNPs that spanned a 400-kb region and included the class II HLA region. Conditional analysis controlling for the effect of rs151043342 found that rs12722051 independently influenced sJIA risk [P = 1.0 × 10(-5), OR 0.7 (0.6, 0.8)]. Meta-analysis of imputed classic HLA-type associations in six study populations of Western European ancestry revealed that HLA-DRB1*11 and its defining amino acid residue, glutamate 58, were strongly associated with sJIA [P = 2.7 × 10(-16), OR 2.3 (1.9, 2.8)], as was the HLA-DRB1*11-HLA-DQA1*05-HLA-DQB1*03 haplotype [6.4 × 10(-17), OR 2.3 (1.9, 2.9)]. By examining the MHC locus in the largest collection of sJIA patients assembled to date, this study solidifies the relationship between the class II HLA region and sJIA, implicating adaptive immune molecules in the pathogenesis of sJIA.

Genetic architecture distinguishes systemic juvenile idiopathic arthritis from other forms of juvenile idiopathic arthritis: clinical and therapeutic implications.

OBJECTIVES: Juvenile idiopathic arthritis (JIA) is a heterogeneous group of conditions unified by the presence of chronic childhood arthritis without an identifiable cause. Systemic JIA (sJIA) is a rare form of JIA characterised by systemic inflammation. sJIA is distinguished from other forms of JIA by unique clinical features and treatment responses that are similar to autoinflammatory diseases. However, approximately half of children with sJIA develop destructive, long-standing arthritis that appears similar to other forms of JIA. Using genomic approaches, we sought to gain novel insights into the pathophysiology of sJIA and its relationship with other forms of JIA. METHODS: We performed a genome-wide association study of 770 children with sJIA collected in nine countries by the International Childhood Arthritis Genetics Consortium. Single nucleotide polymorphisms were tested for association with sJIA. Weighted genetic risk scores were used to compare the genetic architecture of sJIA with other JIA subtypes. RESULTS: The major histocompatibility complex locus and a locus on chromosome 1 each showed association with sJIA exceeding the threshold for genome-wide significance, while 23 other novel loci were suggestive of association with sJIA. Using a combination of genetic and statistical approaches, we found no evidence of shared genetic architecture between sJIA and other common JIA subtypes. CONCLUSIONS: The lack of shared genetic risk factors between sJIA and other JIA subtypes supports the hypothesis that sJIA is a unique disease process and argues for a different classification framework. Research to improve sJIA therapy should target its unique genetics and specific pathophysiological pathways.

A genome-wide association study and biological pathway analysis of epilepsy prognosis in a prospective cohort of newly treated epilepsy.

We present the analysis of a prospective multicentre study to investigate genetic effects on the prognosis of newly treated epilepsy. Patients with a new clinical diagnosis of epilepsy requiring medication were recruited and followed up prospectively. The clinical outcome was defined as freedom from seizures for a minimum of 12 months in accordance with the consensus statement from the International League Against Epilepsy (ILAE). Genetic effects on remission of seizures after starting treatment were analysed with and without adjustment for significant clinical prognostic factors, and the results from each cohort were combined using a fixed-effects meta-analysis. After quality control (QC), we analysed 889 newly treated epilepsy patients using 472 450 genotyped and 6.9 × 10(6) imputed single-nucleotide polymorphisms. Suggestive evidence for association (defined as Pmeta < 5.0 × 10(-7)) with remission of seizures after starting treatment was observed at three loci: 6p12.2 (rs492146, Pmeta = 2.1 × 10(-7), OR[G] = 0.57), 9p23 (rs72700966, Pmeta = 3.1 × 10(-7), OR[C] = 2.70) and 15q13.2 (rs143536437, Pmeta = 3.2 × 10(-7), OR[C] = 1.92). Genes of biological interest at these loci include PTPRD and ARHGAP11B (encoding functions implicated in neuronal development) and GSTA4 (a phase II biotransformation enzyme). Pathway analysis using two independent methods implicated a number of pathways in the prognosis of epilepsy, including KEGG categories 'calcium signaling pathway' and 'phosphatidylinositol signaling pathway'. Through a series of power curves, we conclude that it is unlikely any single common variant explains >4.4% of the variation in the outcome of newly treated epilepsy.

Estimating genome-wide significance for whole-genome sequencing studies.

Although a standard genome-wide significance level has been accepted for the testing of association between common genetic variants and disease, the era of whole-genome sequencing (WGS) requires a new threshold. The allele frequency spectrum of sequence-identified variants is very different from common variants, and the identified rare genetic variation is usually jointly analyzed in a series of genomic windows or regions. In nearby or overlapping windows, these test statistics will be correlated, and the degree of correlation is likely to depend on the choice of window size, overlap, and the test statistic. Furthermore, multiple analyses may be performed using different windows or test statistics. Here we propose an empirical approach for estimating genome-wide significance thresholds for data arising from WGS studies, and we demonstrate that the empirical threshold can be efficiently estimated by extrapolating from calculations performed on a small genomic region. Because analysis of WGS may need to be repeated with different choices of test statistics or windows, this prediction approach makes it computationally feasible to estimate genome-wide significance thresholds for different analysis choices. Based on UK10K whole-genome sequence data, we derive genome-wide significance thresholds ranging between 2.5 × 10(-8) and 8 × 10(-8) for our analytic choices in window-based testing, and thresholds of 0.6 × 10(-8) -1.5 × 10(-8) for a combined analytic strategy of testing common variants using single-SNP tests together with rare variants analyzed with our sliding-window test strategy.

In search of low-frequency and rare variants affecting complex traits.

The allelic architecture of complex traits is likely to be underpinned by a combination of multiple common frequency and rare variants. Targeted genotyping arrays and next-generation sequencing technologies at the whole-genome sequencing (WGS) and whole-exome scales (WES) are increasingly employed to access sequence variation across the full minor allele frequency (MAF) spectrum. Different study design strategies that make use of diverse technologies, imputation and sample selection approaches are an active target of development and evaluation efforts. Initial insights into the contribution of rare variants in common diseases and medically relevant quantitative traits point to low-frequency and rare alleles acting either independently or in aggregate and in several cases alongside common variants. Studies conducted in population isolates have been successful in detecting rare variant associations with complex phenotypes. Statistical methodologies that enable the joint analysis of rare variants across regions of the genome continue to evolve with current efforts focusing on incorporating information such as functional annotation, and on the meta-analysis of these burden tests. In addition, population stratification, defining genome-wide statistical significance thresholds and the design of appropriate replication experiments constitute important considerations for the powerful analysis and interpretation of rare variant association studies. Progress in addressing these emerging challenges and the accrual of sufficiently large data sets are poised to help the field of complex trait genetics enter a promising era of discovery.

Candidate genes for obesity-susceptibility show enriched association within a large genome-wide association study for BMI.

Before the advent of genome-wide association studies (GWASs), hundreds of candidate genes for obesity-susceptibility had been identified through a variety of approaches. We examined whether those obesity candidate genes are enriched for associations with body mass index (BMI) compared with non-candidate genes by using data from a large-scale GWAS. A thorough literature search identified 547 candidate genes for obesity-susceptibility based on evidence from animal studies, Mendelian syndromes, linkage studies, genetic association studies and expression studies. Genomic regions were defined to include the genes ±10 kb of flanking sequence around candidate and non-candidate genes. We used summary statistics publicly available from the discovery stage of the genome-wide meta-analysis for BMI performed by the genetic investigation of anthropometric traits consortium in 123 564 individuals. Hypergeometric, rank tail-strength and gene-set enrichment analysis tests were used to test for the enrichment of association in candidate compared with non-candidate genes. The hypergeometric test of enrichment was not significant at the 5% P-value quantile (P = 0.35), but was nominally significant at the 25% quantile (P = 0.015). The rank tail-strength and gene-set enrichment tests were nominally significant for the full set of genes and borderline significant for the subset without SNPs at P < 10(-7). Taken together, the observed evidence for enrichment suggests that the candidate gene approach retains some value. However, the degree of enrichment is small despite the extensive number of candidate genes and the large sample size. Studies that focus on candidate genes have only slightly increased chances of detecting associations, and are likely to miss many true effects in non-candidate genes, at least for obesity-related traits.

Rare variant association testing for next-generation sequencing data via hierarchical clustering.

OBJECTIVES: It is thought that a proportion of the genetic susceptibility to complex diseases is due to low-frequency and rare variants. Next-generation sequencing in large populations facilitates the detection of rare variant associations to disease risk. In order to achieve adequate power to detect association at low-frequency and rare variants, locus-specific statistical methods are being developed that combine information across variants within a functional unit and test for association with this enriched signal through so-called burden tests. METHODS: We propose a hierarchical clustering approach and a similarity kernel-based association test for continuous phenotypes. This method clusters individuals into groups, within which samples are assumed to be genetically similar, and subsequently tests the group effects among the different clusters. RESULTS: The power of this approach is comparable to that of collapsing methods when causal variants have the same direction of effect, but its power is significantly higher compared to burden tests when both protective and risk variants are present in the region of interest. Overall, we observe that the Sequence Kernel Association Test (SKAT) is the most powerful approach under the allelic architectures considered. CONCLUSIONS: In our overall comparison, we find the analytical framework within which SKAT operates to yield higher power and to control type I error appropriately.

Bayesian semiparametric meta-analysis for genetic association studies.

We present a Bayesian semiparametric model for the meta-analysis of candidate gene studies with a binary outcome. Such studies often report results from association tests for different, possibly study-specific and non-overlapping genetic markers in the same genetic region. Meta-analyses of the results at each marker in isolation are seldom appropriate as they ignore the correlation that may exist between markers due to linkage disequilibrium (LD) and cannot assess the relative importance of variants at each marker. Also such marker-wise meta-analyses are restricted to only those studies that have typed the marker in question, with a potential loss of power. A better strategy is one which incorporates information about the LD between markers so that any combined estimate of the effect of each variant is corrected for the effect of other variants, as in multiple regression. Here we develop a Bayesian semiparametric model which models the observed genotype group frequencies conditional to the case/control status and uses pairwise LD measurements between markers as prior information to make posterior inference on adjusted effects. The approach allows borrowing of strength across studies and across markers. The analysis is based on a mixture of Dirichlet processes model as the underlying semiparametric model. Full posterior inference is performed through Markov chain Monte Carlo algorithms. The approach is demonstrated on simulated and real data.

Application of the optimal discovery procedure to genetic case-control studies: comparison with p values and asymptotic Bayes factors.

OBJECTIVES: The power of genetic association studies is limited by stringent levels of statistical significance. To improve power, Bayes factors (BFs) have been suggested as an alternative measure to the p value, and Storey recently introduced an optimal discovery procedure (ODP) for multiple testing. We aimed to adapt the ODP to genetic case-control studies and to compare its power to p values and asymptotic BFs (ABFs). METHODS: We propose estimators of the ODP based on prospective and retrospective likelihoods. We performed simulations based on independent common SNPs and on sequence data including rare variants. Effects of causal SNPs were simulated under various distributions of effect size. RESULTS: The true ODP is never outperformed, but the estimated ODP has similar power to p values and ABFs. For common SNPs the ODP offers power advantages only in extreme scenarios. However, for rare variants the ODP and ABF detect more associations at low false-positive rates than do p values. CONCLUSIONS: The ODP can provide higher power than p values for genetic case-control studies of common variants. However, as the ABF has similar power to the ODP and is computed more rapidly, it is our currently preferred method.

Bayesian variable selection for survival regression in genetics.

Variable selection in regression with very big numbers of variables is challenging both in terms of model specification and computation. We focus on genetic studies in the field of survival, and we present a Bayesian-inspired penalized maximum likelihood approach appropriate for high-dimensional problems. In particular, we employ a simple, efficient algorithm that seeks maximum a posteriori (MAP) estimates of regression coefficients. The latter are assigned a Laplace prior with a sharp mode at zero, and non-zero posterior mode estimates correspond to significant single nucleotide polymorphisms (SNPs). Using the Laplace prior reflects a prior belief that only a small proportion of the SNPs significantly influence the response. The method is fast and can handle datasets arising from imputation or resequencing. We demonstrate the localization performance, power and false-positive rates of our method in large simulation studies of dense-SNP datasets and sequence data, and we compare the performance of our method to the univariate Cox regression and to a recently proposed stochastic search approach. In general, we find that our approach improves localization and power slightly, while the biggest advantage is in false-positive counts and computing times. We also apply our method to a real prospective study, and we observe potential association between candidate ABC transporter genes and epilepsy treatment outcomes.

Genetic association mapping via evolution-based clustering of haplotypes.

Multilocus analysis of single nucleotide polymorphism haplotypes is a promising approach to dissecting the genetic basis of complex diseases. We propose a coalescent-based model for association mapping that potentially increases the power to detect disease-susceptibility variants in genetic association studies. The approach uses Bayesian partition modelling to cluster haplotypes with similar disease risks by exploiting evolutionary information. We focus on candidate gene regions with densely spaced markers and model chromosomal segments in high linkage disequilibrium therein assuming a perfect phylogeny. To make this assumption more realistic, we split the chromosomal region of interest into sub-regions or windows of high linkage disequilibrium. The haplotype space is then partitioned into disjoint clusters, within which the phenotype-haplotype association is assumed to be the same. For example, in case-control studies, we expect chromosomal segments bearing the causal variant on a common ancestral background to be more frequent among cases than controls, giving rise to two separate haplotype clusters. The novelty of our approach arises from the fact that the distance used for clustering haplotypes has an evolutionary interpretation, as haplotypes are clustered according to the time to their most recent common ancestor. Our approach is fully Bayesian and we develop a Markov Chain Monte Carlo algorithm to sample efficiently over the space of possible partitions. We compare the proposed approach to both single-marker analyses and recently proposed multi-marker methods and show that the Bayesian partition modelling performs similarly in localizing the causal allele while yielding lower false-positive rates. Also, the method is computationally quicker than other multi-marker approaches. We present an application to real genotype data from the CYP2D6 gene region, which has a confirmed role in drug metabolism, where we succeed in mapping the location of the susceptibility variant within a small error.

Bayesian survival analysis in genetic association studies.

MOTIVATION: Large-scale genetic association studies are carried out with the hope of discovering single nucleotide polymorphisms involved in the etiology of complex diseases. There are several existing methods in the literature for performing this kind of analysis for case-control studies, but less work has been done for prospective cohort studies. We present a Bayesian method for linking markers to censored survival outcome by clustering haplotypes using gene trees. Coalescent-based approaches are promising for LD mapping, as the coalescent offers a good approximation to the evolutionary history of mutations. RESULTS: We compare the performance of the proposed method in simulation studies to the univariate Cox regression and to dimension reduction methods, and we observe that it performs similarly in localizing the causal site, while offering a clear advantage in terms of false positive associations. Moreover, it offers computational advantages. Applying our method to a real prospective study, we observe potential association between candidate ABC transporter genes and epilepsy treatment outcomes. AVAILABILITY: R codes are available upon request. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.