r763: fine control long join flank len (#128)

This commit is contained in:
Heng Li 2018-03-29 14:16:58 -04:00
parent 2d7ec75d50
commit ee4cd089f7
5 changed files with 10 additions and 4 deletions

7
hit.c
View File

@ -304,7 +304,7 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs_, mm_r
for (i = n_aux - 1; i >= 1; --i) {
mm_reg1_t *r0 = &regs[(int32_t)aux[i-1]], *r1 = &regs[(int32_t)aux[i]];
mm128_t *a0e, *a1s;
int max_gap, min_gap, sc_thres;
int max_gap, min_gap, sc_thres, min_flank_len;
// test
if (r0->as + r0->cnt != r1->as) continue; // not adjacent in a[]
@ -318,8 +318,9 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs_, mm_r
if (max_gap > opt->max_join_long || min_gap > opt->max_join_short) continue;
sc_thres = (int)((float)opt->min_join_flank_sc / opt->max_join_long * max_gap + .499);
if (r0->score < sc_thres || r1->score < sc_thres) continue; // require good flanking chains
if (r0->re - r0->rs < max_gap>>1 || r0->qe - r0->qs < max_gap>>1) continue; // require enough flanking length
if (r1->re - r1->rs < max_gap>>1 || r1->qe - r1->qs < max_gap>>1) continue;
min_flank_len = (int)(max_gap * opt->min_join_flank_ratio);
if (r0->re - r0->rs < min_flank_len || r0->qe - r0->qs < min_flank_len) continue; // require enough flanking length
if (r1->re - r1->rs < min_flank_len || r1->qe - r1->qs < min_flank_len) continue;
// all conditions satisfied; join
a[r1->as].y |= MM_SEED_LONG_JOIN;

4
main.c
View File

@ -10,7 +10,7 @@
#include "getopt.h"
#endif
#define MM_VERSION "2.10-r761"
#define MM_VERSION "2.10-r763-dirty"
#ifdef __linux__
#include <sys/resource.h>
@ -57,6 +57,7 @@ static struct option long_options[] = {
{ "max-clip-ratio", required_argument, 0, 0 }, // 27
{ "min-occ-floor", required_argument, 0, 0 }, // 28
{ "MD", no_argument, 0, 0 }, // 29
{ "lj-min-ratio", required_argument, 0, 0 }, // 30
{ "help", no_argument, 0, 'h' },
{ "max-intron-len", required_argument, 0, 'G' },
{ "version", no_argument, 0, 'V' },
@ -173,6 +174,7 @@ int main(int argc, char *argv[])
else if (c == 0 && long_idx ==27) opt.max_clip_ratio = atof(optarg); // --max-clip-ratio
else if (c == 0 && long_idx ==28) opt.min_mid_occ = atoi(optarg); // --min-occ-floor
else if (c == 0 && long_idx ==29) opt.flag |= MM_F_OUT_MD; // --MD
else if (c == 0 && long_idx ==30) opt.min_join_flank_ratio = atof(optarg); // --lj-min-ratio
else if (c == 0 && long_idx == 14) { // --frag
yes_or_no(&opt, MM_F_FRAG_MODE, long_idx, optarg, 1);
} else if (c == 0 && long_idx == 15) { // --secondary

View File

@ -115,6 +115,7 @@ typedef struct {
int max_join_long, max_join_short;
int min_join_flank_sc;
float min_join_flank_ratio;
int a, b, q, e, q2, e2; // matching score, mismatch, gap-open and gap-ext penalties
int noncan; // cost of non-canonical splicing sites

View File

@ -31,6 +31,7 @@ void mm_mapopt_init(mm_mapopt_t *opt)
opt->max_join_long = 20000;
opt->max_join_short = 2000;
opt->min_join_flank_sc = 1000;
opt->min_join_flank_ratio = 0.5f;
opt->a = 2, opt->b = 4, opt->q = 4, opt->e = 2, opt->q2 = 24, opt->e2 = 1;
opt->zdrop = 400, opt->zdrop_inv = 200;

View File

@ -24,6 +24,7 @@ cdef extern from "minimap.h":
int best_n
int max_join_long, max_join_short
int min_join_flank_sc
float min_join_flank_ratio;
int a, b, q, e, q2, e2
int noncan
int zdrop, zdrop_inv