Bioinformatics Zen

Modelling genome assembler efficacy

The aim of generating the ‘nucleotid.es’ data was to estimate the efficacy of each assembler. As each assembler was benchmarked across multiple different genomes, any estimate for assembler efficacy needed to account for some genomes being more difficult to assemble than others. Any estimate should be an aggregate across all assessed genomes while holding any genome-specific differences constant.

To this end, linear models were used to estimate two different coefficients, one for each assembler, and one for each genome. As an example the NA50 value of any given genome assembly was split into two components: the first component measured the overall difficulty of assembling the genome, such as might be limited by genomic repeats. The second component measured the assembler efficacy, where given the same input data, each assembler produces more or less contiguous assemblies.

Equation 1 shows an example where NA50 is the sum of the estimated coefficients for two components plus any error unaccounted for in the model, using SPAdes and E. coli. In this equation, the model estimated value for \(\beta_1\) represents the efficacy of an assembler holding the differences in genome \(\beta_2\) constant. This value, along with the associated confidence intervals, was used as a measure for assembler efficacy. The natural logarithm of the genome length (\(ln. length\)) was included to account for differences in genome length. Using a log-link function for the linear model means that changes in the estimated explanatory coefficients have a log. normal change in the response variable.

We wanted to use model fitting to measure genome assembler efficacy, accounting for noise and differences in genomes, however it was not clear a priori which model to use. Selecting a model ourselves could introduce biases based on our own intuitions about what the most accurate model should be. The ideal scenario was to allow the data to inform us as to which model to use. Therefore, the first step in this analysis was to test different mixed-effect linear models¹ to determine which best describes genome assembler efficacy. The aim was to compare different hypotheses, including the null, and select the most parsimonious given the data.

The Bayesian information criterion (BIC) was used to rank to each model, where BIC is the log-likelihood of the model penalised by the number of free parameters, multiplied by number of data points. The model with the lowest BIC was the one selected as the most parsimonious - that which best fits the data without overfitting by using too many parameters.² The models tested to explain the variance in assembles are listed in Table 6.

Random slopes were not considered as all the dependent variables were categorical. Each model had a random effect intercept fitted. Each random intercept was fitted for the terms described below.

• Null - No random intercept terms were fitted.

• Genome - A different random intercept coefficient was fitted for each source genome term.

• FASTQ file - A different random intercept coefficient was fitted for each FASTQ file term.

These random intercepts allowed for each model to account for genome or input file differences in the data that were not specific to the assembler. Each model also had four variations in fixed-effect explanatory variable. These were:

• Null - assembly quality is random regardless of which assembler is used.

• Assembler - the explanatory effect for assembly quality is the assembler used, regardless of the version or command line arguments.

• Assembler Version - the explanatory effect for assembly quality is the version of the assembler used, regardless of command line arguments.

• Assembler Args - the explanatory effect for assembly quality is the command line arguments for that version of the assembler used.

The “Fixed Effect” and “Random Effect” components allowed for fitting different coefficients (\(\beta_1, \beta_2\)) averaged across the data. This was suitable for measuring single point-estimates, however we also wanted to estimate 95% confidence intervals for each assembler. Confidence intervals allowed comparing each assembler conservatively using the expected worst-case result. I.e. instead of comparing by the mean performance of each assembler, compare by the lower 5% confidence interval.

In addition to estimating mean coefficients for each genome or assembler, different variances were considered too. This modelled whether the efficacy of one assembler had a wider range of results than another. These variance components are listed in the ‘Variance Model’ column of Table 6. Each row lists which term was used to estimate a coefficient for the variance. Null indicates a single global term was estimated for the entire model, i.e. a mixed-effect model with no specific modelling of variance.

Example model fitting for contiguity

As an example, model fitting and comparison was performed by fitting each these model using NA50 as the explanatory variable. Figure 1 shows the distribution of NA50 values scaled by the length of the source genome. NA50, a metric used to assess assembly quality based on contiguity, is defined as the contig length where the cumulative sum of the ordered contigs lengths in the assembly was at least half the total assembly length, removing any contigs that don’t align to the reference sequence. Scaling NA50 by the reference length controls for the size of the genome when comparing assemblies from two different source organisms. This becomes a unitless metric between 0-1, where 1 represents an assembly without any gaps.

Figure 1 shows that scaled NA50 is distributed across four orders of magnitude from \(10^{-3} - 10^{0}\), and is variable both between organisms and within the same organism. The figure highlights that assembly quality appears to be highly dependent on the organism being sequenced, even when accounting for the differences in the length of the source organism using the \(ln. length\) term in Equation 1. This effect is modelled by the \(\beta_2\) component in Equation 1, the per-organism variability in assembly quality. In the case of contiguity, this is likely due to the repeat content of the genome and short read length [28]. The effect of GC content or number of repeats in the genome were not considered, all these factors would instead be lumped into the single genome-specific term.

All the models described (Table 6) were fitted with NA50 as the response variable. The best five models according to lowest BIC are listed in Table 1. This table shows the dBIC, residual deviance, and residual degrees of freedom for each model.³ The dBIC value shows the delta between each model and the model with the lowest BIC. Therefore the most parsimonious model for explaining variance in NA50, according to BIC, was assembler command line arguments with a random effect for source genome and variance components for both assembler arguments and genome. The next most parsimonious model had a dBIC much higher than 10, indicating there is no support for it or any other models [29].

As an aside, the null model with only the random effect for source organism was more parsimonious compared with the subset of models using only software as an explanatory effect without a random effect for source genome(not shown in the table for clarity). This emphasises the importance considering the biological source when benchmarking bioinformatics software, and any performance metrics must account for inherent random effects in the benchmark data used. The model comparison suggests not accounting for this is likely to conflate the efficacy of the software with the inherent differences in the data used for benchmarking.

Examining assembly metrics

Contiguity metrics such as NA50, count only the contigs that align to the reference. In addition to NA50, misassemblies, indels, mismatches, and Ns are regularly used to examine the quality of an assembly. These other metrics are useful to describe an assembly because incorrect regions in assemblies are likely to be annotated incorrectly, which can then subsequently change downstream analysis of the genomic features in these erroneous regions.

The above modelling approach used to examine NA50 was also applied to the number of misassemblies, and the number of incorrect bases in the same assemblies determined by aligning to the reference. The number of incorrect bases was defined as the sum of indels, mismatches and Ns. The estimated assembler coefficients (\(\beta_1\) in Equation 1) from the most parsimonious model in each case are shown in Figure 2. These assembler coefficients are point estimates with 95% confidence intervals from a log-normal distribution.⁴ The assemblers are ranked by the worst bound on the confidence intervals, which is the lower bound in the case of NA50 and the upper bound for number of misassemblies, and incorrect bases. In the cases where multiple different command line arguments for the same version were tested, the assemblers were ranked by the command line arguments that gave the “least-worst” coefficient by confidence interval. The least-worst was determined by ranking all the command line parameter sets within the same tool version, and selecting the set with the highest lower bound on the estimated coefficient confidence intervals. Ranking by least-worst on the lower bound of the confidence interval estimates was because assemblers with large confidence intervals could generate both very good and very poor results. Using the lower bound ranks assemblers by the worst-case expected result.

Figure 2A shows the linear-model estimated NA50 coefficient for each assembler, the model coefficients are exponentiated to appear on 0-1 scale as NA50 scaled by reference genome length. The estimated coefficients show that the choice of genome assembler had a relatively minor effect on the contiguity of the assembly, where most assemblers can be expected to produce similar contiguity. The exceptions to this are LightAssembler and MEGAHIT/1.0.6. MEGAHIT versions 1.0.6 and 1.1.1 have different estimates with non-overlapping confidence intervals, highlighting the important of considering and reporting which version of an assembler was used in any genomic analysis. Overall though the choice of assembler appears to have only a minor affect on contiguity.

Figure 2B shows the estimated average misassembly coefficients for each genome assembler, exponentiated to return to the original scale. The estimated coefficients are more uniformly distributed across a range of values than the estimated NA50 coefficients. In this case there are assemblers such as LightAssembler, tadpole and MEGAHIT that generate fewer misassemblies than assemblers such as Ray, Minia, StriDe and the JGI isolate pipeline. The assemblers that have the largest coefficients for NA50, also tend to have the largest coefficients for misassemblies per megabase, indicating an positive relationship between contiguity and misassemblies.

Comparing the estimated coefficients for incorrect bases per megabase in Figure 2C shows that LightAssembler and Minia tend to generate larger numbers of incorrect bases, while all the other assemblers tend to produce much lower values.

In two metric plots, two of the assemblers velvet and LightAssembler had a small estimated coefficients with a large confidence interval range. This may indicate model difficulty in estimating accurate coefficients, which could be due to both of these assemblers failing to produce results for some benchmarks (Table 7). To put this another way, the worst results for both assemblers produced no assemblies at all, which could skew the data. Future analysis could develop an unbiased method of handling missing data.

For contiguity the most parsimonious model indicates that assembler arguments are relevant (Table 1), for misassemblies the most parsimonious model indicates that only the assembler is relevant regardless of version or assembler arguments (Table 3), and for incorrect bases the most parsimonious model indicates only the assembler version is relevant (Table 4).

Effect of assembler choice on annotations

The NA50, misassemblies and incorrect bases metrics all measure different components of assembly quality. Assemblies with better values for these metrics should be more representative of the true genome sequence. As mentioned earlier, an outstanding issue is the lack of a single metric to determine which genome assembly software to use. As an example of this are the results in Figure 2A and 2B indicate that optimising for contiguity is likely to increase misassembly rate, and vice versa. The ideal case is a single metric which may be used as an objective function not only for selecting assembly software, but also for making other optimisations, such as to library and sequencing protocols. The ideal function would be a weighted combination of different metrics drawn from the requirements of the users of the pipeline, and representing different orthogonal axes of assembly quality.

All of the metrics examined so far are likely to affect the annotations in the assembly. A break in contiguity, a misassembly, or an incorrect base in a genome feature could result in no feature being annotated, or possibly worse a misannotation. Therefore using annotated genomic features may be a useful proxy for a single metric, combining multiple different metrics at once. This is an issue specifically raised by Florea et. al [24] where differences in cow genome assemblies result in thousands of predicted gene differences.

Annotation changes between assemblies can be estimated by comparing the difference between the genome feature sequences in the assembly and the reference sets⁵. This distance can be measured as the L1 norm of annotations. The L1 norm is the sum over the following function: for each distinct annotation in the superset of the assembly and the reference, 0 if the annotation appears in both the reference and the assembly, and 1 otherwise. The recently open sourced JGI tool GAET was developed to estimate this.

Figure 3 compares annotation L1 norm distance per megabase for all the genomes used in benchmarking, across all assemblers. The higher the L1 norm, the larger the difference between the reference and assembly annotated feature sets. In the majority of cases, the number of annotations different between the reference and the assembly is fewer than 200 per megabase, not including any genome specific effects. The exceptions to this were MEGAHIT/1.0.6 and LightAssembler which in many cases had a much larger number of differences compared with the other assemblers.

Similarly to the previous sections, the same models listed in Table 6 were fitted to describe the L1 norm distance between each assembly and the corresponding reference. The model for assembler with a random effect for source genome is the most parsimonious (Table 5). The again emphasises the importance of including the source genome when estimating software efficacy. Figure 4 shows the model estimated coefficients for the least-worst set of arguments for each assembler version. The figure shows that MEGAHIT/1.0.6 had higher number of incorrect annotations than the majority of the other assemblers, while this appears to be much improved in MEGAHIT/1.1.1. This again emphasises the importance of considering versioning in bioinformatics software. Light assembler could not be shown on the same scale as the other assemblers due to the estimated coefficient being much greater than the other assemblers. Outside of MEGAHIT/1.0.6 and LightAssembler the difference between the best and worst performing assembler was approximately 25 incorrect annotations per megabase. Therefore the choice of genome assembler is expected to have a non-trivial effect on downstream annotated features.

Table 2 lists the top five assembler biobox Docker images according the upper 95% confidence interval on the linear model estimated efficacy of each assembler. This table illustrates how linear models may be used to generate a single metric for each assembler across a range of benchmark genomes, which can then be used to automatically rank and select the assembly pipeline providing the likely greatest number of correct annotations. This analysis suggests at this time version v0.11.0 of the JGI isolate pipeline biobox, running the ‘default’ task was the best choice of assembler for microbial isolate Illumina sequencing data. This improves over running SPAdes alone by approximately 0.5 correct features per megabase, highlighting how preprocessing raw sequence data before assembly can result in improvements over using the raw data alone. The high ranking of shovill and StriDe assemblers, both of which also contain trimming and error correction, also points to how preprocessing can increase the expected number of correct annotations.

These conclusions are based on the data used, and different data may produce different conclusions. Specifically, earlier analysis with a smaller data set showed that the most parsimonious model was only the assembler used, regardless of the version or command line arguments. This highlights how this analysis is data-dependent on the source genomes and assemblers, and greater numbers of benchmark genomes, or including additional models terms would likely lead to smaller confidence intervals in the estimated coefficients.

Acknowledgements

The work was conducted at the U.S. Department of Energy Joint Genome Institute, a DOE Office of Science User Facility, and is supported under Contract No. DE-AC02-05CH11231. Multiple instances of advice and suggestions provided on Cross Validated to make the model fitting part of this analysis possible.

Data and Model Fitting Notes

Each model was fit with the natural log of the reference genome length as an constant term offset as shown in Equation 1. This accounts for the different sizes in genomes, where in the context of a Poisson-like model the size of the genome corresponds to the exposure rate for a response effect. For example, when comparing contiguity the response was the NA50 score, and the exposure rate was the size of the genome. Maximum NA50 scores are proportional to the size of the genome examined, and an offset term removes this effect rescaling it to be in the range 0-1.

Each model was regressed through the origin because all explanatory variables are categorical. If an intercept term were included, this would fit as a constant term when all variables are removed. For example this would estimate what the NA50 would be for no genome and no genome assembler. Therefore by removing the intercept from a model estimated using only categorical variables, the estimated coefficient will be the mean of the response for each term. No interaction terms between assembler and genome were considered when fitting each model. Adding interaction terms would likely produce a better model fit, however the aim of the analysis was to determine the average assembler efficacy without adding individual genome-specific effects for each assembler. Future analysis including model terms for genome-entropy such as percent GC or genome repeat content and compare with the models examined here to see if there would be an improvement in model fit and explanatory power.

A consideration when estimating confidence intervals was whether the model residuals were over- or under-dispersed, or whether the model fitted values trend with residual size (heteroscedasticity). This was accounted for by fitting a negative binomial model when the data was only over-dispersed, or a double Poisson model when the data was both under- and over-dispersed, [30], [31]. A negative binomial fit a parameter to allow for over-dispersion, a double Poisson model allowed for fitting both under- and over-dispersion. This was in contrast to Poisson regression which assumes the variance is equal to the mean. Both the double Poisson and negative binomial models use a log-link function for the \(\beta\) and \(\sigma\) variables, assuming a log-normal distribution of the response variable. The estimated model coefficients were converted back to the original scale by taking the exponent in Figure 2 and Figure 4.

The negative binomial model was favoured over the double Poisson, where possible, because the negative binomial was faster to fit, and with fewer parameters was more parsimonious. These models were fit using the gamlss function from the package of the same name [30]. Each model was run for 75 convergence cycles or until the model converged at a tolerance of 0.001 or less. The model fitting was run in parallel using 36 cores on an Amazon c4.8xlarge instance to allow the model fitting to finish in reasonable time.

Figure 5, Figure 6, Figure 7, and Figure 8 are diagnostics plots for each of the final models from which the assembler coefficients were drawn from. These diagnostic plots were used to ensure that the model residuals were approximately normally distributed, there was no trend in residuals, and that there were no large assembler or genome specific effects in the model residuals. Table 1, Table 3, Table 4, and Table 5 contrast the top five linear models for each of the assembly metrics considered. Table 6 shows all the different linear model parameterisations considered in this analysis.

Table 8 lists all the different genome assembler biobox Docker images used in this analysis. Table 9 lists all the genomes used to benchmark the assemblers. Table 7 lists all the instances when a benchmarked image was not able to produce quantifiable output as a result of the benchmarking.

References