Highly accurate long-read HiFi sequencing data for five complex genomes

1Until recently, DNA sequencing technologies produced either short highly accurate reads (up to 300 bases at 99% accuracy), or less-accurate long reads (10–100 s of kb at 75–90% accuracy),. Highly accurate short reads are appropriate for germline and somatic variant detection, exome sequencing, liquid biopsy, non-invasive prenatal testing, and counting applications such as transcript profiling or single-cell analysis. In contrast, error-prone long reads are more appropriate for de novo genome assembly–, haplotype phasing, structural variant identification–, full-length mRNA sequencing and mRNA isoform discovery.

2To increase the utility of noisy long-read sequencing, several error correction methods have been devised to improve the accuracy of long reads by combining the data from either multiple independent long-read molecules or combining data from long- and short-read technologies,. These error-corrected reads can then be used for assembly or other downstream applications. In general, these error correction methods suffer from mis-mapping induced errors inherent to the multi-molecule approach that hinder downstream applications.

3A third sequencing data type leveraging multiple pass circular consensus sequencing of long (up to ~25 kb) individual molecules produces highly accurate long sequencing reads (HiFi reads). The HiFi sequencing protocol, data generation, and applications are described in Fig. 1. In the initial publication, 28-fold coverage of a human genome was sequenced with average read length of 13.5 kb and an average accuracy of 99.8%. The data has demonstrated superior assembly and haplotyping results for the human genome as measured by contiguity and accuracy when compared to traditional noisy long- or short-read methods. Additionally, single nucleotide variants were called at comparable precision and recall to Illumina® NovaSeq™ data. Since the initial publication, greatly improved assembly results have been observed in other human sequencing projects– using HiFi reads.

4To encourage further application, software development, and interest in the HiFi sequencing data type, we report here the release of five deep coverage data sets spanning a set of complex genomes including Mus musculus, Zea mays, Fragaria × ananassa (Weston) Duchesne ex Rozier, Rana muscosa, and a standard metagenomic collection of 20 microbes formulated at staggered concentrations (ATCC® MSA-1003™). The data released in this study covers a wide breadth of highly complex plant, animal, and microbial organisms and will provide a useful sequence resource, driving the sequencing standards toward higher quality in the future.

5Organisms sequenced in this study include M. musculus, Z. mays, F. × ananassa, and R. muscosa. The strain of each organism, source of the material, ploidy level, inbreeding status, reference genome sequence, and genome sizes are described in Tables 1 and 2. Additionally, we are releasing sequencing reads from a mock metagenomic sample (ATCC MSA-1003) consisting of 20 bacterial DNA samples at staggered concentrations ranging from 0.02% to 18% composition of the sample. The composition of the mock metagenomic sample as well as genome sizes of the individual bacterial species and their reference sequence accessions are listed in Supplementary Table 1.

6Excluding the metagenomic sample, the expected assembly sizes for the genomes sequenced in this study ranged from the 1,600 Mb for the outbred and octoploid F. × ananassa to approximately 18,000 Mb for the outbred and diploid R. muscosa (estimate based on genome sizes of two related species Rana aurora and Rana cascadae). The individual genome sizes of the metagenomic sample range from 1.67 to 6.34 Mb, totaling 67 Mb of bacterial sequence (Supplementary Table 1).

7Genomic DNA extraction methods and details of individual library preparations are described in the sample specific sections below. In general, if the starting genomic DNA sample was larger than 25 kb, the DNA was sheared to between 15 kb and 23 kb using the Megaruptor® 3 (Diagenode). HiFi sequencing libraries were prepared using SMRTbell™ Express Template Prep Kit 2.0 and followed by immediate treatment with the Enzyme Clean Up Kit (PN: 101-843-100). The libraries were further size selected electrophoretically using either the SageELF or BluePippin Systems from SAGE Science. The appropriate fractions for sequencing runs were identified on the Femto Pulse System (Agilent). After pooling the desired size fractions, the final libraries were further cleaned up and concentrated using AMPure PB beads (Pacific Biosciences PN:100-265-900). Finally, all libraries were checked for concentration using Qubit™ 1X dsDNA HS Assay Kit (Thermo Fisher PN: Q33231) and final size distribution was confirmed on the Femto Pulse. All library sizes are described in Table 3.

8C57BL/6 J genomic DNA was obtained from The Jackson Laboratory (PN: GTC4560). The DNA arrived at an appropriate size for HiFi library preparation (~20 kb) and no shearing was required. Library preparation method, kit, and conditions were as described above. In order to tighten the size distribution of the SMRTbell library, the DNA was size fractionated using the SageELF following library preparation. The SMRTbell library was prepared with loading solution/Marker75 then loaded onto a 0.75% agarose 1kb-18 kb gel cassette (PN: ELD7510). Size fractionation was performed electrophoretically with a target size of 3,500 bp set for elution well 12, which allowed for the collection of the appropriately sized library fractions (15–23 kb) in other elution wells of the SageELF device.

9Leaf tissue for the B73 maize inbred was frozen and provided by Matthew Hufford at Iowa State University, Department of Ecology, Evolution, and Organismal Biology. Genomic DNA was isolated from the frozen leaf tissue at the University of Arizona Genomics Institute using methods previously described. The high molecular weight DNA was sheared using the Megaruptor 3 targeting a size distribution between 15 and 20 kb. Library preparation method, kit and conditions were as described above. Library size selection was performed on the Sage BluePippin using the 0.75% Agarose dye-free Gel Cassette (PN: BLF7510) and the S1 Marker. To ensure suitable yields, the 3–10 kb Improved Recovery cassette definition was run for the size selection and high pass elution mode was chosen to target recovery of molecules greater than 15 kb.

10The plant material was obtained from foundation stock of the cultivar ‘Royal Royce’ maintained by the UC Davis Strawberry Breeding Program. DNA was isolated as previously described. The genomic DNA was larger than required for HiFi library production and was sheared using the Megaruptor 3 targeting a size distribution centered around 22 kb. Library preparation method, kit, and conditions were as described above. The SageELF was used for size selection, with similar conditions as described for M. musculus above, in order to generate a library with an appropriately sized distribution.

11R. muscosa, the Mountain Yellow-legged Frog, is an endangered species endemic to California. To prevent sacrificing an individual, DNA was prepared from a fibroblast cell line (KB 21384; ISIS # 916035) originally derived from a 25-day old tadpole of undetermined sex. The cells were grown at room temperature in low O2 from explants in alpha MEM with 1% NEAA. Approximately two million cells were harvested at passage 7 and frozen in a 1X solution of PBS buffer with 10% DMSO and 10% glycerol. Genomic DNA was isolated from these cells using Qiagen’s MagAttract HMW DNA Kit (PN: 67563) following the manufacture’s protocol. The resulting HMW gDNA was sheared to a target size of 22 kb on the MegaRuptor 3 prior to library preparation. Library preparation, kit and conditions were as described above. In order to tighten the size distribution, the SMRTbell library was size fractionated using SageELF System from Sage Science. The DNA was premixed with loading solution/Marker40 and loaded onto a 0.75% Agarose 10–40 kb Cassette (PN: ELD4010). Size fractionation was performed electrophoretically with a target size of 7,000 bp set for elution well 12 in order to achieve the appropriate resolution in size separation. Fractions having the desired size distribution ranges were identified on the Femto Pulse to generate a final size selected library used in the Sequel II sequencing runs. An additional DNA damage repair step was performed using the SMRTbell Damage Repair Kit (PN:100-992-200) as this was found helpful to improve library performance in sequencing runs.

12ATCC offers a mock metagenomic community (MSA 1003) of 20 bacteria species ranging in composition from 0.02% to 18% of the sample. Isolated DNA from this sample arrived with genomic DNA having a broad distribution of sizes and was sheared using the MegaRuptor 3 to a uniform size of 13.7 kb. Library preparation method, kit and condition were described above. Rather than using electrophoretic size selection, the resulting library was size selected using AMPure PB beads (35% v/v) to remove all small fragments.

13SMRTbell libraries were bound to the sequencing polymerase enzyme using the Sequel II Binding Kit 2.0 (PN:101-842-900) with the modification that the Sequencing Primer v2 (PN:101-847–900) was annealed to the template instead of the standard primer which comes with Sequel II Binding Kit 2.0. All incubations were performed per manufacturer’s recommendations. Prior to sequencing, unbound polymerase enzyme was removed using a modified AMPure PB bead method as previously described,. Shotgun genomic DNA sequence data was collected on the Pacific Biosciences Sequel II system using HiFi sequencing protocols and Sequencing kit V2 (PN: 101-820-200). Sequence data collection was standardized to 30 hours for this study to allow ample time for multiple pass sequencing around SMRTbell template molecules of 10–25 kb which yields high quality circular consensus sequencing (HiFi) results. Raw base-called data was moved from the sequencing instrument and the imported into SMRTLink to generate HiFi reads using the CCS algorithm (version 8.0.0.80529) which processed the raw data and generated the HiFi fastq files with the following settings: minimum pass 3, minimum predicted RQ 20.

14Using Jellyfish (v.2.2.10) a k-mer analysis was performed on each of the HiFi data sets individually using a k-mer size of 21. Counting was done using a two-pass method. First, a Bloom counter was created for each HiFi read dataset using the command described in Box 1.

15After generating the Bloom counter, a frequency count of k-mers (size = 21) was run using the command shown in Box 2:

16Finally, a histogram of the k-mer frequency was generated for each dataset by using the command in Box 3.

17These outputs were then used to generate the additional summary analysis and determine genome sizes for each sample where applicable. Genome sizes were estimated from the ratio of total HiFi bases divided by the frequency mode from each k-mer distribution.

18jellyfish bc -m 21 -s <Input Size> -t <nproc> -C -o

19HiFiReadSetFilename.bc HiFiReadSet.fasta

20where Input Size = 100G (M. musculus, Z. mays, F. × ananassa and R. muscosa) and 5G (ATCC MSA-1003).

21jellyfish count -m 21 -s <Input Size> -t <nproc> -C --bc

22HiFiReadSetFilename.bc HiFiReadSet.fasta

23Where Input Size = 20G (R. muscosa), 3G (M. musculus and Z. mays), 2G (F. × ananassa) and 200M (ATCC MSA-1003).

24jellyfish histo HiFiReadSet_21mer counts.jf >

25HiFiReadSet_21mer_Histogram.out

26In the cases where references were available (M. musculus, Z. mays, and the concatenated genomes comprising the ATCC MSA-1003 sample), HiFi reads were mapped to the references using pbmm2 version 1.2.0 (https://github.com/PacificBiosciences/pbmm2) which is a customized wrapper for minimap2 using the command demonstrated in Box 4.

27To extract accuracy metrics from each bam file using Samtools version 1.9, the command shown in Box 5 was used:

28Box 6 shows the command used to extract read length metrics from each bam file using Samtools,

29Finally, coverage metrics were obtained from each bam files using the Samtools with the command listed in Box 7.

30pbmm2 align REF.fasta HiFiReadSet.fastq

31HiFiReadSet.REF.sorted.bam --preset CCS --sort -j 48 -J 16

32(where j + J= nproc=64)

33samtools view HiFiReadSet.REF.sorted.bam | awk '{ mc="";

34for(i=12;i<=NF;i++) { split($i,TAG,":"); if(TAG[1]=="mc") {

35mc=TAG[3]; break; } } if(mc != "") { print $1 "\t" mc; } }' >

36MappedConcordance.HiFiReadSet.Genome.out

37samtools view HiFiReadSet.REF.sorted.bam | head -n <input # of

38HiFi Reads> | cut -f 10 | perl -ne 'chomp;print length($_).

39"\n"' | sort | uniq -c > MappedRL.HiFiReadSet.Genome.out

40samtools depth -a HiFiReadSet.REF.sorted.bam >

41HiFiReadSet.REF.sorted.Depth.out

42All sequencing data presented are available at the Sequencing Read Archive (SRA) under the SRA study accession SRP258341. The HiFi sequencing data is stored as fastq files with one file for each Sequel II sequencing run. Information describing each data record is presented in Table 3 and described below.

43SRR11606870 The M. musculus ‘C57BL/6 J’ data record is composed of two Sequel II runs (total of two SMRT Cell 8 M) containing 4.1 M sequencing reads and 66.5 Gb of sequence which corresponds to 25-fold coverage of the mouse genome. The average read length is 16.4 kb with an average PacBio predicted quality value (RQ) of 31.

44SRR11606869 The Z. mays ‘B73’ data record is composed of two Sequel II runs (total of two SMRT Cell 8 M) containing 3.1 M sequencing reads and 48.1 Gb of sequence which corresponds to 22-fold coverage of the maize genome. The average read length is 15.6 kb with an average PacBio predicted quality value of 30.

45SRR11606867 The F. × ananassa ‘Royal Royce’ data record is composed of one Sequel II run (total of one SMRT Cell 8 M) containing 1.4 M sequencing reads and 29.7 Gb of sequence of the octoploid Royal Royce genome. The average read length is 21.7 kb with an average RQ value of 28.

46SRR11606868 The R. muscosa (cell line KB 21384; ISIS # 916035) data record is composed of 8 Sequel II runs (total of eight SMRT Cell 8 M) containing 12.1 M sequencing reads and 189.1 Gb bases of sequence which corresponds to approximately 20-fold coverage of the R. mucosa genome. The average read length is 15.7 kb with an average RQ of 31.

47SRR11606871 is the data record for the ATTC MSA-1003 mock metagenome community which is composed of 20 bacterial organisms reported to be mixed at relative amounts differing by 900-fold from highest to lowest (Supplementary Table 1). The files in this sequence record span two Sequel II runs (total of two SMRT Cell 8 M) containing 5.6 M sequencing reads with 59.1 Gb of sequence which corresponds to between ~3 and ~5,000-fold coverage of the individual bacterial genomes. The average read length is 10.5 kb with an average RQ value of 35.

48Additionally, the raw base-called subreads from which the HiFi consensus reads were derived have been made available as a resource for developers interested in improving circular consensus sequencing algorithms. These reads have been deposited to the SRA under the following data records SRR12358174 and SRR12371718 (M. musculus); SRR12358173 (Z. mays); SRR12358171 (F. x ananassa); SRR12371721, SRR12371723, SRR12371724, SRR12371725, SRR12371726, SRR12371727, SRR12371722, SRR12358172 (R. muscosa); and SRR12371719 and SRR12358170 (mock metagenome), and are further described in Supplementary Table 2.

49Two of the non-microbial organisms sequenced, M. musculus, and Z. mays, have high quality reference genomes available, allowing for detailed validation of the sequencing data. Additionally, reads from sequencing the mock metagenome sample were aligned to a concatenated file containing all microbial references listed in Supplementary Table 1 and used for validation. All reads were aligned to their corresponding references using pbmm2 resulting in over 98.9% of reads mapping to their respective genomes and 98.4% to 99.3% of the alignments being unique within the respective references (Supplementary Table 3). The mapped read lengths and read accuracies are reported in Tables 3 and 4 and distributions are presented in Fig. 2 with a breakdown of error types and their distribution described in Supplementary Figure 1. In agreement with previously published reports, the accuracy of the HiFi reads exceeds 99.5% with sequencing errors predominantly arising from indels (Supplementary Figure 1b–d). Median accuracies are 99.87%, 99.84%, and 99.99% for the mouse, maize and mock metagenome samples respectively with mean accuracies of 99.18% (mouse), 99.69% (maize) and 99.73% (mock metagenome). Sequencing read lengths (Table 3) ranged from 10.5 kb (mock metagenome) to 21.7 kb (F. × ananassa) and were dependent on the final size distributions of the sequencing libraries.

50The data for all five organisms was used to generate k-mer plots using a k-mer size of 21 (Fig. 3a–e) to estimate the sequencing coverage and complexity of for each sample. K-mer based sequencing coverage was measured at 17 to 25-fold (Table 4) for each of the individual diploid genomes sequenced (M. musculus, Z. mays, and R. muscosa) and as expected produced a multimodal distribution for the octoploid F. × ananassa and a complex curve for the metagenome sample.

51Additionally, the k-mer plots can be used to characterize genome complexity such as ploidy and/or genome duplications as evidenced by multimode distributions within the k-mer plots caused by inherent polymorphism within the respective genomes. As expected, the inbred mouse C57BL/6 J, shows a single k-mer distribution consistent with the single haplotype present in the inbred animal. The inbred B73 maize shows a dominant k-mer coverage peak at 21-fold as one would expect, but also a minor peak at 42-fold which is consistent with an ancient duplication and polyploidization of this inbred sample. The high heterozygocity and ploidy of F. × ananassa contribute to a complex and ill-defined k-mer spectrum which is consistent with previous observations. Major k-mer frequency peaks at 17, 37, and minor peaks at 74 and 109-fold presumably represent the 8n, 4n, and less well defined 2n and 1n components of the genome, respectively. These k-mer peak identifications are consistent with an 8n genome size of 1.7 Gb (29.5 Gb of sequence/17-fold coverage) which agrees with previously published genome size for octoploid strawberry.

52The diploid R. muscosa sample demonstrates a more interesting case with respect to k-mer analysis as the frequency distribution shows one single haplotype at 20-fold coverage. The presence of a single k-mer peak in the genomic reads likely speaks to population bottlenecking which reduced the level of polymorphism in the genome resulting in collapse of the paternal and maternal haplotypes into one frequency peak for a k-mer size of 21. This is further supported by an apparent haploid genome size of 9,000 Mb (as calculated by the total number of sequenced bases / frequency mode of the k-mer histogram) which is equal to one half the size of the measured diploid genome sizes (~18 Gb) of two closely related species (R. aurora and R. cascadae).

53Alternatively, the genome coverage can be measured by mapping the HiFi reads to published references. The genome wide mapping-based coverages are reported in Table 4 and distributions are show in Supplementary Figures 2–4 and agree with the k-mer based estimates for those samples with known references. Minimal impact of GC composition is observed on HiFi sequencing coverage for the mouse and maize samples (Supplementary Figure 2 and 3). The read mapping method for genome coverage also produces coverage distributions and values for each member of the mock metagenome community sample (Supplementary Figure 4) and is consistent with the genomic complexity displayed in the k-mer plot (Fig. 3e), and agrees with the uneven representation of the abundance of each microbe in the mixture (Supplementary Figure 5).

54The data presented in this manuscript should provide ample DNA sequence for genome assembly, variant detection, evaluation of metagenome completeness and metagenome assembly for the samples covered. Additionally, the data should prove useful for bioinformaticians developing, improving, and validating assembly algorithms, developing haplotyping tools, and variant detection algorithms. High contiguity and high-quality genome assemblies should also be possible for the two unpublished genomes presented in this study (F. × ananassa ‘Royal Royce’, and the endangered anuran R. muscosa). Recently, HiFi read based assemblies have reconstructed several centromeres of the human genome, and the HiFi data presented here will be useful for future updates of the reference genomes for both Z. mays ‘B73’ and M. muscosa ‘C57BL/6 J’ by adding previously unresolvable regions, possibly including some complete centromeres of these genomes. The data from the metagenome mock community should prove valuable for metagenome assembly algorithms, and other analytical tool development allowing for the assembly of complete bacterial genomes from metagenomic samples displaying high heterogeneity in individual bacterial species and relative abundance.