FAQs
LexicMap is mainly designed for sequence alignment with a small number of queries (gene/plasmid/virus/phage sequences) longer than 150 bp by default.
If you want to search some short reads, you need to build the index with small values of -D/--seed-max-desert (default 100) and -d/--seed-in-desert-dist (default 50), e.g., -D 60 -d 30 for 125bp reads, or -D 50 -D 25 for 100bp reads. This will increase indexing time and index size. This is less of a concern with a small number of genomes, such as < 10,000.
If you just want to search long (>1kb) queries for highly similar (>95%) targets, you can build an index with a bigger -D/--seed-max-desert (default 100) and -d/--seed-in-desert-dist (default 50), e.g., -D 300 -d 150. Bigger values decrease the search sensitivity for distant targets, speed up the indexing
process, decrease indexing memory usage, and decrease index size. Alignment speed is largely unaffected.
Note that LexicMap is slow for ultra-long (>1Mb) queries, and the alignment might be fragmented.
Yes. LexicMap mainly supports small genomes including prokaryotic, viral, and plasmid genomes.
Fungal genomes are also supported; remember to increase the value of -g/--max-genome when running lexicmap index,
which is used to skip genomes larger than 15Mb by default.
-g, --max-genome int ► Maximum genome size. Extremely large genomes (e.g., non-isolate
assemblies from Genbank) will be skipped. (default 15000000)
Maximum genome size is about 268 Mb (268,435,456). More precisely:
$total_bases + ($num_contigs - 1) * 1000 <= 268,435,456
as we concatenate contigs with 1000-bp intervals of N’s to reduce the number of sequences to index.
For big and complex genomes, like the human genome (chr1 is ~248 Mb) which has many repetitive sequences, LexicMap alignment can be slow.
Since v0.9.0, you can set a small value of -N, --top-n-chains to keep a few matches if you only want to check if sequences match any position in a human genome, which would be faster.
- For index building, see the hardware requirements.
- For searching, see the hardware requirements.
Use lexicmap utils remerge (available since v0.5.0), which reruns the merging step for an unfinished index.
When to use this command?
- Only one thread is used for merging indexes, which happens when there are a lot (>200 batches) of batches (
$inpu_files / --batch-size) and the value of--max-open-filesis not big enough.- The Slurm/PBS job time limit is almost reached and the merging step won’t be finished before that.
- Disk quota is reached in the merging step.
So you can stop the indexing command by pressing Ctrl + C (make sure it is in the merging step, see the example below), and run lexicmap utils remerge -d index.lmi,
where index.lmi is the output index directory in lexicmap index.
Optionally, you might set bigger values of
the flags --max-open-files and -J/--seed-data-threads if you have hundreds of thousands of input genomes or have set
a small batch size with -b/--batch-size. E.g.,
22:54:24.420 [INFO] merging 297 indexes...
22:54:24.455 [INFO] [round 1]
22:54:24.455 [INFO] batch 1/1, merging 297 indexes to xxx.lmi.tmp/r1_b1 with 1 threads...
Only one thread was used for seed data merging, so it would take a long time.
So we can set a larger --max-open-files, e.g., 4096,
and it would allow (4096-9) / (297+2) = 13.7 threads for merging. Let’s set --seed-data-threads 12.
# specify the maximum open files per process
ulimit -n 4096
lexicmap utils remerge -d index.lmi --max-open-files 4096 --seed-data-threads 12
Update: since v0.8.0, it’s simpler with lexicmap utils subseq.
See more examples.
# Extracting similar sequences for a query gene.
# search matches with query coverage >= 90%
lexicmap search -d demo.lmi/ bench/b.gene_E_faecalis_SecY.fasta -o results.tsv \
--min-qcov-per-hsp 90
# extract matched sequences as FASTA format
lexicmap utils subseq -d demo.lmi -f results.tsv -o results.tsv.aligned.fasta
seqkit head -n 1 results.tsv.aligned.fasta | head -n 3
>NZ_KB944588.1:228637-229935:+ sgenome=GCF_000392875.1 sseqid=NZ_KB944588.1 qcovGnm=100.000 cls=1 hsp=1 qcovHSP=100.000 alenHSP=1299 pident=100.000 gaps=0 qstart=1 qend=1299 sstart=228637 send=229935 sstr=+ slen=274762 evalue=0.00e+00 bitscore=2343
TTGTTCAAGCTATTAAAGAACGCCTTTAAAGTCAAAGACATTAGATCAAAAATCTTATTT
ACAGTTTTAATCTTGTTTGTATTTCGCCTAGGTGCGCACATTACTGTGCCCGGGGTGAAT
For v0.7.0 or older versions:
Yes, lexicmap search has a flag
-a, --all ► Output more columns, e.g., matched sequences. Use this if you
want to output blast-style format with "lexicmap utils 2blast".
to output the CIGAR string and aligned query and subject sequences.
21. cigar, CIGAR string of the alignment. (optional with -a/--all)
22. qseq, Aligned part of query sequence. (optional with -a/--all)
23. sseq, Aligned part of subject sequence. (optional with -a/--all)
24. align, Alignment text ("|" and " ") between qseq and sseq. (optional with -a/--all)
An example:
# Extracting similar sequences for a query gene.
# search matches with query coverage >= 90%
lexicmap search -d gtdb_complete.lmi/ b.gene_E_faecalis_SecY.fasta -o results.tsv \
--min-qcov-per-hsp 90 --all
# extract matched sequences as FASTA format
sed 1d results.tsv | awk -F'\t' '{print ">"$5":"$15"-"$16":"$17"\n"$23;}' \
| seqkit seq -g > results.fasta
seqkit head -n 1 results.fasta | head -n 3
>NZ_JALSCK010000007.1:39224-40522:-
TTGTTCAAGCTATTAAAGAACGCCTTTAAAGTCAAAGACATTAGATCAAAAATCTTATTT
ACAGTTTTAATCTTGTTTGTATTTCGCCTAGGTGCGCACATTACTGTGCCCGGGGTGAAT
And lexicmap util 2blast can help to convert the tabular format to Blast-style format,
see examples.
lexicmap utils subseq can extract subsequences via genome ID, sequence ID and positions. So you can use this information from the search result and extend the regions to extract flanking sequences.
Update: since v0.8.0, you can extract the extended aligned region with lexicmap utils subseq.
See more examples.
It happens if there are some degenerate bases (e.g., N) in the query sequence.
In the indexing step, all degenerate bases are converted to their lexicographically first bases. E.g., N is converted to A.
For query sequences, we don’t convert them.
LexicMap is mainly designed for sequence alignment with a small number of queries against a database with a huge number (millions) of genomes.
There are some ways to improve the search speed of lexicmap search.
How can I know if an index is compatible with a LexicMap version? Should I rebuild an existing index?
LexicMap is under active development, but we are striving to preserve index compatibility as we implement new features and improvements. The change history and compatibility information are available here.
Queries with very many seed matches in large indexes can run out of memory (200GB) during seed collection.
- This happens when searching in huge Logan Project indexes, where each sample (treated as a genome in LexicMap) has a large number of short contigs.
- Searching with long queries, such as plasmids or even bacterial genomes in large indexes (GenBank or AllTheBacteria), may consume hundreds of GB of memory.
See how to trade off speed for memory usage.