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

r1276: skip unnecessary reverse spliced alignment

Browse Source 
for short-read RNA-seq only
This commit is contained in:
Heng Li
2025-04-14 14:57:23 -04:00
parent ff312a2957
commit 94d171b01e
2 changed files with 25 additions and 20 deletions
Download Patch File Download Diff File Expand all files Collapse all files

43
align.c
View File

@@ -1012,30 +1012,35 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
n_a = mm_squeeze_a(km, n_regs, regs, a);
memset(&ez, 0, sizeof(ksw_extz_t));
for (i = 0; i < n_regs; ++i) {
mm_reg1_t r2;
mm_reg1_t r2; // only used for inversion
if ((opt->flag&MM_F_SPLICE) && (opt->flag&MM_F_SPLICE_FOR) && (opt->flag&MM_F_SPLICE_REV)) { // then do two rounds of alignments for both strands
mm_reg1_t s[2], s2[2], *r;
int which, trans_strand;
s[0] = s[1] = regs[i];
mm_align1(km, opt, mi, qlen, qseq0, &s[0], &s2[0], n_a, a, &ez, MM_F_SPLICE_FOR);
mm_align1(km, opt, mi, qlen, qseq0, &s[1], &s2[1], n_a, a, &ez, MM_F_SPLICE_REV);
if (s[0].p->dp_score > s[1].p->dp_score) which = 0, trans_strand = 1;
else if (s[0].p->dp_score < s[1].p->dp_score) which = 1, trans_strand = 2;
else trans_strand = 3, which = (qlen + s[0].p->dp_score) & 1; // randomly choose a strand, effectively
if (which == 0) {
mm_align1(km, opt, mi, qlen, qseq0, &s[0], &s2[0], n_a, a, &ez, MM_F_SPLICE_FOR); // assume the transcript is on the + strand of the genome
if ((opt->flag&MM_F_SR_RNA) && regs[i].qe - regs[i].qs == regs[i].re - regs[i].rs && s[0].qe - s[0].qs == s[0].re - s[0].rs && s[0].qs == 0 && s[0].qe == qlen) {
regs[i] = s[0], r2 = s2[0];
free(s[1].p);
regs[i].p->trans_strand = 0;
} else {
regs[i] = s[1], r2 = s2[1];
free(s[0].p);
}
r = &regs[i];
r->p->trans_strand = trans_strand;
if (r->is_spliced) {
if (trans_strand == 1 || trans_strand == 2) // this is an *approximate* way to tell if there are splice signals.
r->p->dp_max += (opt->a + opt->b) + ((opt->a + opt->b) >> 1);
else if (trans_strand == 3)
r->p->dp_max -= opt->a + opt->b;
int which, trans_strand;
mm_align1(km, opt, mi, qlen, qseq0, &s[1], &s2[1], n_a, a, &ez, MM_F_SPLICE_REV); // assume the transcript on the - strand
if (s[0].p->dp_score > s[1].p->dp_score) which = 0, trans_strand = 1;
else if (s[0].p->dp_score < s[1].p->dp_score) which = 1, trans_strand = 2;
else trans_strand = 3, which = (qlen + s[0].p->dp_score) & 1; // randomly choose a strand, effectively
if (which == 0) {
regs[i] = s[0], r2 = s2[0];
free(s[1].p);
} else {
regs[i] = s[1], r2 = s2[1];
free(s[0].p);
}
r = &regs[i];
r->p->trans_strand = trans_strand;
if (r->is_spliced) {
if (trans_strand == 1 || trans_strand == 2) // this is an *approximate* way to tell if there are splice signals.
r->p->dp_max += (opt->a + opt->b) + ((opt->a + opt->b) >> 1);
else if (trans_strand == 3)
r->p->dp_max -= opt->a + opt->b;
}
}
} else { // one round of alignment
mm_align1(km, opt, mi, qlen, qseq0, &regs[i], &r2, n_a, a, &ez, opt->flag);

2
minimap.h
View File

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