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

code backup

This commit is contained in:
Heng Li 2013-02-12 09:22:47 -05:00
parent dcb190069a
commit 13288e2dcd
3 changed files with 33 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
bwamem.c
View File

@ -493,7 +493,7 @@ ret_gen_cigar:
void bwa_hit2sam(kstring_t *str, const int8_t mat[25], int q, int r, int w, const bntseq_t *bns, const uint8_t *pac, bseq1_t *s, bwahit_t *p, int is_hard)
{
int score, n_cigar, is_rev, nn, rid, mid, is_unmapped = 0;
int score, n_cigar, is_rev = 0, nn, rid, mid, is_unmapped = 0;
uint32_t *cigar = 0;
int64_t pos;
@ -652,6 +652,7 @@ static void *worker1(void *data)
static void *worker2(void *data)
{
extern 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]);
worker_t *w = (worker_t*)data;
int i;
if (!w->opt->is_pe) {
@ -661,6 +662,7 @@ static void *worker2(void *data)
}
} else {
for (i = 0; i < w->n>>1; i += w->step) { // not implemented yet
mem_sam_pe(w->opt, w->bns, w->pac, w->pes, &w->seqs[i<<1], &w->regs[i<<1]);
free(w->regs[i<<1|0].a); free(w->regs[i<<1|1].a);
}
}

33
bwamem_pair.c
View File

@ -103,17 +103,40 @@ 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])
{
vec64_t v;
int r, i;
int r, i, y[4]; // y[] keeps the last hit
kv_init(v);
for (r = 0; r < 2; ++r) {
for (i = 0; i < a[r].n; ++i) {
int64_t pos;
uint64_t key;
mem_alnreg_t *e = &a[r].a[i];
pos = (e->rb < bns->l_pac? e->rb<<1 : ((bns->l_pac<<1) - 1 - e->rb)<<1 | 1)<<1 | r;
kv_push(uint64_t, v, pos);
key = ((e->rb < bns->l_pac? e->rb<<1 : ((bns->l_pac<<1) - 1 - e->rb)<<1 | 1)<<1 | r) << 30 | e->score;
kv_push(uint64_t, v, key);
}
}
ks_introsort_uint64_t(v.n, v.a);
y[0] = y[1] = y[2] = y[3] = -1;
printf("**** %ld\n", v.n);
for (i = 0; i < v.n; ++i) {
printf("%lld\t%c\t%lld\t%lld\n", v.a[i]>>32, "+-"[v.a[i]>>31&1], v.a[i]>>30&1, v.a[i]<<34>>34);
for (r = 0; r < 2; ++r) {
int dir = r<<1 | (v.a[i]>>31&1), which, k;
if (pes[dir].failed) continue; // invalid orientation
which = r<<1 | ((v.a[i]>>30&1)^1);
if (y[which] < 0) continue; // no previous hits
for (k = y[which]; k >= 0; --k) { // TODO: this is a O(n^2) solution in the worst case; remember to check if this loop takes a lot of time (I doubt)
int dist;
double ns;
if ((v.a[k]>>30&3) != which) continue;
dist = (v.a[i]>>32) - (v.a[k]>>32);
printf("%d\t%d\t%d\n", r, which, dist);
if (dist > pes[dir].high) break;
if (dist < pes[dir].low) continue;
ns = (dist - pes[dir].avg) / pes[dir].std;
printf("%f\n", ns);
}
}
y[v.a[i]>>30&3] = i;
}
free(v.a);
}
@ -123,8 +146,10 @@ void mem_sam_pe(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, c
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;
*/
}

2
ksort.h
View File

@ -139,7 +139,7 @@ typedef struct {
tmp = *l; *l = l[i]; l[i] = tmp; ks_heapadjust_##name(0, i, l); \
} \
} \
inline void __ks_insertsort_##name(type_t *s, type_t *t) \
static inline void __ks_insertsort_##name(type_t *s, type_t *t) \
{ \
type_t *i, *j, swap_tmp; \
for (i = s + 1; i < t; ++i) \