last-map-probs

Typical usage

These commands map DNA reads to the human genome:

lastdb -m1111110 hu human/chr*.fa
lastal -Q1 -e120 hu reads.fastq | last-map-probs.py -s150 > myalns.maf

Explanation of typical usage

These commands find alignments with mismap probability <= 0.01 and score >= 150 (-s150). The score threshold should be high enough to avoid random, spurious alignments: otherwise, the mismap probabilities will not be reliable. A threshold of 150 is often reasonable. For instance, if we compare 50 bp reads to the human genome, we expect a random alignment with score >= 150 once every few thousand reads.

The lastal command finds alignments with score >= 120 (-e120). This is because the mismap probability of an alignment with score 150 may depend on alignments with score < 150. Setting e = s-30 seems to work well.

(We used to recommend option -d108, but this is no longer necessary.)

Options

-h, --help	Print a help message and exit.
-m M, --mismap=M
	Don't write alignments with mismap probability > M.
-s S, --score=S
	Don't write alignments with score < S.

Details

• This script can read alignments in either of the formats produced by lastal (maf or tabular).
• The script reads one batch of alignments at a time (by looking for lines starting with "# batch"). If the batches are huge (e.g. because there are no lines starting with "# batch"), it might need too much memory.

Using multiple CPUs

With large datasets, it's important to go faster by using multiple CPUs. One way to do that is by using GNU parallel (http://www.gnu.org/software/parallel/):

parallel --pipe -L4 "lastal -Q1 -e120 hu - | last-map-probs.py -s150" < 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. Beware that older versions of GNU parallel were slow when using --pipe -L, so be sure to use a recent version.

Limitations

• It is possible that two or more alignments reflect the origin of one query sequence, for instance if the query arose by splicing. This script makes no allowance for that possibility.

Method

Suppose one query sequence has three alignments, with scores: s1, s2, s3. The probability that the first alignment is the one that reflects the origin of the query, is:

exp(s1/t) / [exp(s1/t) + exp(s2/t) + exp(s3/t)]

Here, t is a parameter that depends on the scoring scheme: it is written in the lastal header.

Reference

For more information, please see this article:

Incorporating sequence quality data into alignment improves DNA read mapping. Frith MC, Wan R, Horton P. Nucleic Acids Research 2010 38:e100.