lastz \ hg18.chr4.fa galGal3.chr4.fa \ --notransition --step=20 \ --nogapped |
lastz \ hg18.chr4.fa galGal3.chr4.fa |
This document describes installation and usage of the LASTZ sequence alignment program. LASTZ is a drop-in replacement for BLASTZ, 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.
| LASTZ: | A tool for (1) aligning two DNA sequences, and (2) inferring appropriate scoring parameters automatically. | |
|---|---|---|
| Platform: | 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. | |
| Author: | Bob Harris <rsharris at bx dot psu dot edu> | |
| Date: | August 2021 | |
| Mailing list: | http://lists.bx.psu.edu/listinfo/lastz-users |
LASTZ is available from github at https://github.com/lastz/lastz.
A packed archive containing source code for LASTZ is available from the Miller Lab at Penn State. Note: 1.04.15 is the last release that will be available via the Miller Lab website.
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 <somepath>/lastz‑distrib‑X.XX.XX that contains
a src subdirectory (and some others). You may find it convenient
to remove the revision number (‑X.XX.XX) from the directory name.
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
$LASTZ_INSTALL, or by editing the make‑include.mak
file to set the definition of installDir. Also, be sure to add
the directory you choose to your $PATH.
Then to build the LASTZ executable, enter the following commands from bash
or a similar command-line shell (Solaris users should substitute
gmake for make). This will build two executables
(lastz and lastz_D) and copy them into your
installDir.
cd <somepath>/lastz-distrib-X.XX.XX/src
make
make install
The two executables are basically the same program; the only difference is
that lastz uses integer scores, while lastz_D uses
floating-point scores.
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:
make -f Makefile.warnings
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.
A simple self test is included so you can test whether the build succeeded. To run it, enter the following command:
make test
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:
make: *** [test] Error 1
An additional executable (lastz_32) can be built, to handle
genomes larger than 2 gigabases. For details, see the section on
aligning to whole genomes.
Any executable can be built to allow adjacent indels (by default, these are not allowed). For details, see the section on adjacent indels.
LASTZ is designed to preprocess one sequence or set of sequences (which we collectively call the target) and then align several query 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.
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.
The stages are:
The usual flow is as follows (though most of these steps are optional,
and some settings like ‑‑anyornone
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 word 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 seeds, in the target. The seeds are extended to
longer matches called HSPs (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 anchors.
The anchors are then extended to local alignments (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.
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.
For those eager to try it out, here are some illustrative examples to get you started. Detailed reference material begins with the next section.
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.
This example compares human chromosome 4 to chicken chromosome 4. These sequences can be found in the downloads section of the UCSC Genome Browser, and are 191 and 94 megabases long, respectively. To run a quick low-sensitivity alignment of these sequences, use a command like this:
lastz hg18.chr4.fa galGal3.chr4.fa \
--notransition --step=20 --nogapped \
--format=maf > hg18_4_vs_galGal3_4.maf
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
‑‑format=rdotplot output option and
the R statistical package.
(When in MAF format, LASTZ output can be browsed with
the GMAJ interactive viewer for multiple alignments, available from the
Miller Lab at Penn State.)
Using ‑‑notransition lowers
seeding sensitivity and reduces runtime (by a factor of about 10 in this case).
‑‑step=20 also lowers seeding
sensitivity, reducing runtime and also reducing memory consumption (by a factor
of about 3.3 in this case).
‑‑nogapped 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.
|
Figure 1(a)
|
Figure 1(b)
|
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 UCSC Genome Browser
(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 test_data/fake_chimp_reads.2bit).
To see where these reads map onto the human chromosome, use this command:
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
Attaching [unmask] 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
‑‑step=10, we will only be looking for
seeds at every 10th base. Instead of the default seed pattern, we use
‑‑seed=match12 and
‑‑notransition so our
seeds will be exact matches of 12 bases. Instead of the default
x-drop extension method we use
‑‑exact=20 so that a 20-base
exact match is required to qualify as an HSP. Because we are aligning short
reads, we specify
‑‑noytrim so the alignment ends will
not be trimmed back to the highest scoring locations during gapped extension.
We replace the default score set, which is for more distant species, with the
stricter ‑‑match=1,5. This scores
matching bases as +1 and mismatches as −5. We also use
‑‑ambiguous=n so that Ns
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 ‑‑filter=coverage:90 and
‑‑filter=identity:95.
For output, we are only interested in where the reads align, so we use the
‑‑format=general option and specify
that we want the position on the chromosome (name1,
start1, length1) and the read name and orientation
(name2, strand2). 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:
cat hg18_21_vs_reads.dat | grep -v "#" | awk '{print $4}' | sort -u | wc
This example demonstrates the primary
alignment processing stages, using the
α-globin regions of cow and human. This data is included in the LASTZ
distribution in test_data/aglobin.2bit, 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.
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.
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.
Figure 2(c) shows the local alignment blocks resulting from gapped extension of the HSPs. There are four alignment blocks.
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.
|
Figure 2(a)
|
Figure 2(b)
| |||
|
Figure 2(c)
|
Figure 2(d)
| |||
|
Figure 2(e)
|
Figure 2(f)
|
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.
‑‑notrivial
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)).
‑‑self 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
‑‑notrivial 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)).
‑‑self in place of the
query, and also add the ‑‑nomirror
option. In this case LASTZ reports only one copy of each mirror-image pair,
as well as omitting the trivial block (Figure 3(d)).
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.
|
Figure 3(a)
|
Figure 3(b)
| |||
|
Figure 3(c)
|
Figure 3(d)
|
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.
The general format of the LASTZ command line is
lastz <target> [<query>] [<options>]
The angle brackets <> indicate meta-syntactic variables that
should be replaced with your values, while the square ones []
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 "". Elements can appear in
any order, the only constraint being that, if present, the
<query> must appear after the <target>.
Output is generally written to stdout, unless specified otherwise
for a particular option.
The <target> and <query> are usually
just the names of files containing the sequences to be aligned, in either
FASTA, Nib,
or 2Bit format. However they can be
HSX 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 Sequence Specifiers for
details). With certain options such as
‑‑self the <query>
is not needed; otherwise if it is left unspecified the query sequences are read
from stdin
(though this does not work with random-access formats
like 2Bit).
As a special case, the <target> is
omitted when the ‑‑targetcapsule
option is used, since the target sequence is embedded within the capsule file.
For options, the general format is ‑‑<keyword> or
‑‑<keyword>=<value>, but for BLASTZ compatibility
some options also have an alternative syntax
<letter>=<number>.
(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.)
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.
Running the command lastz without any arguments prints a help
message with the most commonly used options, while running
lastz --help
lists all of the options.
| Option | BLASTZ equivalent | Meaning |
--strand=both |
B=2 |
Search both strands. |
--strand=plus |
B=0 |
Search the forward strand only (the one corresponding to the query specifier). |
--strand=minus |
Search only the reverse complement of the query specifier. | |
--self |
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. | |
--nomirror |
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 ‑‑self
option is used.
|
|
--queryhsplimit=<n> |
Discard queries that have more than <n> HSPs. Any queries
that exceed this limit are reported as a warning (to stderr), and
no alignments are reported.
This is useful for mapping reads to a reference genome, when some reads align to too many places in the reference. |
|
--queryhsplimit=nowarn:<n> |
Same as ‑‑queryhsplimit=<n> but warnings for queries that
exceed the limit are witheld.
|
|
--queryhsplimit=keep,nowarn:<n> |
Same as ‑‑queryhsplimit=<n> but queries that exceed the
limit are not discarded and warnings are witheld. For such a query, the first
<n> HSPs found are passed along to downstream processing.
Note that the HSPs reported are not the best |
|
--queryhspbest=<n> |
For queries that have more than <n> HSPs, discard any HSPs
that score below the nth best.
This is useful for mapping reads to a reference genome, when some reads align to too many places in the reference. |
|
--querydepth=<n> |
Stop processing gapped alignments for a query/strand if its ratio of aligned
bases to query length exceeds <n>. A warning is written to
stderr, all alignments for the query/strand are discarded, and processing
continues with the next query (or strand).
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. Moreover, back-end filtering options are not 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. |
|
--querydepth=keep:<n> |
Same as ‑‑querydepth=<n> but any alignments discovered for
this query/strand, before it exceeds the threshold, are reported.
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. |
|
--querydepth=nowarn:<n> |
Same as ‑‑querydepth=<n> but warnings for queries that
exceed the limit are witheld.
|
|
--querydepth=keep,nowarn:<n> |
Same as ‑‑querydepth=<n> but any alignments discovered for
this query/strand, before it exceeds the threshold, are reported and warnings
are witheld.
|
|
--anyornone |
Stop processing after the first qualifying alignment has been found and written to the output, and move on to the next query. "Qualifying" means an alignment that meets all of the thresholds, etc. set by other options as usual. See Any-or-None Alignment for more details. This option is not compatible with chaining or interpolation. | |
| Defaults: |
By default both strands are searched, and the target is assumed to be different
from the query.
If |
These are fundamental parameters for alignment scoring, used in several of the stages.
| Option | BLASTZ equivalent | Meaning |
--scores=<scoring_file> |
Q=<file> |
Read the substitution scores and gap penalties (and possibly other options) from a scoring file. This option cannot be used in conjunction with ‑‑match or inference. |
--match=<reward>[,<penalty>] |
Set the score values for a match (+<reward>)
and mismatch (−<penalty>).
These are both specified as positive values; the "+" and "−" are
implicitly assumed. When <penalty> is not specified,
it is the same as <reward>.
Note that specifying |