An Online approach to Overlap Detection for DNA Fragment Assembly

An Online approach to Overlap Detection for DNAFragment AssemblyShantanu Joshi Subramanian ArumugamUniversity of FloridaDepartment of Computer and Information Sciences and EngineeringGainesville,FL,USA{ssjoshi,sa2}@cise.ufl.eduDecember17,2005AbstractWe consider the problem of estimating the number of sequences overlap-ping with a certain length from a large dataset of sequences.Such a prob-lem is central to several challenging bioinformatics applications like DNAfragment assembly.A naive approach of calculating this numberis to per-form an all-pairs comparison over the dataset.We propose an entirely dif-ferent approach to tackle the problem by considering random samples fromthe dataset and performing a ripple-join based comparison between two ran-domized copies of the dataset.One key feature of our approach is that ateach step of the ripple join it produces an estimate of the eventual number offragments having a certain overlap length along with confidence bounds inan online fashion.Since this estimate converges after a sufficient number ofsamples are considered,the algorithm can be terminated much earlier thanthe point where it exhausts processing the entire dataset.We present experi-mental results that indicate that our approach produces sufficiently accurateestimates in a fraction of the time required by a naive approach.1IntroductionOverlap detection in biological sequences is a commonly encountered problem.Anall-pairs comparison between the sequences returns overlap information betweeneach pair of sequences.This information can be processed to compute the fre-quency of each overlap string.Biologists would often be interested infinding outthe number of fragments having overlap length higher than a pre-defined threshold.This information would be very helpful in an application like genome assemblywith shotgun sequencing.In such an application we would typically consider a7to9fold coverage of the genome of interest.Each copy of the genome would thenbe randomly split into several fragments with some average length.The genome assembly problem would then be to reconstruct the original genome by assembling1the various fragments using overlap information among them.Given that the hu-man genome consists about3billion base-pairs and with a typical coverage and an average fragment length of about700base-pairs,a quick calculation shows that we would obtain over38million ing an all-pairs comparison tofind overlaps between such a large number of fragments is extremely computationally expensive and can require several days to complete execution.We formulate overlap computation as a database problem and then propose a novel sampling-based approach to report overlaps among fragment pairs.Consid-ering the collection of fragments to be assembled as a tuples in a database relation R the following SQL query can be used to compute all overlaps of a certain length l.SELECT DISTINCT COUNT(*)FROMFROM R r1,R r2WHERE size(overlap(r1.s,r2.s))=l))We provide a sampling-based solution to overlap computation that is reminescent of an online aggregation system developed by Haas and Hellerstein[HH99].The basic idea is to evaluate the query over a small random sample of tuples from the relation and use that answer to give an estimate about thefinal answer to the query. The process is repeated for a certain number of sampling steps and terminated when the estimate for thefinal answer has converged.In practice,we are interested in not just thefinal answer.We would also like to know at all times how accurate our guess(estimate)of thefinal answer is,since at the point when we are satisfied with the guess,we would terminate our algorithm.The rest of this report is organized as follows.In section2,we provide a brief overview of the Whole Genome Sequencing technique and outline the steps in fragment assembly.In section3,we formulate the fragment assembly problem as the problem offinding the shortest common superstring,highlight the challenges and describe a simple greedy heuristic to solve it.We establish certain properties of overlap computation that make it particularly amenable to a sampling-based approach.In section5,we present an overview of the ripple join algorithm which is central to our online algorithm.Section4describes the database formulation of the problem and provides algorithms for computing overlap using both an online and traditional approach.Section6presents results of our benchmarking experiments. Section7concludes the report with closing remarks.2Whole Genome SequencingThe genome of an organism holds the key to the many biological processes that govern life.Whole Genome Sequencing(WGS)[MS00]is a fundamental problem of determining the genetic sequence of an organism.The common technique for large scale genome sequencing is the shotgun sequencing approach[Ve98].Other2less commonly used techniques include sequencing by hybridization,primer walk-ing and pyrosequencing.In WGS through the shotgun approach,we repeatedly read a short length(<1000bp)of the DNA from the whole genome sequence using basic sequencing techniques(Sanger method,cycle sequencing[SN77]and form a set of DNA fragments.Reads are scheduled in an overlapping manner and the entire DNA is covered many times.Since reads are overlapping it is guaranteed that two consecutive reads will produce a pair of fragments that have a non-zero overlap.Thus,we can piece together overlapping fragments and reconstruct the original whole genome sequence.However,in practice we are not able to control the order of the reads nor are we able to guarantee the overlap conditions necessary to compose the whole DNA sequence.Thus,we rely on techniques that duplicate the DNA sequence many times and produce reads of varying lengths over the set of duplicated DNA sequences.The entire genome sequencing process proceeds in 2distinct steps:1.Short Fragment Sequencing:Take the whole genome and produce a listL of short DNA fragments from the genome.This is entirely a biological process.2.Fragment Assembly:The next step is to assemble the fragments and con-struct the whole genome sequence.This is mainly a computational task.This task can again be divided into subtasks:1.Overlap:Identify all overlapping fragmentsyout:Produce layout l of thefinal genome sequence3.Consensus:If there is more than one layout,use a multiple sequencealignment algorithm to arrive at a consensus genome sequenceIt is important to realize that several things can go wrong with sequence assem-bly.Some DNA sequences contain a huge number of repeats.For example,the human Y chromosome[CT01]has upto90%of repeats.The biological sequenc-ing process(step1)is itself not completely error free.Towards the tail-end of a read,usually there are some errors,which needs to be accounted for when an over-lap is computed.Even if all this works out well,it is a matter of faith that the shortest common superstring is really the representative sequence of the original sequence.For the sake of simplicity,we consider the fragment assembly problem within the framework of a simple model.It is straightforward to incorporate the standard constraints in fragment assembly like error-handling and different criteria for overlap.3SCS Model for Fragment AssemblyWithout loss of generality,we can treat the set S of fragments to be assembled as a a set of strings and reduce the problem of producing a layout into the shortest com-3mon superstring(SCS)problem[KM95].1Given a set S of strings,the objective of the SCS problem is to construct the shortest common superstring that contains each element of S as a substring.Unfortunately,the SCS problem is NP-complete2 and hence heuristic approaches are needed.The greedy approach is the standard heuristic to approximating the shortest common superstring.The idea behind the greedy heuristic is very simple.First enumerate all possible overlap strings.Find a pair of strings with the maximum overlap and replace them with a single merged string.Repeat this process for all remaining overlap strings in a non-increasing order and stop until only one merged superstring remains.The main problem with this heuristic is the time-consuming task offinding the overlap between all pairs of strings.The overlap operation itself is expensive. Given two strings s1and s2of length l,computing the overlap by a brute-force approach requires O(l2)time.And this process needs to be repeated over n×n string pairs.Even for reasonable sizes of n,this is a computationally intensive task.It is not uncommon for n(number of fragments)to be in millions in WGS. Even with the use of parallel-processing and dedicated hardware for acceleration determining the number of unique overlap strings can potentially take several days tofinish.Traditional fragment assembly tools use a number of heuristics to avoid doing the all-pairs comparison needed to compute the overlap strings.In the next section,we propose a novel approach to fast overlap computation that is based on the technique of sampling.Empirical ObservationsGiven a set of n strings,the set O of all possible overlap strings isfinite.Note that|O|≤n(n−1)/2(though in practice|O|≪n(n−1)/2).Observe that for fragment assembly n is large because the original genome may contain a lot of duplicates.The duplicates are needed so that there is some sequence of frag-ments that are overlapping that will correspond to the whole genome sequence. Frequently,in biological sequencing we are interested in only overlaps of a certain length l(choice of l usually ranges between20and30).An interesting empirical observation is that,in general,as l increases,|O|decreases3for standard datasets with fragments that are sampled uniformly over the original genome. Estimating frequencies of highly repeating overlapsOne useful statistic during overlap detection is the count of such overlap strings which repeat very frequently.Such overlap strings represent the fact that severalfragment pairs exist which have identical overlap regions.Moreover,if the number of such fragment pairs exceeds some pre-defined threshold then it may be a good indicator of the fact that their overlap string is a potential repeat region in the original genome.If the set of fragments is denoted by relations R and S(for convenience),then a SQL-like query to denote overlaps of all possible pairs is given bySELECT OVERLAP(R,S)FROM R,SWHERE R.id<>S.idIf the output of the above query is represented as a relation T,then a SQL-like query to count the number of overlap strings occurring more frequently than a threshold will have a nested sub-query with a complex predicate.Estimating the frequencies of such highly occurring overlap strings is a challenging problem as is evident from the literature[FM:83],[CMN:98],[CCMN:00].Some of the techniques suggest using an auxiliary data structure constructed over T or using a random sample from T for estimation.In fact in[CCMN:00],the authors present a lower bound on the error of a sample-based estimation.In a recent paper[JDPJ:05], the authors have shown that if we have the relation T completely materialized and an index over the relation which can accurately return count information for every overlap string in a random sample of T,it can be possible to estimate the number of highly occurring strings.We observe however that if we can reliably estimate all the“distinct”overlap strings for the fragments in our dataset,then the size of such a relation,say T dist will be much smaller than T.This gives us the benefit that processing all records of T dist will not be a time-consuming task and can be feasible even in our online framework.The estimation of frequencies of each of the strings in T dist would then proceed as follows.We draw a random sample from the set of fragments and for each overlap string o i of T dist,we count the number of fragments in the sample which have o i as a prefix or suffix.This quantity gives us the number of fragments from the sample which would have produced o i as an overlap string. Moreover since we have a random sample of fragments,we can scale up this count to estimate the total number of fragments in R which have o i as an overlap string with some other fragment.If this estimate is k,then the estimate of the frequency of o i will be given by k×(k−1)/2since all possible pairs of the k fragments will produce o i as an overlap.4Ripple Join OverviewIn this Section we describe the ripple join algorithm as proposed by Haas and Hellerstein.In the simplest two-table version of the algorithm one new random tuple is retrieved at each sampling step from the two joining relations,R and S. Thus the sample size grows as the algorithm moves from one sampling step to5another.The most important requirement for the algorithm is that the tuples of the two relations are retrieved in a random order.We can thus imagine that a cartesian product R×S of the two randomized relations is swept out as we increase the number of iterations.At each step,if the sample sizes are denoted by n R and n S, the newly-read tuple from R is then joined with all the n S tuples of S.Also,the newly read tuple from S is then joined with all the n R tuples in the sample.At this point of time letµsam denote the sample-based join result.If we are running a SUM(expr)join query over all tuples satisfying predicate P,thenµsam=(r,s)∈(R n,S n)expr P(r,s)(1)where R n and S n are the sets of tuples that have been read from R and S by the end of the n th sampling step and expr P(r,s)evaluates to expr(r,s)if(r,s)satisfies the WHERE clause and0otherwise.Also,a natural estimator for the actual query result over the entire relations R and S is given by|R|×|S|Nni=1v i andσ2=1n.Here z p is a constant obtained from the ta-ble of the standard normal curve and depends only on the desired confidence level. Sinceσis also unknown,we need to estimateσin order to compute the confi-dence bounds.The standard deviation of the n numbers in the sample is a good choice of an estimator forσto calculate confidence bounds.This theory hinges on a well-known theorem in statistics called the Central-Limit theorem.Given this background,we are now ready to describe the overlap computation problem through a sampling based approach.65A Database Formulation of the Overlap Computation Traditional fragment assembly tools like CAP[Hua92],PCAP[HW03],SBH[PF99] e a number of techniques to avoid an all-pairs comparison,but it is still a time-consuming task.We propose an entirely new database approach to comput-ing O through random sampling,a powerful and well-studied statistical technique. Using this technique,we can determine the set of overlap strings O,by performing only a fraction of the comparisons done by traditional algorithms.Let seq=a1a2...a|seq|be the whole sequence of the genome.The short frag-ment sequencing step produces a list of n readsR={R1,R2,..,R n}where read R i=seq.substring(start i,end i)where1≤start≤|seq|∧start< end≤|seq|.Given the set of strings R,the overlap computation problem is to de-termine the list of overlaps over all pairs of strings in the set.Formally,∀i=j com-pute overlap(R i,R j)=max{suffix(R i)=prefix(R j),suffix(R j)=prefix(R i)} Consider the set of strings R to be a database relation.Now,the overlap com-putation can be solved by the following SQL query:SELECT DISTINCT COUNT(*)FROM R,SWHERE SIZE(OVERLAP(R.s,S.s))==lWe propose a sampling-based approach to report overlaps among fragment pairs.Specifically,we use a ripple join algorithm[HH99]to perform pairwise comparison for overlap detection.The collection of all fragments is viewed as a database relation R where each fragment represents a tuple of the database.The relation R is then randomly shuffled so that the strings appear in a random order within the relation.Logically,we can view a different randomized ordering of the strings as another relation S.We then read a block of records from relation R as well as S into memory and perform an all-pairs comparison over them tofind overlaps.The statistics of the overlapping strings is stored in a table and is scaled up by the inverse of the sampling fractions of R and S.The sample is then grown by reading in the next block from R and S and repeating the previous step while updating the statistics accordingly.This process is continued as more and more blocks are included in the ever-growing sample.The key observation is that the number of fragments having overlap length equal to a threshold,starts to converge as the sample size grows.This is important because this information implies that we will notfind any more fragments satisfy-ing the overlap criteria,by performing pairwise comparisons over the rest of the data set.Since we now have a set of unique overlap strings,the next step is to compute frequencies of each of these overlap strings,that is the number of fragment pairs whose overlap will produce such an overlap string.This is done by drawing a ran-dom sample of the fragments.Each overlap string o i is then compared with all the7fragments in the sample to check if it is either a prefix or a suffix of that fragment. This count information is then scaled up by the inverse of the sampling fraction to obtain an estimate of the total number of fragments which would have o i as an overlap string with some other fragment.If this estimate is denoted as k i,then the frequency of o i is k×(k−1)/2.At this point of time,the overlap information that has been collected by the ripple join should be sufficient to start the next step of a“overlap-layout-consensus”approach to sequence assembly.The pseudocode of the algorithm is given below.———————————————————————————————procedure ComputeOverlapOnline(Relation R,Relation S,int l)//N=the#of unique overlap strings of length l1.List table=∅2.Randomize tuples in R and S3.for(int i=1,j=1;i<|R|,j<|S|;i++,j++)4.for(int k=1;k<j;k++)5.if(overlap(R[i].s,S[k].s)|>l)6.Insert(table,overlap(R[i].s,S[k].s)7.for(int k=1;k<i;k++)8.if(overlap(R[i].s,S[k].s)|>l)9.Insert(table,overlap(R[i].s,S[k].s)10.N=table.size();pute estimates of the frequencies of unique overlap strings from the sample seen so farprocedure Insert(List table,string s)1.for(int i=1;i<table.size();i++)2.if(table[i].s==s)3.table[i].cnt++4.table.add(s,cnt=1)——————————————————————————————–As a contrast to the ripple join approach,a simple traditional algorithm to deter-mine the set of overlap strings is given below.——————————————————————————————–procedure ComputeOverlap(List R,int l)1.List table=∅2.for(int i=1;i<R.size();i++)3.for(int j=i;i<R.size();j++)4.if(overlap(R[i].s,S[k].s)|>l)5.Insert(table,overlap(R[i].s,S[k].s)6.N=table.size();——————————————————————————————–86BenchmarkingIn this section we give some experimental results which confirm the accuracy and usefulness of our online approach.We implemented two approaches to estimating overlaps as outlined in this report:an online approach employing a ripple join al-gorithm,and a conventional all-pairs computation algorithm for detecting overlaps with several optimizations.6.1MethodologyWe tested our algorithms on chromosomes obtained from the Human Genome Project[HS04].We used Human Chromosome20and Human Chromosome21 in our experiments.Each chromosome is approximately30MB in size.We simu-lated the fragment-forming process and partitioned the chromosomes into fragment pieces.Since we start out with an already existing genome sequence,we can con-sider our fragments to be error free.The individual length of each fragment ranges between500to1000base pairs and we employed5-fold coverage.For chromo-some20,we obtained approximately400,000fragments requiring about300MB in disk space.For chromosome21,we obtained around230,000fragments requiring about170MB in disk space.The implementation was done in C++programming language and the experiments were carried out in a dual-processor i386Linux sys-tem with2GB of RAM available.6.2ExperimentsWe executed the following query on the dataset.SELECT DISTINCT COUNT(*)FROM R,SWHERE SIZE(OVERLAP(R.s,S.s))==lWe tested the algorithms with two different settings for overlap length l=4and l=7respectively(i.e we are interested in all unique overlap strings of size l along with a frequency count for each of the value).We began our experiments byfirst running the traditional algorithm to completion.During this process,we keep track of every unique value encountered in a list along with its count.If the same value is encountered again,its count is incremented as described in the Insert procedure (section5).At the end of the process,we have an accurate measure of the actual #of unique overlap strings in the relation along with a frequency value for each of the strings.The frequency table for chromosome20is depicted infigures1-4.For chromosome21the frequency table is plotted infigures9-12.Next,we executed our online approach as discussed in the paper and maintain a estimate of the frequency table.The basic idea is the same.Every time we encounter a value that we haven’t seen before it is inserted into a table and its frequency is estimated on-the-fly from the point of encounter.We can now compare9the estimated frequencies against the actual frequency table to see how well we have done.Since the frequencies are computed in a different order in both the algorithms,we need to calibrate the results before comparing them.In order to do this,we do a simple sort-merge operation to identify common overlap strings in both the table and then assign a common ID to them.In the end,the two frequency tables are compared to determine the accuracy of the online approach.Since the goal is to provide estimates very quickly,we terminated the online approach after 5%of the dataset is processed.The online approach produced a speedup of5-10 times in our experiments.This has to be qualified by the fact that we tested it only on two datasets.Though this can vary and depends heavily on the nature of the underlying data distribution.If particularly accurate estimates are needed then the online approach may need to process as much as60-70%of the dataset.6.3Discussion and Future WorkThe results of our experiments is plotted infigures1-12.The results show that our online approach performs extremely well when N,the#of unique overlap strings is low.The reason for this is that,when N is low(figure7),the individual frequency of each of the unique strings is very high.Thus,there is a greater chance for the online estimation algorithm to be able to quickly narrow down the frequency es-timates.Whereas,when N is high(figure8),the online approach performs well only for those strings which have a relatively high actual frequency.The reason being the opposite of the low N case.When N is high,the individual frequency of any string is likely to be low.This implies we will have reasonable convergence for the individual frequencies but a very slow convergence on N itself.Our results show that the online approach is a promising option to pared to the traditional algorithm we can produce reasonably accurate estimates very early on particularly for frequency values(Figure5and6),thus saving valuable compu-tation time.Often,in genome sequencing,if an overlap occurs very commonly (indicated by a high frequency)it is an indication of a repeat.For such heuristic pruning,the online approach is clearly well-suited.As an extension to this work.we would like to consider a more sophisticated model of fragment assembly where we take in to account errors in fragments.Also, we would like to explore the possibility of using our approach as a plug-in com-ponent into traditional fragment assembly packages and test how much of an im-provement in running time is obtained as well as the accuracy of thefinal assembled genome.It would be interesting to observe the results of our online approach on heavily skewed data(like genomes with large#of repeats or very large#of unique overlap strings).107ConclusionIn this paper we have proposed a novel sampling-based approach to the problemof overlap computation in fragment pared to traditional techniqueswhich rely on heuristics our technique is based on a rigorous mathematical modelthat provides confidence bounds to the user about the accuracy of thefinal result.This enables us to terminate computation as soon as the accuracy of thefinal resultis within acceptable limits.Experimental results validate the utility of our tech-nique.8ReferencesCTO1Tilford C.A,Kawaguchi T et al.A Physical Map of the Human Y Chromo-some,Nature.2001HH99Haas P.J and Hellerstein J.M.Ripple Joins for Online Aggregation.SIG-MOD Conference.1999KM95Kececioglu,J.D and Myers binatorial Algorithms for DNA Se-quence Assembly.Algorithmica.Jan1995.Hua92Huang X.A Contig Assembly Program Based on Sensitive Detection ofFragment Overlaps.Genomics.1992HW03Huang X,Wang J,Aluru S,Yang SP,Hillier L.PCAP:A Whole-GenomeAssembly Program.Genome Res.2003MS00Myers E.W,Sutton G et al.A Whole-Genome Assembly of Drosophila.Science.2000SN77Sanger F,Nicklen S,Coulson AR.DNA sequencing with Chain-TerminatingInhibitors.Proc A,1977PF99Preparata F.P,Frieze A.M,Upfal E.On the Power of Universal Bases inSequencing by Hybridization.3rd Intl.Conf.on Computational MolecularBiology.1999Ve98Venter et al.Shotgun Sequencing of the Human Genome.Science.1998FM:83P.Flajolet,G.N.Martin;Probabilistic Counting,FOCS1983CMN:98S.Chaudhuri,R.Motwani,V.Narasayya;Using Random Sampling for His-togram Construction,SIGMOD1998.CCMN:00M.Charikar,S.Chaudhuri,R.Motwani,V.Narasayya;Towards EstimationError Guarantees for Distinct Values.HS04The Human Genome:Complete Sequence Information.Available at ftp:///genomes/ Accessed on Dec4th,2005.111101001000100001000000 50 100150200 250F r e q u e n c yOverlap String IDActual Frequencies of All Possible Overlap Strings (N = 261)Figure 1:Chromosome20:True Frequency Table for overlap length l =4100100010000 1000000 50 100150200 250F r e q u e n c yOverlap String IDEstimated Frequencies of Possible Overlap Strings (N = 261)Figure 2:Chromosome20:Estimated Frequency Table for overlap length l =4121101001000100000 2000 4000 60008000 10000 12000 14000 16000F r e q u e n c yOverlap String IDActual Frequencies of All Possible Overlap Strings (N = 16280)Figure 3:Chromosome20:True Frequency Table for overlap length l =710 100 1000 10000 0 2000 4000 6000 8000 10000 12000 14000 16000Overlap String ID Estimated Frequencies of Possible Overlap Strings (N = 9402) 100100010000 1000002000 4000 60008000 1000012000F r e q u e n c yOverlap String IDEstimated Frequencies of Possible Overlap Strings (N = 261)Figure 4:Chromosome20:Estimated Frequency Table for overlap length l =713-400-2002004000.5 11.52 2.5F r e q u e n c y% of Total Fragments ProcessedConvergence of Frequency value for a particular string (Actual frequency 50)Figure 5:Chromosome20:Convergence of Frequency Values for a random string in the overlap table l =4-400-2002004000.5 11.52 2.5F r e q u e n c y% of Total Fragments ProcessedConvergence of Frequency value for a particular string (Actual frequency 50)Figure 6:Chromosome20:Convergence of Frequency Values for a random string in the overlap table l =714。