segemehl
segemehl
short read mapping with gaps:
Steve Hoffmann, Christian Otto, Stefan Kurtz, David Langenberger, Cynthia Sharma, Joerg Vogel, Philipp Khaitovitch, Peter F. Stadler & Joerg Hackermueller
1. Introduction
segemehl is a software to map short sequencer reads to reference genomes. Unlike other methods, segemehl is able to detect not only mismatches but also insertions and deletions. Furthermore, segemehl is not limited to a specific read length and is able to map primer- or polyadenylation contaminated reads correctly. segemehl implements a matching strategy based on enhanced suffix arrays (ESA). Segemehl now supports the SAM format, reads gzip'ed queries to save both disk and memory space and allows bisulfite sequencing mapping and split read mapping.
2. Matching model
For each suffix of a read, segemehl aims to find the best scoring seed. Seeds might contain insertions, deletions and mismatches (differences). The number of differences allowed within a single seed is user controlled [parameter -D, --differences] and is crucial for the runtime of the program. Subsequently, seeds that undercut the user defined E-value [parameter -E, --evalue] are passed on to an exact semi global alignment procedure. Finally, reads with a minimum accuracy of [-A, --accuracy] percent are reported to the user.
3. Quick Start
To build segemehl from source you simply download the archive and extract it with
> tar -xvzf segemehl_0_*.tar.gz
subsequently go to the new directory and type
> make
First, you have to build the index structures for the reference sequences (i.e. chromosomes, genomes) you want to match against. Please make sure the reference provided in fasta format. To build the index structure, say for some fasta file chr1.fa, you callsegemehl with option -x
If you want to build indices for multiple fasta files, say chr1.fa chr2.fa chr3.fa you simply type
Please note: building and storing of index structures needs memory. The index for the human chromosome 1, e.g., takes a about 4GB of disk space. Building the index needs approximately 6 GB of RAM. Please make sure that your machine has at least 8GB of RAM available. Machines with 32GB of RAM, e.g., allow building a single index for half of the Human Chromosomes at once.
Matching is just as easy. To match your reads (provided in a single multiple fasta file) just type
All results will be written to the file mymap.map
If the index was build from mutliple fasta files, please provide the fasta files in the same order used when building the index. Otherwise you will recieve an error message (m5 keys do not match). If you are sure you have supplied a correct pair of index and reference sequence you can continue and the md5 field will be automatically updated.
4. Threads
Parallel threads [parameter -t, --threads] will make matching much faster on machines with multiple cores. Use them! To start segemehl with, e.g. 8 parallel threads, type
Using threads ten million reads with 35bp can be matched in 30 minutes.
For additional information please type
5. Downloads
Please note: the current version has a dependency on openssl/md5. Please make sure you have the approriate package (most likely a developer package) installed. Indices build with previous versions of should be readable with later versions of the program
Initial version
Bugfix: floppy alignment ends
Optionally now displays the semi-global alignments; header description corrected; parameter -D now defaults to 1 (sufficient in most cases)
Bugfix: Invalid %N$ use detected (reported on Ubuntu systems; -DDIRECTOUTPUT flag); output contains representation of semi-global alignment; sensitivity increased for equally scoring seeds
Bugfix: Unmatched files bug fixed. Several bugs and flaws occuring during generation of large indices fixed.
Bugfix: missing alignments due to missing carry bit for reads of size 32. Thanks go to Jonathan Hoser!
Bugfix: incorrect margins in alignment step. Due to this bug alignments at the ends of reference sequences were not reported. Thanks go to Benoit Nabholz! A problem with chopping fasta files with few sequences into pieces was fixed.
New release:
Critical Bugfix: this version corrects an error that occured in the MD:Z: of some alignments. Important for variant calling:
6. Contact
If you have any further questions, complaints or bug reports please mail to steve (thesign) bioinf (thesign) uni-leipzig (thesign) de
7. Reference
If you use this program in your work you might want to cite:
Hoffmann S, Otto C, Kurtz S, Sharma CM, Khaitovich P, Vogel J, Stadler PF, Hackermueller J: "Fast mapping of short sequences with mismatches, insertions and deletions using index structures", PLoS Comput Biol (2009) vol. 5 (9) pp. e1000502

./segemehl.x -x chr1.idx -d chr1.fa
./segemehl.x -x chr1_2_3.idx -d chr1.fa chr2.fa chr3.fa
./segemehl.x -i chr1.idx -d chr1.fa -q myreads.fa > mymap.map
./segemehl.x -i chr1_2_3.idx -d chr1.fa chr2.fa chr3.fa -q myreads.fa --threads 8 > mymap.map
./segemehl.x --help