r1285: better --spsc support

This commit is contained in:
Heng Li
2025-05-26 15:47:14 -04:00
parent 005c9a1f6b
commit 3187782b1a
4 changed files with 29 additions and 7 deletions
+8 -1
View File
@@ -967,7 +967,7 @@ typedef struct mm_idx_spsc_s {
uint64_t *a; // pos<<56 | score<<1 | acceptor uint64_t *a; // pos<<56 | score<<1 | acceptor
} mm_idx_spsc_t; } mm_idx_spsc_t;
int32_t mm_idx_spsc_read(mm_idx_t *idx, const char *fn, int32_t max_sc) int32_t mm_idx_spsc_read2(mm_idx_t *idx, const char *fn, int32_t max_sc, float scale)
{ {
gzFile fp; gzFile fp;
kstring_t str = {0,0,0}; kstring_t str = {0,0,0};
@@ -1007,6 +1007,8 @@ int32_t mm_idx_spsc_read(mm_idx_t *idx, const char *fn, int32_t max_sc)
} }
} }
if (i < 4) continue; // not enough fields if (i < 4) continue; // not enough fields
if (scale > 0.0f && scale < 1.0f)
score = score > 0.0f? (int)(score * scale + .499) : (int)(score * scale - .499);
if (score > max_sc) score = max_sc; if (score > max_sc) score = max_sc;
if (score < -max_sc) score = -max_sc; if (score < -max_sc) score = -max_sc;
cid = mm_idx_name2id(idx, name); cid = mm_idx_name2id(idx, name);
@@ -1030,6 +1032,11 @@ int32_t mm_idx_spsc_read(mm_idx_t *idx, const char *fn, int32_t max_sc)
return 0; return 0;
} }
int32_t mm_idx_spsc_read(mm_idx_t *idx, const char *fn, int32_t max_sc)
{
return mm_idx_spsc_read2(idx, fn, max_sc, 1.0f);
}
static int32_t mm_idx_find_intv(int32_t n, const uint64_t *a, int64_t x) static int32_t mm_idx_find_intv(int32_t n, const uint64_t *a, int64_t x)
{ {
int32_t s = 0, e = n; int32_t s = 0, e = n;
+6 -2
View File
@@ -85,6 +85,8 @@ static ko_longopt_t long_options[] = {
{ "jump-min-match", ko_required_argument, 360 }, { "jump-min-match", ko_required_argument, 360 },
{ "write-junc", ko_no_argument, 361 }, { "write-junc", ko_no_argument, 361 },
{ "pass1", ko_required_argument, 362 }, { "pass1", ko_required_argument, 362 },
{ "spsc-scale", ko_required_argument, 363 },
{ "spsc0", ko_required_argument, 364 },
{ "dbg-seed-occ", ko_no_argument, 501 }, { "dbg-seed-occ", ko_no_argument, 501 },
{ "help", ko_no_argument, 'h' }, { "help", ko_no_argument, 'h' },
{ "max-intron-len", ko_required_argument, 'G' }, { "max-intron-len", ko_required_argument, 'G' },
@@ -134,6 +136,7 @@ int main(int argc, char *argv[])
mm_mapopt_t opt; mm_mapopt_t opt;
mm_idxopt_t ipt; mm_idxopt_t ipt;
int i, c, n_threads = 3, n_parts, old_best_n = -1; int i, c, n_threads = 3, n_parts, old_best_n = -1;
float spsc_scale = 0.7f;
char *fnw = 0, *rg = 0, *fn_bed_junc = 0, *fn_bed_jump = 0, *fn_bed_pass1 = 0, *fn_spsc = 0, *s, *alt_list = 0; char *fnw = 0, *rg = 0, *fn_bed_junc = 0, *fn_bed_jump = 0, *fn_bed_pass1 = 0, *fn_spsc = 0, *s, *alt_list = 0;
FILE *fp_help = stderr; FILE *fp_help = stderr;
mm_idx_reader_t *idx_rdr; mm_idx_reader_t *idx_rdr;
@@ -240,7 +243,6 @@ int main(int argc, char *argv[])
else if (c == 338) opt.max_qlen = mm_parse_num(o.arg); // --max-qlen else if (c == 338) opt.max_qlen = mm_parse_num(o.arg); // --max-qlen
else if (c == 340) fn_bed_junc = o.arg; // --junc-bed else if (c == 340) fn_bed_junc = o.arg; // --junc-bed
else if (c == 341) opt.junc_bonus = atoi(o.arg); // --junc-bonus else if (c == 341) opt.junc_bonus = atoi(o.arg); // --junc-bonus
else if (c == 358) opt.junc_pen = atoi(o.arg); // --junc-pen
else if (c == 342) opt.flag |= MM_F_SAM_HIT_ONLY; // --sam-hit-only else if (c == 342) opt.flag |= MM_F_SAM_HIT_ONLY; // --sam-hit-only
else if (c == 343) opt.chain_gap_scale = atof(o.arg); // --chain-gap-scale else if (c == 343) opt.chain_gap_scale = atof(o.arg); // --chain-gap-scale
else if (c == 351) opt.chain_skip_scale = atof(o.arg); // --chain-skip-scale else if (c == 351) opt.chain_skip_scale = atof(o.arg); // --chain-skip-scale
@@ -260,6 +262,8 @@ int main(int argc, char *argv[])
else if (c == 361) opt.flag |= MM_F_OUT_JUNC | MM_F_CIGAR; // --write-junc else if (c == 361) opt.flag |= MM_F_OUT_JUNC | MM_F_CIGAR; // --write-junc
else if (c == 362) fn_bed_pass1 = o.arg; // --jump-pass1 else if (c == 362) fn_bed_pass1 = o.arg; // --jump-pass1
else if (c == 501) mm_dbg_flag |= MM_DBG_SEED_FREQ; // --dbg-seed-occ else if (c == 501) mm_dbg_flag |= MM_DBG_SEED_FREQ; // --dbg-seed-occ
else if (c == 363) spsc_scale = atof(o.arg); // --spsc-scale
else if (c == 358 || c == 364) opt.junc_pen = atoi(o.arg); // --junc-pen or --spsc0
else if (c == 330) { else if (c == 330) {
fprintf(stderr, "[WARNING] \033[1;31m --lj-min-ratio has been deprecated.\033[0m\n"); fprintf(stderr, "[WARNING] \033[1;31m --lj-min-ratio has been deprecated.\033[0m\n");
} else if (c == 313) { // --sr } else if (c == 313) { // --sr
@@ -476,7 +480,7 @@ int main(int argc, char *argv[])
fprintf(stderr, "[WARNING] failed to load the pass-1 jump BED file\n"); fprintf(stderr, "[WARNING] failed to load the pass-1 jump BED file\n");
} }
if (fn_spsc) { if (fn_spsc) {
mm_idx_spsc_read(mi, fn_spsc, mm_max_spsc_bonus(&opt)); mm_idx_spsc_read2(mi, fn_spsc, mm_max_spsc_bonus(&opt), spsc_scale);
if (mi->spsc == 0 && mm_verbose >= 2) if (mi->spsc == 0 && mm_verbose >= 2)
fprintf(stderr, "[WARNING] failed to load the splice score file\n"); fprintf(stderr, "[WARNING] failed to load the splice score file\n");
} }
+4 -2
View File
@@ -5,7 +5,7 @@
#include <stdio.h> #include <stdio.h>
#include <sys/types.h> #include <sys/types.h>
#define MM_VERSION "2.29-r1283" #define MM_VERSION "2.29-r1285-dirty"
#define MM_F_NO_DIAG (0x001LL) // no exact diagonal hit #define MM_F_NO_DIAG (0x001LL) // no exact diagonal hit
#define MM_F_NO_DUAL (0x002LL) // skip pairs where query name is lexicographically larger than target name #define MM_F_NO_DUAL (0x002LL) // skip pairs where query name is lexicographically larger than target name
@@ -163,7 +163,8 @@ typedef struct {
int transition; // transition mismatch score (A:G, C:T) int transition; // transition mismatch score (A:G, C:T)
int sc_ambi; // score when one or both bases are "N" int sc_ambi; // score when one or both bases are "N"
int noncan; // cost of non-canonical splicing sites int noncan; // cost of non-canonical splicing sites
int junc_bonus, junc_pen; int junc_bonus; // bonus for a splice site in annotation
int junc_pen; // penalty for GT- or -AG not scored in --spsc
int zdrop, zdrop_inv; // break alignment if alignment score drops too fast along the diagonal int zdrop, zdrop_inv; // break alignment if alignment score drops too fast along the diagonal
int end_bonus; int end_bonus;
int min_dp_max; // drop an alignment if the score of the max scoring segment is below this threshold int min_dp_max; // drop an alignment if the score of the max scoring segment is below this threshold
@@ -420,6 +421,7 @@ int mm_idx_bed_junc(const mm_idx_t *mi, int32_t ctg, int32_t st, int32_t en, uin
int mm_max_spsc_bonus(const mm_mapopt_t *mo); int mm_max_spsc_bonus(const mm_mapopt_t *mo);
int32_t mm_idx_spsc_read(mm_idx_t *idx, const char *fn, int32_t max_sc); int32_t mm_idx_spsc_read(mm_idx_t *idx, const char *fn, int32_t max_sc);
int32_t mm_idx_spsc_read2(mm_idx_t *idx, const char *fn, int32_t max_sc, float scale);
int64_t mm_idx_spsc_get(const mm_idx_t *db, int32_t cid, int64_t st0, int64_t en0, int32_t rev, uint8_t *sc); int64_t mm_idx_spsc_get(const mm_idx_t *db, int32_t cid, int64_t st0, int64_t en0, int32_t rev, uint8_t *sc);
// deprecated APIs for backward compatibility // deprecated APIs for backward compatibility
+11 -2
View File
@@ -443,14 +443,23 @@ line corresponds to a donor site and `A' for an acceptor site.
A positive score suggests the junction is preferred and a negative score A positive score suggests the junction is preferred and a negative score
suggests the junction is not preferred. suggests the junction is not preferred.
.TP .TP
.BR --junc-pen \ INT .BR --spsc0 \ INT
Penalty for a position not in FILE specified by Penalty for positions not in
.I FILE
specified by
.B --spsc .B --spsc
[5]. Effective with [5]. Effective with
.B --spsc .B --spsc
but not but not
.BR --junc-bed . .BR --junc-bed .
.TP .TP
.BR --spsc-scale \ FLOAT
Scale splice scores in
.B --spsc
by
.IR FLOAT
rounded to the nearest integer [0.7].
.TP
.BR --junc-bed \ FILE .BR --junc-bed \ FILE
Junctions to prefer during base alignment []. Junctions to prefer during base alignment [].
Same format as Same format as