r1263: ~5-10% performance improvement

Via larger batches and more short-read heuristics. Identical alignment on 2 million reads. Short DNA-seq read alignment may be improved in corner cases.
2026-09-15 13:07:55 +08:00 · 2025-04-06 20:48:55 -04:00
parent d43f356ef9
commit af094640e5
5 changed files with 36 additions and 26 deletions
--- a/align.c
+++ b/align.c
@@ -603,7 +603,7 @@ static inline void mm_get_junc(const mm_idx_t *mi, int32_t ctg, int32_t st, int3

 static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], mm_reg1_t *r, mm_reg1_t *r2, int n_a, mm128_t *a, ksw_extz_t *ez, int splice_flag)
 {
-	int is_sr = !!(opt->flag & MM_F_SR), is_splice = !!(opt->flag & MM_F_SPLICE);
+	int is_sr = !!(opt->flag & MM_F_SR), is_splice = !!(opt->flag & MM_F_SPLICE), is_sr_rna = (!!(opt->flag & MM_F_SR_RNA) && is_splice);
 	int32_t rid = a[r->as].x<<1>>33, rev = a[r->as].x>>63, as1, cnt1;
 	uint8_t *tseq, *qseq, *junc;
 	int32_t i, l, bw, bw_long, dropped = 0, ksw_flag = 0, rs0, re0, qs0, qe0;
@@ -773,14 +773,18 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 				mm_idx_getseq(mi, rid, rs, re, tseq);
 			}
 			mm_get_junc(mi, rid, rs, re, !!(ksw_flag&KSW_EZ_SPLICE_REV), junc);
-			if (is_sr) { // perform ungapped alignment
+			if (is_sr || (is_sr_rna && qe - qs == re - rs)) { // perform ungapped alignment
+				int32_t max_gapped_score = (qe - qs - 2) * opt->a - 2 * (opt->q + opt->e);
 				assert(qe - qs == re - rs);
 				ksw_reset_extz(ez);
 				for (j = 0, ez->score = 0; j < qe - qs; ++j) {
-					if (qseq[j] >= 4 || tseq[j] >= 4) ez->score += opt->e2;
+					if (qseq[j] >= 4 || tseq[j] >= 4) ez->score += opt->sc_ambi > 0? -opt->sc_ambi : opt->sc_ambi;
 					else ez->score += qseq[j] == tseq[j]? opt->a : -opt->b;
 				}
-				ez->cigar = ksw_push_cigar(km, &ez->n_cigar, &ez->m_cigar, ez->cigar, MM_CIGAR_MATCH, qe - qs);
+				if (ez->score > max_gapped_score)
+					ez->cigar = ksw_push_cigar(km, &ez->n_cigar, &ez->m_cigar, ez->cigar, MM_CIGAR_MATCH, qe - qs);
+				else
+					mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, junc, mat, bw1, -1, opt->zdrop, ksw_flag|KSW_EZ_APPROX_MAX, ez);
 			} else { // perform normal gapped alignment
 				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, junc, mat, bw1, -1, opt->zdrop, ksw_flag|KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop
 			}
--- a/main.c
+++ b/main.c
@@ -35,7 +35,7 @@ static ko_longopt_t long_options[] = {
 	{ "splice",         ko_no_argument,       310 },
 	{ "cost-non-gt-ag", ko_required_argument, 'C' },
 	{ "no-long-join",   ko_no_argument,       312 },
-	{ "sr",             ko_no_argument,       313 },
+	{ "sr",             ko_optional_argument, 313 },
 	{ "frag",           ko_required_argument, 314 },
 	{ "secondary",      ko_required_argument, 315 },
 	{ "cs",             ko_optional_argument, 316 },
@@ -218,7 +218,6 @@ int main(int argc, char *argv[])
 		else if (c == 309) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ, n_threads = 1; // --print-aln-seq
 		else if (c == 310) opt.flag |= MM_F_SPLICE; // --splice
 		else if (c == 312) opt.flag |= MM_F_NO_LJOIN; // --no-long-join
-		else if (c == 313) opt.flag |= MM_F_SR; // --sr
 		else if (c == 317) opt.end_bonus = atoi(o.arg); // --end-bonus
 		else if (c == 318) opt.flag |= MM_F_INDEPEND_SEG; // --no-pairing (deprecated)
 		else if (c == 320) ipt.flag |= MM_I_NO_SEQ; // --idx-no-seq
@@ -257,6 +256,17 @@ int main(int argc, char *argv[])
 		else if (c == 501) mm_dbg_flag |= MM_DBG_SEED_FREQ; // --dbg-seed-occ
 		else if (c == 330) {
 			fprintf(stderr, "[WARNING] \033[1;31m --lj-min-ratio has been deprecated.\033[0m\n");
+		} else if (c == 313) { // --sr
+			if (o.arg == 0 || strcmp(o.arg, "dna") == 0) {
+				opt.flag |= MM_F_SR;
+			} else if (strcmp(o.arg, "rna") == 0) {
+				opt.flag |= MM_F_SR_RNA;
+			} else if (strcmp(o.arg, "no") == 0) {
+				opt.flag &= ~(uint64_t)(MM_F_SR|MM_F_SR_RNA);
+			} else if (mm_verbose >= 2) {
+				opt.flag |= MM_F_SR;
+				fprintf(stderr, "[WARNING]\033[1;31m --sr only takes 'dna' or 'rna'. Invalid values are assumed to be 'dna'.\033[0m\n");
+			}
 		} else if (c == 314) { // --frag
 			yes_or_no(&opt, MM_F_FRAG_MODE, o.longidx, o.arg, 1);
 		} else if (c == 315) { // --secondary
--- a/minimap.h
+++ b/minimap.h
@@ -5,7 +5,7 @@
 #include <stdio.h>
 #include <sys/types.h>

-#define MM_VERSION "2.28-r1261-dirty"
+#define MM_VERSION "2.28-r1263-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
@@ -46,6 +46,7 @@
 #define MM_F_SECONDARY_SEQ (0x1000000000LL)	//output SEQ field for seqondary alignments using hard clipping
 #define MM_F_OUT_DS        (0x2000000000LL)
 #define MM_F_WEAK_PAIRING  (0x4000000000LL)
+#define MM_F_SR_RNA        (0x8000000000LL)

 #define MM_I_HPC          0x1
 #define MM_I_NO_SEQ       0x2
--- a/minimap2.1
+++ b/minimap2.1
@@ -297,11 +297,13 @@ maximum alignment gap is mostly controlled by
 .B --splice
 Enable the splice alignment mode.
 .TP
-.B --sr
-Enable short-read alignment heuristics. In the short-read mode, minimap2
-applies a second round of chaining with a higher minimizer occurrence threshold
-if no good chain is found. In addition, minimap2 attempts to patch gaps between
-seeds with ungapped alignment.
+.BR --sr [= no | dna | rna ]
+Enable short-read alignment heuristics [no]. If this option is used with no argument,
+.RB ` dna '
+is set. In the DNA short-read mode, minimap2 applies a second round of chaining
+with a higher minimizer occurrence threshold if no good chain is found. In
+addition, minimap2 attempts to patch gaps between seeds with ungapped
+alignment.
 .TP
 .BI --split-prefix \ STR
 Prefix to create temporary files. Typically used for a multi-part index.
@@ -520,20 +522,13 @@ Copy input FASTA/Q comments to output.
 .B -c
 Generate CIGAR. In PAF, the CIGAR is written to the `cg' custom tag.
 .TP
-.BI --cs[= STR ]
+.BR --cs [= short | long ]
 Output the
 .B cs
 tag.
-.I STR
-can be either
-.I short
-or
-.IR long .
-If no
-.I STR
-is given,
-.I short
-is assumed. [none]
+If no argument is given,
+.RB ` short '
+is set. [none]
 .TP
 .B --MD
 Output the MD tag (see the SAM spec).
@@ -689,7 +684,7 @@ Spliced alignment for accurate long RNA-seq reads such as PacBio iso-seq
 .B splice:sr
 Spliced alignment for short RNA-seq reads
 .RB ( -xsplice:hq
-.B --frag=yes -m25 -s40 -2K50m --heap-sort=yes --pairing=weak
+.B --frag=yes -m25 -s40 -2K100m --heap-sort=yes --pairing=weak --sr=rna
 .BR --secondary=no ).
 .TP
 .B sr
--- a/options.c
+++ b/options.c
@@ -179,13 +179,13 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		if (strcmp(preset, "splice:hq") == 0) {
 			mo->noncan = 5, mo->b = 4, mo->q = 6, mo->q2 = 24;
 		} else if (strcmp(preset, "splice:sr") == 0) {
-			mo->flag |= MM_F_NO_PRINT_2ND | MM_F_2_IO_THREADS | MM_F_HEAP_SORT | MM_F_FRAG_MODE | MM_F_WEAK_PAIRING;
+			mo->flag |= MM_F_NO_PRINT_2ND | MM_F_2_IO_THREADS | MM_F_HEAP_SORT | MM_F_FRAG_MODE | MM_F_WEAK_PAIRING | MM_F_SR_RNA;
 			mo->noncan = 5, mo->b = 4, mo->q = 6, mo->q2 = 24;
 			mo->min_chain_score = 25;
 			mo->min_dp_max = 40;
 			mo->pe_ori = 0<<1|1; // FR
 			mo->best_n = 10;
-			mo->mini_batch_size = 50000000;
+			mo->mini_batch_size = 100000000;
 		}
 	} else return -1;
 	return 0;