Commit 9ce640c8 authored by Charles Plessy's avatar Charles Plessy

Imported Upstream version 0.6.0

parent a01f6ca1
......@@ -3,14 +3,14 @@ CXX= g++
CFLAGS= -g -Wall -O2
CXXFLAGS= $(CFLAGS)
DFLAGS= -DHAVE_PTHREAD #-D_FILE_OFFSET_BITS=64
OBJS= utils.o bwt.o bwtio.o bwtaln.o bwtgap.o is.o \
bntseq.o bwtmisc.o bwtindex.o stdaln.o simple_dp.o \
OBJS= QSufSort.o bwt_gen.o utils.o bwt.o bwtio.o bwtaln.o bwtgap.o \
is.o bntseq.o bwtmisc.o bwtindex.o ksw.o stdaln.o simple_dp.o \
bwaseqio.o bwase.o bwape.o kstring.o cs2nt.o \
bwtsw2_core.o bwtsw2_main.o bwtsw2_aux.o bwt_lite.o \
bwtsw2_chain.o bamlite.o
bwtsw2_chain.o bamlite.o fastmap.o bwtsw2_pair.o
PROG= bwa
INCLUDES=
LIBS= -lm -lz -lpthread -Lbwt_gen -lbwtgen
LIBS= -lm -lz -lpthread
SUBDIRS= . bwt_gen
.SUFFIXES:.c .o .cc
......@@ -22,21 +22,11 @@ SUBDIRS= . bwt_gen
all:$(PROG)
lib-recur all-recur clean-recur cleanlocal-recur install-recur:
@target=`echo $@ | sed s/-recur//`; \
wdir=`pwd`; \
list='$(SUBDIRS)'; for subdir in $$list; do \
cd $$subdir; \
$(MAKE) CC="$(CC)" CXX="$(CXX)" DFLAGS="$(DFLAGS)" CFLAGS="$(CFLAGS)" \
INCLUDES="$(INCLUDES)" $$target || exit 1; \
cd $$wdir; \
done;
lib:
bwa:lib-recur $(OBJS) main.o
bwa:$(OBJS) main.o
$(CC) $(CFLAGS) $(DFLAGS) $(OBJS) main.o -o $@ $(LIBS)
QSufSort.o:QSufSort.h
bwt.o:bwt.h
bwtio.o:bwt.h
bwtaln.o:bwt.h bwtaln.h kseq.h
......@@ -44,12 +34,11 @@ bwt1away.o:bwt.h bwtaln.h
bwt2fmv.o:bwt.h
bntseq.o:bntseq.h
bwtgap.o:bwtgap.h bwtaln.h bwt.h
fastmap:bwt.h
bwtsw2_core.o:bwtsw2.h bwt.h bwt_lite.h stdaln.h
bwtsw2_aux.o:bwtsw2.h bwt.h bwt_lite.h stdaln.h
bwtsw2_main.o:bwtsw2.h
cleanlocal:
clean:
rm -f gmon.out *.o a.out $(PROG) *~ *.a
clean:cleanlocal-recur
Release 0.5.10 and 0.6.0 (12 November, 2011)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The 0.6.0 release comes with two major changes. Firstly, the index data
structure has been changed to support genomes longer than 4GB. The forward and
reverse backward genome is now integrated in one index. This change speeds up
BWA-short by about 20% and BWA-SW by 90% with the mapping acccuracy largely
unchanged. A tradeoff is BWA requires more memory, but this is the price almost
all mappers that index the genome have to pay.
Secondly, BWA-SW in 0.6.0 now works with paired-end data. It is more accurate
for highly unique reads and more robust to long indels and structural
variations. However, BWA-short still has edges for reads with many suboptimal
hits. It is yet to know which algorithm is the best for variant calling.
0.5.10 is a bugfix release only and is likely to be the last release in the 0.5
branch unless I find critical bugs in future.
Other notable changes:
* Added the `fastmap' command that finds super-maximal exact matches. It does
not give the final alignment, but runs much faster. It can be a building
block for other alignment algorithms. [0.6.0 only]
* Output the timing information before BWA exits. This also tells users that
the task has been finished instead of being killed or aborted. [0.6.0 only]
* Sped up multi-threading when using many (>20) CPU cores.
* Check I/O error.
* Increased the maximum barcode length to 63bp.
* Automatically choose the indexing algorithm.
* Bugfix: very rare segfault due to an uninitialized variable. The bug also
affects the placement of suboptimal alignments. The effect is very minor.
This release involves quite a lot of tricky changes. Although it has been
tested on a few data sets, subtle bugs may be still hidden. It is *NOT*
recommended to use this release in a production pipeline. In future, however,
BWA-SW may be better when reads continue to go longer. I would encourage users
to try the 0.6 release. I would also like to hear the users' experience. Thank
you.
(0.6.0: 12 November 2011, r85)
Beta Release 0.5.9 (24 January, 2011)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
......
This diff is collapsed.
......@@ -29,12 +29,17 @@
#ifndef __QSUFSORT_H__
#define __QSUFSORT_H__
#include <stdint.h>
#define KEY(V, I, p, h) ( V[ I[p] + h ] )
#define INSERT_SORT_NUM_ITEM 16
void QSufSortSuffixSort(int* __restrict V, int* __restrict I, const int numChar, const int largestInputSymbol,
const int smallestInputSymbol, const int skipTransform);
void QSufSortGenerateSaFromInverse(const int *V, int* __restrict I, const int numChar);
typedef int64_t qsint_t;
#define QSINT_MAX INT64_MAX
void QSufSortSuffixSort(qsint_t* __restrict V, qsint_t* __restrict I, const qsint_t numChar, const qsint_t largestInputSymbol,
const qsint_t smallestInputSymbol, const int skipTransform);
void QSufSortGenerateSaFromInverse(const qsint_t *V, qsint_t* __restrict I, const qsint_t numChar);
#endif
......@@ -163,83 +163,94 @@ void bns_destroy(bntseq_t *bns)
}
}
int64_t bns_fasta2bntseq(gzFile fp_fa, const char *prefix)
{
kseq_t *seq;
char name[1024];
bntseq_t *bns;
bntamb1_t *q;
int l_buf;
unsigned char buf[0x10000];
int32_t m_seqs, m_holes, l, i;
int64_t ret = -1;
FILE *fp;
#define _set_pac(pac, l, c) ((pac)[(l)>>2] |= (c)<<((~(l)&3)<<1))
#define _get_pac(pac, l) ((pac)[(l)>>2]>>((~(l)&3)<<1)&3)
// initialization
seq = kseq_init(fp_fa);
bns = (bntseq_t*)calloc(1, sizeof(bntseq_t));
bns->seed = 11; // fixed seed for random generator
srand48(bns->seed);
m_seqs = m_holes = 8;
bns->anns = (bntann1_t*)calloc(m_seqs, sizeof(bntann1_t));
bns->ambs = (bntamb1_t*)calloc(m_holes, sizeof(bntamb1_t));
q = bns->ambs;
l_buf = 0;
strcpy(name, prefix); strcat(name, ".pac");
fp = xopen(name, "wb");
memset(buf, 0, 0x10000);
// read sequences
while ((l = kseq_read(seq)) >= 0) {
static uint8_t *add1(const kseq_t *seq, bntseq_t *bns, uint8_t *pac, int64_t *m_pac, int *m_seqs, int *m_holes, bntamb1_t **q)
{
bntann1_t *p;
int lasts;
if (bns->n_seqs == m_seqs) {
m_seqs <<= 1;
bns->anns = (bntann1_t*)realloc(bns->anns, m_seqs * sizeof(bntann1_t));
int i, lasts;
if (bns->n_seqs == *m_seqs) {
*m_seqs <<= 1;
bns->anns = (bntann1_t*)realloc(bns->anns, *m_seqs * sizeof(bntann1_t));
}
p = bns->anns + bns->n_seqs;
p->name = strdup((char*)seq->name.s);
p->anno = seq->comment.s? strdup((char*)seq->comment.s) : strdup("(null)");
p->gi = 0; p->len = l;
p->gi = 0; p->len = seq->seq.l;
p->offset = (bns->n_seqs == 0)? 0 : (p-1)->offset + (p-1)->len;
p->n_ambs = 0;
for (i = 0, lasts = 0; i < l; ++i) {
for (i = lasts = 0; i < seq->seq.l; ++i) {
int c = nst_nt4_table[(int)seq->seq.s[i]];
if (c >= 4) { // N
if (lasts == seq->seq.s[i]) { // contiguous N
++q->len;
++(*q)->len;
} else {
if (bns->n_holes == m_holes) {
m_holes <<= 1;
bns->ambs = (bntamb1_t*)realloc(bns->ambs, m_holes * sizeof(bntamb1_t));
if (bns->n_holes == *m_holes) {
(*m_holes) <<= 1;
bns->ambs = (bntamb1_t*)realloc(bns->ambs, (*m_holes) * sizeof(bntamb1_t));
}
q = bns->ambs + bns->n_holes;
q->len = 1;
q->offset = p->offset + i;
q->amb = seq->seq.s[i];
*q = bns->ambs + bns->n_holes;
(*q)->len = 1;
(*q)->offset = p->offset + i;
(*q)->amb = seq->seq.s[i];
++p->n_ambs;
++bns->n_holes;
}
}
lasts = seq->seq.s[i];
{ // fill buffer
if (c >= 4) c = lrand48()&0x3;
if (l_buf == 0x40000) {
fwrite(buf, 1, 0x10000, fp);
memset(buf, 0, 0x10000);
l_buf = 0;
if (c >= 4) c = lrand48()&3;
if (bns->l_pac == *m_pac) { // double the pac size
*m_pac <<= 1;
pac = realloc(pac, *m_pac/4);
memset(pac + bns->l_pac/4, 0, (*m_pac - bns->l_pac)/4);
}
buf[l_buf>>2] |= c << ((3 - (l_buf&3)) << 1);
++l_buf;
_set_pac(pac, bns->l_pac, c);
++bns->l_pac;
}
}
++bns->n_seqs;
bns->l_pac += seq->seq.l;
return pac;
}
int64_t bns_fasta2bntseq(gzFile fp_fa, const char *prefix, int for_only)
{
extern void seq_reverse(int len, ubyte_t *seq, int is_comp); // in bwaseqio.c
kseq_t *seq;
char name[1024];
bntseq_t *bns;
uint8_t *pac = 0;
int32_t m_seqs, m_holes;
int64_t ret = -1, m_pac, l;
bntamb1_t *q;
FILE *fp;
// initialization
seq = kseq_init(fp_fa);
bns = (bntseq_t*)calloc(1, sizeof(bntseq_t));
bns->seed = 11; // fixed seed for random generator
srand48(bns->seed);
m_seqs = m_holes = 8; m_pac = 0x10000;
bns->anns = (bntann1_t*)calloc(m_seqs, sizeof(bntann1_t));
bns->ambs = (bntamb1_t*)calloc(m_holes, sizeof(bntamb1_t));
pac = calloc(m_pac/4, 1);
q = bns->ambs;
strcpy(name, prefix); strcat(name, ".pac");
fp = xopen(name, "wb");
// read sequences
while (kseq_read(seq) >= 0) pac = add1(seq, bns, pac, &m_pac, &m_seqs, &m_holes, &q);
if (!for_only) { // add the reverse complemented sequence
m_pac = (bns->l_pac * 2 + 3) / 4 * 4;
pac = realloc(pac, m_pac/4);
memset(pac + (bns->l_pac+3)/4, 0, (m_pac - (bns->l_pac+3)/4*4) / 4);
for (l = bns->l_pac - 1; l >= 0; --l, ++bns->l_pac)
_set_pac(pac, bns->l_pac, 3-_get_pac(pac, l));
}
xassert(bns->l_pac, "zero length sequence.");
ret = bns->l_pac;
{ // finalize .pac file
ubyte_t ct;
fwrite(buf, 1, (l_buf>>2) + ((l_buf&3) == 0? 0 : 1), fp);
fwrite(pac, 1, (bns->l_pac>>2) + ((bns->l_pac&3) == 0? 0 : 1), fp);
// the following codes make the pac file size always (l_pac/4+1+1)
if (bns->l_pac % 4 == 0) {
ct = 0;
......@@ -253,51 +264,56 @@ int64_t bns_fasta2bntseq(gzFile fp_fa, const char *prefix)
bns_dump(bns, prefix);
bns_destroy(bns);
kseq_destroy(seq);
free(pac);
return ret;
}
int bwa_fa2pac(int argc, char *argv[])
{
int c, for_only = 0;
gzFile fp;
if (argc < 2) {
fprintf(stderr, "Usage: bwa fa2pac <in.fasta> [<out.prefix>]\n");
while ((c = getopt(argc, argv, "f")) >= 0) {
switch (c) {
case 'f': for_only = 1; break;
}
}
if (argc == optind) {
fprintf(stderr, "Usage: bwa fa2pac [-f] <in.fasta> [<out.prefix>]\n");
return 1;
}
fp = xzopen(argv[1], "r");
bns_fasta2bntseq(fp, (argc < 3)? argv[1] : argv[2]);
fp = xzopen(argv[optind], "r");
bns_fasta2bntseq(fp, (optind+1 < argc)? argv[optind+1] : argv[optind], for_only);
gzclose(fp);
return 0;
}
int bns_coor_pac2real(const bntseq_t *bns, int64_t pac_coor, int len, int32_t *real_seq)
int bns_cnt_ambi(const bntseq_t *bns, int64_t pos_f, int len, int *ref_id)
{
int left, mid, right, nn;
if (pac_coor >= bns->l_pac)
err_fatal("bns_coor_pac2real", "bug! Coordinate is longer than sequence (%lld>=%lld).", pac_coor, bns->l_pac);
// binary search for the sequence ID. Note that this is a bit different from the following one...
if (ref_id) {
left = 0; mid = 0; right = bns->n_seqs;
while (left < right) {
mid = (left + right) >> 1;
if (pac_coor >= bns->anns[mid].offset) {
if (pos_f >= bns->anns[mid].offset) {
if (mid == bns->n_seqs - 1) break;
if (pac_coor < bns->anns[mid+1].offset) break;
if (pos_f < bns->anns[mid+1].offset) break; // bracketed
left = mid + 1;
} else right = mid;
}
*real_seq = mid;
// binary search for holes
*ref_id = mid;
}
left = 0; right = bns->n_holes; nn = 0;
while (left < right) {
int64_t mid = (left + right) >> 1;
if (pac_coor >= bns->ambs[mid].offset + bns->ambs[mid].len) left = mid + 1;
else if (pac_coor + len <= bns->ambs[mid].offset) right = mid;
mid = (left + right) >> 1;
if (pos_f >= bns->ambs[mid].offset + bns->ambs[mid].len) left = mid + 1;
else if (pos_f + len <= bns->ambs[mid].offset) right = mid;
else { // overlap
if (pac_coor >= bns->ambs[mid].offset) {
nn += bns->ambs[mid].offset + bns->ambs[mid].len < pac_coor + len?
bns->ambs[mid].offset + bns->ambs[mid].len - pac_coor : len;
if (pos_f >= bns->ambs[mid].offset) {
nn += bns->ambs[mid].offset + bns->ambs[mid].len < pos_f + len?
bns->ambs[mid].offset + bns->ambs[mid].len - pos_f : len;
} else {
nn += bns->ambs[mid].offset + bns->ambs[mid].len < pac_coor + len?
bns->ambs[mid].len : len - (bns->ambs[mid].offset - pac_coor);
nn += bns->ambs[mid].offset + bns->ambs[mid].len < pos_f + len?
bns->ambs[mid].len : len - (bns->ambs[mid].offset - pos_f);
}
break;
}
......
......@@ -70,11 +70,16 @@ extern "C" {
bntseq_t *bns_restore(const char *prefix);
bntseq_t *bns_restore_core(const char *ann_filename, const char* amb_filename, const char* pac_filename);
void bns_destroy(bntseq_t *bns);
int64_t bns_fasta2bntseq(gzFile fp_fa, const char *prefix);
int bns_coor_pac2real(const bntseq_t *bns, int64_t pac_coor, int len, int32_t *real_seq);
int64_t bns_fasta2bntseq(gzFile fp_fa, const char *prefix, int for_only);
int bns_cnt_ambi(const bntseq_t *bns, int64_t pos_f, int len, int *ref_id);
#ifdef __cplusplus
}
#endif
static inline int64_t bns_depos(const bntseq_t *bns, int64_t pos, int *is_rev)
{
return (*is_rev = (pos >= bns->l_pac))? (bns->l_pac<<1) - 1 - pos : pos;
}
#endif
.TH bwa 1 "24 January 2011" "bwa-0.5.9" "Bioinformatics tools"
.TH bwa 1 "12 November 2011" "bwa-0.6.0" "Bioinformatics tools"
.SH NAME
.PP
bwa - Burrows-Wheeler Alignment Tool
......@@ -20,19 +20,19 @@ BWA is a fast light-weighted tool that aligns relatively short sequences
(queries) to a sequence database (targe), such as the human reference
genome. It implements two different algorithms, both based on
Burrows-Wheeler Transform (BWT). The first algorithm is designed for
short queries up to ~200bp with low error rate (<3%). It does gapped
short queries up to ~150bp with low error rate (<3%). It does gapped
global alignment w.r.t. queries, supports paired-end reads, and is one
of the fastest short read alignment algorithms to date while also
visiting suboptimal hits. The second algorithm, BWA-SW, is designed for
long reads with more errors. It performs heuristic Smith-Waterman-like
alignment to find high-scoring local hits (and thus chimera). On
low-error short queries, BWA-SW is slower and less accurate than the
reads longer than 100bp with more errors. It performs a heuristic Smith-Waterman-like
alignment to find high-scoring local hits and split hits. On
low-error short queries, BWA-SW is a little slower and less accurate than the
first algorithm, but on long queries, it is better.
.PP
For both algorithms, the database file in the FASTA format must be
first indexed with the
.B `index'
command, which typically takes a few hours. The first algorithm is
command, which typically takes a few hours for a 3GB genome. The first algorithm is
implemented via the
.B `aln'
command, which finds the suffix array (SA) coordinates of good hits of
......@@ -72,8 +72,7 @@ reimplemented by Yuta Mori.
.TP
.B bwtsw
Algorithm implemented in BWT-SW. This method works with the whole human
genome, but it does not work with database smaller than 10MB and it is
usually slower than IS.
genome.
.RE
.RE
......@@ -260,9 +259,17 @@ Specify the read group in a format like `@RG\\tID:foo\\tSM:bar'. [null]
.B bwasw
bwa bwasw [-a matchScore] [-b mmPen] [-q gapOpenPen] [-r gapExtPen] [-t
nThreads] [-w bandWidth] [-T thres] [-s hspIntv] [-z zBest] [-N
nHspRev] [-c thresCoef] <in.db.fasta> <in.fq>
Align query sequences in the <in.fq> file.
nHspRev] [-c thresCoef] <in.db.fasta> <in.fq> [mate.fq]
Align query sequences in the
.I in.fq
file. When
.I mate.fq
is present, perform paired-end alignment. The paired-end mode only works
for reads Illumina short-insert libraries. In the paired-end mode, BWA-SW
may still output split alignments but they are all marked as not properly
paired; the mate positions will not be written if the mate has multiple
local hits.
.B OPTIONS:
.RS
......@@ -413,20 +420,19 @@ subsequence contains no more than
differences.
.PP
When gapped alignment is disabled, BWA is expected to generate the same
alignment as Eland, the Illumina alignment program. However, as BWA
alignment as Eland version 1, the Illumina alignment program. However, as BWA
change `N' in the database sequence to random nucleotides, hits to these
random sequences will also be counted. As a consequence, BWA may mark a
unique hit as a repeat, if the random sequences happen to be identical
to the sequences which should be unqiue in the database. This random
behaviour will be avoided in future releases.
to the sequences which should be unqiue in the database.
.PP
By default, if the best hit is no so repetitive (controlled by -R), BWA
By default, if the best hit is not highly repetitive (controlled by -R), BWA
also finds all hits contains one more mismatch; otherwise, BWA finds all
equally best hits only. Base quality is NOT considered in evaluating
hits. In paired-end alignment, BWA pairs all hits it found. It further
performs Smith-Waterman alignment for unmapped reads with mates mapped
to rescue mapped mates, and for high-quality anomalous pairs to fix
potential alignment errors.
hits. In the paired-end mode, BWA pairs all hits it found. It further
performs Smith-Waterman alignment for unmapped reads to rescue reads with a
high erro rate, and for high-quality anomalous pairs to fix potential alignment
errors.
.SS Estimating Insert Size Distribution
.PP
......@@ -447,20 +453,20 @@ error output.
.SS Memory Requirement
.PP
With bwtsw algorithm, 2.5GB memory is required for indexing the complete
With bwtsw algorithm, 5GB memory is required for indexing the complete
human genome sequences. For short reads, the
.B `aln'
command uses ~2.3GB memory and the
.B `sampe'
command uses ~3.5GB.
.B aln
command uses ~3.2GB memory and the
.B sampe
command uses ~5.4GB.
.SS Speed
.PP
Indexing the human genome sequences takes 3 hours with bwtsw
algorithm. Indexing smaller genomes with IS or divsufsort algorithms is
several times faster, but requires more memory.
algorithm. Indexing smaller genomes with IS algorithms is
faster, but requires more memory.
.PP
Speed of alignment is largely determined by the error rate of the query
The speed of alignment is largely determined by the error rate of the query
sequences (r). Firstly, BWA runs much faster for near perfect hits than
for hits with many differences, and it stops searching for a hit with
l+2 differences if a l-difference hit is found. This means BWA will be
......@@ -475,36 +481,39 @@ r>0.02.
Pairing is slower for shorter reads. This is mainly because shorter
reads have more spurious hits and converting SA coordinates to
chromosomal coordinates are very costly.
.PP
In a practical experiment, BWA is able to map 2 million 32bp reads to a
bacterial genome in several minutes, map the same amount of reads to
human X chromosome in 8-15 minutes and to the human genome in 15-25
minutes. This result implies that the speed of BWA is insensitive to the
size of database and therefore BWA is more efficient when the database
is sufficiently large. On smaller genomes, hash based algorithms are
usually much faster.
.SH NOTES ON LONG-READ ALIGNMENT
.PP
Command
.B `bwasw'
is designed for long-read alignment. The algorithm behind, BWA-SW, is
similar to BWT-SW, but does not guarantee to find all local hits due to
the heuristic acceleration. It tends to be faster and more accurate if
the resultant alignment is supported by more seeds, and therefore
BWA-SW usually performs better on long queries than on short ones.
On 350-1000bp reads, BWA-SW is several to tens of times faster than the
existing programs. Its accuracy is comparable to SSAHA2, more accurate
than BLAT. Like BLAT, BWA-SW also finds chimera which may pose a
challenge to SSAHA2. On 10-100kbp queries where chimera detection is
important, BWA-SW is over 10X faster than BLAT while being more
sensitive.
BWA-SW can also be used to align ~100bp reads, but it is slower than
the short-read algorithm. Its sensitivity and accuracy is lower than
SSAHA2 especially when the sequencing error rate is above 2%. This is
the trade-off of the 30X speed up in comparison to SSAHA2's -454 mode.
.B bwasw
is designed for long-read alignment. BWA-SW essentially aligns the trie
of the reference genome against the directed acyclic word graph (DAWG) of a
read to find seeds not highly repetitive in the genome, and then performs a
standard Smith-Waterman algorithm to extend the seeds. A key heuristic, called
the Z-best heuristic, is that at each vertex in the DAWG, BWA-SW only keeps the
top Z reference suffix intervals that match the vertex. BWA-SW is more accurate
if the resultant alignment is supported by more seeds, and therefore BWA-SW
usually performs better on long queries or queries with low divergence to the
reference genome.
BWA-SW is perhaps a better choice than BWA-short for 100bp single-end HiSeq reads
mainly because it gives better gapped alignment. For paired-end reads, it is yet
to know whether BWA-short or BWA-SW yield overall better results.
.SH CHANGES IN BWA-0.6
.PP
Since version 0.6, BWA has been able to work with a reference genome longer than 4GB.
This feature makes it possible to integrate the forward and reverse complemented
genome in one FM-index, which speeds up both BWA-short and BWA-SW. As a tradeoff,
BWA uses more memory because it has to keep all positions and ranks in 64-bit
integers, twice larger than 32-bit integers used in the previous versions.
The latest BWA-SW also works for paired-end reads longer than 100bp. In
comparison to BWA-short, BWA-SW tends to be more accurate for highly unique
reads and more robust to relative long INDELs and structural variants.
Nonetheless, BWA-short usually has higher power to distinguish the optimal hit
from many suboptimal hits. The choice of the mapping algorithm may depend on
the application.
.SH SEE ALSO
BWA website <http://bio-bwa.sourceforge.net>, Samtools website
......@@ -529,12 +538,12 @@ If you use the short-read alignment component, please cite the following
paper:
.PP
Li H. and Durbin R. (2009) Fast and accurate short read alignment with
Burrows-Wheeler transform. Bioinformatics, 25, 1754-60. [PMID: 19451168]
Burrows-Wheeler transform. Bioinformatics, 25, 1754-1760. [PMID: 19451168]
.PP
If you use the long-read component (BWA-SW), please cite:
.PP
Li H. and Durbin R. (2010) Fast and accurate long-read alignment with
Burrows-Wheeler transform. Bioinformatics. [PMID: 20080505]
Burrows-Wheeler transform. Bioinformatics, 26, 589-595. [PMID: 20080505]
.SH HISTORY
BWA is largely influenced by BWT-SW. It uses source codes from BWT-SW
......
This diff is collapsed.
This diff is collapsed.
......@@ -31,7 +31,7 @@ void bwa_aln2seq_core(int n_aln, const bwt_aln1_t *aln, bwa_seq_t *s, int set_ma
const bwt_aln1_t *p = aln + i;
if (p->score > best) break;
if (drand48() * (p->l - p->k + 1 + cnt) > (double)cnt) {
s->n_mm = p->n_mm; s->n_gapo = p->n_gapo; s->n_gape = p->n_gape; s->strand = p->a;
s->n_mm = p->n_mm; s->n_gapo = p->n_gapo; s->n_gape = p->n_gape;
s->score = p->score;
s->sa = p->k + (bwtint_t)((p->l - p->k + 1) * drand48());
}
......@@ -67,8 +67,7 @@ void bwa_aln2seq_core(int n_aln, const bwt_aln1_t *aln, bwa_seq_t *s, int set_ma
for (l = q->k; l <= q->l; ++l) {
s->multi[z].pos = l;
s->multi[z].gap = q->n_gapo + q->n_gape;
s->multi[z].mm = q->n_mm;
s->multi[z++].strand = q->a;
s->multi[z++].mm = q->n_mm;
}
rest -= q->l - q->k + 1;
} else { // Random sampling (http://code.activestate.com/recipes/272884/). In fact, we never come here.
......@@ -78,18 +77,19 @@ void bwa_aln2seq_core(int n_aln, const bwt_aln1_t *aln, bwa_seq_t *s, int set_ma
while (x < p) p -= p * j / (i--);
s->multi[z].pos = q->l - i;
s->multi[z].gap = q->n_gapo + q->n_gape;
s->multi[z].mm = q->n_mm;
s->multi[z++].strand = q->a;
s->multi[z++].mm = q->n_mm;
}
rest = 0;
break;
}
}
s->n_multi = z;
/*// the following code removes the primary hit, but this leads to a bug in the PE mode
for (k = z = 0; k < s->n_multi; ++k)
if (s->multi[k].pos != s->sa)
s->multi[z++] = s->multi[k];
s->n_multi = z < n_multi? z : n_multi;
*/
}
}
......@@ -109,54 +109,50 @@ int bwa_approx_mapQ(const bwa_seq_t *p, int mm)
return (23 < g_log_n[n])? 0 : 23 - g_log_n[n];
}
bwtint_t bwa_sa2pos(const bntseq_t *bns, const bwt_t *bwt, bwtint_t sapos, int len, int *strand)
{
bwtint_t pos_f;
int is_rev;
pos_f = bns_depos(bns, bwt_sa(bwt, sapos), &is_rev); // pos_f
*strand = !is_rev;
/* NB: For gapped alignment, pacpos may not be correct, which will be fixed
* in refine_gapped_core(). This line also determines the way "x" is
* calculated in refine_gapped_core() when (ext < 0 && is_end == 0). */
if (is_rev) pos_f = pos_f + 1 < len? 0 : pos_f - len + 1; // mapped to the forward strand
return pos_f; // FIXME: it is possible that pos_f < bns->anns[ref_id].offset
}
/**
* Derive the actual position in the read from the given suffix array
* coordinates. Note that the position will be approximate based on
* whether indels appear in the read and whether calculations are
* performed from the start or end of the read.
*/
void bwa_cal_pac_pos_core(const bwt_t *forward_bwt, const bwt_t *reverse_bwt, bwa_seq_t *seq, const int max_mm, const float fnr)
void bwa_cal_pac_pos_core(const bntseq_t *bns, const bwt_t *bwt, bwa_seq_t *seq, const int max_mm, const float fnr)
{
int max_diff;
int max_diff, strand;
if (seq->type != BWA_TYPE_UNIQUE && seq->type != BWA_TYPE_REPEAT) return;
max_diff = fnr > 0.0? bwa_cal_maxdiff(seq->len, BWA_AVG_ERR, fnr) : max_mm;
if (seq->strand) { // reverse strand only
seq->pos = bwt_sa(forward_bwt, seq->sa);
seq->seQ = seq->mapQ = bwa_approx_mapQ(seq, max_diff);
} else { // forward strand only
/* NB: For gapped alignment, p->pos may not be correct, which
* will be fixed in refine_gapped_core(). This line also
* determines the way "x" is calculated in
* refine_gapped_core() when (ext < 0 && is_end == 0). */
seq->pos = reverse_bwt->seq_len - (bwt_sa(reverse_bwt, seq->sa) + seq->len);
seq->pos = bwa_sa2pos(bns, bwt, seq->sa, seq->len, &strand);
seq->strand = strand;
seq->seQ = seq->mapQ = bwa_approx_mapQ(seq, max_diff);
}
}
void bwa_cal_pac_pos(const char *prefix, int n_seqs, bwa_seq_t *seqs, int max_mm, float fnr)
void bwa_cal_pac_pos(const bntseq_t *bns, const char *prefix, int n_seqs, bwa_seq_t *seqs, int max_mm, float fnr)
{
int i, j;
int i, j, strand;
char str[1024];
bwt_t *bwt;
// load forward SA
strcpy(str, prefix); strcat(str, ".bwt"); bwt = bwt_restore_bwt(str);
strcpy(str, prefix); strcat(str, ".sa"); bwt_restore_sa(str, bwt);
for (i = 0; i != n_seqs; ++i) {
if (seqs[i].strand) bwa_cal_pac_pos_core(bwt, 0, &seqs[i], max_mm, fnr);
for (j = 0; j < seqs[i].n_multi; ++j) {
bwt_multi1_t *p = seqs[i].multi + j;
if (p->strand) p->pos = bwt_sa(bwt, p->pos);
}
}
bwt_destroy(bwt);
// load reverse BWT and SA
strcpy(str, prefix); strcat(str, ".rbwt"); bwt = bwt_restore_bwt(str);
strcpy(str, prefix); strcat(str, ".rsa"); bwt_restore_sa(str, bwt);
for (i = 0; i != n_seqs; ++i) {
if (!seqs[i].strand) bwa_cal_pac_pos_core(0, bwt, &seqs[i], max_mm, fnr);
bwa_cal_pac_pos_core(bns, bwt, &seqs[i], max_mm, fnr);
for (j = 0; j < seqs[i].n_multi; ++j) {
bwt_multi1_t *p = seqs[i].multi + j;
if (!p->strand) p->pos = bwt->seq_len - (bwt_sa(bwt, p->pos) + seqs[i].len);
p->pos = bwa_sa2pos(bns, bwt, p->pos, seqs[i].len, &strand);
p->strand = strand;
}
}
bwt_destroy(bwt);
......@@ -174,7 +170,7 @@ static bwa_cigar_t *refine_gapped_core(bwtint_t l_pac, const ubyte_t *pacseq, in
int l = 0, path_len, ref_len;
AlnParam ap = aln_param_bwa;
path_t *path;
int64_t k, __pos = *_pos > l_pac? (int64_t)((int32_t)*_pos) : *_pos;
int64_t k, __pos = *_pos;
ref_len = len + abs(ext);
if (ext > 0) {
......@@ -192,7 +188,7 @@ static bwa_cigar_t *refine_gapped_core(bwtint_t l_pac, const ubyte_t *pacseq, in
aln_global_core(ref_seq, l, (ubyte_t*)seq, len, &ap, path, &path_len);
cigar = bwa_aln_path2cigar(path, path_len, n_cigar);
if (ext < 0 && is_end_correct) { // fix coordinate for reads mapped on the forward strand
if (ext < 0 && is_end_correct) { // fix coordinate for reads mapped to the forward strand
for (l = k = 0; k < *n_cigar; ++k) {
if (__cigar_op(cigar[k]) == FROM_D) l -= __cigar_len(cigar[k]);
else if (__cigar_op(cigar[k]) == FROM_I) l += __cigar_len(cigar[k]);
......@@ -238,7 +234,6 @@ char *bwa_cal_md1(int n_cigar, bwa_cigar_t *cigar, int len, bwtint_t pos, ubyte_
} else ++u;
}
x += l; y += l;
/* } else if (cigar[k]>>14 == FROM_I || cigar[k]>>14 == 3) { */
} else if (__cigar_op(cigar[k]) == FROM_I || __cigar_op(cigar[k])