fork/minimap2
1
0
Fork 0
mirror of https://github.com/lh3/minimap2.git synced 2026-09-15 13:07:55 +08:00
Code Issues Packages Projects Releases Wiki Activity

r1266: fine tune mapq for short RNA-seq reads

Browse Source
This commit is contained in:
Heng Li
2025-04-07 00:57:39 -04:00
parent 75619c7b51
commit 3542a3d153
4 changed files with 19 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files

18
hit.c
View File

@@ -418,16 +418,19 @@ static void mm_set_inv_mapq(void *km, int n_regs, mm_reg1_t *regs)
kfree(km, aux);
}
void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr)
void mm_set_mapq2(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr, int is_splice)
{
static const float q_coef = 40.0f;
int64_t sum_sc = 0;
float uniq_ratio;
int i;
int i, n_2nd_splice = 0;
if (n_regs == 0) return;
for (i = 0; i < n_regs; ++i)
for (i = 0; i < n_regs; ++i) {
if (regs[i].parent == regs[i].id)
sum_sc += regs[i].score;
else if (regs[i].is_spliced)
++n_2nd_splice;
}
uniq_ratio = (float)sum_sc / (sum_sc + rep_len);
for (i = 0; i < n_regs; ++i) {
mm_reg1_t *r = &regs[i];
@@ -440,13 +443,18 @@ void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int ma
pen_cm = pen_s1 < pen_cm? pen_s1 : pen_cm;
subsc = r->subsc > min_chain_sc? r->subsc : min_chain_sc;
if (r->p && r->p->dp_max2 > 0 && r->p->dp_max > 0) {
float identity = (float)r->mlen / r->blen;
float x = (float)r->p->dp_max2 * subsc / r->p->dp_max / r->score0;
float x, identity = (float)r->mlen / r->blen;
if (is_sr && is_splice)
x = (float)r->p->dp_max2 / r->p->dp_max; // ignore chaining score; for short RNA-seq reads, unspliced chaining score tends to be higher
else
x = (float)r->p->dp_max2 * subsc / r->p->dp_max / r->score0;
mapq = (int)(identity * pen_cm * q_coef * (1.0f - x * x) * logf((float)r->p->dp_max / match_sc));
if (!is_sr) {
int mapq_alt = (int)(6.02f * identity * identity * (r->p->dp_max - r->p->dp_max2) / match_sc + .499f); // BWA-MEM like mapQ, mostly for short reads
mapq = mapq < mapq_alt? mapq : mapq_alt; // in case the long-read heuristic fails
}
if (is_splice && is_sr && r->is_spliced && n_2nd_splice == 0)
mapq += 10;
} else {
float x = (float)subsc / r->score0;
if (r->p) {

8
map.c
View File

@@ -227,7 +227,7 @@ static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *k
void mm_map_frag_core(const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, int *n_regs, mm_reg1_t **regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname)
{
int i, j, rep_len, qlen_sum, n_regs0, n_mini_pos;
int max_chain_gap_qry, max_chain_gap_ref, is_splice = !!(opt->flag & MM_F_SPLICE), is_sr = !!(opt->flag & MM_F_SR);
int max_chain_gap_qry, max_chain_gap_ref, is_splice = !!(opt->flag & MM_F_SPLICE), is_sr = !!(opt->flag & MM_F_SR), is_sr_rna = !!(opt->flag & MM_F_SR_RNA);
uint32_t hash;
int64_t n_a;
uint64_t *u, *mini_pos;
@@ -338,7 +338,7 @@ void mm_map_frag_core(const mm_idx_t *mi, int n_segs, const int *qlens, const ch
if (n_segs == 1) { // uni-segment
regs0 = align_regs(opt, mi, b->km, qlens[0], seqs[0], &n_regs0, regs0, a);
regs0 = (mm_reg1_t*)realloc(regs0, sizeof(*regs0) * n_regs0);
mm_set_mapq(b->km, n_regs0, regs0, opt->min_chain_score, opt->a, rep_len, is_sr);
mm_set_mapq2(b->km, n_regs0, regs0, opt->min_chain_score, opt->a, rep_len, is_sr || is_sr_rna, is_splice);
n_regs[0] = n_regs0, regs[0] = regs0;
} else { // multi-segment
mm_seg_t *seg;
@@ -347,7 +347,7 @@ void mm_map_frag_core(const mm_idx_t *mi, int n_segs, const int *qlens, const ch
for (i = 0; i < n_segs; ++i) {
mm_set_parent(b->km, opt->mask_level, opt->mask_len, n_regs[i], regs[i], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop); // update mm_reg1_t::parent
regs[i] = align_regs(opt, mi, b->km, qlens[i], seqs[i], &n_regs[i], regs[i], seg[i].a);
mm_set_mapq(b->km, n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len, is_sr);
mm_set_mapq2(b->km, n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len, is_sr || is_sr_rna, is_splice);
}
mm_seg_free(b->km, n_segs, seg);
if (n_segs == 2 && opt->pe_ori >= 0 && (opt->flag&MM_F_CIGAR))
@@ -525,7 +525,7 @@ static void merge_hits(step_t *s)
mm_select_sub(km, opt->pri_ratio, s->p->mi->k*2, opt->best_n, 0, opt->max_gap * 0.8, &s->n_reg[k], s->reg[k]);
mm_set_sam_pri(s->n_reg[k], s->reg[k]);
}
mm_set_mapq(km, s->n_reg[k], s->reg[k], opt->min_chain_score, opt->a, rep_len, !!(opt->flag & MM_F_SR));
mm_set_mapq2(km, s->n_reg[k], s->reg[k], opt->min_chain_score, opt->a, rep_len, !!(opt->flag & (MM_F_SR|MM_F_SR_RNA)), !!(opt->flag & MM_F_SPLICE));
}
if (s->n_seg[f] == 2 && opt->pe_ori >= 0 && (opt->flag&MM_F_CIGAR))
mm_pair(km, frag_gap_part[0], opt->pe_bonus, opt->a * 2 + opt->b, opt->a, qlens, &s->n_reg[k0], &s->reg[k0]);

2
minimap.h
View File

@@ -5,7 +5,7 @@
#include <stdio.h>
#include <sys/types.h>
#define MM_VERSION "2.28-r1265-dirty"
#define MM_VERSION "2.28-r1266-dirty"
#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

2
mmpriv.h
View File

@@ -103,7 +103,7 @@ void mm_select_sub_multi(void *km, float pri_ratio, float pri1, float pri2, int
int mm_filter_strand_retained(int n_regs, mm_reg1_t *r);
void mm_filter_regs(const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs);
void mm_hit_sort(void *km, int *n_regs, mm_reg1_t *r, float alt_diff_frac);
void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr);
void mm_set_mapq2(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr, int is_splice);
void mm_update_dp_max(int qlen, int n_regs, mm_reg1_t *regs, float frac, int a, int b);
void mm_jump_split(void *km, const mm_idx_t *mi, const mm_mapopt_t *opt, int32_t qlen, const uint8_t *qseq, mm_reg1_t *r, int32_t ts_strand);