{'fraction_non_alphanumeric': 0.030437782718699687, 'fraction_numerical': 0.03423292187380204, 'mean_word_length': 4.762535895736691, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
ConnectedReads: machine-learning optimized long-range genome analysis workflow for next-generation sequencing

Current human genome sequencing assays in both clinical and research settings primarily utilize short-read sequencing and apply resequencing pipelines to detect genetic variants. However, theses mapping-based data analysis pipelines remains a considerable challenge due to an incomplete reference genome, mapping errors and high sequence divergence. To overcome this challenge, we propose an efficient and effective whole-read assembly workflow with unsupervised graph mining algorithms on an Apache Spark large-scale data processing platform called ConnectedReads. By fully utilizing short-read data information, ConnectedReads is able to generate assembled contigs and then benefit downstream pipelines to provide higher-resolution SV discovery than that provided by other methods, especially in high diversity against reference and N-gap regions of reference. Furthermore, we demonstrate a cost-effective approach by leveraging ConnectedReads to investigate all spectra of genetic changes in population-scale studies.

that provided by other methods, especially in high diversity against reference 1 and N-gap regions of reference. Furthermore, we demonstrate a cost-effective 2 approach by leveraging ConnectedReads to investigate all spectra of genetic 3 changes in population-scale studies. 4 Background 5 Whole-genome sequencing (WGS) is increasingly used in biomedical research, 6 clinical, and personalized medicine applications to identify disease-and drug-7 associated genetic variants in humans, all with the goal of advancing precision 8 medicine (1). At present, next-generation sequencing (NGS, also called short-9 read sequencing (SRS)) is a well-established technology used to generate 10 whole-genome data due to its high throughput and low cost (2). Resequencing, 11 especially of human samples, is one of the popular applications of NGS. This 12 process maps raw reads against a reference genome and determines all kinds of 13 genomic variations, including single nucleotide polymorphisms (SNPs) and 14 indels as well as genetic rearrangements and copy-number variants (CNVs) (3). 15 However, a fundamental flaw in the resequencing pipeline is that it ignores the 16 correlation between sequence reads; thus, resequencing does not fully and properly utilize sequence data and may generate inconsistent alignments, 1 which make variant calling, especially structural variant (SV) calling, more 2 complicated (4, 5). Since the human reference genome is incomplete and 3 contains many low-complexity regions, assembling sequence reads without 4 reference bias would be a proper way to overcome the above challenges. 5 Nonetheless, assembly-based approaches for WGS data suffer from several 6 computational challenges, such as high computing resource requirements and 7 long turnaround times. 8 In this article, we propose an efficient whole-read assembly workflow with 9 unsupervised graph mining algorithms on an Apache Spark large-scale big data 10 processing platform called ConnectedReads. By leveraging the in-memory 11 cluster computing framework of Apache Spark (6), ConnectedReads takes less 12 than 20 hours to process 30-fold human WGS data and generates 13 corresponding long assembled contigs for downstream analysis, such as read 14 mapping, variant calling or phasing. The spec. of cluster and the execution time 15 for each step of ConnectedReads are described in the method session. To 16 evaluate the data correctness of ConnectedReads, we use 68 high-confidence 17 insertions in the NA12878 sample detected by svclassify (7). To demonstrate 1 the influence of reference incompleteness, mapping error and high sequence 2 divergence on NGS data analysis, three public available samples from different 3 populations are used. Through ConnectedReads, we are able to investigate 4 unique non-reference insertions (UNIs) and non-repetitive, non-reference 5 (NRNR) sequences from population datasets (8,9). Furthermore, high 6 resolution for SVs are provided by leveraging ConnectedReads, especially on 7 breakpoints of insertions. In conclusion, ConnectedReads optimizes NGS reads 8 to generate long assembled contigs, not only reducing mapping error but also 9 benefit the downstream data analysis pipelines. 10

11
Data preparation 12 To exemplify the challenges of WGS data analysis mentioned above and assess 13 the utility of ConnectedReads, three WGS datasets from three different ethnic 14 groups were selected from publicly available databases, as listed in Table 1, 15 including NA12878 of European ancestry, NA24694 of Asian ancestry, and 16 NA19240 of African ancestry. In addition, the samples were around 30X and 1 sequenced by three different Illumina platforms, namely, the NovaSeq 6000 2 (NA12878), HiSeq 2500 (NA24694), and HiSeq X Ten (NA19240) platforms. 3 Therefore, we believe that these datasets are representative of the majority of 4 WGS data on cohort datasets and national bio-bank projects. 5 based on paired-end or barcode information, any two contigs would not be 6 assembled together without sufficient supports for paired-end information or 7 overlaps. Although ConnectedReads aims to construct more accurate 8 assembled contigs rather than longer ones without haplotypes, it is usually able 9 to construct several contigs of more than 30 Kbps. In addition, there are 10 1,402,511, 1,224,389 and 2,082,886 contigs of >=1 Kbps for NA12878, 11 NA24694 and NA19240, respectively. Furthermore, the length of contigs is 12 strongly correlated with coverage and read length. The deeper the coverage is, 13 the longer the contigs that can be generated are. 14 7 3 ConnectedReads offers two advantages for downstream analysis. One is 4 mapping recovery, and the other is SV detectability. Both of these advantages 5 are described comprehensively by the following experimental results and case 6 Genomes Project and several studies (15,16), the number of SVs decreases as 1 the length of the SVs increases. Therefore, the majority of the SVs are small 2 indels (<50 bps). Then, the trend of the distribution slightly decreases as the 3 length of the SVs increases, except for the peak at 250-300 bps. The reason is 4 due to abundant Alu elements whose body lengths are approximately 280 bases 5 (17). In addition, several studies have reported this phenomenon when using 6 PacBio LRS (13). Therefore, ConnectedReads is able to complement SRS in not 7 only mapping recovery but also SV detection. 8 Comparison against a high-confidence 9 truth set 10 If ConnectedReads generates a mis-assembled contig, then it should be observed via 11 the CIGAR string after applying any read mapping tools. For example, two unrelated 12 unitigs are assembled together, the contig usually is marked as large soft-or hard-13 clipping in its CIGAR string. In addition, insertions and deletions in the CIGAR string are 14 the candidates should be investigated in the evaluation process. To evaluate the data Since insertions are more difficult to be detected by most of existing SV callers 1 12 (PC), and potentially detectable (PD). The CCC-INS CC+PC caller has a strict 1 criterion because it is based on both completely constructed and partially 2 constructed insertions. The CCC-INS CC+PC+PD caller has a lenient criterion 3 because it accommodates all three levels of detection. As shown in Table 3, CCC-4 INS CC+PC and CCC-INS CC+PC+PD achieve 86.8% and 95.6% detection rates, 5 respectively, indicating that using ConnectedReads is able to achieve the same 6 level of performance as PacBio LRS. Therefore, the above experimental results 7 give us confidence to investigate SVs in population-scale data. 8 Table 3.
Long insertion dataset in NA12878 from svclassify 9

Method
Detection rate (%) Therefore, leveraging ConnectedReads is competent at SV detection, especially 14 for insertions. Then, the next question that we are interested in is how many 15 unique SVs exist in each population. According to Kehr's findings (in Table S4 of (9)), there are 372 SVs in all Icelanders. After removing redundant SVs and 1 merging adjacent SVs, 248 distinct SVs are represented as the second set of SV 2 dataset in this paper. We are eager to know whether these 248 SVs are unique 3 to Icelanders or shared by all populations. Therefore, the three samples from 4 different continents shown in Table 1 are processed by ConnectedReads, and 5 then minimap2 is used to generate alignment records for their ConnectedReads 6 contigs. Then, we manually evaluate whether the 248 SVs exist in the three 7 samples and the results are listed in Table 4. More than 95% of the common 8 SVs in all Icelanders are also found in the three different populations, and only 9 one SV cannot be found in any of the three given samples. It is obvious that 10 most common SVs in Icelanders are not unique to Icelanders. Furthermore, 11 approximately 40 SVs should not be classified as SVs in the given samples 12 because they are composed of several adjacent small variants and incorrectly 13 detected as SVs. More details will be discussed in the next section. It is believed 14 that ConnectedReads provides us with higher resolution to observe SVs than 15 some of SV calling tools. Accuracy of insertion length 4 In the above section, ConnectedReads not only was complementary to the 5 resequencing pipeline of SRS in terms of uncovered regions but also was able 6 to benefit SV discovery, especially for long insertions. ConnectedReads and 7 PacBio adopt HG38, but svclassify uses HG19. For comparison, the insertions 8 provided by Spiral Genetics are converted from HG19 to HG38 by using UCSC 9 LiftOver. Then, all of SVs from different source or methods are comparable 10 since they are based on the same genome coordinate. As shown in Table 3, at 11 least 85% of long insertions can be detected by using ConnectedReads followed 12 by minimap2 and CCC-INS, and 32 insertions are completely constructed. As 13 shown in Fig. 4, most of the insertions identified by ConnectedReads are almost the same as those identified by Spiral Genetics and PacBio. 24 of 32 insertions 1 identified from ConnectedReads have identical lengths with PacBio. In addition, 2 the remaining insertions have only slight differences. Therefore, 3 ConnectedReads with SRS can construct complete and accurate insertion 4 sequences as well as PacBio LRS can. 5

Granularity of SV detection 6
Another interesting phenomenon shown in Table 4  Furthermore, ConnectedReads is able to mitigate the impact on downstream 3 analysis due to the incomplete human reference genome. The most recent 4 human genome still has many ambiguous areas (N-gaps), and they are mainly 5 located in centromeres and telomeres. Fig. 6 illustrates that two ambiguous 6 gaps can be assembled by using ConnectedReads and that the sequences are 7 totally identical among the three individuals from different populations. This 8 finding gives us strong confidence that most humans might have the same 9 sequence in the two N-gap regions. By randomly selecting two Chinese adults, 10 the occurrence of identical sequences in the N-gap regions is confirmed by 11 Sanger sequencing, as performed by a Clinical Laboratory Improvement 12 Amendments (CLIA)-certificated laboratory. The length of the ambiguous 13 region in Fig. 6(a) should be corrected from 382 to 223 bps. In addition, the 14 length of ambiguous regions in Fig. 6(b) should be shortened from 45 to 31. 15 is able to provide a cost-effective approach with which to complete the human 1 reference genome. 2 Correctness of reconstructed contigs 3 The next topic for the evaluation of ConnectedReads is correctness. To 4 comprehensively investigate the data correctness of contigs constructed by 5 ConnectedReads, all of imperfect matched sequences should be considered 6 because any wrong constructed contigs should not be perfectly aligned to 7 reference genome. Therefore, insertion sequences, soft-clipped sequences and 8 unmapped contigs are selected and then identified as being from Homo sapiens 9 or just chimeric DNA sequences resulting from false reconstruction. Take 10 NA12878 for example, there are 37 insertions, 384 soft-clipped and 411 11 unmapped contigs by removing sequences with a length < 1,000 bps. By using 12 BLAST to find any homologous sequences in the National Center for 13 Biotechnology Information (NCBI) non-redundant sequence (nr) database, 14 each sequence can be identified as human or not. As shown in Fig. 7, 35.1%, 15 46.6% and 37.7% of the insertions, soft-clipped sequences and unmapped 16 contigs are identical to Homo sapiens DNA sequences in the nr database. As the threshold of similarity is continuously lowered, more evidence can be found 1 to support the sequence assembly correctness of ConnectedReads. However, 2 there are six sequences without any support when the threshold is set to 20%. 3 Two of the contigs have low-complexity content, and two are matched to several 4 entries but with less than 10% support. The last two unmapped contigs 5 (CONTIG-8337086 and CONTIG-15793805) are more than 10 Kbps in length. 6 CONTIG-15793805 covers all of CONTIG-8337086 in reverse-complement 7 mode, so CONTIG-15793805 is represented for evaluation. In Fig. 8, CONTIG-8 15793805 is almost fully covered by several long Homo sapiens sequences, and 9 some of these sequences are highly overlapped (at the position of 3000-3500 10 and 7000-7700). Therefore, the results give us more confident that these non-11 reference contigs constructed by ConnectedReads should be from Homo 12 sapiens. Although more effort on verifying data correctness is emerging and 13 expected, we believe that ConnectedReads achieves high data correctness based 14 on the above findings. 15 1 Another way to evaluate the data correctness of ConnectedReads is to check 2 whether there is any translocation. If any two unrelated sequences are 3 incorrectly assembled together, it will cause a fake translocation event, in which 4 a contig is mapped to multiple chromosomes. Table 5 lists the 13,686 contigs 5 with multiple alignment records on two chromosomes in sample NA12878. 6 Table 5.

Item Number
Number of contigs mapped on at least two chromosomes 13,686 Number of contigs after removing singletons 2,075 Number of qualified contigs 71 Number of high-confidence translocation groups 22 Number of translocation-based insertions 8 Generally speaking, a translocation-based insertion should have at least three 8 breakpoints (two for translocation sequences and one for insertion) on genome 9 reference, so it's difficult to find a single contig with three breakpoints. If a 10 breakpoint for insertion is mapped in only one contig (named as singleton), it 11 won't be a candidate for translocation-based insertions. After removing 12 singletons and filtering out the contigs with low mapping quality (MAPQ < 60), 13 71 qualitied contigs are represented by 22 translocation groups. 14 Interestingly, eight translocation groups have one clear breakpoint on one 1 chromosome but two breakpoints on another chromosome. Those 2 translocation-based insertions can be observed from the VCF of NA12878 3 generated by pbsv. As shown in Fig. 9(a) bp sequence in intron 1 of CLCC1 is inserted at chr3:110694547. Fig. 9(b) also 8 shows that the 697-bp sequence in exon 10 of BTBD7 is inserted into intron 1 of 9 SLC2A5. These translocation groups are called translocation-based insertions. 10 ConnectedReads proposes a naïve way to investigate these translocation-based 11 insertions. In summary, ConnectedReads not only constructs genome 12 sequences precisely but also facilitates SV detection by mitigating reference 13 mapping bias. 14 Conclusions ConnectedReads, SRS is able to be constructed from short reads to long 1 assembled contigs. Using ConnectedReads not only prevents mapping errors 2 but also facilitates SV discovery. In addition, ConnectedReads are able to 3 mitigate the influence on reference incompleteness. In summary, 4 ConnectedReads can serve as an NGS gateway for facilitating downstream data 5 analysis, such as false positive prevention, SV detection, and haplotype 6 identification. 7 Materials and methods 8 Workflow 9 ConnectedReads leverages Apache Spark (6), a distributed in-memory 10 computing framework, to accelerate its whole workflow, as illustrated in Fig. 10. 11 To fully utilize the power of the distributed framework, the preprocessing step 12 involves splitting a large compressed file into several small files and then 13 uploading those files into the Hadoop distributed file system (HDFS) since 14 most of the WGS samples exist in two separate FASTQ files in GZIP format. 15 First, we adopt Apache Adam (19), as shown in Fig. 10(a), to transform data 16 from FASTQ format to column-based Parquet format for data parallel access. 1 To facilitate data processing in the following steps, we extend Adam to not only 2 encode paired-end information and barcodes into the read name column but 3 also place all reads in different subfolders based on their sequencing quality and 4 sequence complexity. Then, we propose a distributed suffix tree algorithm with 5 supervised graph mining on Spark to mitigate the influence of improper string 6 graph construction due to sequencing errors. Using an outlier detection 7 algorithm on a suffix tree, the process in Fig. 10(b) can be configured as a highly 8 sensitive error detector for low-coverage regions. After that, we are able to 9 adopt the parallel string graph construction shown in Fig. 10(c) to represent the 10 relation of each qualified read and the read overlaps by suffix-prefix 11 information. More details are available in our previous paper (20). Based on 12 the string graph, we propose the parallel haplotype-sensitive assembly (HSA) 13 depicted in Fig. 10(d) to construct assembled contigs; the detailed procedure of 14 this module will be described in the next session. For example, these generated 15 contigs (Fig. 10)

include (i) heterozygous SNPs, (ii) heterozygous deletions, (iii)
ConnectedReads transforms noisy short reads with low quality and sequencing 1 errors to long assembled and qualified contigs. 2 Parallel HSA 3 ConnectedReads adopts a string graph, a lossless data representation that is 4 fundamental for many de novo assemblers based on the overlap-layout-5 consensus paradigm (21,22), to represent the overlaps of each read. Here, we 6 propose a Spark-based HSA based on a string graph. By leveraging GraphFrame, 7 HSA is able to perform efficient and scalable graph traversal operations and 8 supervised graph mining algorithms on Apache Spark. Before introducing the 9 detailed algorithms of HSA, we formally define the notations of the sequencing 10 data and string graph that we will use in the following sections. 11

12
Let G(V, E) be a directed graph and S, T Í V. We define V = {v1,v2, …,vk} and 13 E(S, T) to be the set of edges going from S to T, i.e., 14 Definition 1. Let R be the set of short reads and R RC be the reverse complement 1 set of R. G and G RC are the string graphs based on R and R RC , respectively. 2 If eij exists in G and i, j Î R, then emn also exists in G RC and m,n Î R RC such 3 that m and n are the reverse complements of j and i, respectively. 4 We define G EXT (V, E) to be the expanded graph of vertices V Î R È R RC and 5 edges E to be the set of edges in either G or G RC . 6

Definition 2.
Let IN(vi) and OUT(vi) be the number of in-degree and out-7 degree edges of vi, respectively. We define a vertex vi to be

Definition 3.
We define the priority order of the vertex labels to be T > S > 14

B and A. 15
This means that B could be relabeled as S or T. Once a vertex becomes T, it will 16 always be T regardless of what the graph property propagation is.
To keep the depth information in FASTQ format, a new quality encoding 1 function is proposed. In addition, we define the quality-depth (QD)-encoding function 10 System flow 13 In Fig. 11, HSA makes use of five modules, namely, (a) string graph expansion, traversal algorithm, (d) the ploidy-based routing algorithm and (e) the triplet 1 traversal algorithm. To simplify the following graph traversal algorithms, each 2 vertex and its reverse complement should be separated in a string graph. This 3 means that the string graph G generated by Fig. 10(c) should first be expanded 4 by Definition 1 and named G EXT , which contains G and G RC . For that reason, Fig.  5 11(a) generates the expanded string graph with all reads and their reverse 6 complements. After removing all singletons in G EXT , the remaining vertices will 7 be classified by the vertex connectivity classification module shown in Fig. 11 into four classes, which are defined in Definition 2: start (S), bridge (B), 9 terminator (T) and ambiguity (A). To generate assembled contigs, the graph 10 traversal/pairing algorithms shown in Fig. 11(c-e) are proposed based on the 11 above properties of vertices and described comprehensively in the next section. 12 Graph mining algorithms 13 ConnectedReads leverages Apache Spark and its derived packages to propose 14 three efficient and scalable graph algorithms for large graphical datasets, i.e., 15 string graphs. Spark GraphFrame is a powerful tool for performing distributed 16 computations with large graphical datasets. In addition, Spark Dataset is a type-safe interface that provides the benefits of resilient distributed datasets 1 (RDDs) and Spark SQL optimization. By leveraging GraphFrame and Dataset, 2 we propose several Spark-based graph traversal/routing algorithms for HSA 3 with high performance and scalability. 4 From our observations of NGS sequencing data, two challenges must be 5 overcome if we want to enhance the performance of graph operation for HSA. 6 The first challenge in the string graph involves long diameters, and the second 7 challenge is how to properly connect the vertices with multiple in-/out-degrees. 8 Taking NA12878 as an example, more than 90% of the vertices are B, and the 9 longest diameter from S to T is more than 500. This means that the traversal 10 algorithm must perform the propagation operation in at least 500 iterations 11 from one vertex (S) to another vertex (T). For most graph frameworks, the 12 performance of a graph algorithm is strongly related to the number of iterations 13 for its traversal operation. Here, we propose the S-B-T traversal algorithm, 14 shown in Fig. 11(c), to connect all vertices from S to T via all adjacent Bs. By 15 following Definition 3, the data preprocessing flow for the expanded string vertex A should be relabeled as T and S, respectively. To overcome the long-1 diameter problem, the mechanism used to relabel B to T by customized random 2 selection is applied to shorten the diameter of the given graph, and then, the 3 terminator-intensive propagated string graph is acquired. Based on the 4 terminator-intensive graph, the S-B-T traversal algorithm, which integrates a 5 belief propagation algorithm with iterative graph merging, is able to 6 theoretically reduce the time complexity from O(N) to O(√ ) (N: the number 7 of iterations for graph propagation). 8 When each simple routing path (i.e., S to T via Bs) is completely traversed and 9 merged into a single vertex in the new graph, the majority of the vertices have 10 multiple in-/out-degrees. Therefore, we have to solve the second challenge -11 how to select a proper routing path in the multiple-in-/out-degree graph. Many

13
To mitigate false assembly due to sequencing errors, the function de_noise() 14 is used to remove the noisy vertices in Iv by using a naïve clustering algorithm 15 based on the read-pair information in this paper. In addition, the function 16 find_best_matching() adopts a tripartite clustering algorithm based 17 primarily on the number of supports from read-pair information in (u, v), (u,18 Ov) and (v,Ov) to find the best combination for sequence assembly. Using the 19 ploidy-based pairing algorithm, we obtain sufficient information to overcome 20 the second challenge and then apply the triplet traversal algorithm shown in 21 Here, sample NA12878 is processed by ConnectedReads on our Spark cluster. 5 It takes approximately 18 hours for the 30X WGS sample; the detailed 6 performance of ConnectedReads is described in Table 6. 7

QD encoder/decoder
Since ConnectedReads not only connects reads with their overlaps but also 11 aggregates identical reads, the size of ConnectedReads contigs will be reduced 12 by at least 90% in comparison with that of contigs from short-read data.
However, the downstream analysis tools might not work well due to the loss of 1 depth information. We have to take the depth information into consideration 2 and keep the output compatible with FASTQ format. Therefore, the QD encoder 3 and decoder are proposed based on Definition 4 to mitigate the impact on 4 information loss. The mechanism of the QD-encoding function is quite flexible, 5 allowing it to fit most use cases. For example, if the depth information is more 6 critical than quality, (L, N) = (3, 14) is the best choice. Fig. 13 provides an 7 example to demonstrate how the QD encoder transforms short reads into 8 ConnectedReads contigs. The data reduction rate is approximately 77% (from 9 144 characters to 33 characters in both sequence and quality). In addition, Fig.  10 14 shows the result of applying the QD-decoding function. The recovery rate in 11 terms of sequences and quality is 95.8% (138/144) and 47.3% (71/144), 12 respectively. Therefore, ConnectedReads provides an efficient QD-decoding 13 tool for some use cases that heavily leverages depth information. 14 Availability of data and materials  counting is insertion > deletion > soft-/hard-clipped sequence. Any two 3 SVs will be merged if their distance is less than 50 bps. This means that if 4 one insertion and one deletion are encountered in the same location, only 5 one insertion is counted.

Title Text
Quality+Depth Decoder

Sequence
Quality