forked from lastz/lastz
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathREADME.lastz.html
7893 lines (6944 loc) · 306 KB
/
README.lastz.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
"http://www.w3.org/TR/html4/loose.dtd">
<html>
<head>
<title>LASTZ</title>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<meta http-equiv="Content-Style-Type" content="text/css">
<style type="text/css">
/* see resource://gre/res/html.css for default firefox style sheet */
body { margin: 0 4% 0 3%;
color: black } /* ; background-color: #CCCCCC */
.shaded { background-color: #FFFFFF }
p.vvvlarge { margin-top: 8ex; margin-bottom: 0 }
p.vvlarge { margin-top: 6ex; margin-bottom: 0 }
p.vlarge { margin-top: 4ex; margin-bottom: 0 }
p.large { margin-top: 3ex; margin-bottom: 0 }
p { margin-top: 2ex; margin-bottom: 0 }
p.small { margin-top: 1ex; margin-bottom: 0 }
p.tiny { margin-top: 0.6ex; margin-bottom: 0 }
p.scrollspace { margin-top: 100em; margin-bottom: 0 }
/* margins are relative to header's font size */
h1 { font-size: 1.6em; margin: 4ex 0 4ex } /* T LR B */
h2 { font-size: 1.4em; margin: 1.5ex 0 3ex } /* preceded by <hr.hdr> */
h3 { font-size: 1.1em; margin: 4ex 0 2.5ex }
h4 { font-size: 1em; margin: 2.5ex 0 1ex;
font-weight: normal; font-style: italic }
hr.hdr { margin-top: 4ex; margin-bottom: 0; border: 1px inset }
hr { margin-top: 0; margin-bottom: 0; border: 1px inset }
pre { margin-top: 0.7ex; margin-bottom: 0.7ex }
.notop { margin-top: 0 }
.nobottom { margin-bottom: 0 }
code { padding-left: 0.3ex; padding-right: 0.3ex }
.pad { padding-left: 0.3ex; padding-right: 0.3ex }
.nopad { padding-left: 0; padding-right: 0 }
.ttab { padding-left: 0.3ex; padding-right: 0 }
.stab { padding-left: 1ex; padding-right: 0 }
.mtab { padding-left: 4ex; padding-right: 0 }
.wtab { padding-left: 6ex; padding-right: 0 }
dfn, em { font-weight: normal; font-style: italic }
sub { font-size: 0.6em; vertical-align: baseline }
table { margin-top: 0; margin-bottom: 0;
border-collapse: collapse }
tbody { text-align: left; vertical-align: baseline }
td { padding-top: 0.5ex; padding-bottom: 0.5ex;
padding-left: 0.5ex; padding-right: 0.5ex }
td.indent { padding-left: 1.5ex; padding-right: 0.5ex }
table.withlines { margin-top: 1ex; margin-bottom: 1ex;
border: 2px solid black }
table.withlines tr.newsect { border-top: 2px solid black }
table.withlines tr.sectend { border-bottom: 2px solid black }
table.withlines td { border: 1px dotted gray }
table.nolines { border: none }
table.nolines td { border: none }
table.figure { border: none }
table.figure td { padding-top: 2ex; padding-bottom: 2ex;
padding-left: 0; padding-right: 0;
border: none }
table.holder { border: none }
table.holder td { padding: 0; border: none }
table.holder td.indent { padding-left: 5ex }
table.holder pre { margin: 0 }
table.holder table.withlines { margin-top: 0.5ex; margin-bottom: 0.5ex;
border: 1px solid black }
table.holder table.withlines td { padding-top: 0.3ex; padding-bottom: 0.3ex;
padding-left: 0.3ex; padding-right: 0.3ex;
border: 1px dotted gray }
.usernote { font-style: italic; font-weight: bold;
background-color: #FF9900 }
</style>
</head>
<body>
<p class=vvlarge>
<h1>LASTZ Release 1.04.15,
built August 27, 2021</h1>
TABLE OF CONTENTS
<p class=small>
<ul class=notop>
<li><a href="#intro">Introduction</a>
<li><a href="#avail">Availability</a>
<li><a href="#install">Installation</a>
<li><a href="#build_options">Build Options</a>
<li><a href="#overview">Overview of Processing Stages and Terminology</a>
<li><a href="#examples">Examples</a>
<ul class=notop>
<li><a href="#ex_chrom">Comparing a Human Chromosome and a Chicken Chromosome</a>
<li><a href="#ex_shotgun">Aligning Shotgun Reads to a Human Chromosome</a>
<li><a href="#ex_stages">Seeds, HSPs, Gapped Alignments, Chaining</a>
<li><a href="#ex_self">Aligning a Sequence With Itself</a>
</ul>
<li><a href="#syntax">Command-line Syntax</a>
<ul class=notop>
<li><a href="#options_where">Where to Look</a>
<li><a href="#options_scoring">Scoring</a>
<li><a href="#options_indexing">Indexing</a>
<li><a href="#options_seeding">Seeding</a>
<li><a href="#options_hsps">HSPs (Gap-free Extension)</a>
<li><a href="#options_chaining">Chaining</a>
<li><a href="#options_gapped">Gapped Extension</a>
<li><a href="#options_backend">Back-end Filtering</a>
<li><a href="#options_interpolation">Interpolation</a>
<li><a href="#options_output">Output</a>
<li><a href="#options_housekeeping">Housekeeping</a>
<li><a href="#options_yasra">Shortcuts for Yasra</a>
<li><a href="#options_help">Help</a>
<li><a href="#seq_spec">Sequence Specifiers</a>
</ul>
<li><a href="#stages_detail">Processing Stages in Detail</a>
<ul class=notop>
<li><a href="#stage_input">Target Sequence Input</a>
<li><a href="#stage_inference">Scoring Inference</a>
<li><a href="#stage_indexing">Indexing Target Seed Words</a>
<li><a href="#stage_seeding">Seeding</a>
<li><a href="#stage_gapfree">Gap-free Extension</a>
<li><a href="#stage_chaining">HSP Chaining</a>
<li><a href="#stage_gapped">Gapped Extension</a>
<li><a href="#stage_backend">Back-end Filtering</a>
<li><a href="#stage_interpolation">Interpolation</a>
<li><a href="#stage_output">Alignment Output</a>
</ul>
<li><a href="#formats">File Formats</a>
<ul class=notop>
<li><a href="#fmt_fasta">FASTA</a>
<li><a href="#fmt_fastq">FASTQ</a>
<li><a href="#fmt_nib">Nib</a>
<li><a href="#fmt_2bit">2Bit</a>
<li><a href="#fmt_qdna">Quantum DNA</a>
<li><a href="#fmt_qcode">Quantum Code File</a>
<li><a href="#fmt_names">Sequence Name File</a>
<li><a href="#fmt_mask">Sequence Masking File</a>
<li><a href="#fmt_mask_3fields">Sequence Masking File, Three Fields</a>
<li><a href="#fmt_scoring">Scoring File</a>
<li><a href="#fmt_inference">Inference Control File</a>
<li><a href="#fmt_hsx">HSX (Hashed Sequence Index)</a>
<li><a href="#fmt_capsule">Target Capsule File</a>
<li><a href="#fmt_chores">Alignment Chores File</a>
<li><a href="#fmt_segments">Segment File</a>
<li><a href="#fmt_lav">LAV</a>
<li><a href="#fmt_axt">AXT</a>
<li><a href="#fmt_maf">MAF</a>
<li><a href="#fmt_sam">SAM</a>
<li><a href="#fmt_cigar">CIGAR</a>
<li><a href="#fmt_blastn">BLASTN</a>
<li><a href="#fmt_paf">PAF</a>
<li><a href="#fmt_paf">PAF:wfmash</a>
<li><a href="#fmt_diffs">Differences</a>
<li><a href="#fmt_rdotplot">R Dotplot</a>
<li><a href="#fmt_text">Human-Readable Text</a>
<li><a href="#fmt_general">General Output</a>
<li><a href="#fmt_other">Other Output</a>
</ul>
<li><a href="#advanced">Advanced Topics</a>
<ul class=notop>
<li><a href="#adv_whole_genome">Aligning to Whole Genomes</a>
<li><a href="#adv_adjacent_indels">Adjacent Indels</a>
<li><a href="#adv_coords">Interval Coordinates</a>
<li><a href="#adv_non-acgt">Non-ACGT Characters, Splicing, and Separation</a>
<li><a href="#adv_mangling">Sequence Name Mangling</a>
<li><a href="#adv_patterns">Seed Patterns</a>
<li><a href="#adv_anyornone">Any-or-None Alignment</a>
<li><a href="#adv_shadow">Y-drop Mismatch Shadow</a>
<li><a href="#adv_shingle">Shingle Overlap</a>
<li><a href="#adv_capsule">Using Target Capsule Files</a>
<li><a href="#adv_inference">Inferring Score Sets</a>
<li><a href="#adv_dp">Dynamic Programming Matrix</a>
<li><a href="#adv_shell">Filtering With Shell Commands</a>
<li><a href="#adv_selfmasking">Self-Masking a Sequence</a>
<li><a href="#adv_lotta_subs">Aligning Many Subintervals</a>
</ul>
<li><a href="#blastz">Differences from BLASTZ</a>
<li><a href="#history">Change History</a>
<ul class=notop>
<li><a href="#history_recent">Most Recent Changes</a>
</ul>
<li><a href="#references">References</a>
<li><a href="#acknowledgments">Acknowledgments</a>
</ul>
<p class=vlarge>
<!---->
<!-- Introduction -->
<!---->
<hr class=hdr>
<div><a name="intro"></a></div>
<h2>Introduction</h2>
<p>
This document describes installation and usage of the LASTZ sequence alignment
program. LASTZ is a drop-in replacement for
<a href="http://www.bx.psu.edu/miller_lab/">BLASTZ</a>, and is backward
compatible with BLASTZ’s command-line syntax. That is, it supports
all of BLASTZ’s options but also has additional ones, and may produce
slightly different alignment results.
<p>
<table class=nolines><tbody>
<tr><th>LASTZ:</th><td><span class=mtab></span></td>
<td>A tool for (1) aligning two DNA sequences, and
(2) inferring appropriate scoring parameters automatically.
</td></tr>
<tr><th>Platform:</th><td></td>
<td>This package was developed on a Macintosh OS X system, but should work on
other Unix or Linux platforms with little change (if any). LASTZ is written in
C and compiled with gcc; other C compilers can probably be used by adjusting
the Makefile. Some ancillary tools are written in Python, but only use modules
available in typical python installations.
</td></tr>
<tr><th>Author:</th><td></td>
<td>Bob Harris <rsharris at bx dot psu dot edu>
</td></tr>
<tr><th>Date:</th><td></td>
<td>August 2021
</td></tr>
<tr><th>Mailing list:</th><td></td>
<td>
<a href="http://lists.bx.psu.edu/listinfo/lastz-users"
>http://lists.bx.psu.edu/listinfo/lastz-users</a>
</td></tr>
</tbody></table>
<!---->
<!-- Availability -->
<!---->
<hr class=hdr>
<div><a name="avail"></a></div>
<h2>Availability</h2>
<p>
LASTZ is available from github at
<a href="https://github.com/lastz/lastz/">https://github.com/lastz/lastz</a>.
<p>
A packed archive containing source code for LASTZ is available from the
<a href="http://www.bx.psu.edu/~rsharris/lastz/">Miller Lab</a> at Penn State.
Note: 1.04.15 is the last release that will be available via the Miller Lab
website.
<!---->
<!-- Installation -->
<!---->
<hr class=hdr>
<div><a name="install"></a></div>
<h2>Installation</h2>
<p>
If you have received the distribution as a packed archive, unpack it
by whatever means are appropriate for your computer. The result should be
a directory <code><somepath>/lastz‑distrib‑X.XX.XX</code> that contains
a <code>src</code> subdirectory (and some others). You may find it convenient
to remove the revision number (<code>‑X.XX.XX</code>) from the directory name.
<p>
Before building or installing any of the programs, you will need to tell the
installer where to put the executable, either by setting the shell variable
<code>$LASTZ_INSTALL</code>, or by editing the <code>make‑include.mak</code>
file to set the definition of <code>installDir</code>. Also, be sure to add
the directory you choose to your <code>$PATH</code>.
<p>
Then to build the LASTZ executable, enter the following commands from bash
or a similar command-line shell (Solaris users should substitute
<code>gmake</code> for <code>make</code>). This will build two executables
(<code>lastz</code> and <code>lastz_D</code>) and copy them into your
<code>installDir</code>.
<pre>
cd <somepath>/lastz-distrib-X.XX.XX/src
make
make install
</pre>
The two executables are basically the same program; the only difference is
that <code>lastz</code> uses integer scores, while <code>lastz_D</code> uses
floating-point scores.
<p>
The build process should not report any warnings or errors. Because of this,
the Makefile is set up so that warnings are considered errors and will stop the
build. If you encounter this situation, you can use Makefile.warnings instead:
<pre>
make -f Makefile.warnings
</pre>
This should allow the build to complete, while still reporting the warnings.
You'll need to decide whether the warnings indicate something is really wrong.
Usually they don't, but please report them to the author regardless.
<p>
A simple self test is included so you can test whether the build succeeded.
To run it, enter the following command:
<pre>
make test
</pre>
If the test is successful, you will see no output from this command.
Otherwise, you will see the differences between the expected output and the
output of your build, plus a line that looks like this:
<pre>
make: *** [test] Error 1
</pre>
<!---->
<!-- Build Options -->
<!---->
<hr class=hdr>
<div><a name="build_options"></a></div>
<h2>Build Options</h2>
<p>
An additional executable (<code>lastz_32</code>) can be built, to handle
genomes larger than 2 gigabases. For details, see the section on
<a href="#adv_whole_genome">aligning to whole genomes</a>.
<p>
Any executable can be built to allow adjacent indels (by default, these are
not allowed). For details, see the section on
<a href="#adv_adjacent_indels">adjacent indels</a>.
<!---->
<!-- Overview -->
<!---->
<hr class=hdr>
<div><a name="overview"></a></div>
<h2>Overview of Processing Stages and Terminology</h2>
<p>
LASTZ is designed to preprocess one sequence or set of sequences (which we
collectively call the <dfn>target</dfn>) and then align several
<dfn>query</dfn> sequences to it. The general flow of the program is like a
pipeline: the output of one stage is the input to the next. The user can
choose to skip most stages via command-line options; any stages that are
skipped pass their input along to the next stage unchanged. Two of the stages,
scoring inference and interpolation, are special in that they perform a
miniature version of the pipeline within them.
<p>
Note that the following discussion is a generalization, intended to describe
the basic idea of LASTZ’s operation. There are many exceptions that
depend on the particular options specified.
<p>
The stages are:
<ul class=notop>
<li><a href="#stage_input">read sequence(s) from target file</a>
<li><a href="#stage_inference">infer scoring parameters (if necessary)</a>
<li><a href="#stage_indexing">build a seed word position table for the target</a>
<li>for each sequence in query file:
<ul class=notop>
<li><a href="#stage_seeding">find seeds</a>
<li><a href="#stage_gapfree">perform gap-free extension of seeds to HSPs</a>
<li><a href="#stage_chaining">chain HSPs</a>
<li><a href="#stage_gapped">perform gapped extension of HSPs to local alignments</a>
<li><a href="#stage_backend">filter alignments according to specified criteria</a>
<li><a href="#stage_interpolation">interpolate between alignments</a>
<li><a href="#stage_output">print output</a>
<li>repeat with the reverse complement of the query
</ul>
</ul>
<p>
The usual flow is as follows (though most of these steps are optional,
and some settings like <code><a href="#option_anyornone">‑‑anyornone</a></code>
may affect the processing order).
We first read the target sequence(s) into memory, and use that to build a seed
word position table that will allow us to quickly map any word in the target to
all of the positions where it appears. (For the purposes of this discussion
you can think of a <dfn>word</dfn> as a 12-mer of DNA.) Then we read each
query sequence in turn, processing them more or less independently. We examine
the word starting at each base in the query and use the position table to find
matches, called <dfn>seeds</dfn>, in the target. The seeds are extended to
longer matches called <dfn>HSPs</dfn> (high-scoring segment pairs) and filtered
based on score. The HSPs are chained into the highest-scoring set of syntenic
alignments, and then reduced to single locations called <dfn>anchors</dfn>.
The anchors are then extended to <dfn>local alignments</dfn> (which may contain
gaps) and again filtered by score, followed by back-end filtering to discard
alignment blocks that do not meet specified criteria for certain traits. We
then interpolate, repeating the entire process at a higher sensitivity in the
holes between the alignment blocks. And finally, we write out the alignment
information to a file. Then these steps are repeated with the reverse
complement of the query sequence, before moving on to the next sequence in the
query file.
<p>
The scoring inference stage is not usually performed. Typically it is used
only when sequences for new species are acquired, to create scoring files for
subsequent alignments of those species.
<!---->
<!-- Examples -->
<!---->
<hr class=hdr>
<div><a name="examples"></a></div>
<h2>Examples</h2>
<p>
For those eager to try it out, here are some illustrative examples to get you
started. Detailed reference material begins with the
<a href="#syntax">next section</a>.
<!-- Human vs. Chicken -->
<div><a name="ex_chrom"></a></div>
<h3>Comparing a Human Chromosome and a Chicken Chromosome</h3>
<p>
It is often adequate to use a lower sensitivity level than is achieved with
LASTZ’s defaults. For example, to compare two complete chromosomes, even
for species as distant as human and chicken, the alignment landscape is evident
even at very low sensitivity settings. This can speed up the alignment process
considerably.
<p>
This example compares human chromosome 4 to chicken chromosome 4. These
sequences can be found in the downloads section of the
<a href="http://genome.ucsc.edu/">UCSC Genome Browser</a>, and are 191 and 94
megabases long, respectively. To run a quick low-sensitivity alignment of
these sequences, use a command like this:
<pre>
lastz hg18.chr4.fa galGal3.chr4.fa \
--notransition --step=20 --nogapped \
--format=maf > hg18_4_vs_galGal3_4.maf
</pre>
<p>
This runs in about two and a half minutes on a 2-GHz workstation, requiring
only 400 Mb of RAM. Figure 1(a) shows the results, plotted using the
<code><a href="#option_format">‑‑format=rdotplot</a></code> output option and
the <a href="http://www.r-project.org/">R statistical package</a>.
(When in <a href="#fmt_maf">MAF</a> format, LASTZ output can be browsed with
the GMAJ interactive viewer for multiple alignments, available from the
<a href="http://www.bx.psu.edu/miller_lab/">Miller Lab</a> at Penn State.)
<p>
Using <code><a href="#option_notransition">‑‑notransition</a></code> lowers
seeding sensitivity and reduces runtime (by a factor of about 10 in this case).
<code><a href="#option_step">‑‑step=20</a></code> also lowers seeding
sensitivity, reducing runtime and also reducing memory consumption (by a factor
of about 3.3 in this case).
<code><a href="#option_nogapped">‑‑nogapped</a></code> eliminates the
computation of gapped alignments. The complete alignment process using default
settings (shown in Figure 1(b)) uses 1.3 Gb of RAM and takes 4.5 hours on a
machine running at 2.83 GHz.
<p>
<table class=figure><tbody>
<tr><td>
Figure 1(a)
<p class=small>
<img src="images/human_vs_chicken.png" alt="human vs. chicken: low sensitivity">
<p class=small>
<table class=holder align=center><tr><td>
<pre>
lastz \
hg18.chr4.fa galGal3.chr4.fa \
--notransition --step=20 \
--nogapped
</pre>
</td></tr></table>
</td><td><span class=wtab></span></td><td>
Figure 1(b)
<p class=small>
<img src="images/human_vs_chicken_full.png" alt="human vs. chicken: defaults">
<p class=small>
<table class=holder align=center><tr><td>
<pre>
lastz \
hg18.chr4.fa galGal3.chr4.fa
</pre>
</td></tr></table>
</td></tr>
</tbody></table>
<!-- Shotgun Reads -->
<div><a name="ex_shotgun"></a></div>
<h3>Aligning Shotgun Reads to a Human Chromosome</h3>
<p>
Short read mapping for close species requires parameters very different from
LASTZ’s defaults. This example compares a simulated set of primate shotgun
reads to human chromosome 21. The chromosome can be found in the downloads
section of the <a href="http://genome.ucsc.edu/">UCSC Genome Browser</a>
(it is about 47 megabases). Ten thousand simulated reads were generated by
extracting 60-bp intervals from chimp chr21, subjecting them to mild mutation
(including short gaps), and then truncating them to 50 bp (these are included
in the LASTZ distribution, in <code>test_data/fake_chimp_reads.2bit</code>).
<p>
To see where these reads map onto the human chromosome, use this command:
<pre>
lastz hg18.chr21.fa[unmask] fake_chimp_reads.2bit \
--step=10 --seed=match12 --notransition --exact=20 --noytrim \
--match=1,5 --ambiguous=n \
--filter=coverage:90 --filter=identity:95 \
--format=general:name1,start1,length1,name2,strand2 \
> hg18_21_vs_reads.dat
</pre>
<p>
Attaching <code><a href="#action_unmask">[unmask]</a></code> to the chromosome
filename instructs LASTZ to ignore masking information and treat repeats the
same as any other part of the chromosome, in order to accurately assess the
uniqueness of the read mappings. Since we know the two species are close, we
want to reduce sensitivity. Using
<code><a href="#option_step">‑‑step=10</a></code>, we will only be looking for
seeds at every 10th base. Instead of the default seed pattern, we use
<code><a href="#option_seed_match">‑‑seed=match12</a></code> and
<code><a href="#option_notransition">‑‑notransition</a></code> so our
seeds will be exact matches of 12 bases. Instead of the default
<code><a href="#stage_gapfree">x-drop extension method</a></code> we use
<code><a href="#option_exact">‑‑exact=20</a></code> so that a 20-base
exact match is required to qualify as an HSP. Because we are aligning short
reads, we specify
<code><a href="#option_noytrim">‑‑noytrim</a></code> so the alignment ends will
not be trimmed back to the highest scoring locations during gapped extension.
<p>
We replace the default score set, which is for more distant species, with the
stricter <code><a href="#option_match">‑‑match=1,5</a></code>. This scores
matching bases as +1 and mismatches as −5. We also use
<code><a href="#option_ambign">‑‑ambiguous=n</a></code> so that <code>N</code>s
will be scored appropriately.
We are only interested in alignments that involve nearly an entire read, and
since the species are close we don't want alignments with low identity;
therefore we use <code><a href="#option_coverage">‑‑filter=coverage:90</a></code> and
<code><a href="#option_identity">‑‑filter=identity:95</a></code>.
<p>
For output, we are only interested in where the reads align, so we use the
<code><a href="#option_format">‑‑format=general</a></code> option and specify
that we want the position on the chromosome (<code>name1</code>,
<code>start1</code>, <code>length1</code>) and the read name and orientation
(<code>name2</code>, <code>strand2</code>). This creates a tab-delimited
output file with one line per alignment block, a format that is well-suited for
downstream processing by other programs. For example, to count the number of
different reads we've mapped, we can run this Unix shell command:
<pre>
cat hg18_21_vs_reads.dat | grep -v "#" | awk '{print $4}' | sort -u | wc
</pre>
<!-- Seeds et al. -->
<div><a name="ex_stages"></a></div>
<h3>Seeds, HSPs, Gapped Alignments, Chaining</h3>
<p>
This example demonstrates the primary
<a href="#stages_detail">alignment processing stages</a>, using the
α-globin regions of cow and human. This data is included in the LASTZ
distribution in <code>test_data/aglobin.2bit</code>, and consists of a 70K bp
segment of human DNA and a 66K bp segment of cow DNA. We will follow this
example through the major stages of seeding, gap-free extension, chaining, and
gapped extension.
<p>
Figure 2(a) shows the result of default seeding on a small window (3K bp) in the
middle of these segments. Seeds are short near-matches; in this case each seed
is 19 bp and could have as many as 8 mismatches (12-of-19 with one transition).
There are 338 seeds in this window, but regions where there are many seeds are
indistinguishable from line segments.
<p>
Figure 2(b) shows high-scoring segment pairs, the result of gap-free extension
of the seeds. There are 11 HSPs (only 10 are apparent in the figure, but one
of those is split by a 1-bp shift to the next diagonal). Note that many seeds
were discarded because their extensions were low scoring or overlapped.
<p>
Figure 2(c) shows the local alignment blocks resulting from gapped extension of
the HSPs. There are four alignment blocks.
<p>
Then we zoom out and show the results for the full sequences; the red box
indicates the small region shown in the earlier figures. Figure 2(d) shows
the HSPs, 2(e) shows the gapped alignment blocks, and 2(f) illustrates how
chaining reduces the alignment blocks to a single syntenic line (or two lines,
if there were matches on both strands). Note that one can already tell
quite a bit about how the sequences align just from looking at the HSPs.
<p>
<table class=figure><tbody>
<tr><td>
Figure 2(a)
<p class=small>
<img src="images/aglobin_closeup_seeds.png" alt="alpha-globin: seeds (closeup)">
<p class=small>
<table class=holder align=center><tr><td>
<pre>
lastz \
aglobin.2bit/human[34000..37000] \
aglobin.2bit/cow[35000..38000] \
--nogfextend --nochain --nogapped
</pre>
</td></tr></table>
</td><td><span class=wtab></span></td><td>
Figure 2(b)
<p class=small>
<img src="images/aglobin_closeup_hsps.png" alt="alpha-globin: HSPs (closeup)">
<p class=small>
<table class=holder align=center><tr><td>
<pre>
lastz \
aglobin.2bit/human[34000..37000] \
aglobin.2bit/cow[35000..38000] \
--gfextend --nochain --nogapped
</pre>
</td></tr></table>
</td></tr>
<tr><td>
Figure 2(c)
<p class=small>
<img src="images/aglobin_closeup_gapped.png" alt="alpha-globin: gapped blocks (closeup)">
<p class=small>
<table class=holder align=center><tr><td>
<pre>
lastz \
aglobin.2bit/human[34000..37000] \
aglobin.2bit/cow[35000..38000] \
--gfextend --nochain --gapped
</pre>
</td></tr></table>
</td><td><span class=wtab></span></td><td>
Figure 2(d)
<p class=small>
<img src="images/aglobin_hsps.png" alt="alpha-globin: HSPs">
<p class=small>
<table class=holder align=center><tr><td>
<pre>
lastz \
aglobin.2bit/human \
aglobin.2bit/cow \
--gfextend --nochain --nogapped
</pre>
</td></tr></table>
</td></tr>
<tr><td>
Figure 2(e)
<p class=small>
<img src="images/aglobin_unchained.png" alt="alpha-globin: gapped blocks">
<p class=small>
<table class=holder align=center><tr><td>
<pre>
lastz \
aglobin.2bit/human \
aglobin.2bit/cow \
--gfextend --nochain --gapped
</pre>
</td></tr></table>
</td><td><span class=wtab></span></td><td>
Figure 2(f)
<p class=small>
<img src="images/aglobin_chained.png" alt="alpha-globin: gapped blocks with chaining">
<p class=small>
<table class=holder align=center><tr><td>
<pre>
lastz \
aglobin.2bit/human \
aglobin.2bit/cow \
--gfextend --chain --gapped
</pre>
</td></tr></table>
</td></tr>
</tbody></table>
<!-- Self-Alignment -->
<div><a name="ex_self"></a></div>
<h3>Aligning a Sequence With Itself</h3>
<p>
When a sequence is aligned to itself, the full result will contain mirror-image
copies of each alignment block. It is computationally wasteful to process both
copies. LASTZ can handle this situation in four different ways.
<ol>
<li> Simply give LASTZ the same sequence for both the
<a href="#target">target and query</a>. In this case, LASTZ does not know that
it is aligning a sequence to itself, and performs the full computation on both
copies (Figure 3(a)).
<p>
<li> Specify the <code><a href="#option_notrivial">‑‑notrivial</a></code>
option. This performs the full computation on both copies, but doesn't report
the trivial self-alignment block along the main diagonal (Figure 3(b)).
<p>
<li> Specify the <code><a href="#option_self">‑‑self</a></code> option in place
of the query sequence. LASTZ will save work by computing with only one block
of each mirror-image pair, though it still reports both copies in the output by
reconstructing the second copy from the first. It also invokes
<code>‑‑notrivial</code> automatically to omit the trivial self-alignment block
along the main diagonal. This gives the same output as the previous method,
but runs faster (Figure 3(c)).
<p>
<li> Specify <code><a href="#option_self">‑‑self</a></code> in place of the
query, and also add the <code><a href="#option_nomirror">‑‑nomirror</a></code>
option. In this case LASTZ reports only one copy of each mirror-image pair,
as well as omitting the trivial block (Figure 3(d)).
</ol>
<p>
In the following figure, we suppose we have a sequence with repeated motifs,
in the order
α1 β1 γ1 β2 δ1 α2 δ2′ γ2.
That is, α1 and α2 are ancient duplications, as are β1 and
β2, and γ1 and γ2. δ2′ is an inversion, a
reverse-complement duplicate of δ1.
<p>
<table class=figure><tbody>
<tr><td>
Figure 3(a)
<p class=small>
<img src="images/seq_vs_same.png" alt="rearranged sequence: vs. itself, default options">
<p class=small>
<table class=holder align=center><tr><td>
<pre>
lastz target target
</pre>
</td></tr></table>
</td><td><span class=wtab></span></td><td>
Figure 3(b)
<p class=small>
<img src="images/seq_vs_same_notrivial.png" alt="rearranged sequence: vs. itself, --notrivial">
<p class=small>
<table class=holder align=center><tr><td>
<pre>
lastz target target --notrivial
</pre>
</td></tr></table>
</td></tr>
<tr><td>
Figure 3(c)
<p class=small>
<img src="images/seq_vs_self.png" alt="rearranged sequence: --self">
<p class=small>
<table class=holder align=center><tr><td>
<pre>
lastz target --self
</pre>
</td></tr></table>
</td><td><span class=wtab></span></td><td>
Figure 3(d)
<p class=small>
<img src="images/seq_vs_self_no_mirror.png" alt="rearranged sequence: --self --nomirror">
<p class=small>
<table class=holder align=center><tr><td>
<pre>
lastz target --self --nomirror
</pre>
</td></tr></table>
</td></tr>
</tbody></table>
<!---->
<!-- Command-line Syntax -->
<!---->
<hr class=hdr>
<div><a name="syntax"></a></div>
<h2>Command-line Syntax</h2>
<p>
If you are familiar with BLASTZ, you can run LASTZ the same way you ran BLASTZ,
with the same options and input files. In addition to this BLASTZ compatibility,
LASTZ provides other options.
<p>
The general format of the LASTZ command line is
<pre>
lastz <target> [<query>] [<options>]
</pre>
<p>
The angle brackets <code><></code> indicate meta-syntactic variables that
should be replaced with your values, while the square ones <code>[]</code>
indicate elements that are optional. Spaces separate fields on the command
line; a field that needs to contain a space (e.g. within a file name) must be
enclosed in double quotes <code>""</code>. Elements can appear in
any order, the only constraint being that, if present, the
<code><query></code> must appear after the <code><target></code>.
Output is generally written to <code>stdout</code>, unless specified otherwise
for a particular option.
<div><a name="target"></a></div>
<p>
The <code><target></code> and <code><query></code> are usually
just the names of files containing the sequences to be aligned, in either
<a href="#fmt_fasta">FASTA</a>, <a href="#fmt_nib">Nib</a>,
or <a href="#fmt_2bit">2Bit</a> format. However they can be
<a href="#fmt_hsx">HSX</a> index files that refer to the sequences indirectly,
and they also can specify pre-processing actions such as selecting a
subsequence from the file (see <a href="#seq_spec">Sequence Specifiers</a> for
details). With certain options such as
<code><a href="#option_self">‑‑self</a></code> the <code><query></code>
is not needed; otherwise if it is left unspecified the query sequences are read
from <code>stdin</code>
(though this does not work with random-access formats
like <a href="#fmt_2bit">2Bit</a>).
As a special case, the <code><target></code> is
omitted when the <code><a href="#option_capsule">‑‑targetcapsule</a></code>
option is used, since the target sequence is embedded within the capsule file.
<p>
For options, the general format is <code>‑‑<keyword></code> or
<code>‑‑<keyword>=<value></code>, but for BLASTZ compatibility
some options also have an alternative syntax
<code><letter>=<number></code>.
(Be careful when copying options from the tables below, as some of the hyphens
here are special characters to avoid awkward line wrapping in certain web
browsers. If you have trouble, replace the pasted hyphens with ordinary typed
ones on your command line.)
<p>
Please understand that LASTZ is a complex program and its options are not all
independent, i.e., some options are not valid in combination with certain
others. It would be difficult and cumbersome to attempt to list every possible
conflict here; instead we just mention some of the major ones. If you are not
sure about a particular combination, go ahead and try it — LASTZ will
tell you if it’s not allowed.
<p>
Running the command <code>lastz</code> without any arguments prints a help
message with the most commonly used options, while running
<pre>
lastz --help
</pre>
lists all of the options.
<!-- Where to Look -->
<div><a name="options_where"></a></div>
<h3>Where to Look</h3>
<p>
<table class=withlines><tbody>
<tr class=sectend><td>Option</td><td>BLASTZ equivalent</td><td>Meaning</td></tr>
<tr>
<td><a name="option_strand"></a><code>--strand=both</code></td>
<td><code>B=2</code></td>
<td>
Search both strands.
</td>
</tr>
<tr>
<td><code>--strand=plus</code></td>
<td><code>B=0</code></td>
<td>
Search the forward strand only (the one corresponding to the query specifier).
</td>
</tr>
<tr>
<td><code>--strand=minus</code></td>
<td></td>
<td>
Search only the reverse complement of the query specifier.
</td>
</tr>
<tr>
<td><a name="option_self"></a><code>--self</code></td>
<td></td>
<td>
Perform a self-alignment: the target sequence is also the query.
Computation is more efficient than it would be without this option, since only
one of each mirror-image pair of alignment blocks is processed (the other,
redundant one is skipped during processing, but re-created in the output).
Also, the trivial self-alignment block along the main diagonal is omitted from
the output.
This option cannot be used if the target is comprised of multiple sequences.
</td>
</tr>
<tr>
<td><a name="option_nomirror"></a><code>--nomirror</code></td>
<td></td>
<td>
Inhibit the re-creation of mirror-image alignments. Output consists of only
one copy of each meaningful alignment block in a self-alignment. This option
is only applicable when the <code><a href="#option_self">‑‑self</a></code>
option is used.
</td>
</tr>
<tr>
<td><a name="option_queryhsplimit"></a><code>--queryhsplimit=<n></code></td>
<td></td>
<td>
Discard queries that have more than <code><n></code> HSPs. Any queries
that exceed this limit are reported as a warning (to <code>stderr</code>), and
no alignments are reported.
<p class=small>
This is useful for mapping reads to a reference genome, when some reads align
to too many places in the reference.
</td>
</tr>
<tr>
<td><code>--queryhsplimit=nowarn:<n></code></td>
<td></td>
<td>
Same as <code>‑‑queryhsplimit=<n></code> but warnings for queries that
exceed the limit are witheld.
</td>
</tr>
<tr>
<td><code>--queryhsplimit=keep,nowarn:<n></code></td>
<td></td>
<td>
Same as <code>‑‑queryhsplimit=<n></code> but queries that exceed the
limit are not discarded and warnings are witheld. For such a query, the first
<code><n></code> HSPs found are passed along to downstream processing.
<p class=small>
Note that the HSPs reported are not the best <code><n></code> HSPs. They
are simply the first <code><n></code> found; they very likely have a
positional bias.
</td>
</tr>
<tr>
<td><a name="option_queryhspbest"></a><code>--queryhspbest=<n></code></td>
<td></td>
<td>
For queries that have more than <code><n></code> HSPs, discard any HSPs
that score below the <code>n</code>th best.
<p class=small>
This is useful for mapping reads to a reference genome, when some reads align
to too many places in the reference.
</td>
</tr>
<tr>
<td><a name="option_querydepth"></a><code>--querydepth=<n></code></td>
<td></td>
<td>
Stop processing gapped alignments for a query/strand if its ratio of aligned
bases to query length exceeds <code><n></code>. A warning is written to
stderr, all alignments for the query/strand are discarded, and processing
continues with the next query (or strand).
<p class=small>
<code>‑‑querydepth=keep:<n></code> can be used if the preference is to
keep some alignments for such query/strands.
<p class=small>
<code><n></code> is a real number and corresponds to a depth of coverage
threshold. For example, a value of <code>5.0</code> would cause termination
once a query/strand has an average of five alignments for every base in the
query. The numerator is the number of matches or substitutions (but not gaps);
the denominator is the length of the query sequence.
<p class=small>
The purpose of this option is one of saving time. It is useful for
automatically terminating the processing of queries with high repeat content,
for which other methods of dealing with repetitive content fail.
<p class=small>
Moreover, <a href="#options_backend">back-end filtering</a> options are
<em>not</em> considered. In other words, matches are counted for any alignment
that meets the scoring threshold, regardless of whether that alignment would be
reported. The justification is that we are trying to abort the processing of
queries that have too many bounding alignments in the DP matrix, and back-end
filtering occurs later in the process.
</td>
</tr>
<tr>
<td><code>--querydepth=keep:<n></code></td>
<td></td>
<td>
Same as <code>‑‑querydepth=<n></code> but any alignments discovered for
this query/strand, before it exceeds the threshold, are reported.
<p class=small>
Note that the alignments reported are not guaranteed to be the highest scoring
alignments that would achieve the threshold. They are simply the first
alignments found. In other words, the purpose of this option is one of saving
time, not one of finding optimal alignments.
</td>
</tr>
<tr>
<td><code>--querydepth=nowarn:<n></code></td>
<td></td>
<td>
Same as <code>‑‑querydepth=<n></code> but warnings for queries that
exceed the limit are witheld.