zzh
/
fast-bwa
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

code backup; to upgrade ksw.{c,h}

This commit is contained in:
Heng Li 2013-02-12 16:15:26 -05:00
parent 325ba8213b
commit 604e3d8da1
3 changed files with 21 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
bwamem.c
View File

@ -13,8 +13,6 @@
#include "kvec.h"
#include "ksort.h"
#define MAPQ_COEF 40.
int mem_verbose = 3; // 1: error only; 2: error+warning; 3: message+error+warning; >=4: debugging
void mem_fill_scmat(int a, int b, int8_t mat[25])
@ -583,7 +581,7 @@ static inline int approx_mapq_se(const mem_opt_t *opt, const mem_alnreg_t *a)
double identity;
sub = a->csub > sub? a->csub : sub;
l = a->qe - a->qb > a->re - a->rb? a->qe - a->qb : a->re - a->rb;
mapq = a->score? (int)(MAPQ_COEF * (1. - (double)sub / a->score) * log(a->seedcov) + .499) : 0;
mapq = a->score? (int)(MEM_MAPQ_COEF * (1. - (double)sub / a->score) * log(a->seedcov) + .499) : 0;
identity = 1. - (double)(l * opt->a - a->score) / (opt->a + opt->b) / l;
mapq = identity < 0.95? (int)(mapq * identity * identity + .499) : mapq;
if (a->sub_n) mapq -= (int)(4.343 * log(a->sub_n) + .499);

2
bwamem.h
View File

@ -5,6 +5,8 @@
#include "bntseq.h"
#include "utils.h"
#define MEM_MAPQ_COEF 40.0
struct __smem_i;
typedef struct __smem_i smem_i;

29
bwamem_pair.c
View File

@ -97,11 +97,12 @@ void mem_pestat(const mem_opt_t *opt, int64_t l_pac, int n, const mem_alnreg_v *
}
}
void mem_pair(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, const mem_pestat_t pes[4], bseq1_t s[2], mem_alnreg_v a[2], bwahit_t h[2])
int mem_pair(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, const mem_pestat_t pes[4], bseq1_t s[2], mem_alnreg_v a[2], bwahit_t h[2])
{
extern void mem_alnreg2hit(const mem_alnreg_t *a, bwahit_t *h);
pair64_v v;
pair64_t o, subo; // score<<32 | raw_score<<8 | hash
int r, i, y[4]; // y[] keeps the last hit
int r, i, k, y[4]; // y[] keeps the last hit
kv_init(v);
for (r = 0; r < 2; ++r) { // loop through read number
for (i = 0; i < a[r].n; ++i) {
@ -117,7 +118,7 @@ void mem_pair(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, con
o.x = o.y = subo.x = subo.y = 0;
for (i = 0; i < v.n; ++i) {
for (r = 0; r < 2; ++r) { // loop through direction
int dir = r<<1 | (v.a[i].y>>1&1), which, k;
int dir = r<<1 | (v.a[i].y>>1&1), which;
if (pes[dir].failed) continue; // invalid orientation
which = r<<1 | ((v.a[i].y&1)^1);
if (y[which] < 0) continue; // no previous hits
@ -133,7 +134,7 @@ void mem_pair(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, con
raw_score = (v.a[i].y>>32) + (v.a[i].y>>32);
if (raw_score + 20 * opt->a < (subo.x>>8&0xffffff)) continue; // skip the following if the score is too small
ns = (dist - pes[dir].avg) / pes[dir].std;
score = (int)(23. * raw_score / (opt->a + opt->b) - 4.343 * log(.5 * erfc(fabs(ns) * M_SQRT1_2)) + .499);
score = (int)(raw_score - 4.343 / 23. * (opt->a + opt->b) * log(erfc(fabs(ns) * M_SQRT1_2)) + .499);
pair = (uint64_t)k<<32 | i;
x = (uint64_t)score<<32 | (int64_t)raw_score<<8 | (hash_64(pair)&0xff);
if (x > o.x) subo = o, o.x = x, o.y = pair;
@ -142,7 +143,13 @@ void mem_pair(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, con
}
y[v.a[i].y&3] = i;
}
if (o.x > 0) {
i = o.y >> 32; k = o.y << 32 >> 32;
mem_alnreg2hit(&a[v.a[i].y&1].a[v.a[i].y<<32>>34], &h[v.a[i].y&1]);
mem_alnreg2hit(&a[v.a[k].y&1].a[v.a[k].y<<32>>34], &h[v.a[k].y&1]);
}
free(v.a);
return o.x == 0? -1 : 0;
}
void mem_sam_pe(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, const mem_pestat_t pes[4], bseq1_t s[2], mem_alnreg_v a[2])
@ -150,11 +157,11 @@ void mem_sam_pe(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, c
kstring_t str;
bwahit_t h[2];
str.l = str.m = 0; str.s = 0;
mem_pair(opt, bns, pac, pes, s, a, h);
/*
bwa_hit2sam(&str, opt->mat, opt->q, opt->r, opt->w, bns, pac, &s[0], &h[0], opt->is_hard);
s[0].sam = strdup(str.s); str.l = 0;
bwa_hit2sam(&str, opt->mat, opt->q, opt->r, opt->w, bns, pac, &s[1], &h[1], opt->is_hard);
s[1].sam = str.s;
*/
if (mem_pair(opt, bns, pac, pes, s, a, h) == 0) { // successful
bwa_hit2sam(&str, opt->mat, opt->q, opt->r, opt->w, bns, pac, &s[0], &h[0], opt->is_hard);
s[0].sam = strdup(str.s); str.l = 0;
bwa_hit2sam(&str, opt->mat, opt->q, opt->r, opt->w, bns, pac, &s[1], &h[1], opt->is_hard);
s[1].sam = str.s;
} else {
}
}