Physlr: Construct a physical map from linked reads
Creators: name: Jackman, Shaun
name: Chu, Justin
name: Nikolić, Vladimir
name: Amirhossein Afshinfard
name: Coombe, Lauren
name: Gokce Dilek
name: Sauparna Palchowdhury
name: Mohamadi, Hamid
name: Warren, René
name: Joerg Bohlmann
name: İnanç Birol
Publisher: name: F1000 Research Limited
Year: 2019

description: <i>Background</i><br><br>Sequencing large molecules of DNA has drastically improved the contiguity of genome sequence assemblies. Long read sequencing has reduced sequence fidelity compared to short read sequencing, and it is currently more expensive. With linked read sequencing from 10x Genomics Chromium, reads originating from the same large molecule all have the same 16 nucleotide barcode identifier. Linked read sequencing combines the benefits of large DNA molecules with the sequence fidelity and cost of short read sequencing. Two key challenges present themselves when assembling linked read sequencing. First, each of the roughly one million barcodes contain approximately ten large molecules of DNA. Second, each individual molecule is sequenced to less than one fold depth of coverage. For this reason, the assembly of linked reads typically follows a bottom-up approach, starting with an assembly of all the reads, at first ignoring the linked-read barcodes. These initial contigs are then scaffolded using the long-range information of the linked reads. Scaffolding may be difficult if, for any reason, the contiguity of this initial assembly is low. A top-down approach could be modeled after the overlap-layout-consensus (OLC) assembly paradigm, first determining the order of the large DNA molecules, before determining the fine-scale sequence. Our tool, Physlr, implements this top-down approach to construct a physical map of large molecules sequenced using linked reads.<br><br><i>Methods</i><br><br>Each barcode is a set of reads, and Physlr treats each barcode as a bag of <i>k</i>-mers. To reduce memory requirements and improve computational efficiency, only those <i>k</i>-mers that are minimizers of each read are retained. Minimizers that occur in only a single barcode are discarded, as likely resulting from a sequencing error. The most frequent minimizers are also removed, as likely derived from repetitive sequence. All pairs of barcodes that share a substantial number of minimizers are then identified. A weighted, undirected barcode-overlap graph is constructed, wherein each vertex represents a barcode, each edge represents the intersection of the minimizers of two barcodes, and the weight of the edge is the number of common minimizers. If each barcode contained a single molecule of DNA, we would have the desired molecule-overlap graph, and we could move on to the layout step of OLC. Because each barcode is composed of multiple molecules, each vertex must first be separated into its component molecules.<br><br>For each vertex we identify a vertex-induced subgraph of its adjacent vertices and the edges between those vertices. This subgraph is composed of communities of vertices, where each community is born from a DNA molecule. These communities are identified using <i>k</i>-clique community detection. This barcode vertex is then removed and replaced by one molecule vertex for each of its detected communities, and each of its incident edges is assigned to one of these newly-created molecule vertices. This method reconstructs the underlying molecule-overlap graph from the barcode-overlap graph.<br><br>With the molecule-overlap graph constructed, we proceed to the layout stage of OLC assembly. A maximum spanning tree (MST) is used to compute a fast and approximate solution to the traveling salesman problem. The MST is computed for each connected component of the graph, identifying the longest path through the MST. This path naturally does not visit every vertex, but it is not necessary at this point to include every molecule in the physical map. All vertices in a path are removed from the graph, and this process is repeated to complete the construction of the physical map.<br><br>The physical map is a set of contigs, where each contig is an ordered list of barcodes. With the physical map constructed, two applications readily present themselves. First, a subset of reads taken from a set of barcodes from one physical map contig (or a portion of one) may be assembled to yield a local assembly of a particular genomic region. This process may be repeated to assemble the entire genome. Second, an existing assembly may be mapped to the contigs of the physical map to order and orient the sequences of that existing assembly.<br><br><i>Results</i><br><br>We constructed a physical map of the 137 Mbp fruit fly (<i>Drosophila melanogaster</i>) genome from linked reads using Physlr. The physical map is composed of one large contig for each chromosome, with an NG50 of 21.5 Mbp. Comparing the physical map to the reference genome indicates no interchromosomal misassemblies. An ABySS 2.0 assembly of the linked reads was scaffolded by mapping it to this physical map. The resulting assembly has an NG50 of 16.9 Mbp, approaching the 20.8 Mbp NG50 of Supernova. When the Supernova contigs (derived by splitting its scaffolds at Ns) are scaffolded using the physical map, the resulting assembly has an NG50 of 20.5 Mbp.<br><br>We constructed a physical map of the 1.34 Gbp zebrafish (<i>Danio rerio</i>) genome from linked reads using Physlr. The physical map has an NG50 of 12.7 Mbp. Comparing the physical map to the reference genome indicates one interchromosomal misassembly. The scaffolds of the Supernova assembly of the linked reads was further scaffolded by mapping it to this physical map. The NG50 of the Supernova assembly improved from 4.8 Mbp to 13.0 Mbp by scaffolding with Physlr.<br><br>We constructed a physical map of the 3.09 Gbp human (<i>Homo sapiens</i>) genome from linked reads using Physlr. The physical map has an NG50 of 40.9 Mbp, and twelve chromosome arms are each assembled into a single contig per arm. Comparing the physical map to the reference genome indicates one interchromosomal misassembly. We use this physical map to scaffold a Supernova assembly, which has a scaffold NG50 of 38.5 Mbp and contig NG50 of 150 kbp, derived by splitting its scaffolds at Ns. Mapping the Supernova contigs to the physical map produces scaffolds with an NG50 of 33.0 Mbp. Mapping the Supernova scaffolds to the physical map to produce super-scaffolds yields an assembly with an NG50 of 59.2 Mbp.<br><br><i>Conclusions</i><br><br>Physlr constructs a physical map of large DNA molecules sequenced using linked reads from 10x Genomics Chromium without first assembling those reads. We plan to scale Physlr to assemble conifer genomes, which can exceed 20 Gbp.