IterCluster: a barcode clustering algorithm for long fragment read analysis

IterCluster algorithm

Given a set of long-fragment reads, IterCluster selects appropriate barcodes as the seed and iteratively cluster target barcodes overlapping with the seed barcode. In order to minimize the overlap between clusters in the final clustering result, the target barcode that has become the seed barcode cannot be used as seed in future clusters. The size of each cluster depends on the number of iterations (parameter -k, default 3). The IterCluster algorithm proceeds in three steps: (1) constructing the adjacency matrix between barcodes; (2) generating barcode clusters; (3) extracting the reads of the target regions from each barcode cluster.

In the first step, we first extract the unique kmer in the genomic data set. The distance between the barcodes is measured by the number of common unique k-mers between the barcodes. In terms of time efficiency, IterCluster uses the sparsity of the overlap relationship between barcodes to generate the adjacency matrix.

In the second step, IterCluster evaluates the quality of the barcode by the unique kmer number and unique k-mers proportion, than adds the high quality barcode to the candidate seed list. IterCluster extracts the seed from the candidate seed list and clusters the seed barcode according to the adjacency matrix and the k-mers frequency diversity selection model to obtain the preliminary barcode cluster. Because most barcode carry fragments of multiple regions, the preliminary barcode cluster will also carry barcodes from multiple regions. IterCluster uses the Markov clustering model to disassemble each preliminary barcode cluster and obtain a single barcode cluster.

Finally, IterCluster will count the kmer coverage depth of the reads in each barcode cluster, identify the low coverage depth reads and filter out, and then get the reads in the barcode clustering area.

Using unique k-mer as a similarity measurement

In the stLFR data, if there are overlaps between a long fragment contained in two barcodes, the read overlap in the two barcodes can be found, as these signify the similarity between barcodes. In measuring the similarity between barcodes, we did not use time-consuming methods such as read alignment between barcodes. Instead, we built a k-mer set for each barcode and measured the relation of barcodes based on the similarity of its k-mer set. Due to the fact that repeats are widely distributed in the genome, the k-mer set generated from a repeat would introduce a false-positive barcode. Therefore, we only select unique k-mers as the feature of a barcode based on the k-mer frequency distribution. Let K (b) denote the unique k-mer set of barcode, the distance between barcodes can be described as: ${\displaystyle \mathrm{Dist}\left(b_1,b_2\right)=K\left(b_1\right)\cap K\left(b_2\right).}$

A time-efficient way to generate adjacency matrix

IterCluster implements barcode data clustering based on an adjacency matrix between barcodes. Many solutions to calculate the collective similarity exist, such as a local sensitive hash (Charikar, 2002), which has been applied to the de novo assembly of single-molecule sequencing data (Berlin et al., 2015) and the construction of phylogenetic trees (species trees) (Ondov et al., 2016). Although these techniques can efficiently detect the overlap between barcodes, a probabilistic method inevitably produces a certain error rate, and an absolute common unique k-mer number between barcodes cannot be calculated. To calculate the common unique k-mer number between barcodes, a simple, time-consuming method is to compute the intersection independently between every two barcodes and the computational complexity of this method is given as: O(n²).

IterCluster, however, utilizes an innovative algorithm to calculate the absolute common unique k-mer number between barcodes. The core logic behind this approach was that the adjacency matrix is sparse and each unique k-mer can determine which barcodes are related to each other. IterCluster extracts the unique k-mer of a barcode and builds a hash from the unique k-mer from each barcode. Barcodes that contain the same unique k-mer mean that the number of co-unique k-mers between them increases by one. Since there are few barcodes that simultaneously contain the same unique k-mer (theoretically less than or equal to the sequencing depth), it takes a constant time to accumulate the number of co-unique k-mers between the same unique k-mer numbers of supported barcodes. By traversing all of the unique k-mers, the absolute value of the co-unique k-mers between all of the barcodes can be calculated, and the time complexity of the program is given as: O (n) (Fig. 2).

Barcode enrichment based on a k-mer frequency diversity selection model

We developed an algorithm to iteratively cluster the barcodes from a target area in the genome. The reasoning behind this approach was that the barcodes covering a target area tend to contain unique k-mers derived from the target area while the content of unique k-mers derived from other areas should be random (Fig. 1C), since these barcode cover multiple long-fragment molecule. After each iteration of clustering, the unique k-mer frequency of target and non-target areas are significantly different. Using this difference to choose the unique k-mers belonging to a target area as a feature of the next round of clustering effectively controls the false-positive rate. IterCluster uses common unique k-mers to measure the similarity between barcodes. The number of common unique k-mers greater than the parameter -g between barcodes is considered to contain the overlapping long fragment molecule. In each iterative clustering process, IterCluster first captures the target barcodes directly overlapping the seed barcode and obtains the target barcode set T. Then it statistically determines the frequency of all of the unique k-mers of the barcode in T. The high-frequency unique k-mer set is then selected as the seed unique k-mer set for the next iteration and added to the feature k-mer set F (Table 1).

MCL model

Since each barcode covers multiple long fragment molecules (Fig. 1A), using a seed clustering strategy may result in a barcode set representing more than one target area in the genome. IterCluster uses a Markov clustering model (Van Dongen, 2000) to ensure that the result of seed clustering is a set of barcodes for individual areas (Fig. 1B). Specifically, assuming that a seed contains two long fragment molecules from different areas, the target barcode set will contain two clusters from different areas in the genome. However, the barcodes generated from the same area are closely related and barcodes generated from different area are distantly related. Therefore, IterCluster use a Markov clustering model to detect and classify closely-coupled barcode groups after the first clustering iteration (Table 2).

MCL algorithm has two main processes, expansion and inflation, which operate on the state transition matrix. When a state transition matrix is M, the dimension of M is the number of points in the graph. Each list in M represents the probability of starting from a certain point at a certain time and arriving at the remaining points at the next time.

The extended process is to simulate the random walk process. Taking positive integer e and multiply the current state transition matrix by e times to get a new state transition matrix, which is equivalent to a random walk on the original state transition matrix. This process can be described as: ${\displaystyle M^{\prime }=M^e.}$

The expansion process is a matrix regularization process, which regularizes the columns of the state transition matrix. The processing formula is shown as follow: ${\displaystyle {\left(M^{\ast }\right)}_{pq}=\frac{{\left(M_{pq}\right)}^T}{\sum _{i=1}^K{\left(M_{iq}\right)}^T}.}$

Where M is a state transition matrix; M* is a regularized matrix; T is a relaxation coefficient; K is the number of rows of M; P is a row subscript; q is a column subscript. The function above is to regularize the columns of the transfer matrix to get the regularized matrix M*

For the N target barcodes obtained after the seed barcode’s first clustering iteration, the N*N adjacency matrix M was built and Markov clustering occurred. Finally, the sub-cluster can be found after MCL process.

Data sources

We select 2 human genome data sets to evaluate IterCluster’s performance: BGI stLFR data sets and 10× Genomics Chromium linked reads data sets (Table 3). BGI stLFR data were obtained for sample NA12878 from the CNGB Nucleotide Sequence Archive (CNSA) with accession ID: CNS0007594. The data were also available from the European Nucleotide Archive with accession ID: PRJEB27414.

Human Chromium linked reads data from 10× Genomics was obtained for a HGP sample, and were downloaded from the 10× Genomics company website (https://support.10Xgenomics.com/de-novo-assembly/datasets/2.1.0/hgp).

Data analysis

Evaluation

In our IterCluster algorithm evaluation, we build actual positive dataset by aligning reads to genome reference with bwa (Li & Durbin, 2009) to know each barcode’s position on the genomes. According to the range of LFR length and reads count in each barcode (Fig. 3), there are so many barcodes have little reads that haven’t enough unique kmer to find other barcodes together. we defined initial seed barcodes that must cover areas where the reads spanned larger than 500 bp in length and had more than 30 paired-ends to represent an LFR fragment and the target barcodes that belong to an LFR fragment need to have more than 5 paired-ends, which can ensure the aggregation ability of the barcode.

After getting the cluster set, compare the IterCluster results with actual positive dataset, count the true positive barcode number that included by cluster set and actual target region barcode set, false positive barcode number that included by cluster set but actual target region barcode set in every random initial seeds’ cluster set target region. And the span region for every initial seed may extend during IterCluster, we decided extend 80 kp∼400 kbp both left and right for the initial seed span region on reference as target region to evaluate. Using a P-R curve to evaluate, determining the Precision and Recall as follows: ${\displaystyle Precision=\frac{TruePositive}{TruePositive+FalsePositive}.}$ ${\displaystyle Recall=\frac{TruePositive}{TruePositive+FalseNegative}.}$

IterCluster for de novo assembly

In the IterCluster result, there were enough reads data in an LFR cluster set to support de novo assembly. To improve the effect of the assembly, cluster false positive reads were removed using the FalseRemove module according to an average k-mer frequency for each LFR, as we assumed that the non-target LFR k-mer frequency was lower than a targets LFR. Then, we used SOAPdenovov2 (Luo et al., 2012) to assemble each IterCluster result set. For each subassembly result, we filtered contigs less than 1 kb because these were most likely sequences assembled by false positive reads. After that, we merged each assembly result together, and evaluated the results with QUAST (Gurevich et al., 2013).

IterCluster performance

To evaluate the cluster size after IterCluster in two dataset. We show the barcodes number distribution (Fig. 4) and reads number distribution (Fig. 5) of BGI stLFR datasets and 10× Chromium reads datasets.

To evaluate the performance of IterCluster on BGI stLFR datasets and 10× Chromium reads datasets, we plotted R-P to analyze the accuracy and recall of IterCluster. As can be seen from Figs. 6A, 6B IterCluster’s accuracy and recall rate were high compared to baseline (k = 1). IterCluster performed better on BGI stLFR data than on 10× Chromium reads data. Figures 6C and 6E shows the performance of different parameters using the BGI stLFR data, where the f parameter represents the lowest frequency threshold in the frequency selection model and the g parameter represents the unique k-mer number threshold between barcodes with overlap. Increasing the -f and -g values did not significantly improve the accuracy of IterCluster, but it greatly reduced the recall rate. However, from Fig. 6F, we found that a slight increase in the -g value can improve the accuracy and recall rate using 10× Genomics Chromium reads data.

An important factor affecting the precision of IterCluster was the repeat sequences in genome. That is why we choose unique kmers as features. The K value (the base number of each kmer) is an important parameter because it determines the specificity of a kmer. We thus studied IterCluster’s performance with different K values using BGI’s stLFR data. Figure 7 shows that IterCluster can achieve higher precision and recall with larger K values.

k-mer frequency diversity selection model improve clusting

IterCluster needs to deal with two key issues in clusting: (1) how to determine the overlap between barcodes; (2) how to avoid false positives overlap caused by multiple fragments contained in the same barcode. These two issues are described in detail below.

In regard to the first issue, because overlap between different barcodes means that these barcodes will contain common unique k-mer, we use the number of common unique k-mer to measure the distance between barcodes. In the IterCluster, a user-specified common unique k-mer threshold (-g parameter) is used to determine whether two barcodes contain overlap. This threshold must be high enough that repeats and sequencing error do not result in false overlaps yet low enough that slight overlaps overlaps between fragments can be detected. This threshold can be estimated by the the distance distribution between all barcodes (Fig. 8A).

Another issue is caused by multiple fragment contained in barcode. Most of target barcodes carry a fragment from seed barcode area but carry fragments from other areas of the genome randomly. Therefore, after each cluster iteration, the feature of the seed barcode area (unique k-mer) need to be updated to prevent non-target area’s unique k-mer from being used as feature for the next iteration of clusting. When choosing feature unique k-mer, IterCluster uses a unique k-mer selection model based on k-mer frequency diversity. The key idea is that after each iteration of the cluster, the unique k-mer frequency from seed barcode area is higher than the unique k-mer frequency from non-seed barcode area, because target barcode always has an overlap with seed barcode. Base on frequency diversity, IterCluster provides a user-specified minimum frequency threshold (-f parameter) to filter the unique k-mer from non-seed barcode area. This threshold can be estimated from the unique k-mer frequency distribution after each clusting iteration (Fig. 8B).

To evaluate the effect of the frequency diversity selection model in controlling false positives, we compared the performance of IterCluster with or without this model (Fig. 9). When iteration k = 2, the number of unique k-mers in non-seed fragment area is still very small, and the frequency selection model has slight effect on IterCluster. However, when the number of iteration increased (k = 3), the false positive rate of IterCluster without frequency diversity selection model is significantly increased. This shows that the frequency selection model can control false positive rate with the number of iteration increased.

Runtime and memory cost

IterCluster’s runtime and memory cost largely depends on the parameter K, the kmer size used to calculate the adjacency matrix. The runtime and RAM usage for IterCluster with different K values is shown in Table 4. The time and memory cost of IterCluster mainly stem from the adjacency matrix construction step of this algorithm. The larger the K value, the more kmer numbers, and the more time is needed to calculate the relation between all kmers. However, changes in K do not greatly affect memory consumption. Although a larger K value means a larger hash is needed to map all of the kmers, for each single kmer, the number of barcodes containing the kmer decreases. Therefore, the consumption of memory is stable and only decided by the amount of data. IterCluster is a highly parallel tool, both in the matrix build step and in its clustering step. We demonstrated that IterCluster can finish a 50-fold human whole genome dataset in 4 days with 50 threads.

IterCluster improves de novo assembly using a divide-and-conquer strategy

We used IterCluster to improve the de novo assembly of human genome stLFR data from BGI via a divide-and-conquer strategy. First, we ran IterCluster on whole human genome stLFR data with default parameters. For each barcode cluster, we filtered the false positive reads with the FalseRemove module which was also added to our pipeline and assembly tasks were performed independently using SOAPdenovo2 (Luo et al., 2012). After the subassembly of each cluster, the assembly results from different clusters were integrated together. We call this assembly strategy “divide-and-conquer” or subassembly. In comparison to a baseline SOAPdenovo2 assembly using all reads (Table 5), assembly with this divide-and-conquer strategy achieved a longer contig length, as measured by the N50 length metric, of 11.9k, representing a 67% increase over the initial baseline achieved using SOAPdenovo2 assembly from all of the reads. Assembly with our divide-and-conquer strategy achieved a 30% increase in scaffold creation (N50 = 17.1K) and a slight improvement in genome coverage. This suggests that IterCluster uses more of the long-range information present in stLFR data to break down a complex chromosomal genome into multiple simple sub-sections.