forked from PacificBiosciences/DALIGNER
-
Notifications
You must be signed in to change notification settings - Fork 0
Find all significant local alignments between reads
License
Fyphia/DALIGNER
Folders and files
Name | Name | Last commit message | Last commit date | |
---|---|---|---|---|
Repository files navigation
*** PLEASE GO TO THE DAZZLER BLOG (https://dazzlerblog.wordpress.com) FOR TYPESET *** DOCUMENTATION, EXAMPLES OF USE, AND DESIGN PHILOSOPHY. /************************************************************************************\ UPGRADE & DEVELOPER NOTES ! ! ! If you have already performed a big comparison and don't want to recompute all your local alignments in .las files, but do want to use a more recent version of the software that entails a change to the data structures (currently the update on December 31, 2014), please note the routine LAupgrade.Dec.31.2014. This take a .las file, say X.las, as an argument, and writes to standard output the .las file in the new format. The program can be made with "make" but is not by default created when make is called without an argument. For those interested in the details, on December 30, the "alen" and "blen" fields were dropped to save space as they can always be gotten from the underlying DB. \************************************************************************************/ The Daligner Overlap Library Author: Gene Myers First: July 17, 2013 Current: December 31, 2014 The commands below permit one to find all significant local alignments between reads encoded in Dazzler database. The assumption is that the reads are from a PACBIO RS II long read sequencer. That is the reads are long and noisy, up to 15% on average. Recall that a database has a current partition that divides it into blocks of a size that can conveniently be handled by calling the "dalign" overlapper on all the pairs of blocks producing a collection of .las local alignment files that can then be sorted and merged into an ordered sequence of sorted files containing all alignments between reads in the data set. The alignment records are parsimonious in that they do not record an alignment but simply a set of trace points, typically every 100bp or so, that allow the efficient reconstruction of alignments on demand. 1. daligner [-vbAI] [-k<int(14)>] [-w<int(6)>] [-h<int(35)>] [-t<int>] [-M<int>] [-e<double(.70)] [-l<int(1000)] [-s<int(100)>] [-H<int>] [-m<track>]+ <subject:db|dam> <target:db|dam> ... Compare sequences in the trimmed <subject> block against those in the list of <target> blocks searching for local alignments involving at least -l base pairs (default 1000) or more, that have an average correlation rate of -e (default 70%). The local alignments found will be output in a sparse encoding where a trace point on the alignment is recorded every -s base pairs of the a-read (default 100bp). Reads are compared in both orientations and local alignments meeting the criteria are output to one of several created files described below. The -v option turns on a verbose reporting mode that gives statistics on each major step of the computation. The options -k, -h, and -w control the initial filtration search for possible matches between reads. Specifically, our search code looks for a pair of diagonal bands of width 2^w (default 2^6 = 64) that contain a collection of exact matching k-mers (default 14) between the two reads, such that the total number of bases covered by the k-mer hits is h (default 35). k cannot be larger than 32 in the current implementation. If the -b option is set, then the daligner assumes the data has a strong compositional bias (e.g. >65% AT rich), and at the cost of a bit more time, dynamically adjusts k-mer sizes depending on compositional bias, so that the mers used have an effective specificity of 4^k. If there are one or more interval tracks specified with the -m option, then the reads of the DB or DB's to which the mask applies are soft masked with the union of the intervals of all the interval tracks that apply, that is any k-mers that contain any bases in any of the masked intervals are ignored for the purposes of seeding a match. An interval track is a track, such as the "dust" track created by DBdust, that encodes a set of intervals over either the untrimmed or trimmed DB. Invariably, some k-mers are significantly over-represented (e.g. homopolymer runs). These k-mers create an excessive number of matching k-mer pairs and left unaddressed would cause daligner to overflow the available physical memory. One way to deal with this is to explicitly set the -t parameter which suppresses the use of any k-mer that occurs more than t times in either the subject or target block. However, a better way to handle the situation is to let the program automatically select a value of t that meets a given memory usage limit specified (in Gb) by the -M parameter. By default daligner will use the amount of physical memory as the choice for -M. If you want to use less, say only 8Gb on a 24Gb HPC cluster node because you want to run 3 daligner jobs on the node, then specify -M8. Specifying -M0 basically indicates that you do not want daligner to self adjust k-mer suppression to fit within a given amount of memory. For each subject, target pair of blocks, say X and Y, the program reports alignments where the a-read is in X and the b-read is in Y, and vice versa. However, if the -A option is set ("A" for "asymmetric") then just overlaps where the a-read is in X and the b-read is in Y are reported, and if X = Y, then it further reports only those overlaps where the a-read index is less than the b-read index. In either case, if the -I option is set ("I" for "identity") then when X = Y, overlaps between different portions of the same read will also be found and reported. Each found alignment is recorded as -- a[ab,ae] x bo[bb,be] -- where a and b are the indices (in the trimmed DB) of the reads that overlap, o indicates whether the b-read is from the same or opposite strand, and [ab,ae] and [bb,be] are the intervals of a and bo, respectively, that align. The program places these alignment records in files whose name is of the form X.Y.[C|N]#.las where C indicates that the b-reads are complemented and N indicates they are not (both comparisons are performed) and # is the thread that detected and wrote out the collection of alignments contained in the file. That is the file X.Y.O#.las contains the alignments produced by thread # for which the a-read is from X and the b-read is from Y and in orientation O. The command "daligner -A X Y" produces 2*NTHREAD thread files X.Y.?.las and "daligner X Y" produces 4*NTHREAD files X.Y.?.las and Y.X.?.las (unless X=Y in which case only NTHREAD files, X.X.?.las, are produced). By default daligner compares all overlaps between reads in the database that are greater than the minimum cutoff set when the DB or DBs were split, typically 1 or 2 Kbp. However, the HGAP assembly pipeline only wants to correct large reads, say 8Kbp or over, and so needs only the overlaps where the a-read is one of the large reads. By setting the -H parameter to say N, one alters daligner so that it only reports overlaps where the a-read is over N base-pairs long. While the default parameter settings are good for raw Pacbio data, daligner can be used for efficiently finding alignments in corrected reads or other less noisy reads. For example, for mapping applications against .dams we run "daligner -k20 -h60 -e.85" and on corrected reads, we typically run "daligner -k25 -w5 -h60 -e.95 -s500" and at these settings it is very fast. 2. LAsort [-v] <align:las> ... Sort each .las alignment file specified on the command line. For each file it reads in all the overlaps in the file and sorts them in lexicographical order of (a,b,o,ab) assuming each alignment is recorded as a[ab,ae] x b^o[bb,be]. It then writes them all to a file named <align>.S.las (assuming that the input file was <align>.las). With the -v option set then the program reports the number of records read and written. 3. LAmerge [-v] <merge:las> <parts:las> ... Merge the .las files <parts> into a singled sorted file <merge>, where it is assumed that the input <parts> files are sorted. Due to operating system limits, the number of <parts> files must be <= 252. With the -v option set the program reports the # of records read and written. Used correctly, LAmerge and LAsort together allow one to perform an "external" sort that produces a collection of sorted files containing in aggregate all the local alignments found by the daligner, such that their concatenation is sorted in order of (a,b,o,ab). In particular, this means that all the alignments for a given a-read will be found consecutively in one of the files. So computations that need to look at all the alignments for a given read can operate in simple sequential scans of these sorted files. 4. LAshow [-caroUF] [-i<int(4)>] [-w<int(100)>] [-b<int(10)>] <src1:db|dam> [ <src2:db|dam> ] <align:las> [ <reads:FILE> | <reads:range> ... ] LAshow produces a printed listing of the local alignments contained in the specified .las file, where the a- and b-reads come from src1 or from src1 and scr2, respectively. If a file or list of read ranges is given then only the overlaps for which the a-read is in the set specified by the file or list are displayed. See DBshow for an explanation of how the file and list of read ranges are interpreted. If the -F option is set then the roles of the a- and b- reads are reversed in the display. If the -c option is given then a cartoon rendering is displayed, and if -a or -r option is set then an alignment of the local alignment is displayed. The -a option puts exactly -w columns per segment of the display, whereas the -r option puts exactly -w a-read symbols in each segment of the display. The -r display mode is useful when one wants to visually compare two alignments involving the same a-read. If a combination of the -c, -a, and -r flags is set, then the cartoon comes first, then the -a alignment, and lastly the -r alignment. The -i option sets the indent for the cartoon and/or alignment displays, if they are requested. The -b option sets the number of symbols on either side of the aligned segments in an alignment display, and -U specifies that uppercase should be used for DNA sequence instead of the default lowercase. If the -o option is set then only alignments that are proper overlaps (a sequence end occurs at the each end of the alignment) are displayed. 5. LAdump [-cdt] [-o] <src1:db|dam> [ <src2:db|dam> ] <align:las> [ <reads:FILE> | <reads:range> ... ] Like LAshow, LAdump allows one to display the local alignments (LAs) of a subset of the piles in an .las file and select which information to show about them. The difference is that the information is written in a very simple "1-code" ASCII format that makes it easy for one to read and parse the information for further use. For each LA the pair of reads is output on a line. -c requests that one further output the coordinates of the LA segments be output. The -d option requests that the number of difference in the LA be output, and -t requests that the tracepoint information be output. Finally, -o requests that only LAs that are proper overlaps be output. The format is very simple. Each requested piece of information occurs on a line. The first character of every line is a "1-code" character that tells you what information to expect on the line. The rest of the line contains information where each item is separated by a single blank space. The trace point line gives the number of trace point intervals in the LA and is immediately followed by that many lines containing a pair of integers giving the # of differences and b-displacement in each successive trace point interval. P #a #b - (#a,#b) have an LA between them C #ab #ae #bb #be - [#ab,#ae] aligns with [#bb,#be] D # - there are # differences in the LA T #n - there are #n trace point intervals for the LA (#d #y )^#n - there are #d difference aligning the #y bp's of B with the next fixed-size interval of A + X # - Total amount of X (X = P or T) % X # - Maximum amount of X in any pile (X = P or T) @ T # - Maximum number of trace points in any trace 1-code lines that begin with +, %, or @ are always the first lines in the output. They give size information about what is contained in the output. Specifically, '+ X #' gives the total number of LAs (X=P), or the total number of trace point intervals (X=T) in the file . '% X #' gives the maximum number of LAs (X=P) or the maximum number of trace point intervals (X=T) in a given *pile* (collection of LAs all with the same a-read (applies only to sorted .las files). Finally @ T # gives the maximum # of trace point intervals in any trace within the file. 6. LAindex -v <source:las> ... LAindex takes a series of one or more sorted .las files and produces a "pile index" for each one. If the input file has name "X.las", then the name of its index file is ".X.las.idx". For each A-read pile encoded in the .las file, the index contains the offset to the first local alignment with A in the file. The index starts with four 64-bit integers that encode the numbers % P, + T, % T, and @ T described for LAdump above, and then an offset for each pile beginning with the first A-read in the file (which may not be read 0). The index is meant to allow programs that process piles to more efficiently read just the piles they need at any momment int time, as opposed to having to sequentially scan through the .las file. 7. LAcat <source:las> > <target>.las Given argument <source>, find all files <source>.1.las, <source>.2.las, ... <source>.n.<las> where <source>.i.las exists for every i in [1,n]. Then concatenate these files in order into a single .las file and pipe the result to the standard output. 8. LAsplit <target:las> (<parts:int> | <path:db|dam>) < <source>.las If the second argument is an integer n, then divide the alignment file <source>, piped in through the standard input, as evenly as possible into n alignment files with the name <target>.i.las for i in [1,n], subject to the restriction that all alignment records for a given a-read are in the same file. If the second argument refers to a database <path>.db that has been partitioned, then divide the input alignment file into block .las files where all records whose a-read is in <path>.i.db are in <align>.i.las. 9. LAcheck [-vS] <src1:db|dam> [ <src2:db|dam> ] <align:las> ... LAcheck checks each .las file for structural integrity, where the a- and b-sequences come from src1 or from src1 and scr2, respectively. That is, it makes sure each file makes sense as a plausible .las file, e.g. values are not out of bound, the number of records is correct, the number of trace points for a record is correct, and so on. If the -S option is set then it further checks that the alignments are in sorted order. If the -v option is set then a line is output for each .las file saying either the file is OK or reporting the first error. If the -v option is not set then the program runs silently. The exit status is 0 if every file is deemed good, and 1 if at least one of the files looks corrupted. 10. HPCdaligner [-vbAI] [-k<int(14)>] [-w<int(6)>] [-h<int(35)>] [-t<int>] [-M<int>] [-e<double(.70)] [-l<int(1000)] [-s<int(100)>] [-H<int>] [-m<track>]+ [-dal<int(4)>] [-deg<int(25)>] <path:db|dam> [<first:int>[-<last:int>]] HPCdaligner writes a UNIX shell script to the standard output that consists of a sequence of commands that effectively run daligner on all pairs of blocks of a split database and then externally sorts and merges them using LAsort and LAmerge into a collection of alignment files with names <path>.#.las where # ranges from 1 to the number of blocks the data base is split into. These sorted files if concatenated by say LAcat would contain all the alignments in sorted order (of a-read, then b-read, ...). Moreover, all overlaps for a given a-read are guaranteed to not be split across files, so one can run artifact analyzers or error correction on each sorted file in parallel. The data base must have been previously split by DBsplit and all the parameters, except -v, -dal, and -deg, are passed through to the calls to daligner. The defaults for these parameters are as for daligner. The -v flag, for verbose-mode, is also passed to all calls to LAsort and LAmerge. -dal and -deg options are described later. For a database divided into N sub-blocks, the calls to daligner will produce in total 2TN^2 .las files assuming daligner runs with T threads. These will then be sorted and merged into N^2 sorted .las files, one for each block pair. These are then merged in ceil(log_deg N) phases where the number of files decreases geometrically in -deg until there is 1 file per row of the N x N block matrix. So at the end one has N sorted .las files that when concatenated would give a single large sorted overlap file. The -dal option (default 4) gives the desired number of block comparisons per call to daligner. Some must contain dal-1 comparisons, and the first dal-2 block comparisons even less, but the HPCdaligner "planner" does the best it can to give an average load of dal block comparisons per command. The -deg option (default 25) gives the maximum number of files that will be merged in a single LAmerge command. The planner makes the most even k-ary tree of merges, where the number of levels is ceil(log_deg N). If the integers <first> and <last> are missing then the script produced is for every block in the database. If <first> is present then HPCdaligner produces an incremental script that compares blocks <first> through <last> (<last> = <first> if not present) against each other and all previous blocks 1 through <first>-1, and then incrementally updates the .las files for blocks 1 through <first>-1, and creates the .las files for blocks <first> through <last>. Each UNIX command line output by the HPCdaligner can be a batch job (we use the && operator to combine several commands into one line to make this so). Dependencies between jobs can be maintained simply by first running all the daligner jobs, then all the initial sort jobs, and then all the jobs in each phase of the external merge sort. Each of these phases is separated by an informative comment line for your scripting convenience. 9. HPCmapper [-vb] [-k<int(20)>] [-w<int(6)>] [-h<int(50)>] [-t<int>] [-M<int>] [-e<double(.85)] [-l<int(1000)] [-s<int(100)>] [-H<int>] [-m<track>]+ [-dal<int(4)>] [-deg<int(25)>] <ref:db|dam> <reads:db|dam> [<first:int>[-<last:int>]] HPCmapper writes a UNIX shell script to the standard output that consists of a sequence of commands that effectively "maps" every read in the DB <reads> against a reference set of sequences in the DB <ref>, recording all the found local alignments in the sequence of files <ref>.<reads>.1.las, <ref>.<reads>.2.las, ... where <ref>.<reads>.k.las contains the alignments between all of <ref> and the k'th block of <reads>. The parameters are exactly the same as for HPCdaligner save that the -k, -h, and -e defaults are set appropriately for mapping, and the -A and -I options make no sense as <ref> and <reads> are expected to be distinct data sets. If the integers <first> and <last> are missing then the script produced is for every block in the database <reads>. If <first> is present then HPCmapper produces an script that compares blocks <first> through <last> (<last> = <first> if not present) against DAM <ref>. Example: // Recall G.db from the example in DAZZ_DB/README > cat G.db files = 1 1862 G Sim blocks = 2 size = 11 cutoff = 0 all = 0 0 0 1024 1024 1862 1862 > HPCdaligner -mdust -t5 G | csh -v // Run the HPCdaligner script # Dazzler jobs (2) dazzler -d -t5 -mdust G.1 G.1 dazzler -d -t5 -mdust G.2 G.1 G.2 # Initial sort jobs (4) LAsort G.1.G.1.*.las && LAmerge G.L1.1.1 G.1.G.1.*.S.las && rm G.1.G.1.*.S.las LAsort G.1.G.2.*.las && LAmerge G.L1.1.2 G.1.G.2.*.S.las && rm G.1.G.2.*.S.las LAsort G.2.G.1.*.las && LAmerge G.L1.2.1 G.2.G.1.*.S.las && rm G.2.G.1.*.S.las LAsort G.2.G.2.*.las && LAmerge G.L1.2.2 G.2.G.2.*.S.las && rm G.2.G.2.*.S.las # Level 1 jobs (2) LAmerge G.1 G.L1.1.1 G.L1.1.2 && rm G.L1.1.1.las G.L1.1.2.las LAmerge G.2 G.L1.2.1 G.L1.2.2 && rm G.L1.2.1.las G.L1.2.2.las > LAshow -c -a:G -w50 G.1 | more // Take a look at the result ! G.1: 34,510 records 1 9 c [ 0.. 1,876] x [ 9,017..10,825] ( 18 trace pts) 12645 A ---------+====> dif/(len1+len2) = 398/(1876+1808) = 21.61% B <====+--------- 9017 1 ..........gtg-cggt--caggggtgcctgc-t-t-atcgcaatgtta |||*||||**||||||||*||||*|*|*||**|*|*|||| 9008 gagaggccaagtggcggtggcaggggtg-ctgcgtcttatatccaggtta 27.5% 35 ta-ctgggtggttaaacttagccaggaaacctgttgaaataa-acggtgg ||*|||||||||||||*|**|*||*|*||||||*|**|||||*|*||||| 9057 tagctgggtggttaaa-tctg-ca-g-aacctg-t--aataacatggtgg 24.0% 83 -ctagtggcttgccgtttacccaacagaagcataatgaaa-tttgaaagt *||||||||*||||||||*||**||||*|||**|||||||*||||*|||| 9100 gctagtggc-tgccgttt-ccgcacag-agc--aatgaaaatttg-aagt 20.0% 131 ggtaggttcctgctgtct-acatacagaacgacggagcgaaaaggtaccg ||*|||||||||||||*|*||||*|*|*||||||||||*||||||||||* 9144 gg-aggttcctgctgt-tcacat-c-ggacgacggagc-aaaaggtacc- 16.0% ... > LAcat G >G.las // Combine G.1.las & G.2.las into a single .las file > LAshow G G | more // Take another look, now at G.las G: 62,654 records 1 9 c [ 0.. 1,876] x [ 9,017..10,825] : < 398 diffs ( 18 trace pts) 1 38 c [ 0.. 7,107] x [ 5,381..12,330] : < 1,614 diffs ( 71 trace pts) 1 49 n [ 5,493..14,521] x [ 0.. 9,065] : < 2,028 diffs ( 91 trace pts) 1 68 n [12,809..14,521] x [ 0.. 1,758] : < 373 diffs ( 17 trace pts) 1 147 c [ 0..13,352] x [ 854..14,069] : < 2,993 diffs (133 trace pts) 1 231 n [10,892..14,521] x [ 0.. 3,735] : < 816 diffs ( 37 trace pts) 1 292 c [ 3,835..14,521] x [ 0..10,702] : < 2,353 diffs (107 trace pts) 1 335 n [ 7,569..14,521] x [ 0.. 7,033] : < 1,544 diffs ( 70 trace pts) 1 377 c [ 9,602..14,521] x [ 0.. 5,009] : < 1,104 diffs ( 49 trace pts) 1 414 c [ 6,804..14,521] x [ 0.. 7,812] : < 1,745 diffs ( 77 trace pts) 1 415 c [ 0.. 3,613] x [ 7,685..11,224] : < 840 diffs ( 36 trace pts) 1 445 c [ 9,828..14,521] x [ 0.. 4,789] : < 1,036 diffs ( 47 trace pts) 1 464 n [ 0.. 1,942] x [12,416..14,281] : < 411 diffs ( 19 trace pts) ...
About
Find all significant local alignments between reads
Resources
License
Stars
Watchers
Forks
Releases
No releases published
Packages 0
No packages published
Languages
- C 94.5%
- Python 4.4%
- Other 1.1%