r1261: code clean up; renamed --jump-bed to -j

Also added --pairing to replace --no-pairing and --pe-ind-chain

README.md

@@ -20,8 +20,9 @@ cd minimap2 && make
 ./minimap2 -ax splice ref.fa rna-reads.fa > aln.sam       # spliced long reads (strand unknown)
 ./minimap2 -ax splice -uf -k14 ref.fa reads.fa > aln.sam  # noisy Nanopore direct RNA-seq
 ./minimap2 -ax splice:hq -uf ref.fa query.fa > aln.sam    # PacBio Kinnex/Iso-seq (RNA-seq)
-./minimap2 -ax splice --junc-bed anno.bed12 ref.fa query.fa > aln.sam  # use annotated junctions
+./minimap2 -ax splice --junc-bed=anno.bed12 ref.fa query.fa > aln.sam  # use annotated junctions
 ./minimap2 -ax splice:sr ref.fa r1.fq r2.fq > aln.sam     # short-read RNA-seq (r1236+; experimental)
+./minimap2 -ax splice:sr -j anno.bed12 ref.fa r1.fq r2.fq > aln.sam
 ./minimap2 -cx asm5 asm1.fa asm2.fa > aln.paf             # intra-species asm-to-asm alignment
 ./minimap2 -x ava-pb reads.fa reads.fa > overlaps.paf     # PacBio read overlap
 ./minimap2 -x ava-ont reads.fa reads.fa > overlaps.paf    # Nanopore read overlap
@@ -206,6 +207,10 @@ bonus score (tuned by `--junc-bonus`) if an aligned junction matches a junction
 in the annotation. Option `--junc-bed` also takes 5-column BED, including the
 strand field. In this case, each line indicates an oriented junction.
 
+**Note:** that `--junc-bed` is intended for long noisy RNA-seq reads only.
+Applying the option to short RNA-seq reads increase run time with little
+improvement to junction accuracy.
+
 #### <a name="long-overlap"></a>Find overlaps between long reads
 
 ```sh
@@ -234,8 +239,9 @@ mixed.
 #### <a name="short-rna-seq"></a>Map short RNA-seq reads (experimental & evolving)
 
 ```sh
-minimap2 -ax splice:sr ref.fa reads-se.fq > aln.sam     # single-end
+minimap2 -ax splice:sr ref.fa reads-se.fq.gz > aln.sam           # single-end
-minimap2 -ax splice:sr ref.fa r1.fq r2.fq > aln.sam     # paired-end
+minimap2 -ax splice:sr ref.fa r1.fq.gz r2.fq.gz > aln.sam        # paired-end
+minimap2 -ax splice:sr -j anno.bed ref.fa r1.fq r2.fq > aln.sam  # use annotation
 ```
 The new preset `splice:sr` was added between v2.28 and v2.29. It functions
 similarly to `sr` except that it performs spliced alignment. Note that this

main.c

@@ -40,7 +40,7 @@ static ko_longopt_t long_options[] = {
 	{ "secondary",      ko_required_argument, 315 },
 	{ "cs",             ko_optional_argument, 316 },
 	{ "end-bonus",      ko_required_argument, 317 },
-	{ "no-pairing",     ko_no_argument,       318 },
+	{ "no-pairing",     ko_no_argument,       318 }, // deprecated but reserved for backward compatibility
 	{ "splice-flank",   ko_required_argument, 319 },
 	{ "idx-no-seq",     ko_no_argument,       320 },
 	{ "end-seed-pen",   ko_required_argument, 321 },
@@ -81,8 +81,7 @@ static ko_longopt_t long_options[] = {
 	{ "rmq-inner",      ko_required_argument, 356 },
 	{ "spsc",           ko_required_argument, 357 },
 	{ "junc-pen",       ko_required_argument, 358 },
-	{ "pe-ind-chain",   ko_no_argument,       359 },
+	{ "pairing",        ko_required_argument, 359 },
-	{ "jump-bed",       ko_required_argument, 360 },
 	{ "dbg-seed-occ",   ko_no_argument,       501 },
 	{ "help",           ko_no_argument,       'h' },
 	{ "max-intron-len", ko_required_argument, 'G' },
@@ -127,7 +126,7 @@ static inline void yes_or_no(mm_mapopt_t *opt, int64_t flag, int long_idx, const
 
 int main(int argc, char *argv[])
 {
-	const char *opt_str = "2aSDw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:b:O:E:m:N:Qu:R:hF:LC:yYPo:e:U:J:";
+	const char *opt_str = "2aSDw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:b:O:E:m:N:Qu:R:hF:LC:yYPo:e:U:J:j:";
 	ketopt_t o = KETOPT_INIT;
 	mm_mapopt_t opt;
 	mm_idxopt_t ipt;
@@ -194,6 +193,7 @@ int main(int argc, char *argv[])
 		else if (c == 'R') rg = o.arg;
 		else if (c == 'h') fp_help = stdout;
 		else if (c == '2') opt.flag |= MM_F_2_IO_THREADS;
+		else if (c == 'j') jump_bed = o.arg;
 		else if (c == 'J') {
 			int t;
 			t = atoi(o.arg);
@@ -220,7 +220,7 @@ int main(int argc, char *argv[])
 		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
+		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
 		else if (c == 321) opt.anchor_ext_shift = atoi(o.arg); // --end-seed-pen
 		else if (c == 322) opt.flag |= MM_F_FOR_ONLY; // --for-only
@@ -254,8 +254,6 @@ int main(int argc, char *argv[])
 		else if (c == 355) opt.flag |= MM_F_OUT_DS; // --ds
 		else if (c == 356) opt.rmq_inner_dist = mm_parse_num(o.arg); // --rmq-inner
 		else if (c == 357) fn_spsc = o.arg; // --spsc
-		else if (c == 359) opt.flag |= MM_F_PE_IND; // --pe-ind-chain
-		else if (c == 360) jump_bed = o.arg; // --jump-bed
 		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");
@@ -283,6 +281,14 @@ int main(int argc, char *argv[])
 		} else if (c == 347) { // --rmq
 			if (o.arg) yes_or_no(&opt, MM_F_RMQ, o.longidx, o.arg, 1);
 			else opt.flag |= MM_F_RMQ;
+		} else if (c == 359) { // --pairing
+			if (strcmp(o.arg, "no") == 0) opt.flag |= MM_F_INDEPEND_SEG;
+			else if (strcmp(o.arg, "weak") == 0) opt.flag |= MM_F_WEAK_PAIRING, opt.flag &= ~(uint64_t)MM_F_INDEPEND_SEG;
+			else {
+				if (strcmp(o.arg, "strong") != 0 && mm_verbose >= 2)
+					fprintf(stderr, "[WARNING]\033[1;31m unrecognized argument for --pairing; assuming 'strong'.\033[0m\n");
+				opt.flag &= ~(uint64_t)(MM_F_INDEPEND_SEG|MM_F_WEAK_PAIRING);
+			}
 		} else if (c == 'S') {
 			opt.flag |= MM_F_OUT_CS | MM_F_CIGAR | MM_F_OUT_CS_LONG;
 			if (mm_verbose >= 2)
@@ -362,6 +368,7 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -s INT       minimal peak DP alignment score [%d]\n", opt.min_dp_max);
 		fprintf(fp_help, "    -u CHAR      how to find GT-AG. f:transcript strand, b:both strands, n:don't match GT-AG [n]\n");
 		fprintf(fp_help, "    -J INT       splice mode. 0: original minimap2 model; 1: miniprot model [1]\n");
+		fprintf(fp_help, "    -j FILE      junctions in BED12 to extend *short* RNA-seq alignment []\n");
 		fprintf(fp_help, "  Input/Output:\n");
 		fprintf(fp_help, "    -a           output in the SAM format (PAF by default)\n");
 		fprintf(fp_help, "    -o FILE      output alignments to FILE [stdout]\n");

map.c

@@ -379,7 +379,7 @@ void mm_map_frag_core(const mm_idx_t *mi, int n_segs, const int *qlens, const ch
 
 void mm_map_frag(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)
 {
-	if ((opt->flag & MM_F_PE_IND) && n_segs == 2 && opt->pe_ori >= 0 && (opt->flag&MM_F_CIGAR)) {
+	if ((opt->flag & MM_F_WEAK_PAIRING) && n_segs == 2 && opt->pe_ori >= 0 && (opt->flag&MM_F_CIGAR)) {
 		int i;
 		for (i = 0; i < n_segs; ++i)
 			mm_map_frag_core(mi, 1, &qlens[i], &seqs[i], &n_regs[i], &regs[i], b, opt, qname);

minimap.h

@@ -5,7 +5,7 @@
 #include <stdio.h>
 #include <sys/types.h>
 
-#define MM_VERSION "2.28-r1259-dirty"
+#define MM_VERSION "2.28-r1261-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
@@ -45,7 +45,7 @@
 #define MM_F_SPLICE_OLD    (0x800000000LL)
 #define MM_F_SECONDARY_SEQ (0x1000000000LL)	//output SEQ field for seqondary alignments using hard clipping
 #define MM_F_OUT_DS        (0x2000000000LL)
-#define MM_F_PE_IND        (0x4000000000LL)
+#define MM_F_WEAK_PAIRING  (0x4000000000LL)
 
 #define MM_I_HPC          0x1
 #define MM_I_NO_SEQ       0x2

minimap2.1

@@ -79,19 +79,6 @@ Minimizer k-mer length [15]
 .BI -w \ INT
 Minimizer window size [10]. A minimizer is the smallest k-mer
 in a window of w consecutive k-mers.
-.TP
-.BI -j \ INT
-Syncmer submer size [10]. Option
-.B -j
-and
-.B -w
-will override each: if
-.B -w
-is applied after
-.BR -j ,
-.B -j
-will have no effect, and vice versa.
-
 .TP
 .B -H
 Use homopolymer-compressed (HPC) minimizers. An HPC sequence is constructed by
@@ -334,10 +321,6 @@ Only map to the reverse complement strand of the reference sequences.
 If yes, sort anchors with heap merge, instead of radix sort. Heap merge is
 faster for short reads, but slower for long reads. [no]
 .TP
-.B --no-pairing
-Treat two reads in a pair as independent reads. The mate related fields in SAM
-are still properly populated.
-.TP
 .B --no-hash-name
 Produce the same alignment for identical sequences regardless of their sequence names.
 .SS Alignment options
@@ -371,7 +354,16 @@ Splice model [1]. 0 for the original minimap2 splice model that always penalizes
 .B -C
 has no effect with the default
 .BR -J1 .
-.BR -J0 .
+.TP
+.BR -j \ FILE
+Junctions used to extend alignment towards ends of reads [].
+.I FILE
+can be gene annotations in the BED12 format (aka 12-column BED), or intron
+positions in 5-column BED. BED12 file can be converted from GTF/GFF3 with
+`paftools.js gff2bed anno.gtf'.
+This option is intended for short RNA-seq reads, while
+.B --junc-bed
+for long noisy RNA-seq reads.
 .TP
 .BI -C \ INT
 Cost for a non-canonical GT-AG splicing (effective with
@@ -416,9 +408,14 @@ Score bonus when alignment extends to the end of the query sequence [0].
 .BI --score-N \ INT
 Penalty of a mismatch involving ambiguous bases [1].
 .TP
-.BI --pe-ind-chain
+.BR --pairing = strong | weak | no
-For paired-end short reads, perform chaining for each end independently.
+How to pair paired-end reads [strong].
-By default, minimap2 jointly chains the two ends.
+.RB ` no '
+for aligning the two ends in a pair independently with no `properly paired' set.
+.RB ` weak '
+for aligning the two ends independently and then pairing the hits.
+.RB ` strong '
+for jointly aligning and pairing the two ends.
 .TP
 .BR --splice-flank = yes | no
 Assume the next base to a
@@ -453,13 +450,13 @@ but not
 .BR --junc-bed .
 .TP
 .BR --junc-bed \ FILE
-Junctions to prefer during base alignment.
+Junctions to prefer during base alignment [].
-.I FILE
+Same format as
-can be gene annotations in the BED12 format (aka 12-column BED), or intron
+.BR -j .
-positions in 5-column BED. BED12 file can be converted from GTF/GFF3 with
+It is
-`paftools.js gff2bed anno.gtf'. It is
 .I NOT
-recommended to apply this option for short RNA-seq reads. []
+recommended to apply this option to short RNA-seq reads. This would increase
+run time with little improvement to junction accuracy.
 .TP
 .BR --junc-bonus \ INT
 Score bonus for a splice donor or acceptor found in annotation [9]. Effective with
@@ -467,13 +464,6 @@ Score bonus for a splice donor or acceptor found in annotation [9]. Effective wi
 but not
 .BR --spsc .
 .TP
-.BR --jump-bed \ FILE
-Junctions used to extend alignment towards ends of reads. Same format as with
-.BR --junc-bed .
-This option is intended for short RNA-seq reads, while
-.B --junc-bed
-for long noisy RNA-seq reads. []
-.TP
 .BI --end-seed-pen \ INT
 Drop a terminal anchor if
 .IR s <log( g )+ INT ,
@@ -699,7 +689,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 --pe-ind-chain
+.B --frag=yes -m25 -s40 -2K50m --heap-sort=yes --pairing=weak
 .BR --secondary=no ).
 .TP
 .B sr

options.c

@@ -179,7 +179,7 @@ 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_PE_IND;
+			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->noncan = 5, mo->b = 4, mo->q = 6, mo->q2 = 24;
 			mo->min_chain_score = 25;
 			mo->min_dp_max = 40;

setup.py

@@ -33,7 +33,7 @@ setup(
 	keywords = 'sequence-alignment',
 	scripts = ['python/minimap2.py'],
 	ext_modules = [Extension('mappy',
-		sources = ['python/mappy.pyx', 'align.c', 'bseq.c', 'lchain.c', 'seed.c', 'format.c', 'hit.c', 'index.c', 'pe.c', 'options.c',
+		sources = ['python/mappy.pyx', 'align.c', 'bseq.c', 'lchain.c', 'seed.c', 'format.c', 'hit.c', 'index.c', 'pe.c', 'jump.c', 'options.c',
 				   'ksw2_extd2_sse.c', 'ksw2_exts2_sse.c', 'ksw2_extz2_sse.c', 'ksw2_ll_sse.c',
 				   'kalloc.c', 'kthread.c', 'map.c', 'misc.c', 'sdust.c', 'sketch.c', 'esterr.c', 'splitidx.c'],
 		depends = ['minimap.h', 'bseq.h', 'kalloc.h', 'kdq.h', 'khash.h', 'kseq.h', 'ksort.h',