Duplications: a challenge to genome assembly

Segmental duplications have been the focus of much of the biological analysis of the mammalian genomes [1], [2]. They play an important role in the evolution of many species, often providing a substrate for the development of new gene functions. Such duplications present a challenge to genome assembly, particularly when the duplications are recent and the copied sequences are near-identical. Assembly programs sometimes collapse nearby duplications into a single copy, or erroneously incorporate multiple copies of a unique sequence into the assembly. The creation of erroneous duplications can be caused by divergent regions in a diploid genome, in which the two haplotypes are sufficiently different that the assembler fails to merge them together. Identifying such mis-assemblies is critically important for downstream biological analysis.

Two assemblies of the Bos taurus genome

In April 2009, two assemblies of the Bos taurus genome were published simultaneously: Btau 4.0 by the Baylor College of Medicine [3] and UMD Bos taurus 2.0 (UMD 2.0) by the University of Maryland [4]. These assemblies have since been updated, and the current versions are Btau 4.2, available from the sequencing center's website ftp://ftp.hgsc.bcm.tmc.edu/pub/data/Btaurus/fasta/Btau20080815/, and UMD Bos taurus 3.1 (UMD 3.1), available from Genbank as accession DAAA00000000.2. We note that Btau 4.2 has only minor differences from the published Btau 4.0 assembly; the primary update was the replacement of selected contigs by finished BAC sequences. In this paper we analyze the latest available versions of both assemblies.

High Fidelity Duplications (HFDs) in the two assemblies of Bos taurus genome

One striking difference between the assemblies is the disparity in the number of large regions of sequence that are duplicated within the chromosomes with high fidelity between copies. We defined a High Fidelity Duplication (HFD) as any region >5 Kbp in length occurring in two copies in the assembly, such that the copies are >99% identical to each other and reside on the same chromosome. To find the HFDs we used the Nucmer software [5] to map each assembly to itself and looked for non-overlapping self matches longer than 5 kbp with at least 99% identity. Btau 4.2 has 3,111 HFDs, while UMD 3.1 has 69. More surprisingly, only 2 of these HFDs appear in both assemblies. The Btau 4.2 regions cover 83 Mbp of sequence, while the UMD 3.1 duplications cover 1.3 Mbp.

In this paper we present analysis that shows that almost all HFDs in the Btau 4.2 and some in UMD 3.1 are assembly artifacts and therefore should be ignored for biological analysis.

Figure 1 shows the histograms of coverage for all HFDs in which the two assemblies disagree about copy number; i.e., at least one of the assemblies is incorrect. We created the set B1U2 containing the regions with exactly one copy in Btau 4.2 and two copies in the UMD 3.1 assembly; conversely, we created the set B2U1 containing the regions with two copies in Btau 4.2 and one copy in UMD 3.1. We show the distributions of read coverage for regions in B1U2 (dashed line) and B2U1 (solid line) as percentages of all regions. (Note that B2U1 is a much larger set, with 3,111 regions versus just 69 regions in B1U2.) Based on this WGS coverage statistic, 47 of the 69 regions (68%) in B1U2 are more likely to be true segmental duplications, suggesting that the UMD3.1 assembly is correct for these regions. In contrast, only 187 out of 3,111 regions (6%) in B2U1 appear to be true duplications, indicating that Btau 4.2 has a large number of erroneously duplicated sequences.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Histogram of the percentage of HFDs that belong to (i) set B2U1, duplicated in Btau 4.2 and single copy in UMD Bos taurus 3.1 (solid line), and (ii) set B1U2, single copy in Btau 4.2 and duplicated in UMD Bos taurus 3.1 (dashed line).

The area under each curve integrates to 100%. The histograms were computed by mapping the WGS reads to both assemblies. The average WGS read coverage of the assemblies is 5.9. The solid vertical line is placed at 5.9/ln(2), the coverage at which it is equally likely that a region occurs in two copies versus one. 47 of the 69 regions (68%) in B1U2 are on the right hand side of the line and thus they are more likely to be true segmental duplications. 94% of the 3,111 HFDs in Btau 4.2 (set B2U1) are more likely to be unique in the genome and thus probably represent assembly errors in Btau 4.2.

https://doi.org/10.1371/journal.pone.0042680.g001

Independent validation of false duplications in Btau 4.0

The BGSAC authors devote part of their paper to discussing the biological implications of the segmental duplications in their Btau 4.0 assembly. However, in the online supplement, they remark that many of their duplications are likely a product of mis-assembly: “A total of 1,860 pairwise alignments (>20 kbp, >94% identity) corresponding to 92.45 Mbp of apparent duplicated sequence in Btau 4.0 could not be substantiated by WSSD.” Note that these duplicated sequences were omitted from the main analysis, but they are still present in the 4.0 and 4.2 assemblies. Our analysis suggests that the problem is even more extensive since 84% of the regions that we analyzed are shorter than 20 Kb (but longer than 5 kb, see the definition of the HFD above), and therefore they had to be included in the main analysis.

These indications of erroneous duplications in the Btau 4.0 assembly are supported by a recent independent study by [6], which examined intra-chromosomal duplication patterns in the Bos taurus genome using fluorescent in situ hybridization (FISH). They compared Btau 4.0 and UMD 2.0 by analyzing 13 segments of the genome that were duplicated in only one of the assemblies. The FISH results were consistent with the UMD 2.0 assembly at 10 of 13 sites, while only 2 of 13 were consistent with Btau 4.0.

Evaluating the read coverage of duplicated regions in Btau 4.2 and UMD3.1

To determine which HFD regions are likely to be actually duplicated in the genome, we examined their whole-genome shotgun (WGS) read coverage. We aligned the WGS read sequences to each region that was a HFD in either assembly and calculated that region's read coverage, shown in Figure 1. We used Nucmer software to align the reads to the HFD sequences. We used a single copy of each HFD sequence because duplicate copies differed by less than 1% by definition. We accepted all matches with >94% identity over 90% of the read length. Mate pair information was not used. We then computed the coverage for each copy of the HFD as the total number of bases in the reads that match the HFD sequence divided by the length of the HFD. Next, we compared the individual read coverage of each HFD to the mean WGS read coverage over the entire genome, which was approximately 5.9×. In our analysis we assumed that WGS reads, which provided two thirds of the sequence data, were distributed nearly uniformly over the genome. Under this assumption, the coverage of the HFDs should have a Poisson distribution with a mean coverage of 5.9 in unique regions. If a sequence is truly duplicated in the genome and all reads are aligned to a single copy of that sequence, then the expected coverage would be twice the normal coverage, or about 11.8.