For the vast majority of species – including many economically or ecologically important organisms, progress in biological research is hampered due to the lack of a reference genome sequence. Despite recent advances in sequencing technologies, several factors still limit the availability of such a critical resource. At the same time, many research groups and international consortia have already produced BAC libraries and physical maps and now are in a position to proceed with the development of whole-genome sequences organized around a physical map anchored to a genetic map. We propose a BAC-by-BAC sequencing protocol that combines combinatorial pooling design and second-generation sequencing technology to efficiently approach denovo selective genome sequencing. We show that combinatorial pooling is a cost-effective and practical alternative to exhaustive DNA barcoding when preparing sequencing libraries for hundreds or thousands of DNA samples, such as in this case gene-bearing minimum-tiling-path BAC clones. The novelty of the protocol hinges on the computational ability to efficiently compare hundred millions of short reads and assign them to the correct BAC clones (deconvolution) so that the assembly can be carried out clone-by-clone. Experimental results on simulated data for the rice genome show that the deconvolution is very accurate, and the resulting BAC assemblies have high quality. Results on real data for a gene-rich subset of the barley genome confirm that the deconvolution is accurate and the BAC assemblies have good quality. While our method cannot provide the level of completeness that one would achieve with a comprehensive whole-genome sequencing project, we show that it is quite successful in reconstructing the gene sequences within BACs. In the case of plants such as barley, this level of sequence knowledge is sufficient to support critical end-point objectives such as map-based cloning and marker-assisted breeding.
The second generation of DNA sequencing instruments is revolutionizing the way molecular biologists design and carry out investigations in genomics and genetics. These new sequencing technologies (e. g. , Illumina, ABI SOLiD) can produce a significantly greater number of reads at a fraction of the cost of Sanger-based technologies, but with the exception of Roche/454 and Ion Torrent (ABI) read lengths are only 50–150 bases. While the number (and to a lesser extent the length) of reads keeps increasing at each update of these instruments, the number of samples that can be run has remained small (e. g. , two sets of eight independent lanes on the Illumina HiSeq). Since the number of reads produced by the instrument is essentially fixed, when DNA samples to be sequenced are relatively “short” (e. g. , BAC clones) and the correspondence between reads and their source has to be maintained, several samples must be “multiplexed” within a single lane to optimize the trade-off between cost and sequencing depth. Multiplexing is traditionally achieved by adding a DNA barcode to each sample in the form of an additional (oligo) adapter, but this does not scale readily to thousands of samples. Although it is possible to exhaustively barcode such a number of objects [1], the procedure of preparing (and balancing in multiplexes) thousands to ten of thousands of barcoded libraries for sequencing is very labor-intensive and can be quite expensive. Additionally, the resulting distribution of reads for each barcoded sample can be severely skewed (see, e. g. , [2], [3]), necessitating rounds of selective follow-up. Here, we demonstrate that multiplexing can be achieved without exhaustive barcoding by taking advantage of recent advances in combinatorial pooling design (also known as group testing). The essence of this work is to significantly reduce the burden of library production, without severely compromising on the sequencing coverage of each BAC. Combinatorial pooling has been used previously in the context of genome analysis, but here we attempt to use it for de novo genome sequencing. Earlier works use a simple grid design that can be very vulnerable to noise and behaves poorly when several objects are positive; a simple grid design is also far from optimal in terms of the number of pools it produces [4]–[7]. Recent works use more sophisticated pooling designs in combination with second-generation sequencing technology [8]–[12]. The application domain of “DNA Sudoku” is the detection of microRNA targets in Arabidopsis and human genes [9], whereas [8], [10] are focused on targeted resequencing (i. e. , when a reference genome is available). Pooling designs have also been used to recover novel or known rare alleles in groups of individuals [11], [12]. In our approach to de novo sequencing, subsets of non-redundant but overlapping genome-tiling BACs are chosen to form intersecting pools. Each pool is then sequenced individually on a fraction of a flowcell via standard multiplexing. Due to the short length of a BAC (typically 130 kb), cost-effective sequencing requires each sequenced sample to contain thousands of BACs. Assembling short reads originating from a mix of hundreds to thousands of BACs is likely to produce low-quality assemblies, as the assembler is unable to partition the reads according to individual BACs. Moreover, resulting contigs would not be assigned to a specific BAC. If instead reads could be assigned (or deconvoluted) to individual BACs, then the assembly could proceed clone-by-clone. We demonstrate that this objective can be achieved by choosing a pooling strategy in which each BAC is present in a carefully designed set of pools such that the identity of each BAC is encoded within the pooling pattern (rather than by its association with a particular barcode). We report experimental results on simulated data on the genome of Oryza sativa (rice) and real sequencing data on the genome of Hordeum vulgare L. (barley).
While our method can in general be applied to any set of clones that cover a genome or a portion thereof, the protocol we describe here for selective genome sequencing uses a physical map of (gene-bearing) bacterial artificial chromosomes (BACs) to identify a set of minimally redundant clones. A physical map is a partial ordering of a set of genomic clones (usually BACs) encompassing one or more chromosomes. A physical map can be represented by a set of unordered contigs, where each contig is a set of overlapping clones. A physical map is usually obtained by digesting BAC clones via restriction enzymes into DNA fragments and then measuring the length of the resulting fragments (restriction fingerprints) on an agarose gel. The smallest set of clones that spans the region represented by the physical map is called minimum tiling path (MTP). The construction of a physical library and the selection of a MTP from a physical map are well-known procedures, and many organisms now have these resources available. More details can be found in, e. g. [14], [16]–[19], and references therein. Once the set of clones to be sequenced has been identified, they must be pooled according to a scheme that allows the deconvolution of the sequenced reads back to their corresponding BACs. In Combinatorics, the design of a pooling method reduces to the problem of constructing a disjunctive matrix [20]. Each row of the disjunctive matrix corresponds to a BAC to be pooled and each column corresponds to a pool. Consider a subset of the rows (BAC clones) in the disjunctive matrix, and let be the set of pools that contain at least one BAC in. A design (or a matrix) is said to be -decodable if when, , and. The construction of -decodable pooling designs has been extensively studied [20]. The popular 2D grid design is simple to implement but cannot be used for the purposes of this work because it is only one-decodable. Recently, a new family of “smart” pooling methods has generated considerable attention [8]–[10], [15], [21], [22]. Among these, we selected the shifted transversal design [15] due to its ability to handle multiple positives and its robustness to noise. The parameters of a shifted transversal design pooling are defined by three integers, where is a prime number, defines the number of layers, and is a small integer. A layer is one of the classes in the partition of BACs and consists of exactly pools: the larger the number of layers, the higher is the decodability. By construction the total number of pools is. If we set to be the smallest integer such that where is the number of BACs that need to be pooled, then the decodability of the design is. In [15], the shifted transversal design is defined by parameters where is the number of samples to be pooled, is the prime corresponding to the number of pools in each layer, and is the number of layers. In this paper, we use instead of, and instead of. An important property of this pooling design is that any two BACs share at most only pools. By choosing a small value for one can make pooling extremely robust to deconvolution errors. In our experiments, we use, so that at least ten errors are needed to mistakenly assign a read to the wrong BAC. In contrast, two errors are sufficient to draw an erroneous conclusion with the 2D grid-design. Barley BAC pools were obtained as follows. Escherichia coli strain DH10B BAC cultures were grown individually in 96-well plates covered by a porous membrane for 36 hr in 2YT medium with 0. 05% glucose and 30 g/ml chloramphenicol at 37°C in a shaking incubator. Following combinatorial pooling of 50 l aliquots from each of 2197 BAC cultures, each of 91 collected pools (169 BACs, 8. 3 ml each) was centrifuged to create cell pellets. The pellets were frozen and then used for extraction of BAC DNA using Qiagen plasmid DNA isolation reagents. Each BAC pool DNA sample was then dissolved in 225 l of TE buffer at an estimated final concentration of 20 ng/l. For gene-BAC assignment using the Golden Gate assays, a total of 10 l (200 ng) of this DNA was then digested for 1 hour at 37°C by using 2 units of NotI enzyme with 100 g/ml BSA in a volume of 100 l. The NotI enzyme was then heat inactivated at 65°C for 20 min. DNA for three sets of BACs (HV4, HV5, HV6) were prepared using a procedure that yields on average about 60% BAC DNA and 40% E. coli DNA. DNA for set HV3 was prepared using a procedure that was expected to yield about 90% BAC DNA and 10% E. coli DNA. Although these DNA samples performed well for SNP locus detection in the GoldenGate assay, we were unaware of the extent of E. coli in the samples until we began BAC pool sequencing, after all BAC pool DNAs had been prepared. The HV3 set reported here replaced an earlier set containing more E. coli DNA. Attempts were made to remove E. coli DNA from the BAC DNA samples through selective digestion by using exonucleases, and to reduce highly repetitive DNA using a denaturation/renaturation and double strand nuclease method. These procedures provided little or no reduction of the proportion of E. coli DNA in the samples. A cost-benefit analysis determined that the cost of replacing all of the BAC pools by applying an alternative BAC DNA purification procedure yielding an average of 90–92% BAC DNA and 8–10% E. coli DNA would be no more advantageous than simply repeating the sequencing of samples for which more DNA sequence information was needed to support the sequence-to-BAC deconvolution. A video showing 76 seconds of the pooling process is available as Video S1. Sequencing of the barley BAC pools was carried out on an Illumina HiSeq 2000 at UC Riverside. Paired-end reads from each pool were quality-trimmed using a sliding window and a minimum Phred quality of 23. Next, Illumina PCR adapters were removed with Far (Flexible Adapter Remover, http: //sourceforge. net), and the remaining sequence discarded either if shorter than 36 bases or if containing any ‘N’. Finally, reads were cleaned of E. coli (DH10B) and vector contamination (pBeloBAC11) using BWA [23] and custom scripts. According to our simulations, the sequencing depth of each BAC after deconvolution is required to be at least 50x to obtain good BAC assemblies. The parameters of the pooling design should be chosen so that the coverage pre-deconvolution is at least 150x–200x to compensate for non-uniformity in the molar concentrations of individual BACs within each pool, and for losses due to sequencing errors. To understand how deconvolution is achieved, let us make the simplifying assumption that clones in the MTP do not overlap. i. e. , that the MTP BACs form a non-redundant tiling for the genome under study, or a fraction thereof. Let us pool the MTP BACs according to a shifted transversal design with layers and obtain a set of reads from them. Now, consider a read occurring only once in the portion of the genome covered by the BACs. If there are no sequencing errors and sequencing depth is sufficient, will appear in the sequenced output of exactly pools (see Figure 2, case 1). To determine the BAC to which read should be assigned, a search is made for the BAC signature that matches the list of positive pools for. For the most realistic scenario where at most MTP clones overlap, the pooling must be at least -decodable for the deconvolution to work. We expect each non-repetitive read to belong to at most two BACs if the MTP has been computed perfectly, or rarely three BACs when considering imperfections, so we set. When a read belongs to the overlap between two clones (again assuming no sequencing errors), it will appear in the sequenced output for pools (see Figure 2, case 2). The case for three overlapping clones (see Figure 2, case 3) is analogous. Recall that in step E each BAC is assigned to pools, thus the signature of a BAC is a set of numbers in the range, where the first number belongs to the range, the second belongs to, …, and the last one belongs to. In our pooling design two BAC signatures cannot share more than numbers (see Theorem I in [15]). One can think of BAC signatures as -dimensional vectors which are rather “far” from each other. Our deconvolution method proceeds as follows. First, recall the notion of k-mer of a string (read) as a contiguous substring of of length. If we denote by the length of string, observe that has -mers, not necessarily distinct. Let us call the set of reads obtained by sequencing pool, for all, where (and in our pooling design). For each set, we first compute the number of occurrences of each of its distinct -mers. Specifically, for each -mer (i. e. , for each appearing in a read), if or its reverse complement occurs exactly times in. These counts are stored in a hash table such that, given a -mer, we can efficiently retrieve a count vector of numbers, namely. Once the table is built, we process each read as follows. Given a read in pool, we fetch the count vectors for each of its -mers. Given a -mer, where, let be the number of positive (non-zero) entries in its count vector, i. e. , the number of pools where occurs at least once. Several scenarios are possible: At the end of this process, we consider the subset of -mers that have been assigned to one, two or three BACs and compute the union of their signatures, which becomes the signature of read. If a perfect match between the read signature and BAC signatures (either one, or the union of two or three) is found, then the read is assigned to the corresponding BAC (s). Reads for which no such match is found are declared non-deconvolutable and saved in a separate file. This algorithm is implemented in the tool HashFilter, which has been extensively tested under Linux and MacOS. Source code and manual can be downloaded from http: //www. cs. ucr. edu/stelo/hashfilter/, under the GNU General Public License. Once the reads have been assigned to individual BACs, sets of single and paired-end reads are assembled clone-by-clone using Velvet [24]. Velvet requires an expected coverage, which can be computed using the amount of sequenced bases assigned to each BAC and the estimated BAC size. For barley, BAC sizes were estimated from the number of bands in the restriction fingerprinting data. First, we computed the average number of bands in the 72,055 BACs fingerprinted using high-information-content fingerprinting [16], [17] (see also http: //phymap. ucdavis. edu/barley/). Assuming that the average BAC length in this set was 106 kb, we computed the multiplier to apply to the number of bands to obtain the estimated BAC length, which turned out to be 1175 bases. We used that constant to obtain estimated sizes for all BACs in HV3, HV4, HV5 and HV6 (see Dataset S4). Note that the average size is over 125 kb, much larger than the library average size of 106 kb; this indicates that the MTP selection favors larger BACs. We also tested SOAPdenovo [25] and Abyss [26] on simulated data, but there were no obvious performance benefits compared to Velvet in terms of assembly quality (data not shown). We evaluated the assembly for several choices of the -mer (hash) size, but have reported only the assembly that maximized the N50 (N50 indicates the minimum length of all contig/scaffolds that together account for at least 50% of the genome). We recorded the number of contigs, their N50/median/max/sum statistics, and the number of reads used in the assembly. For rice assemblies, we Blast-ed the BAC contigs to the rice genome. We computed the fraction of the original (source) BAC covered by at least one contig, and the number of gaps and overlaps in the assembly. The parameters used for Blast are reported in the legend of Dataset S4. For barley BAC assemblies, we carried out a validation based on the known BAC-unigene associations from the Illumina GoldenGate assay described in the next section. The validation involved Blast-ing EST-derived unigenes (Harvest: Barley assembly #35 unigenes, http: //harvest. ucr. edu) against the BAC assemblies. To reduce spurious hits, we applied three filters. First, we masked highly repetitive regions by computing the frequency of all distinct 26-mers in the cleaned/trimmed HV5 data, then masking any -mers that occurred at least 11,000 times in the reads used for the assembly (80 copies in the genome) from the assembled contigs, by replacing the occurrences of those repetitive -mers with Xs. Second, we ignored any BAC contig that covered a unigene for less than 50% of its length. Third, we excluded from the hit count any unigene that hit more than ten individual BACs overall. We recorded the number of unigenes hitting a BAC, and compared them with the expected unigenes according to the Illumina assay. Samples for the GoldenGate assay were prepared by combining 5 l of NotI-digested BAC pool DNA (10 ng) with 4 l of sonicated E. coli DNA pre-dialyzed into TE buffer at a concentration of 500 ng/l (2000 ng) and 16 l of TE buffer. The final volume of each sample was thus 25 l, composed of 0. 4 ng/l of digested BAC pool DNA and 80 ng/l of additional E. coli DNA. These DNA samples were provided to Joe DeYoung at the University of California, Los Angeles, California, or to Shiaoman Chao at the US Department of Agriculture genotyping facility in Fargo, North Dakota. The DNA concentrations were then readjusted to 50 ng/l and a total of 5 l of each DNA sample was used for each GoldenGate assay. Each Illumina GoldenGate oligonucleotide pool assay (OPA) allows interrogation of a DNA sample for the presence of 1536 SNP loci. In [27], five OPAs were designed from approximately 22,000 SNPs from EST and PCR amplicon sequence alignments. Details of the development of three test phase (POPA1, POPA2, and POPA3) and two production scale (BOPA1 and BOPA2) can be found in [27]. We genotyped the barley BAC pools described in Section “The gene space of barley” on BOPA1 and BOPA2. The output from the Illumina GoldenGate assay was first converted to binary data by visual inspection of the theta/R space in BeadStudio. A positive reading meant that the SNP locus (and its corresponding unigene) is present in at least one BAC within the pool (refer to Figure S1 for an example). Given the genotyping data for all unigene-pool pairs, we designed an algorithm that computes the optimal assignment of unigenes to BACs so that the number of errors is minimized. For a particular unigene under consideration, let be the signature set of corresponding positive pools. Let be an arbitrary set of BACs, where, and be the union of the pools that contain at least a BAC clone in. The number of errors associated with this particular choice of is defined to be the number of extra observations (equal to) plus the number of missing observation (equal to). Among all possible choices of, we chose such that the value of is minimized. When the number of errors associated with the final solution was too large (more than three), we declared that unigene to be non-decodable. This procedure resulted in 1849 unigenes mapped to one, two, or three BACs. As a verification step, when a unigene was mapped to more than one BAC, we checked whether all those BACs belonged to the same contig in the barley physical map [18], [19]. Using the genetic map developed in [27], [28] we were also able to assign these unigene-anchored BACs to genetic map positions (Dataset S6) and to check whether a BAC was associated with more than one genetic map position. The unigene-BAC error rate from these cross-checking methods appeared to be about 5%. The whole genome shotgun sequencing of barley was carried out at several locations: Ambry Genetics (Aliso Viejo, California) sequenced five (277 bases) paired-end lanes and four long-insert paired-end (LIPE) lanes (insert size of 2,3 and 5 kb); University of Minnesota (courtesy of G. Muehlbauer) sequenced two (2100 bases) paired-end lanes; UC Riverside sequenced seven (2100 bases) paired-end lanes. The number of usable paired-end bases after quality-based trimming was 159. 31 Gb and 4. 92 Gb of LIPE, for an overall 31x sequencing depth of the 5. 3 Gb barley genome. An -mer () analysis showed that the effective depth of coverage of the data was about 24x [29].
The physical map for Oryza sativa was assembled from 22,474 BACs fingerprinted at AGCoL, and contained 1,937 contigs and 1,290 singletons. From this map, we selected only BACs whose sequence could be uniquely mapped to the rice genome. We computed an MTP of this smaller map using our tool FMTP [14]. The resulting MTP contained 3,827 BACs with an average length of kb, and spanned 91% of the rice genome (which is Mb). The overlap between rice BACs is significant: 1555 BACs overlap another BAC by at least 50 Kb, and 421 BACs overlap another BAC by at least 100 Kb (see Figure S2). In general, our method makes no assumption on the shared sequence content for pooled BACs. We pooled in silico a subset of 2,197 BACs from the set above according to the shifted transversal design [15]. This pooling design is defined by three parameters (see Materials and Methods for a detailed description of the properties of the pooling design). First observe that if the MTP was truly a set of minimally overlapping clones, a two-decodable pooling design would be sufficient. We decided that a three-decodable pooling scheme would give additional protection against errors and imperfections in the MTP. Taking into consideration the format of the standard 96-well plate and the need for a 3-decodable design, we chose parameters, and, so that and. Each of the layers consisted of pools, for a total of 91 BAC pools, which left some space for a few control DNA samples on a 96-well plate. In this pooling design, each BAC is contained in pools and each pool contains BACs. We call the set of pools to which a BAC is assigned, the BAC signature. Any two BAC signatures can share at most pools, and any triplet of BAC signatures can share at most pools. Specifically, 57. 9% of any BAC signature pairs have no pool in common, 30. 6% share one pool, and 11. 5% share two pools. For triplets of BAC signatures, 18. 5% have no pool in common, 32% share one pool, 29. 6% share two pools, 14. 8% share three pools, 4. 5% share four pools, 0. 6% share five pools, and 0. 01% share six. The 91 resulting rice BAC pools were “sequenced” in silico by generating paired-end reads of 104 bases with an insert size of 327 bases, and 1% sequencing error distributed uniformly along the read. A total of 208 M usable bases gave an expected x sequencing depth for a BAC in a pool. As each BAC is present in seven pools, this is an expected x combined coverage. The 91 read pools were processed for deconvolution using the -mer based algorithm presented in the Materials and Methods section. We set because we wanted to detect an overlap between two reads of bases with a length of at least 75% (78 bases) and at most two mismatches (observe that if the two errors are equally spaced along the 78 overlapping bases, a perfect match of length 26 must occur). The computation was relatively quick, but required a significant amount of memory. The construction of the hash table required about 120 GB of RAM and 164 minutes running on one core of a Dell PowerEdge T710 server (dual Intel Xeon X5660 2. 8 Ghz, 12 cores, 169 GB RAM). The deconvolution phase took 33 minutes on 10 cores; sorting the reads into 2,197 files took 22 minutes on one core. Figure 3- (a) illustrates the distribution of signature sizes for all the distinct -mers in the rice dataset. Observe that the distribution has clear peaks around, around the interval and the interval. These peaks correspond to -mers originating from one, two, and three overlapping BACs, respectively (see Figure 2). We also have a rather large number of -mers appearing in 1–6 pools. Observe that if the sequencing depth was sufficient, and in the absence of technical errors with BACs for a long -mer to have fewer than occurrences, sequencing errors must have occurred. Figure 3- (b) shows the distribution of signature sizes for all the reads in the rice dataset at the outset of our deconvolution algorithm (presented in the Materials and Methods section). Observe that the vast majority of reads now have a signature size in the expected ranges, with the exception of reads that appear in more than 80 pools. This latter set of reads cannot be deconvoluted and is discarded. The set of reads with a signature of size 7,12–14 or 15–21 that could be deconvoluted was % of the total (see Table S2 and Dataset S1). Since we knew the BAC from which each read was generated, we determined that % of the deconvoluted reads were assigned to either the correct BAC or to a BAC overlapping the correct BAC (see Table S2 and Dataset S4). After deconvolution, the average sequencing depth for each BAC was x, about 50% higher than the expected x. Even if we are losing about % of the reads due to invalid signatures, deconvoluted reads are frequently assigned to multiple BACs, thereby amplifying the sequencing depth. Part of this inflation can be attributed to the overlap between BACs in the rice MTP (see Figure S2). In the final step of the protocol, we independently assembled the set of reads assigned to each BAC. We carried out this step with Velvet [24] for each of the 2,197 BACs, for a variety of choices of -mer size (hash length) and reported only the assembly that maximized the N50. This is an arbitrary choice that does not guarantee the “best” overall assembly. Sheet 1 in Dataset S4 summarizes the results. If we average assembly statistics over all the 2,197 BACs, the percentage of reads used in the assembly was 82. 3%, the average number of contigs was 41, the average N50 was 47,551 bp (31. 4% of the average BAC length), the average largest contig was 57,258 bp (37. 8% of the average BAC length), the average sum of all contig sizes was 137,050 bp (90. 7% of the average BAC length). The N50 is quite high, and so is the percentage of reads used by the assembler. While these numbers already indicate high quality assemblies, we determined whether BACs were correctly assembled by Blast-ing BAC contigs against the rice genome. Sheet 1 in Dataset S4 reports the results of this analysis. Considering these statistics over all the 2,197 BACs, the average BAC coverage was 76. 8%, the average gap size was 263 bp, the average number of gaps was 138, the average overlap size was 107 bp, and the average number of overlaps was 75. To establish a comparison “baseline” for these assembly statistics, we considered the most optimistic scenario of a “perfect deconvolution”, which entails using the provenance annotation of each read to assign it back to the correct BAC with 100% accuracy. Sheet 2 in Dataset S4 reports the same statistics for all 2,197 BAC assemblies in this best-case scenario. If we compute the average over all the 2,197 BACs, the average fraction of the reads used by Velvet was 82. 7% and the average N50 was 132,865 bp (88% of the average BAC length). The Blast statistics showed an average BAC coverage of 96. 3%, an average gap size of 52 bp, an average number of gaps of 97, an average overlap size of 29 bp, and an average number of overlaps of 54. While this latter BAC coverage is about 20% higher, the results following deconvolution compare quite favorably with what would be possible sequencing each BAC separately. Barley' s diploid genome size is estimated at 5,300 Mb and is composed of at least 80% highly repetitive DNA, predominantly LTR retrotransposons [30]. We started with a 6. 3x genome equivalent barley BAC library which contains 313,344 BACs with an average insert size of 106 kb [31]. Nearly 84,000 gene-enriched BACs were identified, mainly by the overgo probing method [32]. Gene-enriched BACs were fingerprinted using high-information-content fingerprinting [16], [17]. From the fingerprinting data a physical map was produced [18], [19] and a MTP of about 15,000 clones was derived [14]. Seven sets of 2,197 clones were chosen to be pooled according to the shifted transversal design [15], which we internally call HV3, HV4, …, HV9 (HV1 and HV2 were pilot experiments). We used the same pooling parameters discussed in the previous section (, and). Here we are reporting on four out of seven sets, namely HV3, HV4, HV5, and HV6. Each is comprised of pools with a total of 2,197 MTP gene-rich barley BAC clones. Given the estimated 129. 5 kb size of a BAC in the barley MTP (see Materials and Methods), the total complexity of each pool of 169 BACs is Mb and of the 2,197 BACs is Mb. To take advantage of the high density of sequencing of the Illumina HiSeq2000,13–16 pools were multiplexed on each lane, using custom multiplexing adapters. After each sample was sequenced and the reads demultiplexed, we obtained an average of 22. 9 M, 11. 3 M, 11. 5 M, and 10. 1 M reads per pool in HV3, HV4, HV5, and HV6, respectively, with a read length of 92 bases. Reads were quality-trimmed and cleaned of spurious sequencing adaptors, and then reads derived from E. coli contamination or the BAC vector were discarded (see Dataset S2). The percentage of E. coli contamination was rather high, averaging around 41%, 40%, 51%, 65% in HV3, HV4, HV5, and HV6, respectively. An alternative DNA purification method we used for HV3 showed the potential to lower amount to 8–15% if properly executed (see some of HV3 pools in Dataset S2, column I). The average number of usable reads after trimming and cleaning was about 13. 5 M, 6. 8 M, 5. 5 M, and 3. 6 M per pool in HV3, HV4, HV5, and HV6, respectively, with an average high quality read length of 87–89 bases. The number of reads in the set of 91 pools ranged between about 4. 2 M–27. 7 M in HV3,2. 9 M–17 M in HV4,2. 5 M–11. 3 M in HV5, and 1. 6 M–15 M in HV6. The total number of reads was about 1229 M, 620 M, 503 M, and 327 M in HV3, HV4, HV5, and HV6, respectively, for a total of about 109. 1T, 54. 9T, 44. 8T, and 28. 5T usable bases, respectively. The 91 read pools in the barley datasets were processed using the -mer based algorithm (HashFilter with = 26) described in the Materials and Methods section. The computation took slightly longer than on the analogous rice dataset (i. e. , about 363 minutes on one core of a Dell PowerEdge T710 server to build the hash table for HV5), but used less memory (i. e. , about 43 Gb of RAM for HV5). For HV5, the deconvolution phase took 169 minutes on 10 cores, and the sorting of reads into 2,197 BAC files took 37 minutes on one core. Due to the higher repeat content of the barley genome compared to rice, we were able to deconvolute a smaller fraction of the barley reads, about 68. 14% for HV3,59. 9% for HV4,71. 3% for HV5, and 58% for HV6 (see Tables S3, S4, S5, S6 and Dataset S3). Figure 3- (c) and (d) illustrate the distribution of signature sizes for all the distinct -mers and for all reads in HV5. As expected, the number of reads occurring in over 80 pools is much higher for barley than for rice. After deconvolution, the total number of usable bases was about 97. 9T for HV3 (about 90% of the bases before deconvolution), 34. 6T for HV4 (about 63%), 38. 9T for HV5 (about 87%), and 19. 3T for HV6 (about 68%), which translated to an average sequencing depth of coverage for each BAC of about 431x in HV3,134x in HV4,137x in HV5, and 72x in HV6 (see Dataset S4). On the set HV5, we tested HashFilter for different choices of the -mer. For, , and, the memory used by HashFilter was about 35 GB, 40 GB, and 48 GB, respectively (compared to 43 Gb of RAM for). On one CPU core, the time to build the hash table was about 323 minutes for, 338 minutes for, and 915 minutes for, compared to 363 minutes for. On ten CPU cores, the overlap phase took 108 minutes for, 165 minutes for, and 244 minutes for, compared to 169 minutes for. In terms of the number of reads deconvoluted, for, , and, we deconvoluted 62. 2%, 68. 2%, and 73. 6% of the reads in HV5, respectively (compared to 71. 3% when – see Dataset S3 for more details). We have collected strong evidence that the deconvolution accuracy was quite high in barley. For instance, six BACs that were assigned less than 20 reads in HV5 were noted as not growing during the pooling (for a video of the pooling see Video S1). For two HV sets, we realized that we erroneously swapped two adjacent pools after noticing that the percentage of deconvoluted reads for those pools was significant lower than the average. Later we confirmed the swap by mapping the reads in those pools to the unigenes that were known to be present in the pooled BACs. During the same investigation, we also assessed how the overall percentage of deconvolution would be affected if one pool was missing. We removed all the reads in a pool of HV3, and re-executed the deconvolution algorithm: the number of deconvoluted reads decreased by less than 0. 5%, a very small loss considering that an entire pool was removed. We also carried out an analysis of deconvoluted paired-end reads for HV5 to determine to what extent the left and the right mate agreed on their BAC (s) assignment. The deconvolution treated paired-end reads as two separate single-end reads, which were processed independently. For each paired-end read, we collected in the set of BACs assigned to the left mate, and in the set of BACs assigned to the right mate. Unless, we declared the paired-end read to be concordant when or. For barley HV5,68. 7% of the deconvoluted paired-end reads were concordant, which indicates that the deconvolution was quite accurate (see the second sheet in Dataset S3). We assembled each set of reads assigned to a BAC using Velvet [24] for a variety of choices of -mer size. From the assemblies obtained for different choices of, we decided to report in Dataset S4 the assembly that maximized the N50. If we average the assembly statistics over the 2,197 BACs, the number of deconvoluted reads used in the assemblies was 84% in HV3,86% in HV4,87. 6% in HV5, and 83. 2% in HV6 indicating that Velvet took advantage of most of the data; the average N50 was 8,190 bp in HV3 (7. 0% of the average BAC length in that set), 5,883 bp in HV4 (4. 65% of the average BAC length), 7,210 bp in HV5 (5. 6% of the average BAC length), and 6,032 in HV6 (4. 67% of the average BAC length); the average longest contig was 18,958 bp in HV3,15,674 bp in HV4,19,222 bp in HV5, and 16,018 in HV6; the average sum of all the contigs in each assembly was 104,578 bp in HV3 (89. 8% of the average BAC length in that set), 102,502 bp in HV4 (81. 5% of the average BAC length), 113,678 bp in HV5 (87. 8% of the average BAC length), and 98,087 bp in HV6 (75. 9% of the average BAC length). Barley BAC assemblies were compared against BAC-unigene lists obtained using the Illumina GoldenGate oligonucleotide pool assay (OPA) [33] developed for barley [27]. We used the Illumina OPAs on the same four, and three additional, sets of barley pools described above (637 pools in total) and determined which BAC clones were positive for two sets of 1,536 marker loci/unigenes (see Materials and Methods for details). The Illumina OPAs allowed us to map a total of 1,849 unique unigenes to BACs (estimated error rate of 5%, see Materials and Methods for details). Table S1 summarizes the results of unigene-BAC and BAC-unigene assignment broken down by chromosome and chromosome arms, whereas Dataset S6 contains all the solved BAC-unigene relationships along with their chromosomal location. Analysis of the assembly of the 2,197 barley BACs in each of the HV3–HV6 sets was carried out by assuming the results of the OPA as the “ground truth”, although the Illumina OPA assay has an estimated error rate of 5%. For instance in HV5 we extracted a total of 221 marker loci/unigenes that were mapped to a total of 202 distinct BACs. We obtained the sequence of these 221 unigenes from Harvest (http: //harvest. ucr. edu) and Blast-ed them against the HV5 BAC contigs. Out of 202 BACs that were expected to contain those genes, only 20 BAC assemblies (10%) missed the expected marker loci/unigenes (see Dataset S4). For the other 90% of the assemblies which contained the expected unigenes, the average coverage of those unigenes was about 90% of their length. Similar results were obtained from the HV3, HV4, and HV6 sets (see Dataset S4). This analysis suggests that these BAC assemblies contain the majority of the barley genes, which is the main objective of this work.
Barley raw sequencing data for the barley BAC set can be obtained from NCBI Sequence Read Archive accession numbers SRA051771 and SRA051780 (HV3), SRA051535 (HV4), SRA047913 (HV5) and SRA050074 (HV6). When all the BAC assemblies will be complete, we will make them available in Harvest: Barley (http: //harvest. ucr. edu) and GenBank (http: //www. ncbi. nlm. nih. gov/genbank/). The current set of HV3, HV4, HV5, and HV6 BAC assembly as well as the 31x shotgun genome assembly of barley can be accessed via our Blast server hosted at the address http: //www. harvest-blast. org/, by selecting “Morex Barley BACs” or “Barley Genome” from the database menu. These assemblies can also be downloaded from http: //www. harvest-web. org/utilmenu. wc. The source code of HashFilter is available from http: //www. cs. ucr. edu/stelo/hashfilter/under the GNU General Public License. HashFilter runs under Linux or MacOS. Additional data are available with the online version of this paper.
The problem of obtaining the full genomic sequence of an organism has been solved either via a global brute-force approach (called whole-genome shotgun) or by a divide-and-conquer strategy (called clone-by-clone). Both approaches have advantages and disadvantages in terms of cost, manual labor, and the ability to deal with sequencing errors and highly repetitive regions of the genome. With the advent of second-generation sequencing instruments, the whole-genome shotgun approach has been the preferred choice. The clone-by-clone strategy is, however, still very relevant for large complex genomes. In fact, several research groups and international consortia have produced clone libraries and physical maps for many economically or ecologically important organisms and now are in a position to proceed with sequencing. In this manuscript, we demonstrate the feasibility of this approach on the gene-space of a large, very repetitive plant genome. The novelty of our approach is that, in order to take advantage of the throughput of the current generation of sequencing instruments, we pool hundreds of clones using a special type of “smart” pooling design that allows one to establish with high accuracy the source clone from the sequenced reads in a pool. Extensive simulations and experimental results support our claims.
algorithms
genome sequencing
computer science
genomics
mathematics
biology
computational biology
discrete mathematics
combinatorics
Combinatorial Pooling Enables Selective Sequencing of the Barley Gene Space