LAST finds similar regions between sequences, and aligns them.
For our first example, we wish to find and align similar regions between the human and fugu mitochondrial genomes. You can find these sequences in the examples directory: humanMito.fa and fuguMito.fa. We can compare them like this:
lastdb -c humdb humanMito.fa lastal humdb fuguMito.fa > myalns.maf
The lastdb command creates several files whose names begin with "humdb". The lastal command then compares fuguMito.fa to the humdb files, and writes the alignments to a file called "myalns.maf".
The -c option causes lowercase letters to be soft-masked. Lowercase is often used to indicate repetitive regions, and soft-masking avoids getting uninteresting repetitive alignments.
The output has very long lines, so you need to view it without line-wrapping. For example, with a Unix/Linux/MacOsX command line, you can use:
less -S myalns.maf
Each alignment looks like this:
a score=85 s humanMito 1742 289 + 16571 AGTATAGGCGATAGAAATTGAAACCTGGCGCAAT... s fuguMito 1182 300 + 16447 AGTATAGGAGATAGAAAAGGAA-CTAGGAGCTAT...
The score is a measure of how strong the similarity is. Lines starting with "s" contain: the sequence name, the start coordinate of the alignment, the number of bases spanned by the alignment, the strand, the sequence length, and the aligned bases.
The start coordinates are zero-based. This means that, if the alignment begins right at the start of a sequence, the coordinate is 0. If the strand is "-", the start coordinate is in the reverse strand.
This alignment format is called MAF (multiple alignment format), and it is described in the UCSC Genome FAQ. You can convert it to several other formats using maf-convert, which accompanies LAST.
Although these sequences are much larger, we can compare them in the same way as the mitochondrial genomes:
lastdb -v -c mousedb mouse/chr*.fa lastal -v mousedb cat/chr*.fa > myalns.maf
This might take half a day, and use perhaps 9 gigabytes of memory. If you do not have enough memory, investigate lastdb's -w and -s options. The -v (verbose) option just makes it write progress messages on the screen.
Use the lastdb -p option to indicate that the sequences are proteins:
lastdb -p -c invdb invertebrate.fa lastal invdb vertebrate.fa
Here we use the -F15 option, to specify translated alignment with a score penalty of 15 for frameshifts:
lastdb -p -c protdb proteins.fa lastal -F15 protdb dnas.fa
lastal reports alignments whose score is at least some minimum value, e.g. 40. If this value is too high we may miss meaningful alignments, but if it is too low we may find meaningless alignments.
To solve this dilemma, it is useful to know what alignment scores are likely between completely random sequences. For example, let us find what alignment scores are likely between two random sequences with the same lengths and base frequencies as the human and fugu mitochondrial genomes:
lastdb -x humdb humanMito.fa lastdb -x fugdb fuguMito.fa lastex humdb.prj fugdb.prj
The lastdb commands count bases, and write them in files called humdb.prj and fugdb.prj. The -x option tells it to only count bases and skip its usual preparation steps. The lastex command prints a table of scores and expected numbers of alignments. Here is an abbreviated version:
Score Expected number of alignments 39 8.44e-11 22 0.00805 20 0.0699 12 398
This tells us, for example, that there will be on average 398 alignments of score 12 or more between random sequences with these lengths and base frequencies. Also, 22 is the minimum score such that the average number of alignments is no more than 0.01.
Finally, we can find alignments with score at least 22 like this:
lastdb -c humdb humanMito.fa lastal -e22 humdb fuguMito.fa > myalns.maf
Let's assume we have DNA reads in a file called reads.fastq, in fastq-sanger format. We can align them to the human genome like this:
lastdb -m1111110 humandb human/chr*.fa lastal -Q1 humandb reads.fastq > myalns.maf
The funny-looking -m1111110 option makes it better at finding short, strong alignments. (The default settings are tuned for long, weak alignments.) The -Q1 option indicates that the reads are in fastq-sanger format.
Often, one read will align to more than one genome location. You can use last-map-probs to help judge which location reflects the origin of the read. Please see the "typical usage" recipe in last-map-probs.html.
If you have paired end reads, then last-pair-probs may be useful (see last-pair-probs.html).
Unfortunately, there is more than one fastq format (see http://nar.oxfordjournals.org/content/38/6/1767.long). Recently (2013) fastq-sanger seems to be dominant, but if you have another variant you need to change the -Q option (see lastal.txt).
The default DNA scoring scheme used by lastal is tuned for finding long, weak alignments. It is:
match score = 1, mismatch cost = 1, gap cost = 7 + 1 * (gap length), minimum alignment score = 40
However, if you use option -Q1, it uses a different scoring scheme tuned for finding short, strong alignments:
match score = 6, mismatch cost = 18, gap cost = 21 + 9 * (gap length), minimum alignment score = 180
In the next two examples, we set the scoring scheme by hand.
Suppose we have DNA reads in fasta format (without quality data) instead of fastq. We need to omit the -Q option, but we wish to use the same scoring scheme as -Q1:
lastdb -m1111110 humandb human/chr*.fa lastal -r6 -q18 -a21 -b9 -e180 genomedb reads.fa
In this case we need to use the -Q option, but we wish to find weak alignments:
lastdb -c humandb human/chr*.fa lastal -Q1 -r5 -q5 -a35 -b5 humandb reads.fastq > myalns.maf
This example uses a scaled version of the default alignment scores (1:1:7:1 -> 5:5:35:5). The reason for this is to put them on roughly the same scale as the fastq quality scores.
lastal uses the quality scores to modify the alignment scores, and then rounds the modified scores to integers. By using scaled alignment scores, we reduce the information loss caused by rounding.
WARNING! The default score parameters do not align very short reads. This is because the match score is 6 and the score threshold is 180, so at least 30 high-quality matches are required (or a greater number of low-quality matches). To align very short reads, reduce the score threshold (lastal -e).
If the score threshold is too low, you will get meaningless, random alignments.
You can make LAST more sensitive, at the expense of speed, by increasing lastal's m parameter. The default value is 10. So -m100 makes it more slow and sensitive, and -m1000 makes it much more slow and sensitive.
If you have more than one query sequence, you can go faster by aligning them in parallel. One way to do that is by using GNU parallel (http://www.gnu.org/software/parallel/). (Beware that GNU parallel had some efficiency bugs that were fixed in late 2012 / early 2013, so be sure to use a recent version.)
If you have fasta queries in separate files (e.g. chr*.fa), then instead of this:
lastal mydb chr*.fa > myalns.maf
Try this:
parallel lastal mydb ::: chr*.fa > myalns.maf
If you have fasta queries in one file, then instead of this:
lastal mydb queries.fa > myalns.maf
Try this:
parallel --pipe --recstart '>' lastal mydb - < queries.fa > myalns.maf
If you have fastq queries in one file, then instead of this:
lastal -Q1 mydb reads.fastq > myalns.maf
Try this:
parallel --pipe -L4 lastal -Q1 mydb - < reads.fastq > myalns.maf
(The "-L4" tells it that each fastq record is 4 lines, so there should be no line wrapping or blank lines.)
Consider this alignment:
TGAAGTTAAAGGTATATGAATTCCAATTCTTAACCCCCCTATTAAACGAATATCTTG |||||||| |||||| | || | | | || |||||| ||||||||||| TGAAGTTAGAGGTAT--GGTTTTGAGTAGT----CCTCCTATTTTTCGAATATCTTG
The middle section has such weak similarity that its precise alignment cannot be confidently inferred.
It is sometimes useful to estimate the ambiguity of each column in an alignment. We can do that using lastal option -j4:
lastdb -c humdb humanMito.fa lastal -j4 humdb fuguMito.fa > myalns.maf
The output looks like this:
a score=17 s seqX 0 57 + 57 TGAAGTTAAAGGTATATGAATTCCAATTCTTAACCCCCCTATTAAACGAATATCTTG s seqY 0 51 + 51 TGAAGTTAGAGGTAT--GGTTTTGAGTAGT----CCTCCTATTTTTCGAATATCTTG p %*.14442011.(%##"%$$$$###""!!!""""&'(*,340.,,.~~~~~~~~~~~
The "p" line indicates the probability that each column is wrongly aligned, using a compact code (the same as fastq-sanger format):
Symbol Error probability Symbol Error probability ! 0.79 -- 1 0 0.025 -- 0.032 " 0.63 -- 0.79 1 0.02 -- 0.025 # 0.5 -- 0.63 2 0.016 -- 0.02 $ 0.4 -- 0.5 3 0.013 -- 0.016 % 0.32 -- 0.4 4 0.01 -- 0.013 & 0.25 -- 0.32 5 0.0079 -- 0.01 ' 0.2 -- 0.25 6 0.0063 -- 0.0079 ( 0.16 -- 0.2 7 0.005 -- 0.0063 ) 0.13 -- 0.16 8 0.004 -- 0.005 * 0.1 -- 0.13 9 0.0032 -- 0.004 + 0.079 -- 0.1 : 0.0025 -- 0.0032 , 0.063 -- 0.079 ; 0.002 -- 0.0025 - 0.05 -- 0.063 < 0.0016 -- 0.002 . 0.04 -- 0.05 = 0.0013 -- 0.0016 / 0.032 -- 0.04 > 0.001 -- 0.0013
Note that each alignment is grown from a "core" region, and the ambiguity estimates assume that the core is correctly aligned. The core is indicated by "~" symbols, and it contains exact matches only.
LAST has options to find alignments with optimal column probabilities, instead of optimal score: see lastal.txt.