diff --git a/bwa.1 b/bwa.1 index 1edbf12..9d45a3d 100644 --- a/bwa.1 +++ b/bwa.1 @@ -1,4 +1,4 @@ -.TH bwa 1 "15 March 2013" "bwa-0.7.3a" "Bioinformatics tools" +.TH bwa 1 "15 March 2013" "bwa-0.7.4" "Bioinformatics tools" .SH NAME .PP bwa - Burrows-Wheeler Alignment Tool @@ -220,7 +220,7 @@ deducted. [5] Penalty for an unpaired read pair. BWA-MEM scores an unpaired read pair as .RI scoreRead1+scoreRead2- INT and scores a paired as scoreRead1+scoreRead2-insertPenalty. It compares these -two scores to determine whether we should force pairing. [9] +two scores to determine whether we should force pairing. [17] .TP .B -p Assume the first input query file is interleaved paired-end FASTA/Q. See the command description for details. diff --git a/bwamem.c b/bwamem.c index e188da4..8c2ba9c 100644 --- a/bwamem.c +++ b/bwamem.c @@ -46,7 +46,7 @@ mem_opt_t *mem_opt_init() o->a = 1; o->b = 4; o->q = 6; o->r = 1; o->w = 100; o->T = 30; o->zdrop = 100; - o->pen_unpaired = 9; + o->pen_unpaired = 17; o->pen_clip = 5; o->min_seed_len = 19; o->split_width = 10; diff --git a/bwape.c b/bwape.c index c413b34..b219527 100644 --- a/bwape.c +++ b/bwape.c @@ -106,6 +106,11 @@ static int infer_isize(int n_seqs, bwa_seq_t *seqs[2], isize_info_t *ii, double tmp = (int)(p25 - OUTLIER_BOUND * (p75 - p25) + .499); ii->low = tmp > max_len? tmp : max_len; // ii->low is unsigned ii->high = (int)(p75 + OUTLIER_BOUND * (p75 - p25) + .499); + if (ii->low > ii->high) { + fprintf(stderr, "[infer_isize] fail to infer insert size: upper bound is smaller than read length\n"); + free(isizes); + return -1; + } for (i = 0, x = n = 0; i < tot; ++i) if (isizes[i] >= ii->low && isizes[i] <= ii->high) ++n, x += isizes[i]; @@ -404,7 +409,7 @@ bwa_cigar_t *bwa_sw_core(bwtint_t l_pac, const ubyte_t *pacseq, int len, const u bwa_cigar_t *cigar = 0; ubyte_t *ref_seq; bwtint_t k, x, y, l; - int xtra; + int xtra, gscore; int8_t mat[25]; bwa_fill_scmat(1, 3, mat); @@ -422,12 +427,12 @@ bwa_cigar_t *bwa_sw_core(bwtint_t l_pac, const ubyte_t *pacseq, int len, const u // do alignment xtra = KSW_XSUBO | KSW_XSTART | (len < 250? KSW_XBYTE : 0); r = ksw_align(len, (uint8_t*)seq, l, ref_seq, 5, mat, 5, 1, xtra, 0); - ksw_global(r.qe - r.qb + 1, &seq[r.qb], r.te - r.tb + 1, &ref_seq[r.tb], 5, mat, 5, 1, 50, n_cigar, &cigar32); + gscore = ksw_global(r.qe - r.qb + 1, &seq[r.qb], r.te - r.tb + 1, &ref_seq[r.tb], 5, mat, 5, 1, 50, n_cigar, &cigar32); cigar = (bwa_cigar_t*)cigar32; for (k = 0; k < *n_cigar; ++k) cigar[k] = __cigar_create((cigar32[k]&0xf), (cigar32[k]>>4)); - if (r.score < SW_MIN_MATCH_LEN || r.score2 == r.score) { // poor hit or tandem hits + if (r.score < SW_MIN_MATCH_LEN || r.score2 == r.score || gscore != r.score) { // poor hit or tandem hits or weird alignment free(cigar); free(ref_seq); *n_cigar = 0; return 0; } diff --git a/bwtsw2_pair.c b/bwtsw2_pair.c index f0ced11..df62e3f 100644 --- a/bwtsw2_pair.c +++ b/bwtsw2_pair.c @@ -47,7 +47,7 @@ bsw2pestat_t bsw2_stat(int n, bwtsw2_t **buf, kstring_t *msg, int max_ins) p75 = isize[(int)(.75 * k + .499)]; ksprintf(msg, "[%s] infer the insert size distribution from %d high-quality pairs.\n", __func__, k); if (k < 8) { - ksprintf(msg, "[%s] fail to infer the insert size distribution.\n", __func__); + ksprintf(msg, "[%s] fail to infer the insert size distribution: too few good pairs.\n", __func__); free(isize); r.failed = 1; return r; @@ -56,6 +56,12 @@ bsw2pestat_t bsw2_stat(int n, bwtsw2_t **buf, kstring_t *msg, int max_ins) r.low = tmp > max_len? tmp : max_len; if (r.low < 1) r.low = 1; r.high = (int)(p75 + OUTLIER_BOUND * (p75 - p25) + .499); + if (r.low > r.high) { + ksprintf(msg, "[%s] fail to infer the insert size distribution: upper bound is smaller than max read length.\n", __func__); + free(isize); + r.failed = 1; + return r; + } ksprintf(msg, "[%s] (25, 50, 75) percentile: (%d, %d, %d)\n", __func__, p25, p50, p75); ksprintf(msg, "[%s] low and high boundaries for computing mean and std.dev: (%d, %d)\n", __func__, r.low, r.high); for (i = x = 0, r.avg = 0; i < k; ++i) diff --git a/ksw.c b/ksw.c index f5b24dc..c1ff5a6 100644 --- a/ksw.c +++ b/ksw.c @@ -202,10 +202,11 @@ end_loop16: r.score = gmax + q->shift < 255? gmax : 255; r.te = te; if (r.score != 255) { // get a->qe, the end of query match; find the 2nd best score - int max = -1, low, high, qlen = slen * 16; + int max = -1, tmp, low, high, qlen = slen * 16; uint8_t *t = (uint8_t*)Hmax; for (i = 0; i < qlen; ++i, ++t) if ((int)*t > max) max = *t, r.qe = i / 16 + i % 16 * slen; + else if ((int)*t == max && (tmp = i / 16 + i % 16 * slen) < r.qe) r.qe = tmp; //printf("%d,%d\n", max, gmax); if (b) { i = (r.score + q->max - 1) / q->max; @@ -303,10 +304,11 @@ end_loop8: } r.score = gmax; r.te = te; { - int max = -1, low, high, qlen = slen * 8; + int max = -1, tmp, low, high, qlen = slen * 8; uint16_t *t = (uint16_t*)Hmax; for (i = 0, r.qe = -1; i < qlen; ++i, ++t) if ((int)*t > max) max = *t, r.qe = i / 8 + i % 8 * slen; + else if ((int)*t == max && (tmp = i / 8 + i % 8 * slen) < r.qe) r.qe = tmp; if (b) { i = (r.score + q->max - 1) / q->max; low = te - i; high = te + i; diff --git a/main.c b/main.c index aa60e33..cbd695a 100644 --- a/main.c +++ b/main.c @@ -3,7 +3,7 @@ #include "utils.h" #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.7.3-r376-beta" +#define PACKAGE_VERSION "0.7.3-r382-beta" #endif int bwa_fa2pac(int argc, char *argv[]);