Jasmijn A Baaijens

Alexander Schönhuth

Overlap graph-based generation of haplotigs for diploids and polyploids

Motivation: Haplotype aware genome assembly plays an important role in genetics, medicine, and various other disciplines, yet generation of haplotype-resolved de novo assemblies remains a major challenge. Beyond distinguishing between errors and true sequential variants, one needs to assign the true variants to the different genome copies. Recent work has pointed out that the enormous quantities of traditional NGS read data have been greatly underexploited in terms of haplotig computation so far, which reflects that methodology for reference independent haplotig computation has not yet reached maturity.

Results:

We present POLYTE (POLYploid genome fitTEr) as a new approach to de novo generation of haplotigs for diploid and polyploid genomes of known ploidy. Our method follows an iterative scheme where in each iteration reads or contigs are joined, based on their interplay in terms of an underlying haplotype-aware overlap graph. Along the iterations, contigs grow while preserving their haplotype identity. Benchmarking experiments on both real and simulated data demonstrate that POLYTE establishes new standards in terms of error-free reconstruction of haplotype-specific sequence. As a consequence, POLYTE outperforms state-of-the-art approaches in various relevant aspects, where advantages become particularly distinct in polyploid settings.

Introduction

In most eukaryotic organisms genomes come in copies, where each copy stems from one of the ancestors. The number of copies determines the ploidy of the organism: while diploid relates to two copies, polyploid refers to more than two copies1 . The copy-specific sequences are referred to as haplotypes, which generally differ in terms of the genetic variants affecting them. Distinguishing the two haplotypes in diploid organisms (such as in most vertebrates) or more than two in polyploid organisms (such as many plants and some funghi) plays an important role in various disciplines. Prominent examples are genetics, where assigning variants to ancestors is key [START_REF] Tewhey | The importance of phase informationn for human genomics[END_REF], and medicine, because very often haplotype-specific combinations of variants establish clinically relevant effects, for example, when disease risks have been inherited [START_REF] Glusman | Whole-genome haplotyping approaches and genomic medicine[END_REF]. In general, determining haplotypic sequence, that is, in other words keeping track of ancestry based dependencies is instrumental in many biomedical settings.

Assembling the two (diploid) or more (polyploid) haplotypes from sequencing reads is known as haplotype aware genome assembly, and the resulting assembled pieces of sequence are haplotigs, as a shorthand for haplotype aware contigs. The advent of next-generation sequencing (NGS) has brought about a plethora of NGS read compatible assembly programs. The vast majority of these programs, however, do not yield haplotigs, but consensus genome sequence, as a summary across all haplotypes involved. Even then, sequencing errors, read length and hardware limitations already pose fundamental challenges during the assembly process. Baaijens and Schönhuth Generating haplotigs from NGS reads-which is the challenge that we tackle here-comes with additional obstacles. Beyond distinguishing between errors and true sequential variants, one needs to assign the true sequential variants to the different copies. This requires keeping track of information that allows to link the true sequential variants stemming from identical copies. However, NGS reads in general are rather short: techniques are needed that can link haplotype-specific variants across read boundaries. Despite the many recent advances, this is not (yet) a standard procedure in genome assembly: haplotype aware assembly can still be considered in its early stages of development which explains that further advances are desirable.

Motivation. The majority of sequencing machines installed worldwide perform traditional NGS, such as Illumina sequencing. A plethora of population-scale sequencing studies (e.g. [START_REF] Besenbacher | Novel variation and de novo mutation rates in population-wide de novo assembled Danish trios[END_REF][START_REF] Sudmant | An integrated map of structural variation in 2504 human genomes[END_REF] The Genome of the Netherlands Consortium, 2014; The UK10K Consortium, 2015)) have filled up databases with traditional, short NGS reads. In terms of quantities, traditional short NGS reads exceed the amount of reads stemming from more recent third-generation-sequencing (TGS) protocols by at least one order of magnitude. The increase in read length due to TGS has considerably spurred the development of methods for haplotype-aware assembly (see Related work). While the increase in length is beneficial, the increase in sequencing error rates is also a major obstacle when distinguishing between haplotypes, usually leaving applicants with ambiguities that are hard to resolve.

Recent work has pointed out that targeted examination of NGS (Illumina type) reads can have significant positive effects in haplotype aware assembly [START_REF] Berger | HapTree: A novel bayesian framework for single individual polyplotyping using ngs data[END_REF][START_REF] Patterson | WhatsHap: weighted haplotype assembly for future generation sequencing reads[END_REF]. Seemingly, the enormous quantities of traditional NGS read data have been underexploited in terms of haplotig computation so far. This establishes our major motivation.

To better understand where serious progress can be made, one needs to realize that existing methods for haplotype computation from traditional NGS (Illumina) reads fall into two classes: the first (and arguably more popular) choice of approaches are referred to as haplotype assembly programs. These approaches make use of a reference genome to call variants from aligned reads, which are subsequently phased into separate haplotypes. The advantage of haplotype assembly programs is their stability and their resource-friendly usage. Examples for diploid haplotype assembly are WhatsHap [START_REF] Patterson | WhatsHap: weighted haplotype assembly for future generation sequencing reads[END_REF], Phaser [START_REF] Castel | Rare variant phasing and haplotypic expression from rna sequencing with phaser[END_REF], HapCut2 [START_REF] Edge | HapCUT2: robust and accurate haplotype assembly for diverse sequencing technologies[END_REF], ProbHap [START_REF] Kuleshov | Probabilistic single-individual haplotyping[END_REF] and HapCol [START_REF] Pirola | Hapcol: accurate and memory-efficient haplotype assembly from long reads[END_REF]. Examples for polyploid haplotype assembly are HapCompass [START_REF] Aguiar | HapCompass: A fast cycle basis algorithm for accurate haplotype assembly of sequence data[END_REF], HapTree [START_REF] Berger | HapTree: A novel bayesian framework for single individual polyplotyping using ngs data[END_REF], SDhaP [START_REF] Das | SDhaP: haplotype assembly for diploids and polyploids via semi-definite programming[END_REF], and H-PoP [START_REF] Xie | H-PoP and H-PoPG: heuristic partitioning algorithms for single individual haplotyping of polyploids[END_REF]. The disadvantage of haplotype assembly programs is that they depend on high-quality reference sequence as a backbone. In addition, they depend on external variant call sets. These two factors can introduce non-negligible biases.

The second class of methods are de novo haplotype aware assembly approaches that can deal with traditional NGS (in particular Illumina) reads. The advantage of such approaches is that they are independent of reference genomes and external call sets, which eliminates the externally induced biases. There are only little such approaches available however; to the best of our knowledge, only ALLPATHS-LG [START_REF] Ribeiro | Finished bacterial genomes from shotgun sequence data[END_REF], Platanus [START_REF] Kajitani | Efficient de novo assembly of highly heterozygous genomes from whole-genome shotgun short read[END_REF], and dipSPAdes [START_REF] Safonova | dipSPAdes: Assembler for Highly Polymorphic Diploid Genomes[END_REF] explicitly aim at computation of haplotigs from (diploid) NGS data. However, ALLPATHS-LG and Platanus require particularly tailored libraries, which renders their general application difficult, and the dipSPAdes software is no longer maintained. In results of ours, we further noted that SPAdes [START_REF] Bankevich | SPAdes: A new genome assembly algorithm and its applications to single-cell sequencing[END_REF] can be run in diploid mode (which is not to be confused with the no longer maintained dipSPAdes), and is able to compute haplotigs (surprisingly not only in diploid, but also in conventional mode), thereby likely establishing the only tool among the (myriad of) approaches for consensus oriented genome assembly (see [START_REF] Bradnam | Assemblathon 2: evaluating de novo method of genome assembly in three vertebrate species[END_REF]; [START_REF] Salzberg | GAGE: A critical evaluation of genome assemblies and assembly algorithms[END_REF] for references) that one can use for computation of haplotigs from short NGS reads. In summary, there are no approaches that 1) specialize in the generation of (high-quality) haplotigs, but 2) do not depend on high quality reference sequence as a backbone, 3) do not depend on external variant call sets and 4) do not require particularly tailored sequencing libraries.

Contribution. The contribution of this paper is to close this gap in the landscape of approaches. We present POLYTE (POLYploid genome fitTEr), as an approach to do this for genomes of known ploidy. Our results indicate that POLYTE outperforms state-of-the-art approaches of the two classes-haplotype assembly and de novo assembly approacheswith significant advantages in a variety of relevant aspects. As an example of an application scenario, POLYTE outperforms the other approaches in reconstructing individual haplotypes of the human Major Histocompatibility Complex (MHC). This region of 6Mb on chromosome 6 is essential to the acquired immune system and shows very high genetic variability; haplotype-aware reconstruction of the MHC region therefore usually is particularly challenging during the assembly process. Note finally that the majority of approaches focuses on diploid genomes. Therefore, the lack of approaches that can compute haplotigs for organisms of ploidy larger than two is even more striking. For ploidy larger than two, POLYTE achieves performance rates that are nearly on a par with those achieved for diploid organisms. To the best of our understanding, because of the lack of competitors, one might perceive POLYTE's achievements for polyploid organisms as a novelty in its own right.

Related work.

In terms of assembly paradigms, POLYTE is an overlap graph based approach. It adopts ideas from earlier work that either focused on variant discovery [START_REF] Marschall | CLEVER: clique-enumerating variant finder[END_REF], viral quasispecies assembly [START_REF] Baaijens | De novo assembly of viral quasispecies using overlap graphs[END_REF][START_REF] Töpfer | Viral quasispecies assembly via maximal clique enumeration[END_REF] or metagenome gene assembly [START_REF] Gregor | Snowball: strain aware gene assembly of metagenomes[END_REF] and unites the virtues of [START_REF] Marschall | CLEVER: clique-enumerating variant finder[END_REF]the ability to handle low coverage-on the one hand, and [START_REF] Baaijens | De novo assembly of viral quasispecies using overlap graphs[END_REF]; [START_REF] Töpfer | Viral quasispecies assembly via maximal clique enumeration[END_REF]-dealing with real overlap graphs and contig computation-on the other hand. That is, POLYTE brings forth an iterative overlap graph based scheme for contig generation that reliably works in low coverage settings, requiring coverage of only as low as 5x per haplotype.

Note finally that our approach also draws motivation from the recent technology shifts, such as the advent of third-generation sequencing (TGS) and explicitly haplotype-aware sequencing protocols like StrandSeq [START_REF] Porubsky | Dense and accurate whole-chromosome haplotyping of individual genomes[END_REF], which have put the computation of haplotigs into the focus of current attention. [START_REF] Chin | Phased diploid genome assembly with single molecule realtime sequencing[END_REF]; [START_REF] Jain | Nanopore sequencing and assembly of a human genome with ultra-long reads[END_REF]; [START_REF] Weisenfeld | Direct determination of diploid genome sequences[END_REF] describe approaches that aim to exploit the respective advances in sequencing technology and protocol design. Although there are similarities between these approaches and POLYTE, we focus on NGS data and hence our method fully exploits paired-end read information. We consider the adaptation of POLYTE to TGS data most interesting future work: the framework of POLYTE is generic in terms of choosing reads, such that this is a matter of adapting parameters, more than anything else. We recall, however, that our motivation was to bring forward a method that exploits (the abundantly available) traditional NGS reads in the first place. This, for example, enables to reconstruct MHC region haplotypes in various population-scale studies (e.g. [START_REF] Besenbacher | Novel variation and de novo mutation rates in population-wide de novo assembled Danish trios[END_REF]; [START_REF] Sudmant | An integrated map of structural variation in 2504 human genomes[END_REF]; The Genome of the Netherlands Consortium (2014); The UK10K Consortium (2015)), which has been a major challenge so far.

Methods

We present POLYTE, an algorithm to assemble individual haplotypes of diploid and polyploid genomes from short read sequencing data; see Figure 1 for the complete workflow. POLYTE follows the overlap-layoutconsensus (OLC) paradigm, where consensus refers to removing errors within haplotypes (instead of the common interpretation of reaching consensus across different haplotypes). Our method starts by constructing a read overlap graph which is used for error correction of the input sequences. Subsequently, we make use of an iterative OLC scheme, where in each iteration a contig overlap graph is constructed. This graph is further reduced by applying transitive edge removal and read-based branch reduction. Then, contigs are clustered and merged according to their interplay within the overlap graph, resulting in a collection of extended contigs ('contig extension' in Figure 1). These extended contigs establish the nodes of the contig overlap graph of the next iteration, which is achieved by an updating procedure. When contigs can not be merged any further, POLYTE outputs the final set of contigs. When dealing with diploid organisms, an additional assembly stage can be activated which consists of two additional steps ('diploid branch reduction' and 'contig extension' in Figure 1), creating an optional output that is refined for diploid organisms.

Given that we are dealing with data of relatively low sequencing depth, we need to exploit the information present in the sequencing reads as much as possible. The initial error correction procedure is particularly crucial, as sequencing errors can heavily disturb the process of distinguishing between different haplotypes. For this error correction step, approximate suffix-prefix overlaps are computed to establish an initial read overlap graph. Inspired by [START_REF] Baaijens | De novo assembly of viral quasispecies using overlap graphs[END_REF] and [START_REF] Töpfer | Viral quasispecies assembly via maximal clique enumeration[END_REF], maximal cliques are enumerated in the non-oriented graph and errors are corrected by inspecting the read overlaps within the cliques. By design of the overlap graph-edges indicate that two reads stem from identical haplotypes-every clique only contains reads from identical haplotypes, which allows to eliminate errors based on majority votes. Note that this procedure is particularly tailored to low coverage settings with known ploidy: admissible clique sizes and minimal sequence overlap lengths can heavily vary in comparison to earlier approaches. However, with edge criteria that are much less restrictive than in other approaches, we obtain a larger number of spurious edges. We have developed a procedure for read-based branch reduction to reduce the number of spurious edges in the overlap graph, which is of great importance for accurate reconstruction of haplotigs.

In the following sections we will discuss each of the steps involved in POLYTE, following the workflow depicted in Figure 1.

Read overlap graph construction

The steps outlined in this section refer to the initial step 'approximate suffix-prefix overlaps' that leads to the establishment of the 'read overlap graph' in Figure 1.

Read overlap graph: definition. The read overlap graph follows the idea that nodes are reads and edges indicate that a pair of reads stem from identical haplotypes. Given the input consisting of paired-end sequencing reads (Illumina), let R be the collection of single end sequences from all paired-end reads. The read overlap graph G = (V, E) is a directed graph where V corresponds to the collection of input sequences R. That is, for every paired-end read we have two vertices v, v 0 2 V , one for each single end sequence R 2 R. Directed edges v i ! v j 2 E connect sequences R i , R j whenever the suffix of R i overlaps the prefix of R j for at least 50% of the average sequence length of all reads. Furthermore, for each edge

v i ! v j , we require QS(R i , R j )
, where QS : R ⇥ R ! R is a quality score and is an appropriate threshold. This threshold is determined based on empirical statistics so as to maximize the chances that the edge (v i , v j ) indeed indicates that the corresponding sequences R i and R j stem from identical haplotypes; increasing the -threshold would lead to a higher accuracy but possible loss of low abundance haplotypes. In this, we largely follow ideas presented in earlier work [START_REF] Baaijens | De novo assembly of viral quasispecies using overlap graphs[END_REF][START_REF] Marschall | CLEVER: clique-enumerating variant finder[END_REF][START_REF] Töpfer | Viral quasispecies assembly via maximal clique enumeration[END_REF].

The difference with respect to these prior approaches is that only single ends are considered, whereas in the earlier approaches nodes represent the entire paired-end reads. Also note that here overlap graphs are twice as large in comparison to the earlier approaches, because each pairedend read is represented by two nodes, instead of only one. While this difference imposes substantial methodical and technical challenges, it is key to dealing with low coverage because it decisively increases the recall in terms of recovering reads that stem from identical haplotypes.However, it also implies follow-up complications, because the information that read ends come in pairs is temporarily lost. In POLYTE, paired-end information is stored and used in later steps; see Section 2.4 below.

Construction. Computation of the edges for the read overlap graph requires enumeration of all pairwise approximate suffix-prefix overlaps (of sufficient length) between the single read ends R 2 R and evaluation of a quality score QS(R i , R j ) for each pair of sequences for which a sufficiently good overlap was established during the approximate suffixprefix overlap computation. We further orient the edges (which is necessary because reads can stem from either the forward or the reverse stand) and systematically remove double transitive edges, which ensures that one can enumerate maximal cliques in an efficient manner (see Section 2.2). Each of these graph construction steps is described in detail in Section 1 of the Supplementary Material.

The computation of approximate suffix-prefix overlaps for vertebrate genome sized input read sets is a serious issue, currently hardly conceivable without external auxiliary means (see also [START_REF] Simpson | Efficient de novo assembly of large genomes using compressed data structures[END_REF]).

Here, we suggest a method that aims to suppress externally introduced biases to a maximum degree. We make use of a reference genome for binning reads in an initial step and, after binning, we discard the reference genome and any related information entirely such that POLYTE operates Baaijens and Schönhuth in full de novo mode. This binning step does not require a high-quality reference genome, as long as reads get mapped; any unaligned reads are discarded (see also Supplementary Section 9).

Correction of sequencing errors

After the establishment of the read overlap graph, we cluster its nodes by enumerating the maximal cliques contained in the non-oriented graph. The idea is to collect groups of reads belonging to the same haplotype and produce error-free sequences for subsequent assembly steps ('corrected sequences', Figure 1). By definition of a maximal clique-a maximal group of nodes all of which are connected by edges-maximal cliques represent maximum-sized groups of reads all of which belong to the same haplotype. Once all maximal cliques are determined, it is therefore reasonable to merge the reads within a maximal clique into a single contig. Note that this contig is longer than the individual reads participating in the contig and that sequencing errors can be eliminated by raising majority votes among the reads participating in the maximal clique. While this reflects an approved procedure in its generic form [START_REF] Baaijens | De novo assembly of viral quasispecies using overlap graphs[END_REF][START_REF] Gregor | Snowball: strain aware gene assembly of metagenomes[END_REF][START_REF] Töpfer | Viral quasispecies assembly via maximal clique enumeration[END_REF], accounting for the particular setting we are facing here-namely low coverage in combination with sequence-based edge definition-requires particular care.

The minimum clique size depends on the coverage per haplotype; in all settings considered we are dealing with known ploidy, while the overall coverage of reads can be determined by usual considerations, which yields per-haplotype-coverage estimates. To determine the optimal minimum size of a clique for a given per-haplotype coverage, we compute the probability p c,k that, due to unfortunate fragmentation of sequencing reads, at a perhaplotype coverage of c there is no clique of size k that extends a given sequencing read R to the right when requiring at least 50% read overlap. In other words, we compute the probability that there are at most k 1 reads extending R to the right; the exact same analysis applies to extensions to the left.

For determining p c,k , we assume that sequencing reads are fragmented randomly, which implies that reads are generated independently of one another. Let R be a read and S be a set containing reads from the haplotype of R at exactly 1x coverage, further assuming that all reads R 0 2 S have the same length as R (which reflects that all single read ends have the same length). It is straightforward to see that the probability that there is R 0 2 S that overlaps R at at least 50% of its length (into one direction, left or right) as 0.5. When dealing with a per-haplotype coverage of c, we assume the existence of c sets of reads S i , i = 1, .., c all of which contain reads that cover the haplotype c at 1x. For computing p c,k we consider that for only k 2 of the c sets S i , i = 1, ..., c we have that there is R 0 2 S i that overlaps R at at least 50% of its length (resulting in a clique of size at most k 1), which evaluates as

p c,k = k 2 X i=0 ⇣ c 1 i ⌘ 0.5 i 0.5 c 1 i = k 2 X i=0 ⇣ c 1 i ⌘ 0.5 c 1 . ( 1
)
We aim to have p c,k low to be able to deal with sufficiently many cliques, hence for every choice of c we compute k such that p c,k < 0.001. In this regard, we obtain that for up to 10x per haplotype an appropriate choice for the minimum clique size is 2, for coverages between 10x and 15x a minimum clique size of 3 is required, while for c 15x an optimal choice for the minimum clique size is 4. Note that in practice cliques do not grow larger than size 4 because of double transitive edge removal (see Supplementary Material).

Contig overlap graph construction

Given the corrected sequences obtained by merging maximal cliques, we build a new graph: the contig overlap graph (see Figure 1). , E 0 ) is very similar to the read overlap graph, except that we construct it from a set of contigs assumed to be free of sequencing errors. Therefore, every node v 2 V 0 corresponds to a contig and we add an edge between a pair of nodes whenever they have an exact (i.e. error-free) overlap of sufficient length.

Construction. The contig overlap graph can be constructed very efficiently by making use of the FM-index-based algorithm from Section 2.1 while allowing only exact overlaps. This gives us the complete edge set E 0 without any further computations, since we do not need to compute the overlap quality score for exact overlaps. Note that the minimal overlap length in the contig overlap graph does not need to be as high as before error correction and it is independent of the read length: all experiments were performed using a minimal contig overlap of 50bp.

Remark. Allowing approximate overlaps in this stage of the algorithm, for example by allowing some substitutions, would slow down the contig overlap graph construction considerably. Although the additional edges could lead to improved recovery of true haplotypes, it would also bring the risk of collapsing highly similar sequences and thus missing haplotypes.

Branch reduction in the contig overlap graph

Before using the contig overlap graph to extend our contigs, we trim the graph by removing redundant vertices and edges and resolving branches based on read evidence where possible, now also exploiting the paired-end information. After completing this step, we have a 'reduced graph' (see Figure 1) that is ready for contig extension.

Transitive edge removal. An edge u ! w 2 E 0 is called transitive if there exists a vertex v 2 V 0 and edges u ! v, v ! w 2 E 0 . Now that sequences (contigs) are assumed to be error-free, transitive edges have become fully redundant, hence we remove all transitive edges from the graph before further processing.

Branching edges and nodes. The indegree (resp. outdegree) of a node v 2 V 0 is defined as the total number of incoming (resp. outgoing) edges in G 0 . If v has indegree greater than one, we say v has an in-branch; analogously, if v has outdegree greater than one, we say that v has an outbranch. We refer to the corresponding edges as branching edges and to v as a branching node. Since we did not use any read pairing information during construction of our overlap graphs, we observe many branches in the contig overlap graph. We now use the information how ends are paired to remove any branching edges in the contig overlap graph that do not correspond to a true haplotype.

Merging simple paths. Following the above definition, any edges that are not branching edges constitute simple paths through the contig overlap graph. For such paths, there is only one possible way to combine the corresponding contigs; hence, before processing the graph any further, we merge every simple path into a single contig. Since edges in the graph represent exact overlaps, this is a straightforward procedure.

Branching components. After merging simple paths, all remaining edges are branching edges. We define a branching component as subgraph H of the contig overlap graph, such that (1) H is an induced subgraph, (2) H is connected as an undirected graph, and (3) within H, any vertex has only incoming or outgoing edges in H, but not both. A branching component is defined to be maximal with respect to these three properties; see Figure 2. Intuitively, a branching component reflects all possible haplotypes within a small region of the genome. Note that different components may intersect across their vertex sets, but cannot have any edges in common. In other words, the maximal branching components partition the set of all branching edges, as illustrated in Figure 2. This partition can be found in time linear in the number of branching edges by alternatingly traversing in-branch edges and outbranch edges until every edge has been seen exactly once; see Section 2 in the Supplementary Material for further details. After enumerating all maximal branching components, we evaluate read evidence per component.

Read evidence. The main idea of read-based branch reduction is to remove all branching edges for which there is insufficient read evidence in the input data. For this purpose, we keep track of all original sequencing reads ('subreads') that were used to build a contig; each of these subreads may provide evidence for a branching edge. Within a branching component, we first list all variant positions, i.e., the positions at which the sequences corresponding to the different neighbors differ from each other. Intuitively, these are the positions where we may find sequencing reads supporting a given branching edge. A paired-end sequencing read R = (R1, R2) is marked as evidence for the branching edge u ! v if it satisfies the following conditions:

(i)R spans the branching edge, meaning that at least one of the sequences R1, R2 is a subread of u and at least one of the sequences R1, R2 is a subread of v;

(ii)The sequence spanning the edge is identical to the contig sequence of the corresponding node for all variant positions it covers; (iii)R is unique for this edge: it does not satisfy conditions (i) and (ii) for any other edge involved in this branching component.

Figure 3 shows two examples of contigs creating branches in the overlap graph, along with the sequencing reads ('subreads') that were used to build these contigs; the subreads providing read evidence are highlighted in yellow. Observe that in panel A, in order to satisfy condition (iii) a subread has to cover at least one variant position on either contig. In panel B, we illustrate that also a single read end can provide evidence: the rightmost subread covers a variant position and satisfies all conditions listed above.

Note that condition (ii) ensures that erroneous contigs do not find evidence in correct reads: if a sequencing error accidentally ends up in a contig, it will cause a branch in the overlap graph which can only be supported by reads containing exactly this sequencing error. Whenever such a branch occurs, there will be insufficient evidence and hence the A B erroneous contigs will never be merged. Eventually, these contigs can be filtered out based on their short length. In the Supplementary Material we discuss how an appropriate evidence threshold is determined (using similar considerations as for determining the optimal clique size, Section 2.2).

Increasing the evidence threshold would lead to a higher accuracy but also potential loss of low abundance haplotypes.

Branching edge removal. For every branching component, we count the read evidence per branching edge and remove any edges with evidence count below the evidence threshold.

Contig extension and graph updating.

After applying the read-based branch reduction techniques described above, all branches have been either resolved or removed from the contig overlap graph. Contig extension has become an easy task: any contigs which are connected by an edge in the graph must belong to the same haplotype, and, therefore, we merge each such pair of contigs into a new, longer contig. Then, we update the overlap graph: the extended contigs become the new nodes and the edges are updated accordingly. The resulting updated graph is used for further assembly in an iterative manner, as described in Section 2.6.

Iterative procedure and diploid mode

Our workflow consists of iteratively performing the steps described in Sections 2.3-2.5, as illustrated in Figure 1. The number of edges in the contig overlap graph decreases with every iteration, since contigs connected by an edge in the graph are merged (Section 2.5). The algorithm terminates when the edge set E 0 of the updated contig overlap graph becomes empty, either upon construction or after branch reduction. Thus, our algorithm is guaranteed to converge, and once it does we remove any remaining inclusions from the final contig set. Also any contigs shorter than the fragment size of the original reads are removed from the output.

Diploid mode. Knowing that a given sample is diploid is a very strong piece of information when performing haplotype assembly. We have developed a special module which can be activated for diploid samples. It extends the POLYTE pipeline by two additional steps after the standard algorithm has terminated: construction of a diploid contig graph, followed by contig extension (see Figure 1). In these additional steps, we use the knowledge that the sample is diploid to resolve additional branches (for which there was insufficient evidence in the read set to resolve them during the read-based branch reduction step; see Section 2.4).

In overlap graphs from diploid samples we typically see two types of branching components; Figure 4 illustrates both types (Panel A and B) and gives an example of a possible collection of contigs giving rise to the corresponding branching component. In both situations we have four contigs, two from each haplotype, which have identical sequence where the contigs overlap. In diploid mode, a single read of evidence may already be considered sufficient, depending on the amount of evidence found for the other edges (Supplementary Material, Section 3). Baaijens and Schönhuth This procedure is more risky than default branch reduction (Section 2.4), since it does not require such stringent read evidence. Therefore, we always run the main POLYTE algorithm until convergence before turning to diploid mode (Figure 1). This ensures that all evidence in the original reads has been exploited first.

Results

In this section we show results for POLYTE on both simulated and real Illumina data sets and evaluate the assembly quality in terms of haplotype coverage, N50, NGA50, error rate, and the number of misassembled contigs relative to the total number of contigs. We also compare our method against alternative haplotype reconstruction tools: SPAdes [START_REF] Bankevich | SPAdes: A new genome assembly algorithm and its applications to single-cell sequencing[END_REF], Phaser [START_REF] Castel | Rare variant phasing and haplotypic expression from rna sequencing with phaser[END_REF], HapCut2 [START_REF] Edge | HapCUT2: robust and accurate haplotype assembly for diverse sequencing technologies[END_REF], WhatsHap [START_REF] Patterson | WhatsHap: weighted haplotype assembly for future generation sequencing reads[END_REF], SGA [START_REF] Simpson | Efficient de novo assembly of large genomes using compressed data structures[END_REF], and H-PoP [START_REF] Xie | H-PoP and H-PoPG: heuristic partitioning algorithms for single individual haplotyping of polyploids[END_REF]. Other polyploid assemblers [START_REF] Aguiar | HapCompass: A fast cycle basis algorithm for accurate haplotype assembly of sequence data[END_REF][START_REF] Berger | HapTree: A novel bayesian framework for single individual polyplotyping using ngs data[END_REF][START_REF] Das | SDhaP: haplotype assembly for diploids and polyploids via semi-definite programming[END_REF] were unable to process our benchmarking data due to issues with the available software. All methods were run with default settings and assembly statistics were obtained with QUAST [START_REF] Gurevich | QUAST: quality assessment tool for genome assemblies[END_REF].

Data sets

Simulated data. We generated a collection of simulated data sets of varying ploidy and sequencing depth to evaluate the effect of these characteristics. We selected four human MHC haplotypes from the Vega Genome Browser 2 : COX, DBB, MANN, and SSTO. Subsequently, we used SimSeq 3 to simulate Illumina MiSeq reads of length 2x250bp for each of those haplotypes at a coverage of 5x, 10x, 20x, 30x, 40x, and 50x, respectively, and combined the resulting read sets to form data sets of ploidy 1 (only COX haplotype, a sanity check), ploidy 2 (COX and DBB), ploidy 3 (COX, DBB, and MANN), and ploidy 4 (all).

Real data. For evaluation on real sequencing data, we considered a data set from phase 3 of the 1000 Genomes project (1000 Genomes Project Consortium et al., 2012;[START_REF] Sudmant | An integrated map of structural variation in 2504 human genomes[END_REF] for individuals NA19240. This data set was obtained from a 2x250 bp PCR free Illumina protocol, sequenced to a coverage of 28-68x. Full haplotypes have been reconstructed for this individual as part of a recent study [START_REF] Chaisson | Multi-platform discovery of haplotype-resolved structural variation in human genomes[END_REF] using various specialized sequencing techniques and reconstruction algorithms; we use the resulting haplotypes as a ground truth for a wholechromosome benchmarking experiment on chromosome 22.

Alignments and variant call sets for reference-guided methods. Reference-guided methods Phaser, HapCut2, WhatsHap, and H-PoP require as input a reference genome, read alignments to the reference genome, and a pre-computed set of genomic variants. For the simulated data we performed read alignment to the GRCh38 reference genome using BWA MEM [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF]. The real data was already provided as alignments to the GRCh37 reference genome, also obtained with BWA MEM. We extracted the sequencing reads corresponding to chromosome 22 from the provided BAM files. Finally, we performed variant calling on all data sets with FreeBayes 4 .

Assembly performance criteria

We evaluate assembly performance in terms of several statistics commonly used for de novo assembly evaluation, as reported by QUAST. Haplotype coverage (HC). The completeness of the assembly is measured by the fraction of nucleotides in the target haplotypes (ground truth) covered by haplotigs, referred to as the haplotype coverage. N50 and NGA50. Assembly contiguity is measured using the N50 value, which is defined as the length for which the collection of all contigs of that length or longer covers at least half the assembly. The NGA50 measure is computed in a similar fashion, but only aligned blocks are considered (obtained by breaking contigs at misassembly events and removing all unaligned bases). This measure reports the length for which the total size of all aligned blocks of this length or longer equals at least 50% of the total length of the true haplotypes.

Error rate (ER) and N-rate (NR). We evaluate error rate as the sum of mismatch rate and indel rate when comparing to the ground truth haplotype sequences. In addition, we report the relative number of ambiguous bases ('N's), referred to as the N-rate.

Misassembled contigs (MC).

A contig or haplotig is called misassembled if it contains at least one misassembly, meaning that left and right flanking sequences align to the true haplotypes with a gap or overlap of more than 1kbp, or align to different strands, or even align to different haplotypes. We report the proportion of misassembled contigs.

Benchmarking results

We performed benchmarking experiments on one of the simulated MHC data sets described above (ploidy 2, 20x coverage per haplotype) to compare a variety of haplotype reconstruction tools. In addition, we ran all methods on the chromosome 22 data of the 1000 Genomes individual NA19240. The assembly statistics on both data sets are shown in Table 1. Since both data sets are diploid, we present results for SPAdes in regular mode and in diploid mode, referred to as SPAdes-dip.

In both experiments, we observe that across all methods POLYTE has the largest haplotype coverage (HC,92.4% and 78.2% for MHC and chr22,respectively). In other words, it reconstructs the largest fraction of the true haplotype sequences. In comparison, the other methods are all more or less on a par , respectively). On the real data the haplotype coverage achieved by all methods is rather low; this can be explained by only 86.4% of the target haplotypes being covered by sequencing reads. After normalizing the haplotype coverage values by 86.4, POLYTE achieves a haplotype coverage of 90.5%.

In terms of assembly contiguity, indicated by high N50 and NGA50 values, reference-guided methods (HapCut2, Phaser, WhatsHap, H-PoP) perform better than de novo assemblers (POLYTE, SGA, SPAdes). This reflects a common advantage of reference-guided approaches, which can make use of the external information to bridge regions only poorly covered with informative reads, if appropriate. The increase in length, however, is offset by a substantial decrease in terms of haplotig quality: referenceguided approaches exhibit both substantially more misassemblies (which in particular can lead to severe issues in downsteam interpretations) and increased error rates, here larger by one to two orders of magnitude. Note that several NGA50 values are undefined ('-'), because the aligned blocks are unable to cover at least 50% of the total reference length.

Another important difference between reference-guided methods and de novo approaches is reflected in the N-rates on the real data: the reference genome contains several stretches of ambiguous nucleotides ('N's), which the reference-guided methods cannot correct. De novo approaches, on the other hand, can potentially uncover the true sequence behind these ambiguous regions and show an N-rate of 0% (versus 0.9% for the reference-guided methods).

Between de novo approaches, we compare POLYTE with SGA and SPAdes and observe that POLYTE reconstructs a substantially larger fraction of the true haplotypes. Although SPAdes achieves better N50 values, this comes at the expense of a decrease in terms of error rate and misassemblies, also reflected in a low NGA50 value on the simulated data and the NGA50 being undefined on the real data (see explanation above). On the simulated data set, POLYTE and SPAdes achieve comparable error rates of 0.035% and 0.031%, respectively. On the real data we notice an advantage for POLYTE, with an error rate of only 0.090% compared to 0.131% for SPAdes. In addition, POLYTE is less vulnerable to misassemblies than SPAdes on real data, with 0.2% versus 0.6% MC. SGA is able to reconstruct highly accurate contigs with slightly lower error rates than POLYTE (0.025 vs 0.035% [MHC] and 0.069 vs 0.090% [Chr22], respectively), but covers a significantly lower fraction of the ground truth haplotypes (73.4 vs 92.4% [MHC] and 57.7 vs 78.2% [Chr22], respectively).

In an overall account, we believe that, arguably, the major advantage of POLYTE is established by the increase of 10-15% over the other approaches in terms of haplotype-specific coverage, in combination with the error rates, which are clearly lower than those of the other tools.

In terms of runtime and memory usage, de novo approaches are in general more expensive than reference-guided methods. We also observe this when comparing CPU time and peak memory usage (Supplementary Tables 345). Reference-guided methods have CPU times that are orders of magnitude less compared to de novo methods (where POLYTE requires 9-15 times more (resp. 3-6 times more) runtime and 3 times less (resp. 12-17 times more) memory than SPAdes and SGA, respectively). It is important to notice, however, that these de novo assemblers are highly parallelizable; we demonstrate the effect of increasing the number of available CPU's on the effective runtime in Supplementary Table 6. This leads to feasible runtimes on multi-core computing facilities in practice.

Effect of ploidy and sequencing depth

To study the effect of genome ploidy and sequencing depth on the assembly quality and completeness, we ran POLYTE, SPAdes, SGA, and H-PoP on all simulated data sets described in Section 3.1 (other tools were unsuitable for polyploid genomes). Figure 5 shows the results for the 5x, 10x, and 20x data sets in terms of haplotype coverage (HC), N50, NGA50, error rate (ER), and misassembled contigs (MC). For additional result tables we refer the reader to the Supplementary Tables 891011.

We observe that POLYTE excels regarding haplotype coverage, with advantages becoming more distinct as the ploidy increases. SPAdes and H-PoP achieve more contiguous assemblies (higher N50 values) but, as we already observed on diploid data, this comes at the cost of significantly higher error rates and misassemblies. SGA performs very similar to POLYTE when considering N50, ER, and MC, but obtains much lower HC values. The NGA50 values highlight the improved assembly quality of POLYTE over SGA and SPAdes: while POLYTE achieves NGA50 values comparable to the N50, SGA and SPAdes are unable to cover at least 50% of the ground truth with alignments (hence NGA50 is undefined). Overall, we conclude that in polyploid settings the same advantages of POLYTE apply as in diploid settings -increased haplotype-specific coverage in combination with low error rates -and become even more pronounced.

All other methods evaluated (HapCut2, Phaser, Whatshap, SPAdesdip) are designed for diploid data, so for those we could only assess the effect of sequencing depth. Results indicate that each of the referenceguided methods already performs optimally at a coverage per haplotype of 5x. Moreover, these methods are unaffected by a further increase in sequencing depth (see Supplementary Table 12). SPAdes in diploid mode (SPAdes-dip) performs optimally at a per-haplotype coverage 20x.

Discussion

Assembling the individual haplotypes of an organism from sequencing reads is known as haplotype aware genome assembly and plays a major role in various disciplines, including genetics and medicine [START_REF] Glusman | Whole-genome haplotyping approaches and genomic medicine[END_REF][START_REF] Tewhey | The importance of phase informationn for human genomics[END_REF]. Computing haplotype-specific pieces of sequence, also known as haplotigs, is a difficult task. Algorithms addressing this task do not only need to distinguish between sequencing errors and true variants, but also need to assign the true variants to the individual haplotypes. Enormous quantities of next-generation sequencing (NGS) reads generated worldwide have not been fully exploited in terms of haplotig computation, because methodology for de novo haplotig computation from NGS reads has been in a rather immature state.

We have presented POLYTE (POLYploid genome fitTEr) as a new approach to de novo assembly of haplotigs from NGS data, suitable for diploid genomes as well as genomes of higher ploidy. Unlike the majority of NGS based de novo assemblers, our method follows the overlaplayout-consensus (OLC) paradigm to achieve enhanced performance rates in terms of haplotype-specific computation of contigs. In order to appropriately distinguish between errors and true variants to be assigned to haplotypes, it employs an iterative OLC scheme. Along the iterations, contigs grow in length while preserving their uniqueness in terms of haplotype identity. As a result, POLYTE outperforms the currently available state-of-the-art approaches for haplotig computation, where it performs particularly favorable in terms of quantities that refer to haplotype-specific reconstruction of the genomes.

Experimental results showed that POLYTE can build accurate assemblies from Illumina MiSeq reads (2x250bp) for data sets of varying ploidy (di-, tri-and tetraploid), with results for tetraploid organisms almost on a par with those for diploid organisms. Although building overlap graphs for larger genomes remains a challenge, we provide a read binning step that allows efficient assembly by splitting the work over multiple cores. The typical use-case for POLYTE consists of Illumina NGS reads for a specific gene, region or genome that is highly polymorphic (of any ploidy greater than one).

We showed that POLYTE succeeds in accurate reconstruction of individual haplotypes of the human MHC region. Future work may Baaijens and Schönhuth therefore be to apply POLYTE to NGS data in population-scale human genome projects, where individual genomes are still lacking proper annotation of their MHC region, which applies in the majority of cases.

Advantages become particularly distinct on data of higher ploidy, leading to plant genome assembly as another interesting future application of POLYTE.

Our algorithm is, in its essence, generic in the choice of input reads, so applying it for TGS reads essentially is a matter of adapting parameters, which we will explore in the short-term future.

Fig. 1 .

 1 Fig. 1. Algorithm overview.

Fig. 2 .

 2 Fig. 2. Illustration of branching components in a contig overlap graph. Edges of the same color belong to the same branching component.

Fig. 3 .

 3 Fig. 3. Two examples of contigs creating branches in the overlap graph. Edges corresponding to true haplotypes are highlighted in yellow. The corresponding subreads are aligned below, those providing read evidence are again highlighted.

Fig. 4 .

 4 Fig. 4. Typical branching components in diploid assemblies: four contigs, two from each haplotype, having identical sequence in their overlap. Depending on the contig lengths, all contigs overlap (panel A) or only a subset of the contigs overlap (panel B).

2

  http://vega.

Fig. 5 .

 5 Fig. 5.Assembly results per method for simulated data of increasing ploidy (k=2,3,4) and per-haplotype coverage (5x,10x,20x). N50 and NGA50 values are plotted on a log-scale for increased readability. For SGA (k=2,3,4) and SPAdes (k=3,4) the NGA50 values are undefined.

Table 1 .

 1 Benchmarking results, HC = Haplotype Coverage, ER = Error Rate (mismatches + indels), NR = N-Rate (ambiguous bases), MC = Misassembled Contigs. Top: simulated diploid data for the MHC region. Bottom: real data for chromosome 22 of 1000 Genomes individual NA19240. HC values within parentheses indicate haplotype coverage relative to the amount of bases covered by sequencing reads.

		archive.ensembl.org/info/data/MHC_
	Homo_sapiens.html
	3	https://github.com/jstjohn/SimSeq
	4	https://github.com/ekg/freebayes

Depending on the context, polyploid includes diploid, but here we refer to polyploid as more than two copies.

Funding

This work was supported by the Netherlands Organisation for Scientific Research (NWO) through Vidi grant 679.072.309.

POLYTE is freely available as part of the HaploConduct package at https://github.com/HaploConduct/HaploConduct, implemented in Python and C++.


Kamer Kaya
email: kaya@sabanciuniv.edu

Johannes Langguth
email: langguth@simula.no

Ioannis Panagiotas
email: ioannis.panagiotas@ens-lyon.fr

Bora Uçar
email: bora.ucar@ens-lyon.fr

Karp-Sipser based kernels for bipartite graph matching

come   L'archive ouverte pluridisciplinaire

Introduction

Data reduction, or preprocessing, is applied in different problems to shrink the size of the instance before actually solving the problem. The aim is to significantly reduce the size with fast methods so that the exact algorithm whose run time complexity might be high can benefit from the reduction. In this paper, we focus on data reduction algorithms in the context of finding maximum cardinality matchings in bipartite graphs. Karp and Sipser [START_REF] Karp | Maximum matching in sparse random graphs[END_REF] describe two reduction rules for the maximum cardinality matching problem in unweighted, simple bipartite graphs. The first rule matches a vertex with a single edge to its unique neighbor. The second rule takes a vertex v with two neighbors u and w, removes v and merges u and w into a single vertex. The first rule is easy to implement efficiently, while a straightforward implementation of the second reduction rule can take quadratic time, which is potentially higher than the O(m √ n) time the exact matching algorithm [START_REF] Hopcroft | An n 5/2 algorithm for maximum matchings in bipartite graphs[END_REF] takes for sparse bipartite graphs with m edges and 2n vertices equally distributed to each part. When used as a heuristic for the matching problem, a random edge is added to the matching when the reduction rules do not apply. When an edge u, v is added to the matching, other edges incident on u and v are discarded. This reduces the degree of the vertices adjacent to u and v and hence creates opportunities to apply the reduction rules.

Our focus in this paper is to design and analyze efficient algorithms for implementing both reduction rules for the matching problem. To see the impact of data reduction we conducted several experiments using only the first rule and also, both rules. We also investigate these alternatives when Karp-Sipser is used as a cheap matching heuristic. Note that in the former model, the matching decisions are concordant with at least one maximum cardinality matching. However in the latter case, there may not be a maximum cardinality matching obeying all these decisions.

The paper is organized as follows: Section 2 introduces the notation and the background on the Karp-Sipser heuristic. The original heuristic with two rules as proposed is discussed in Section 3 with various approaches for implementation. Section 4 presents the experimental results and Section 5 summarizes the literature focusing on the original Karp-Sipser heuristic. Section 6 concludes the paper.

Notation and background

Let G = (V A ∪ V B , E) be a simple bipartite graph with |V A | = |V B | = n and |E| = m. Two vertices are neighbors in G, if they share an edge. A vertex with k neighbors is called a degree-k vertex. The degree of a vertex v is shown as d v . A matching in a graph G is a set of disjoint edges, such that no two of them share a common vertex. We use M(G) to refer to the maximum cardinality of a matching in G. If M(G) = n, a maximum cardinality matching is called perfect.

For a graph G, there exist polynomial time algorithms to find a matching with cardinality M(G). In practice, one very popular strategy to find M(G) is to use a two-step process. In the first step, a cheap initialization heuristic is used in order to quickly find a matching of large cardinality which is usually maximal i.e. cannot be extended via simple edge additions. This matching is given as an input to the second step and improved via a (relatively) more expensive exact algorithm. It has been empirically shown that this strategy is much faster than running an exact algorithm from the beginning-see for example Langguth et al. [START_REF] Langguth | Heuristic initialization for bipartite matching problems[END_REF].

While there exist different initialization heuristics for the first step, the Karp-Sipser (KS) heuristic empirically performs better than many, if not all, on aver-age [START_REF] Duff | Design, implementation, and analysis of maximum transversal algorithms[END_REF][START_REF] Langguth | Heuristic initialization for bipartite matching problems[END_REF]. In this paper, we will focus on this heuristic and its implementation. Karp-Sipser is based on performing reductions on a graph with no degree-0 vertices (they are discarded throughout); when a degree-1 or degree-2 vertex appears, KS reduces the problem to a smaller one via the following rules:

• Rule-1: At any time, if a degree-1 vertex u with neighbor v appears the edge {u, v} is added to the matching and both vertices are removed from the graph. This decision is optimal in the sense that there exists at least one maximum cardinality matching in the current graph containing {u, v}.

• Rule-2: At any time, if there are no degree-1 vertices, and a degree-2 vertex u with neighbors v and w appears, u and its edges are removed from the current graph, and v and w are merged to create a new vertex vw whose set of neighbors is the union of those of v and w (excluding u). Karp and Sipser showed that a maximum cardinality matching for the reduced graph can be extended to obtain maximum cardinality matching for the original graph by matching u with either v or w depending on vw's match.

• When none of these rules can be applied, a random edge from the current graph is added to the matching. This decision may not be optimal. Therefore, this step is not performed in the data reduction setting.

Both Rule-1 and Rule-2 have the property that they preserve M(G). We will use the notation G (0) = G to denote the initial graph, and G (t) to denote the graph after t random or rule-based decisions. Let G (k) be the first graph where neither Rule-1 nor Rule-2 is applicable. We call G (k) a kernel of G for the maximum cardinality matching problem, i.e., a reduced, smaller graph where one can obtain a maximum cardinality matching for G given a maximum cardinality matching for G (k) by following the reductions in reverse order. Thus, in addition to the well known initialization strategy, another use of the Karp-Sipser algorithm is to obtain G (k) initially, and then apply an exact algorithm on the smaller subgraph G (k) , rather than G. However, obtaining G (k) can be computationally expensive.

Rule-1 is simple to implement. However, Rule-2 is more complicated and requires more effort to be used in a matching heuristic. That is why in practice, the KS heuristic has usually been associated with its simpler variant KS R1 , which only applies Rule-1, and when necessary, random choices. Although it does not exploit the second rule, KS R1 has been shown to obtain large Figure 1: A toy bipartite graph G = (VA ∪ VB, E) with vertex sets VA and VB. Except for the leftmost vertices from each vertex set, each vertex has two neighbors and will be removed by Rule-2 of KS. For every G (t) , the merged vertices will have degrees (n + 1 -t) and 2. Hence, the total time for the naive KS implementation on this instance is Θ(n 2 ).

cardinality maximal matchings on real life graphs, see e.g. [START_REF] Duff | Design, implementation, and analysis of maximum transversal algorithms[END_REF][START_REF] Kaya | Pushrelabel based algorithms for the maximum transversal problem[END_REF].

3 Karp-Sipser: The Heuristic with two rules KS R1 can be implemented in O(n + m) time by keeping an up-to-date list of all degree-1 vertices. When the edge {u, v} is added to the matching, if u is a degree-1 vertex, one needs to visit only v's neighbors to reduce their degrees. Otherwise, u's and v's neighbors must be visited and their degrees are updated. Hence, with only Rule-1, each adjacency list is accessed at most once and KS R1 runs in linear time. This implies that the complexity of KS depends on the the cost of applying Rule-2.

For any graph G (t) that appears during the execution of the heuristic, let u be a degree-2 vertex with neighbors v and w. According to Rule-2, v and w must be merged to vw ∈ G (t+1) . With a naive approach, this operation takes Θ(d v + d w ) time. Since there can be at most n merge operations, the time complexity of this strategy is O(n 2 ). This bound is tight as one can create highly sparse graphs for which KS requires quadratic time as shown in Figure 1.

3.1 An expected O(m log n) time algorithm. Let u be a degree-2 vertex with neighbors v and w. For KS, a merge operation due to u performs three steps:

1. Merging v's adjacency list with that of w, 2. Performing degree reductions for the vertices that are both v's and w's neighbors, 3. Updating the adjacency lists of the vertices which are neighbors of v and w.

Assume that it is possible to perform this merge operation in O(min(d v , d w )) time. We will label the KS variant with this merge complexity as KS-min and show in the following theorem that its worst-case time complexity is O(m log n) which is better than O(n 2 ) for sparse graphs where m n 2 .

Theorem 3.1. KS-min runs in O(m log n) time.

Proof. Consider a hypothetical algorithm KS * -min which operates on a multi-graph where multiple/parallel edges are allowed between two vertices while merging the vertices via Rule-2. Hence the adjacency lists contain parallel edges. This means that KS * -min keeps all the edges in the graph at any time step and simply combines the adjacency lists of the two vertices during a merge operation. Similar to KS-min, we assume that the cost of each merge is equal to the size of the smaller adjacency list. Hence, for a single merge, the amortized cost per edge in the smaller list becomes O(1).

Let v and w be the vertices in G (t) and d v , d w with d v ≤ d w be the number of edges (parallel edges are allowed) incident on v and w. Let vw be the merged vertex in G (t+1) . Note that d vw is at most m, the total number of edges. We also have 2d v ≤ d vw , since we allow parallel edges. Therefore, a single edge can be in the smaller adjacency list in at most log m merge operations. Having O(1) amortized complexity for each such edge, the total time complexity of KS

* -min becomes O(m log m). Since m is O(n 2 ), the complexity is O(m log n).
To analyze the complexity of KS-min, assume the same merge operations are followed by KS * -min and KS-min. Since the sizes of the merged adjacency lists in KS-min are always smaller than or equal to the corresponding lists in KS * -min, the time complexity of KS-min is also O(m log n).

It remains to discuss how to implement a merge operation on v and w in KS-min. Assume d v ≤ d w . First, the larger list, i.e., w's list is kept intact and w becomes the merged vertex vw in the reduced graph. Then each neighbor x of v is processed one after another. To keep the graph simple and to correctly update the degree reductions for the common neighbors of v and w, x is first searched in the adjacency list of w. If x is already in w's list its degree is decremented and v is removed from x's adjacency list. Otherwise, x is inserted to w's list. Furthermore, v is also replaced by w in x's list.

To perform these operations in expected constant time per edge, we can use a hash table to store all the edges in the graph. This hash table is used to query the existence of x in w's list by looking whether {x, w} exists or not in the hash table. With this data structure, we have O(m log n) expected time complexity for KSmin and linear space. Instead of a hash table, one can use a data structure to store the edges with O(log m) insertion, update and lookup time such as a binary search tree. Such a data structure yields an O(m log 2 n) worst-case time complexity for KS-min in linear space.

3.2 An implementation with list caching. In order to merge the adjacency lists of two vertices v and w, one can use a dense, 0-1 array L of size n. This array represents v's adjacency list with L[x] = 1 iff {v, x} exists in the current graph. Once such a list is created for one of the merging vertices, the edges of the other vertex can be looked up in L. Note that this array needs to be created every time there is a merge involving v and each time, one needs to re-iterate over v's adjacency list.

A natural optimization to this approach would be to allocate such arrays for some vertices and persistently keep them in memory, i.e., cache them. With caching, we do not have to re-iterate over these vertices' adjacency lists each time they participate in a merge. Indeed, if we cache two arrays for the two leftmost vertices in Figure 1, the complexity drops from Θ(n 2 ) to Θ(n). We refer to this variant as KS cache .

There are many different strategies to decide on which vertices to keep in the cache. For example, one can cache the lists for the k most recently merged vertices, or k highest degree vertices. In Theorem 3.2, we show a negative result holding for any arbitrary caching strategy. Theorem 3.2. KS cache using k arrays and applying Rule 1, Rule 2 and random decisions has an instance requiring Ω n 2 /k time for all possible caching policies.

Proof. Assume that we have enough memory to cache the adjacency lists of k vertices in dense form. Let G be the n × n bipartite graph shown in Fig. 2 in which the vertices in V A and V B are shown in black and red, respectively. The graph contains 2k identical subgraphs, each having

( n 2k -1) 2
two-by-two complete bipartite structures and an additional V A (V B ) vertex that is connected to every other V B (V A ) vertex in its subgraph. These extra vertices are labeled as a i and b i for the ith subgraph (in fact, they correspond to the left most vertices in Fig. 1). As the figure shows, the minimum degree in the initial graph is three and a random decision is required.

When a random decision hits into a two-by-two part in the ith subgraph, the corresponding matching and the removal of the matched vertices yield two possible merges involving either a i or b i . These potential merges eliminate each other, i.e., only one of them is possible. Furthermore, after the merge operation, the minimum degree in the graph will go back to three. Hence, assuming all the random edges are from the two-bytwo bipartite subgraphs, each random choice is followed by a merge which is then followed by another random decision. Since the decisions are random, we cannot build a caching strategy based on them.

Hence, regardless of caching policy, it is always possible that the dense form of both a i and b i is not in the cache until at least k components are reduced completely. Doing so for only one component without touching to others takes

( n 2k -1)/2 j=0 n 2k -2j = Θ n 2 k 2
time since the degrees of a i (as well as b i ) vertices are reduced by two after each step. Since at least k components must be consumed, the complexity until this point is Ω n 2 /k which implies a lower bound on the execution time of KS cache for the instance in Fig 2

assuming that the random choices always hit two-bytwo subgraphs.

An alternating component approach.

Here, we consider and extend an idea that was first introduced in [START_REF] Langguth | Heuristic initialization for bipartite matching problems[END_REF]. There, the authors used the term alternating component to refer to a connected subgraph consisting of a set of boundary vertices connected by paths. The boundary vertices are vertices of arbitrary degree greater than one. The paths are required to contain an even number of edges, have a boundary vertex as beginning and endpoint, and all other vertices in each path are required to be of degree 2. An example can be seen in Fig. 3. In addition, we only consider components whose paths are maximal, i.e., they cannot be extended without violating the above definition.

If an alternating component contains a cycle, then there is no need to explicitly merge vertices in the component. Since every edge has at least one incident vertex of degree 2, at least half of the vertices in such a cycle have degree 2. Thus, for any edge of the cycle, there is a maximum cardinality matching containing that edge. Thus, we can pick an arbitrary edge and match it. The remaining component will be matched via the application of Rule-1.

If a component is acyclic, it can be immediately merged into a single vertex via Rule-2. To see this, consider two vertices that are both connected to a degree-2 vertex in a simple alternating component. Thus, each application of Rule-2 will reduce the length of a path between boundary vertices in an alternating component by two, and because the length is even, they will be merged when it reaches zero. The same applies to any such path in an alternating component. The advantage of doing so is that we do not have to perform merges immediately. Instead, we maintain components whose merges can then be performed together more efficiently.

In [START_REF] Langguth | Heuristic initialization for bipartite matching problems[END_REF], only the first operation was performed. Acyclic components were maintained, but not explicitly merged. This resulted in a useful heuristic, but it is not equivalent to an implementation of Rule-2. Here, we show how to use alternating components to implement the full KS algorithm. We will call this version KS comp . Similar to the basic implementation, the algorithm keeps stacks of degree-1 and degree-2 vertices. If a degree-1 vertex exists at any time in the remaining graph, we add its edge to the matching. For a degree-2 vertex u, we create P (u), a maximal even-length path containing u and only degree-2 vertices except for the first and last vertex in P (u). Let l, r be these boundary vertices in P (u). We use a union-find data structure to keep track of membership in components. At the start of the algorithm, every vertex is its own component. Now, if l and r already belong to the same component C, then adding the additional path must close a cycle in C. In that case we arbitrarily add one of the edges of u to the matching and the component is then cleared out by Rule-1 reductions. If no circle exists, then we simply update the labels of l, r to signal that their components have been joined together. No merge is performed at this point.

When there are no degree-2 vertices left, we merge the remaining acyclic components one by one. We refer to this as a component shrink. Let C be component with boundary vertices b 1 , . . . , b k . To perform the shrink we only require O( bi∈C d bi ) since we only need to access the adjacency list of each b i once.

If during these merges, a new degree-2 vertex becomes available (this can happen if b i and b j are both adjacent to some vertex w with d(w) = 3 before the shrink) the merging stops and the degree-2 operation is applied, which can connect two components. If all components have been shrunk, and no degree-1 or degree-2 vertices exists, KS comp selects a random edge to be matched.

Using alternating components can prove useful in the sense that they provide a structured way to perform reductions such that they avoid re-reading adjacency lists of vertices (which was also the motivation of KS cache ). For example, it takes only linear time for the worst case instance of Figure 1. However, it can degrade to the worst case performance of the basic implementation if there are no nontrivial alternating components in the graph:

Theorem 3.3. The KS comp variant requires O(n 2 ) in the worst case.
Proof. The proof is similar to that of Thm. 3.2 hence we skip the details. In short, if the algorithm is given an input instance as in Fig. 2, there will be no nontrivial components. Assuming all random choices hit into the two-by-two complete bipartite subgraphs, there will always be a random choice following a single merge. Hence, on this instance, the component approach behaves just like the basic implementation.

We note that the approach of using alternating components bears some similarity with the approach used in [START_REF] Bartha | A depth-first algorithm to reduce graphs in linear time[END_REF], which performs all currently existing Rule-2 operations in linear time. However, the authors consider solely Rule-2 and do not allow neither Rule-1 nor random decisions. They propose a complicated static linear-time algorithm, operating on the DFS-tree of a given graph and avoiding to do merges explicitly by doing implied reductions. Let us refer to a phase, as a series of consecutive Rule-1 or random operations done by KS. Hence, phases are separated from each other by one or more Rule-2 applications. Their O(n + m) time algorithm would need to be rerun between the different phases and thus can similarly potentially yield a complexity of O(n 2 ) .

3.4 Fast recovery of the matching. As described above, a merge operation creates a new vertex that may be involved in upcoming merges. Hence, the matched vertices returned by KS do not necessarily appear in G (0) . The standard way to recover the actual matching is by keeping a stack S containing all the information related to each merge [START_REF] Magun | Greeding matching algorithms, an experimental study[END_REF]. After the heuristic, the stored merges are popped one after another from S.

Assume that during the heuristic a degree-2 vertex u appears with neighbors v and w; hence, a new vertex vw is created. While recovering the actual matching, if vw is unmatched we can match either v or w with u and continue. Otherwise, we need to check whether vw is matched with one of v's or w's neighbors stored in S. Using similar reasoning as in Theorem 3.1, this approach takes O(m log n) time and space in the worst case, assuming the smaller neighborhood (among v's and w's) is stored in the stack. Here, we propose an algorithm with O(n) time complexity.

To recover the matching efficiently, we use a graph F ⊆ G. Initially, F has no edges but contains all vertices of G. During the heuristic, assume a merge is being performed due to a vertex having two edges e and e . Let {u, v} and {u , w} be the original endpoints of e and e in G, respectively. We add the edges e 1 = {u, v} and e 2 = {u , w} to F and set twin(e 1 ) = e 2 and twin(e 2 ) = e 1 . The association of e 1 and e 2 as twins exploits the implicit Rule-2 property that either e 1 or e 2 can be in a maximum cardinality matching, but not both.

Proposition 3.1. F is a forest.

Proof. Assume that F contains a cycle. Let C ⊆ F ⊆ G be a simple cycle with length . Let C be first initiated with a merge in G (t) , and let u be the degree-2 vertex for that merge with neighbors v and w. Both {u, v} and its twin {u, w} must be in C, since u will not exist in G (t+1) . Let C be the reduced cycle obtained by removing u from C and merging v and w. The above steps can be repeated by using the edge that initiates C . At each step, the cycle will be further reduced and after l -1 steps, only two vertices and a single edge will remain. As parallel edges are not allowed in the graph G (•) , such a merge cannot exist to complete the cycle.

Let M be the set of matched edges found by KS and let M be the the matching constructed with the original versions of M 's edges from G (0) . Let e = {u, w} be an edge in M . Then for all edges e = e of u and e = e of w in F , we add e * = twin(e ) and e * * = twin(e ) to M . We then remove u and w from F (i.e., F = F \ {u, w}) and consider the next matched edge in M . If no more nodes in the current tree are matched, we select a vertex of degree-1 in F , add its unique edge to M and repeat.

Lemma 3.1. The algorithm described above to recover the matching is correct.

Proof. Let e and e be two twin edges which have not been included in the matching M . If one of the endpoints of e is matched in M , then it is no longer possible to include e in M . By Rule-2, one of the twins must necessarily participate in a maximum cardinality matching. Hence, we can extend M by using e's twin edge. This means that as long as there exist matched vertices in F , the algorithm correctly extends the matching.

If there are no matched vertices in F , for any twins e and e , we know that either can be in a maximum cardinality matching. Hence, we can arbitrarily add either to M . For simplicity, we chose the edge of a degree-1 node, which always exists since F is always a forest.

Experiments

We implemented the algorithms discussed in the last section. Except where noted, the Karp-Sipser algorithm variant being used is the default one. When KS is used as a heuristic, we first generate a random permutation of the edges with uniform probability, which we store in a list. Then, in the event of a random decision, the first available edge (that is without any matched endpoints) is returned from this list.

All the codes are written in C++ and compiled with gcc 8.3.0 with -03 optimization flag. They are tested on two different machines; The first machine had 2 x AMD EPYC 7551 CPUs and 256 GB RAM (Arch 1). The other had 4 x Intel Xeon E7-8890 CPUs and 1.5 TB RAM (Arch 2).

Experiments with real graphs:

For the reallife dataset, we used bipartite graphs with n = |V A | = |V B | and 10 6 ≤ n ≤ 10 7 from the SuiteSparse collection [START_REF] Davis | The University of Florida sparse matrix collection[END_REF]. There were 93 such graphs with 3 × 10 6 ≤ m ≤ 3.1 × 10 8 . For each of these matrices, we performed a total of sixteen runs. More specifically, we randomly permuted each matrix four times, and then measured the timings for both KS R1 and the default KS algorithm four times in the permuted matrix. Permuting the matrix between runs can potentially affect the run-time for both KS heuristics as well as the exact algorithm used afterwards. To analyze the impact of adding Rule-2 for data reduction, two exact matching algorithms pushrelabel (PR) [START_REF] Goldberg | A new approach to the maximum-flow problem[END_REF] and Pothen-Fan+ (PF) [START_REF] Duff | Design, implementation, and analysis of maximum transversal algorithms[END_REF][START_REF] Pothen | Computing the block triangular form of a sparse matrix[END_REF] are used. In previous studies [START_REF] Kaya | Pushrelabel based algorithms for the maximum transversal problem[END_REF], these two algorithms are shown to be the best performing algorithms from a set of alternatives, including the asymptotically fastest one [START_REF] Hopcroft | An n 5/2 algorithm for maximum matchings in bipartite graphs[END_REF].

Using KS for kernelization:

We first consider the use of KS as a kernelization tool and run experiments on Arch 1: we first run KS and KS R1 , with no random decisions, to find their kernels G (k) and G (k ) where k ≥ k . After that, an exact algorithm, PR or PF, is executed to find a maximum cardinality match-ing. As KS requires a minimum degree ≤ 2, for this experiment, we only use the graphs having min degree one or two. There are 59 such graphs. The results are summarized in Figure 4.

To measure the impact of the second rule, we measured the kernel quality for each heuristic. The quality is computed as the number of times a rule is applied during KS and KS R1 normalized with respect to size of the maximum cardinality matching. Figure 4(a) shows that in terms of quality, for around 30 matrices, the second rule has an impact and for one matrix, it increased the quality from almost 0 to 0.65. Hence, KS obtains a kernel of smaller size compared to KS R1 . Furthermore, as Fig. 4(b) shows, KS is slightly slower compared to KS R1 differing in most graphs by less than a second.

To evaluate the impact of additional size reduction obtained via Rule 2 on the exact matching process, we measured the execution times of PR and PF given the kernels G (k) and G (k ) obtained by KS and KS R1 . Figure 4(c) shows the speedup obtained, i.e., the execution time of PR/PF on G (k ) divided by the execution time of the same algorithm on G (k) . In the figure, if both PR and PF speedups are positive (negative) for the same matrix the larger (smaller) one is shown on top of the other. The first observation is that the speedups are mostly in the positive side and a good kernel can significantly improve the execution time of the matching phase. The second observation is that the impact of kernelization heavily depends on the algorithm used since PR and PF can behave in a different way. For instance, in one matrix, although G (k) is smaller than G (k ) , PF runs slower on G (k) . Yet, even with such fluctuations, for both algorithms, smaller kernels usually yield better performance.

Using KS as a heuristic:

To check if KS with Rule-2 is also useful as an initial, cheap matching heuristic, we let KS and KS R1 run in their entirety and input the returned matchings to the exact algorithms. The experiments are performed on Arch 2. Figure 5 summarizes the results of these experiments. Both KS R1 and KS obtain excellent results with matching quality almost always larger than 0.97 of M(G). KS is able to match 1% more vertices than KS R1 , which corresponds to more than ten thousand vertices due to the size of our instances.

As in the previous experiments, the exact algorithms tend to be faster when initialized with KS. However, since the impact, i.e., the difference in the quality, is higher for kernelization compared to matching, the speedups are not as large as the previous set of experiments. If we look at the total run time to find a maxi- The kernel quality of KS and KSR1 for the 59/93 graphs having at least a single vertex with degree strictly less than three. The quality is measured as the number of rules applied during a heuristic normalized with respect to size of the maximum cardinality matching. (b) The execution times (in sec.) of these two heuristics are also shown. (c) The speedups for the exact matching algorithms PR and PF when the KS kernel is used instead of KSR1. A value higher than one is in favor of Rule 2 implementation. In all figures, the x-axis, i.e., matrices, is organized with respect to non-decreasing kernel quality for KS. The matching quality of KS and KSR1 for all 93 graphs. The quality is measured as the number of decisions, including random ones, applied during a heuristic normalized with respect to size of the maximum cardinality matching. (b) The execution times (in sec.) of these two heuristics are also shown. (c) The speedups for the exact matching algorithms PR and PF when KS is used instead of KSR1. A value higher than one is in favor of Rule 2 implementation. In all figures, the x-axis, i.e., matrices, is organized with respect to non-decreasing matching quality for KS.

mum cardinality matching, we see that on average using KS rather than KS R1 can yield a 20% slow-down using PF and 40% with PR. This is understandable since KS itself does significantly more work than its simpler variant. On the other hand there exist instances where using KS rather than KS R1 can be crucial. For example, consider the graph com-LiveJournal for which PF, initialized with KS R1 's matching, results in a run time of about 1060 seconds on average, whereas PF requires 65 seconds on average, if initialized with KS. That is why we are also interested in parallel implementations of Rule-1 and Rule-2 which will potentially bring the performance of KS closer to that of KS R1 and hence reap more benefits from the speedups seen in Figure 5(c).

Experiments with synthetic datasets:

Here we focus on special, synthetic instances, for which Rule-2 increases the probability of finding a maximum cardinality matching. Similar examples have also been described in [START_REF] Anastos | Finding perfect matchings in random cubic graphs in linear time[END_REF].

Optimal instances for KS:

The first dataset is a family of synthetic graphs where Rule-2 is sufficient to find a maximum cardinality matching.

Let G = (V A ∪ V B , E) be a bipartite graph with |V A | = |V B | = n. For any 1 ≤ i ≤ j ≤ n, there is the edge {a i , b j } in G where a i ∈ V A and b j ∈ V B . The edges {a 2 , b 1 } and {a n , b n-1 } also exist.
In the graph, the vertices a n-1 , a n ∈ V A and b 1 , b 2 ∈ V B have degree equal to two. Assume without loss of generality that we apply Rule-2 on b 1 , to merge a 1 and a 2 . Then in the reduced graph, a 2 a 1 will be a degree-1 vertex hence Rule-1 can now be applied. We continue with Rule-1 until four vertices remain; applying Rule-2 followed by Rule-1 yields a perfect matching. Hence KS with Rule-2 always find a perfect matching for these synthetic graphs. Furthermore, since Rule-2 is applied only twice, KS requires linear time. If however, only Rule-1 reductions are allowed, no reduction is possible at first and KS R1 will immediately resort to random decisions. In Table 1, we provide the average quality returned by KS R1 for different values of n as well as run times. As can be seen, there is 25% quality difference in favor of KS. Interestingly, KS R1 requires slightly more time than KS. This can be explained by the fact that due to the random decisions, during most steps KS R1 needs to iterate over two adjacency lists, while KS needs to iterate over only one, due to it applying Rule 1. Furthermore, the cost of both exact algorithms increases with n. The average performance on the family of graphs referenced in Section 4.2.1 for n ∈ {7500, 10000, 15000, 20000}. For these graphs, KS finds a perfect matching. The PF and PR timings do not include the time needed for KSR1.

n KS R1 KS 10, 000 0 0.68 25, 000 0 0.64 50, 000 0 0.56 graph. Then each vertex of the graph randomly picks two other vertices with uniform probability from the the other side. A 2-out graph therefore has approximately 4n edges (in case of parallel edges, only one is kept). Walkup [START_REF] Walkup | Matchings in random regular bipartite digraphs[END_REF] showed that it is almost certain that such graphs contain a perfect matching. By construction, each vertex in a 2-out graph has degree at least two and KS R1 will resort to random decisions immediately. This is not the case however for KS. In expectation, a random 2-out graph will initially have about 2n e 2 degree-2 vertices. Hence Rule-2 can potentially be applied a significant amount of times before the graph runs out of degree-2 vertices. Consequently, because it has a smaller kernel and can apply its rules easier, KS will have to make random decisions less frequently than KS R1 and therefore has less of a chance to make an erroneous decision.

In Table 2, we compare KS R1 and KS on varying sizes of 2-out graphs and show the percentage of tests where the two algorithms found a perfect matching. For each n ∈ {10000, 25000, 50000}, we generate 20 different random 2-out graphs. Each distinct random 2-out graph is also used four times to allow alternative random decisions on different runs. As the table shows, KS R1 is never able to find a perfect matching on such graphs. On the other hand, KS is able to find one in a significant percentage of the trials. Nonetheless, both algorithms output matchings of high quality. i.e., greater than 0.99. In addition, both KS and KS R1 require less than a second in all of the instances.

Experiments with worst-case inputs:

In the results with the real graphs for Section 4. KS cache in the worst-case instance shown in Figure 1 for n ∈ {40000, 80000, 160000, 320000}.

observed that KS behaves in practice similar to a linear time algorithm. For this reason, the various modifications discussed in Section 3 in most cases did not have a noticeable positive impact on the performance. However, we do note that our modifications can indeed be useful and we demonstrate in Table 3 the performance for various values of n for the instance in Figure 1. For simplicity, we provide only results with KS cache but an implementation using a hash table, and KS comp both had equivalent performance. As can be seen, while KS R1 and KS cache require in all instances less than a second, KS struggles as n increases. For example, it requires over four minutes in the last case. This shows that in certain situations the sub-quadratic algorithms of Subsection 3.1 might be crucial due to their reduced worst-case behavior. Similarly the other algorithms in Section 3 might prove useful as they can avoid some pitfalls that harm the performance of KS.

Experiments with another implementation:

Korenwein et al. [START_REF] Korenwein | Data reduction for maximum matching on real-world graphs: theory and experiments[END_REF] provide an implementation in C++ which finds the kernel G (k) of a given undirected graph G. We use KS TUB to denote this code. The code KS TUB works along the same lines as our kernelization code. Initially, Rule-1 or Rule-2 are applied for as long as it is possible, and the kernel G (k) is created. Then, the vertices in G (k) are renumbered to be from 1 until |V G (k) |. This smaller representation of the kernel is then given as input for the exact matching algorithm. The code KS TUB was compiled according to the instructions.

We now list some of the differences between our implementation and KS TUB before giving experimental results. One important difference between our implementation and KS TUB is that we prioritize Rule-1 over Rule-2. In KS TUB , vertices of degree one and two are stored together in a container, and the appropriate reduction rule is executed when a vertex is accessed from this container, depending on its degree. As a consequence, KS TUB might apply Rule-2 more frequently, which can harm the performance. Indeed, consider an extension of the graph shown in Figure 1, in which the two leftmost vertices are connected with an additional degree-1 vertex each. Then, by applying Rule-1 on one of the new vertices, it becomes possible to find the kernel in linear time. If, however, the application of Rule-1 is postponed in favor of Rule-2, then the worst case run time becomes quadratic. The code KS TUB also uses linked lists which are not as cache-friendly as our vector-based implementation. Furthermore, our implementation also returns a matching for G, whereas KS TUB returns only the kernel and does not return any information about which edges or vertices of G have been involved in the Rule-2 reductions. Hence, it is not possible to extend the matching found in the kernel G (k) to a maximum cardinality matching for G using KS TUB as is.

We now give the run time comparison of the two codes. The experiments were conducted on Arch 1. For the sake of controlled experimentation, we made the two codes have the same permutation of the input graphs. KS TUB normally prints output; we removed this part for fairness. In Figure 6, we compare the two codes on the bipartite graphs used in Section 4.1.1. For each implementation, we report the minimum time observed over four runs. In both implementations, we measured the time to apply Rule-1 and Rule-2 and the time needed to renumber the vertices in G (k) . In addition, for our implementation, we included additional costs such as the initial cost required to create the adjacency list representation of G and as well as the time needed by the algorithm of Section 3.4 to find a maximum cardinality matching of the original graph (by using and expanding the one on the kernel to comply with the Rule-2 applications). As can be seen, our implementation is almost in all cases faster than KS TUB and in some situations significantly so. In regards to the time required solely for the purpose of kernelization, shown in Figure 6(a), we see that our implementation is about 3.19 times faster on average than KS TUB . The maximum time our code requires is 5.8 seconds, whereas KS TUB exceeds 10 seconds on six instances. The difference in speed becomes more noticeable in instances such as kron g500-logn21, where our implementation takes 3 seconds, while KS TUB takes 34 seconds.

The difference needed to renumber the vertices in the kernel can similarly be as striking as can be observed in Figure 6(b). For example, on the bipartite graph Com-Orkut, our code to renumber vertices takes 5.6 seconds, whereas KS TUB requires 41 seconds. We estimate that this significant difference derives from the fact that KS TUB creates an entirely new class object, for which it then has to initialize a number of linked lists (one per vertex). In contrast, our code needs to initialize only four arrays as G (k) is recreated in the standard CSR/CSC formats of sparse matrices to be given as input to the exact algorithms.

Finally, looking on the overall time, again in Fig- [START_REF] Korenwein | Data reduction for maximum matching on real-world graphs: theory and experiments[END_REF] on the graphs of Section 4.1.1. Figure 6(a) gives only the time needed to find the kernel by applying Rule-1 and Rule-2. Figure 6(b) gives the overall time by adding the time needed to renumber the vertices of the kernel G (k) . In addition, for our implementation Figure 6(b) also includes the time required to recover the maximum matching using the algorithm of Section 3.4, and the time to create the adjacency list representation for G. ure 6(b), our implementation is about 2.17 times faster than KS TUB even with recovery of a maximum cardinality matching for the initial graph. Our implementation never exceeds 15 seconds in overall time, whereas there exist six graphs where KS TUB exceeds this time bound.

Related work

Early work by Möhring and Müller-Hannemann [START_REF] Möhring | Cardinality matching: Heuristic search for augmenting paths[END_REF] popularized the use of Karp-Sipser as the name for KS R1 , and later work [START_REF] Duff | Design, implementation, and analysis of maximum transversal algorithms[END_REF][START_REF] Kaya | Pushrelabel based algorithms for the maximum transversal problem[END_REF][START_REF] Langguth | Heuristic initialization for bipartite matching problems[END_REF] followed suit. To the best of our knowledge, there are only a few studies focusing on implementing KS as in its original, proposed form. Magun [START_REF] Magun | Greeding matching algorithms, an experimental study[END_REF] studied the cardinalities of matchings produced by KS, KS R1 , and other extensions in random graphs. The paper claims linear running time for all heuristics, but from the pseudocode and description in the paper, it is not clear whether the reductions are implemented in a way that is equivalent to ours. Analysis of the source code however does suggests that the implementation is equivalent and hence can similarly require quadratic time. The paper also does not study practical running time. Mertzios et al. [START_REF] Mertzios | The power of linear-time data reduction for maximum matching[END_REF] apply the first reduction rule and a restricted version of the second rule described before in linear time to obtain a linearsize kernel with respect to a graph parameter (known as feedback edge number) for computing maximum cardinality matchings in undirected graphs. A cubic kernel on bipartite graphs with the same reduction rules is also shown. Korenwein et al. [START_REF] Korenwein | Data reduction for maximum matching on real-world graphs: theory and experiments[END_REF] discuss extensions of Karp-Sipser reduction rules for the maximum weighted matching problem. Their results indicate that kernelization is less successful in the weighted case in comparison to the unweighted case. Bartha and Krezs [START_REF] Bartha | A depth-first algorithm to reduce graphs in linear time[END_REF] discuss a linear time algorithm when only the second rule of Karp-Sipser is considered. However, using this algorithm does not result in a linear time matching heuristic. Furthermore, several theoretical papers [START_REF] Aronson | Maximum matchings in sparse random graphs: Karp-sipser revisited[END_REF][START_REF] Chebolu | Finding a maximum matching in a sparse random graph in O(n) expected time[END_REF] analyze the expected matching quality of the Karp-Sipser algorithm and derived variants for sparse random graphs.

Conclusion

We have investigated two data reduction rules for the maximum cardinality matching problem proposed by Karp and Sipser [START_REF] Karp | Maximum matching in sparse random graphs[END_REF]. While the first rule has a simple, worst case linear time implementation, the second rule can take quadratic time. We have focused on the second rule which merges the two neighbors of a degree-2 vertex. We have considered and analyzed three different algorithms with different levels of sophistication, and proposed an efficient algorithm to recover the matching in the original graph. For two of the these algorithms, we showed that their worst-case performance can still be quadratic in terms of n, whereas the third approach has sub-quadratic complexity in sparse graphs. On a set of experiments with real life, random, and constructed problem instances we showed that the second rule indeed increases the cardinality of the matching found by Karp-Sipser and can lead to a drop in the run time of the exact algorithm afterwards.

One open question is whether a linear time implementation for the second rule is actually possible. As discussed in the text, it seems that the approach of Bartha and Krezs [START_REF] Bartha | A depth-first algorithm to reduce graphs in linear time[END_REF] does not achieve a linear time complexity when degree-2 reductions are interleaved with degree-1 reductions or random choices, made in the matching context. With these, we conjecture that the answer to our question is negative. We are also interested in a parallel version of the Karp-Sipser algorithm.

Figure 2 :

 2 Figure 2: A toy bipartite graph G = (VA ∪ VB, E) where |VA| = |VB| = n and the vertices in VA and VB are colored with black and red, respectively.

Figure 3 :

 3 Figure 3: An example of a component: diamond shaped vertices correspond to the degree-2 vertices and circle shaped vertices correspond to boundary vertices. The dashed lines correspond to additional edges, unrelated to the displayed component.

Figure 4 :

 4 Figure 4: (a) The kernel quality of KS and KSR1 for the 59/93 graphs having at least a single vertex with degree strictly

Figure 5 :

 5 Figure 5: (a) The matching quality of KS and KSR1 for all 93 graphs. The quality is measured as the number of

  Comparison only for kernelization between KS and KSTUB.

  Overall time comparison between KS and KSTUB.

Figure 6 :

 6 Figure6: Comparison of our implementation with KS TUB[START_REF] Korenwein | Data reduction for maximum matching on real-world graphs: theory and experiments[END_REF] on the graphs of Section 4.1.1. Figure6(a)gives only the time needed to find the kernel by applying Rule-1 and Rule-2. Figure6(b) gives the overall time by adding the time needed to renumber the vertices of the kernel G (k) . In addition, for our implementation Figure6(b) also includes the time required to recover the maximum matching using the algorithm of Section 3.4, and the time to create the adjacency list representation for G.

Table 1 :

 1 4.2.2 2-out graphs:The random 2-out graphs are created by initially considering an empty bipartite

			KS R1			KS
	n quality	time	PF	PR quality time
	7,500	0.75	1.01	1.29	1.30	1 0.73
	10,000	0.75	1.95	2.86	2.57	1 1.46
	15,000	0.74	5.61	7.48	6.37	1 4.48
	20,000	0.74 10.34 17.08 12.37	1 8.82

Table 2 :

 2 The ratio of tests for which KSR1 and KS find a perfect matching in random 2-out graphs.

Table 3 :

 3 1, it was The run time in seconds for KSR1, KS, and

	n KS R1	KS KS cache
	40,000	0.00	3.44	0.04
	80,000	0.01	15.32	0.11
	160,000	0.03	64.53	0.22
	320,000	0.10 265.89	0.48