Release minimap2-2.5 (r572)

Fix typo in strand property

NEWS.md

@@ -1,3 +1,26 @@
+Release 2.5-r572 (11 November 2017)
+-----------------------------------
+
+This release fixes several bugs and brings a couple of minor improvements:
+
+ * Fixed a severe bug that leads to incorrect mapping coordinates in rare
+   corner cases.
+
+ * Fixed underestimated mapping quality for chimeric alignments when the whole
+   query sequence contain many repetitive minimizers, and for chimeric
+   alignments caused by Z-drop.
+
+ * Fixed two bugs in Python binding: incorrect strand field (#57) and incorrect
+   sequence names for Python3 (#55).
+
+ * Improved mapping accuracy for highly overlapping paired ends.
+
+ * Added option -Y to use soft clipping for supplementary alignments (#56).
+
+(2.5: 11 November 2017, r572)
+
+
+
 Release 2.4-r555 (6 November 2017)
 ----------------------------------
 

README.md

@@ -1,4 +1,4 @@
-[![Release](https://img.shields.io/badge/Release-v2.4-blue.svg?style=flat)](https://github.com/lh3/minimap2/releases)
+[![Release](https://img.shields.io/badge/Release-v2.5-blue.svg?style=flat)](https://github.com/lh3/minimap2/releases)
 [![BioConda](https://img.shields.io/conda/vn/bioconda/minimap2.svg?style=flat)](https://anaconda.org/bioconda/minimap2)
 [![PyPI](https://img.shields.io/pypi/v/mappy.svg?style=flat)](https://pypi.python.org/pypi/mappy)
 [![Python Version](https://img.shields.io/pypi/pyversions/mappy.svg?style=flat)](https://pypi.python.org/pypi/mappy)
@@ -22,7 +22,7 @@ cd minimap2 && make
 ./minimap2 -ax splice -k14 -uf ref.fa reads.fa > aln.sam  # Nanopore Direct RNA-seq
 ./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-one reads.fa reads.fa > overlaps.paf    # Nanopore read overlap
+./minimap2 -x ava-ont reads.fa reads.fa > overlaps.paf    # Nanopore read overlap
 # man page for detailed command line options
 man ./minimap2.1
 ```
@@ -71,9 +71,9 @@ Detailed evaluations are available from the [minimap2 preprint][preprint].
 Minimap2 only works on x86-64 CPUs. You can acquire precompiled binaries from
 the [release page][release] with:
 ```sh
-wget --no-check-certificate -O- https://github.com/lh3/minimap2/releases/download/v2.4/minimap2-2.4_x64-linux.tar.bz2 \
+curl -L https://github.com/lh3/minimap2/releases/download/v2.5/minimap2-2.5_x64-linux.tar.bz2 \
   | tar -jxvf -
-./minimap2-2.4_x64-linux/minimap2
+./minimap2-2.5_x64-linux/minimap2
 ```
 If you want to compile from the source, you need to have a C compiler, GNU make
 and zlib development files installed. Then type `make` in the source code

align.c

@@ -110,8 +110,11 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 	}
 	if ((p->cigar[0]&0xf) == 1 || (p->cigar[0]&0xf) == 2) { // get rid of leading I or D
 		int32_t l = p->cigar[0] >> 4;
-		if ((p->cigar[0]&0xf) == 1) r->qs += l, *qshift = l;
-		else r->rs += l, *tshift = l;
+		if ((p->cigar[0]&0xf) == 1) {
+			if (r->rev) r->qe -= l;
+			else r->qs += l;
+			*qshift = l;
+		} else r->rs += l, *tshift = l;
 		--p->n_cigar;
 		memmove(p->cigar, p->cigar + 1, p->n_cigar * 4);
 	}

chain.c

@@ -49,7 +49,7 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m
 			int64_t dr = ri - a[j].x;
 			int32_t dq = qi - (int32_t)a[j].y, dd, sc, log_dd;
 			int32_t sidj = (a[j].y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT;
-			if (dr == 0 || dq <= 0) continue;
+			if ((sidi == sidj && dr == 0) || dq <= 0) continue; // don't skip if an anchor is used by multiple segments; see below
 			if ((sidi == sidj && dq > max_dist_y) || dq > max_dist_x) continue;
 			dd = dr > dq? dr - dq : dq - dr;
 			if (sidi == sidj && dd > bw) continue;
@@ -61,7 +61,8 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m
 				int c_log, c_lin;
 				c_lin = (int)(dd * .01 * avg_qspan);
 				c_log = log_dd;
-				if (dr > dq || sidi != sidj) sc -= c_lin < c_log? c_lin : c_log;
+				if (sidi != sidj && dr == 0) ++sc; // possibly due to overlapping paired ends; give a minor bonus
+				else if (dr > dq || sidi != sidj) sc -= c_lin < c_log? c_lin : c_log;
 				else sc -= c_lin + (c_log>>1);
 			} else sc -= (int)(dd * .01 * avg_qspan) + (log_dd>>1);
 			sc += f[j];

format.c

@@ -203,7 +203,6 @@ static void write_cs(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_
 static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 {
 	int type = r->inv? 'I' : r->id == r->parent? 'P' : 'S';
-	if (r->iden_flt) mm_sprintf_lite(s, "\tom:i:%d", r->mapq);
 	if (r->p) {
 		mm_sprintf_lite(s, "\tNM:i:%d\tms:i:%d\tAS:i:%d\tnn:i:%d", r->blen - r->mlen + r->p->n_ambi, r->p->dp_max, r->p->dp_score, r->p->n_ambi);
 		if (r->p->trans_strand == 1 || r->p->trans_strand == 2)
@@ -258,7 +257,7 @@ static inline const mm_reg1_t *get_sam_pri(int n_regs, const mm_reg1_t *regs)
 	return NULL;
 }
 
-static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, const mm_reg1_t *r)
+static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, const mm_reg1_t *r, int opt_flag)
 {
 	if (r->p == 0) {
 		mm_sprintf_lite(s, "*");
@@ -267,14 +266,14 @@ static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, co
 		clip_len[0] = r->rev? qlen - r->qe : r->qs;
 		clip_len[1] = r->rev? r->qs : qlen - r->qe;
 		if (in_tag) {
-			int clip_char = (sam_flag&0x800)? 5 : 4;
+			int clip_char = (sam_flag&0x800) && !(opt_flag&MM_F_SOFTCLIP)? 5 : 4;
 			mm_sprintf_lite(s, "\tCG:B:I");
 			if (clip_len[0]) mm_sprintf_lite(s, ",%u", clip_len[0]<<4|clip_char);
 			for (k = 0; k < r->p->n_cigar; ++k)
 				mm_sprintf_lite(s, ",%u", r->p->cigar[k]);
 			if (clip_len[1]) mm_sprintf_lite(s, ",%u", clip_len[1]<<4|clip_char);
 		} else {
-			int clip_char = (sam_flag&0x800)? 'H' : 'S';
+			int clip_char = (sam_flag&0x800) && !(opt_flag&MM_F_SOFTCLIP)? 'H' : 'S';
 			if (clip_len[0]) mm_sprintf_lite(s, "%d%c", clip_len[0], clip_char);
 			for (k = 0; k < r->p->n_cigar; ++k)
 				mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, "MIDN"[r->p->cigar[k]&0xf]);
@@ -336,9 +335,8 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 			mm_sprintf_lite(s, "\t%s\t%d\t0\t*", mi->seq[this_rid].name, this_pos+1);
 		} else mm_sprintf_lite(s, "\t*\t0\t0\t*");
 	} else {
-		int mapq = !r->iden_flt? r->mapq : r->mapq < 3? r->mapq : 3;
 		this_rid = r->rid, this_pos = r->rs, this_rev = r->rev;
-		mm_sprintf_lite(s, "\t%s\t%d\t%d\t", mi->seq[r->rid].name, r->rs+1, mapq);
+		mm_sprintf_lite(s, "\t%s\t%d\t%d\t", mi->seq[r->rid].name, r->rs+1, r->mapq);
 		if ((opt_flag & MM_F_LONG_CIGAR) && r->p && r->p->n_cigar > max_bam_cigar_op - 2) {
 			int n_cigar = r->p->n_cigar;
 			if (r->qs != 0) ++n_cigar;
@@ -350,7 +348,7 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 			if (flag & 0x100) mm_sprintf_lite(s, "0S"); // secondary alignment
 			else if (flag & 0x800) mm_sprintf_lite(s, "%dS", r->re - r->rs); // supplementary alignment
 			else mm_sprintf_lite(s, "%dS", t->l_seq);
-		} else write_sam_cigar(s, flag, 0, t->l_seq, r);
+		} else write_sam_cigar(s, flag, 0, t->l_seq, r, opt_flag);
 	}
 
 	// write mate positions
@@ -383,7 +381,7 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 		if (t->qual) sam_write_sq(s, t->qual, t->l_seq, 0, 0);
 		else mm_sprintf_lite(s, "*");
 	} else {
-		if ((flag & 0x900) == 0) {
+		if ((flag & 0x900) == 0 || (opt_flag & MM_F_SOFTCLIP)) {
 			sam_write_sq(s, t->seq, t->l_seq, r->rev, r->rev);
 			mm_sprintf_lite(s, "\t");
 			if (t->qual) sam_write_sq(s, t->qual, t->l_seq, r->rev, 0);
@@ -431,7 +429,7 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 		if (r->p && (opt_flag & MM_F_OUT_CS))
 			write_cs(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG));
 		if (cigar_in_tag)
-			write_sam_cigar(s, flag, 1, t->l_seq, r);
+			write_sam_cigar(s, flag, 1, t->l_seq, r, opt_flag);
 	}
 
 	s->s[s->l] = 0; // we always have room for an extra byte (see str_enlarge)

hit.c

@@ -76,7 +76,7 @@ mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u,
 		mm_reg1_t *ri = &r[i];
 		ri->id = i;
 		ri->parent = MM_PARENT_UNSET;
-		ri->score = z[i].x >> 32;
+		ri->score = ri->score0 = z[i].x >> 32;
 		ri->hash = (uint32_t)z[i].x;
 		ri->cnt = (int32_t)z[i].y;
 		ri->as = z[i].y >> 32;
@@ -406,27 +406,33 @@ void mm_seg_free(void *km, int n_segs, mm_seg_t *segs)
 void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr)
 {
 	static const float q_coef = 40.0f;
+	int64_t sum_sc = 0;
+	float uniq_ratio;
 	int i;
+	for (i = 0; i < n_regs; ++i)
+		if (regs[i].parent == regs[i].id)
+			sum_sc += regs[i].score;
+	uniq_ratio = (float)sum_sc / (sum_sc + rep_len);
 	for (i = 0; i < n_regs; ++i) {
 		mm_reg1_t *r = &regs[i];
 		if (r->inv) {
 			r->mapq = 0;
 		} else if (r->parent == r->id) {
 			int mapq, subsc;
-			float pen_s1 = (r->score > 100? 1.0f : 0.01f * r->score) * ((float)r->score / (r->score + rep_len));
+			float pen_s1 = (r->score > 100? 1.0f : 0.01f * r->score) * uniq_ratio;
 			float pen_cm = r->cnt > 10? 1.0f : 0.1f * r->cnt;
 			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->score;
+				float 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
 				}
 			} else {
-				float x = (float)subsc / r->score;
+				float x = (float)subsc / r->score0;
 				if (r->p) {
 					float identity = (float)r->mlen / r->blen;
 					mapq = (int)(identity * pen_cm * q_coef * (1.0f - x) * logf((float)r->p->dp_max / match_sc));

main.c

@@ -6,7 +6,7 @@
 #include "mmpriv.h"
 #include "getopt.h"
 
-#define MM_VERSION "2.4-r555"
+#define MM_VERSION "2.5-r572"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -109,6 +109,7 @@ int main(int argc, char *argv[])
 		else if (c == 'X') opt.flag |= MM_F_AVA | MM_F_NO_SELF;
 		else if (c == 'a') opt.flag |= MM_F_OUT_SAM | MM_F_CIGAR;
 		else if (c == 'Q') opt.flag |= MM_F_NO_QUAL;
+		else if (c == 'Y') opt.flag |= MM_F_SOFTCLIP;
 		else if (c == 'L') opt.flag |= MM_F_LONG_CIGAR;
 		else if (c == 'T') opt.sdust_thres = atoi(optarg);
 		else if (c == 'n') opt.min_cnt = atoi(optarg);
@@ -232,6 +233,7 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -R STR       SAM read group line in a format like '@RG\\tID:foo\\tSM:bar' []\n");
 		fprintf(fp_help, "    -c           output CIGAR in PAF\n");
 		fprintf(fp_help, "    --cs[=STR]   output the cs tag; STR is 'short' (if absent) or 'long' [none]\n");
+		fprintf(fp_help, "    -Y           use soft clipping for supplementary alignments\n");
 		fprintf(fp_help, "    -t INT       number of threads [%d]\n", n_threads);
 		fprintf(fp_help, "    -K NUM       minibatch size for mapping [500M]\n");
 //		fprintf(fp_help, "    -v INT       verbose level [%d]\n", mm_verbose);

minimap.h

@@ -24,6 +24,7 @@
 #define MM_F_LONG_CIGAR    0x10000
 #define MM_F_INDEPEND_SEG  0x20000
 #define MM_F_SPLICE_FLANK  0x40000
+#define MM_F_SOFTCLIP      0x80000
 
 #define MM_IDX_MAGIC   "MMI\2"
 
@@ -64,7 +65,7 @@ typedef struct {
 
 typedef struct {
 	int32_t id;                     // ID for internal uses (see also parent below)
-	uint32_t cnt:30, rev:1, seg_split:1; // number of minimizers; if on the reverse strand
+	uint32_t cnt:28, rev:1, seg_split:1, sam_pri:1, proper_frag:1; // number of minimizers; if on the reverse strand
 	uint32_t rid:31, inv:1;         // reference index; if this is an alignment from inversion rescue
 	int32_t score;                  // DP alignment score
 	int32_t qs, qe, rs, re;         // query start and end; reference start and end
@@ -72,7 +73,7 @@ typedef struct {
 	int32_t as;                     // offset in the a[] array (for internal uses only)
 	int32_t mlen, blen;             // seeded exact match length; seeded alignment block length
 	uint32_t mapq:8, split:2, n_sub:22; // mapQ; split pattern; number of suboptimal mappings
-	uint32_t sam_pri:1, proper_frag:1, iden_flt:1, pe_thru:1, dummy:28;
+	uint32_t pe_thru:1, score0:31;
 	uint32_t hash;
 	mm_extra_t *p;
 } mm_reg1_t;

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "6 November 2017" "minimap2-2.4 (r555)" "Bioinformatics tools"
+.TH minimap2 1 "11 November 2017" "minimap2-2.5 (r572)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -302,6 +302,9 @@ is given,
 .I short
 is assumed. [none]
 .TP
+.B -Y
+In SAM output, use soft clipping for supplementary alignments.
+.TP
 .BI --seed \ INT
 Integer seed for randomizing equally best hits. Minimap2 hashes
 .I INT

python/mappy.pyx

@@ -13,7 +13,8 @@ cdef class Alignment:
 	cdef _ctg, _cigar # these are python objects
 
 	def __cinit__(self, ctg, cl, cs, ce, strand, qs, qe, mapq, cigar, is_primary, mlen, blen, NM, trans_strand):
-		self._ctg, self._ctg_len, self._r_st, self._r_en = str(ctg), cl, cs, ce
+		self._ctg = ctg if isinstance(ctg, str) else ctg.decode()
+		self._ctg_len, self._r_st, self._r_en = cl, cs, ce
 		self._strand, self._q_st, self._q_en = strand, qs, qe
 		self._NM, self._mlen, self._blen = NM, mlen, blen
 		self._mapq = mapq
@@ -34,7 +35,7 @@ cdef class Alignment:
 	def r_en(self): return self._r_en
 
 	@property
-	def strand(self): return self.strand
+	def strand(self): return self._strand
 
 	@property
 	def trans_strand(self): return self._trans_strand
@@ -150,9 +151,11 @@ def fastx_read(fn):
 	ks = cmappy.mm_fastx_open(str.encode(fn))
 	if ks is NULL: return None
 	while cmappy.kseq_read(ks) >= 0:
-		if ks.qual.l > 0: qual = str(ks.qual.s)
+		if ks.qual.l > 0: qual = ks.qual.s if isinstance(ks.qual.s, str) else ks.qual.s.decode()
 		else: qual = None
-		yield str(ks.name.s), str(ks.seq.s), qual
+		name = ks.name.s if isinstance(ks.name.s, str) else ks.name.s.decode()
+		seq = ks.seq.s if isinstance(ks.seq.s, str) else ks.seq.s.decode()
+		yield name, seq, qual
 	cmappy.mm_fastx_close(ks)
 
 def verbose(v=None):

setup.py

@@ -23,7 +23,7 @@ def readme():
 
 setup(
 	name = 'mappy',
-	version = '2.4',
+	version = '2.5',
 	url = 'https://github.com/lh3/minimap2',
 	description = 'Minimap2 python binding',
 	long_description = readme(),