Release minimap2-2.10 (r761)

Fix typo in README.md

Makefile

@@ -6,16 +6,16 @@ PROG=		minimap2
 PROG_EXTRA=	sdust minimap2-lite
 LIBS=		-lm -lz -lpthread
 
-ifeq ($(arm_neon),)
-ifeq ($(sse2only),)
+ifeq ($(arm_neon),) # if arm_neon is not defined
+ifeq ($(sse2only),) # if sse2only is not defined
 	OBJS+=ksw2_extz2_sse41.o ksw2_extd2_sse41.o ksw2_exts2_sse41.o ksw2_extz2_sse2.o ksw2_extd2_sse2.o ksw2_exts2_sse2.o ksw2_dispatch.o
-else
+else                # if sse2only is defined
 	OBJS+=ksw2_extz2_sse.o ksw2_extd2_sse.o ksw2_exts2_sse.o
 endif
-else
+else                # if arm_neon is defined
     OBJS+=ksw2_extz2_neon.o ksw2_extd2_neon.o ksw2_exts2_neon.o
     CFLAGS+=-D_FILE_OFFSET_BITS=64 -mfpu=neon -fsigned-char
-    INCLUDES+=-I sse2neon
+    INCLUDES+=-Isse2neon
 endif
 
 .PHONY:all extra clean depend
@@ -42,26 +42,31 @@ sdust:sdust.c getopt.o kalloc.o kalloc.h kdq.h kvec.h kseq.h sdust.h
 
 # SSE-specific targets on x86/x86_64
 
+ifeq ($(arm_neon),)   # if arm_neon is defined, compile this target with the default setting (i.e. no -msse2)
+ksw2_ll_sse.o:ksw2_ll_sse.c ksw2.h kalloc.h
+		$(CC) -c $(CFLAGS) -msse2 $(CPPFLAGS) $(INCLUDES) $< -o $@
+endif
+
 ksw2_extz2_sse41.o:ksw2_extz2_sse.c ksw2.h kalloc.h
-		$(CC) -c -msse4 $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
+		$(CC) -c $(CFLAGS) -msse4.1 $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
 
 ksw2_extz2_sse2.o:ksw2_extz2_sse.c ksw2.h kalloc.h
-		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH -DKSW_SSE2_ONLY $(INCLUDES) $< -o $@
+		$(CC) -c $(CFLAGS) -msse2 -mno-sse4.1 $(CPPFLAGS) -DKSW_CPU_DISPATCH -DKSW_SSE2_ONLY $(INCLUDES) $< -o $@
 
 ksw2_extd2_sse41.o:ksw2_extd2_sse.c ksw2.h kalloc.h
-		$(CC) -c -msse4 $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
+		$(CC) -c $(CFLAGS) -msse4.1 $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
 
 ksw2_extd2_sse2.o:ksw2_extd2_sse.c ksw2.h kalloc.h
-		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH -DKSW_SSE2_ONLY $(INCLUDES) $< -o $@
+		$(CC) -c $(CFLAGS) -msse2 -mno-sse4.1 $(CPPFLAGS) -DKSW_CPU_DISPATCH -DKSW_SSE2_ONLY $(INCLUDES) $< -o $@
 
 ksw2_exts2_sse41.o:ksw2_exts2_sse.c ksw2.h kalloc.h
-		$(CC) -c -msse4 $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
+		$(CC) -c $(CFLAGS) -msse4.1 $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
 
 ksw2_exts2_sse2.o:ksw2_exts2_sse.c ksw2.h kalloc.h
-		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH -DKSW_SSE2_ONLY $(INCLUDES) $< -o $@
+		$(CC) -c $(CFLAGS) -msse2 -mno-sse4.1 $(CPPFLAGS) -DKSW_CPU_DISPATCH -DKSW_SSE2_ONLY $(INCLUDES) $< -o $@
 
 ksw2_dispatch.o:ksw2_dispatch.c ksw2.h
-		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
+		$(CC) -c $(CFLAGS) -msse4.1 $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
 
 # NEON-specific targets on ARM
 

NEWS.md

@@ -1,3 +1,50 @@
+Release 2.10-r761 (27 March 2018)
+---------------------------------
+
+Changes to minimap2:
+
+ * Optionally output the MD tag for compatibility with existing tools (#63,
+   #118 and #137).
+
+ * Use SSE compiler flags more precisely to prevent compiling errors on certain
+   machines (#127).
+
+ * Added option --min-occ-floor to set a minimum occurrence threshold. Presets
+   intended for assembly-to-reference alignment set this option to 100. This
+   option alleviates issues with regions having high copy numbers (#107).
+
+ * Exit with non-zero code on file writing errors (e.g. disk full; #103 and
+   #132).
+
+ * Added option -y to copy FASTA/FASTQ comments in query sequences to the
+   output (#136).
+
+ * Added the asm20 preset for alignments between genomes at 5-10% sequence
+   divergence.
+
+ * Changed the band-width in the ava-ont preset from 500 to 2000. Oxford
+   Nanopore reads may contain long deletion sequencing errors that break
+   chaining.
+
+Changes to mappy, the Python binding:
+
+ * Fixed a typo in Align.seq() (#126).
+
+Changes to paftools.js, the companion script:
+
+ * Command sam2paf now converts the MD tag to cs.
+
+ * Support VCF output for assembly-to-reference variant calling (#109).
+
+This version should produce identical alignment for read overlapping, RNA-seq
+read mapping, and genomic read mapping. We have also added a cook book to show
+the variety uses of minimap2 on real datasets. Please see cookbook.md in the
+minimap2 source code directory.
+
+(2.10: 27 March 2017, r761)
+
+
+
 Release 2.9-r720 (23 February 2018)
 -----------------------------------
 

README.md

@@ -68,9 +68,8 @@ Detailed evaluations are available from the [minimap2 preprint][preprint].
 Minimap2 is optimized for x86-64 CPUs. You can acquire precompiled binaries from
 the [release page][release] with:
 ```sh
-curl -L https://github.com/lh3/minimap2/releases/download/v2.9/minimap2-2.9_x64-linux.tar.bz2 \
-  | tar -jxvf -
-./minimap2-2.9_x64-linux/minimap2
+curl -L https://github.com/lh3/minimap2/releases/download/v2.10/minimap2-2.10_x64-linux.tar.bz2 | tar -jxvf -
+./minimap2-2.10_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
@@ -137,7 +136,7 @@ Nanopore reads.
 #### <a name="map-long-splice"></a>Map long mRNA/cDNA reads
 
 ```sh
-minimap2 -ax splice -uf ref.fa iso-seq.fq > aln.sam          # PacBio Iso-seq/traditional cDNA
+minimap2 -ax splice -uf -C5 ref.fa iso-seq.fq > aln.sam      # PacBio Iso-seq/traditional cDNA
 minimap2 -ax splice ref.fa nanopore-cdna.fa > aln.sam        # Nanopore 2D cDNA-seq
 minimap2 -ax splice -uf -k14 ref.fa direct-rna.fq > aln.sam  # Nanopore Direct RNA-seq
 minimap2 -ax splice --splice-flank=no SIRV.fa SIRV-seq.fa    # mapping against SIRV control
@@ -229,7 +228,7 @@ sorting) still work with such BAM records; tools that read CIGAR will
 effectively ignore these records. It has been decided that future tools will
 will seamlessly recognize long-cigar records generated by option `-L`.
 
-**TD;DR**: if you work with ultra-long reads and use tools that only process
+**TL;DR**: if you work with ultra-long reads and use tools that only process
 BAM files, please add option `-L`.
 
 #### <a name="cs"></a>The cs optional tag
@@ -346,6 +345,10 @@ mappy` or [from BioConda][mappyconda] via `conda install -c bioconda mappy`.
   possible to add non-SIMD support, but it would make minimap2 slower by
   several times.
 
+* Minimap2 does not work with a single query or database sequence ~2
+  billion bases or longer (2,147,483,647 to be exact). The total length of all
+  sequences can well exceed this threshold.
+
 
 
 [paf]: https://github.com/lh3/miniasm/blob/master/PAF.md

bseq.c

@@ -62,7 +62,7 @@ static inline char *kstrdup(const kstring_t *s)
 	return t;
 }
 
-static inline void kseq2bseq(kseq_t *ks, mm_bseq1_t *s, int with_qual)
+static inline void kseq2bseq(kseq_t *ks, mm_bseq1_t *s, int with_qual, int with_comment)
 {
 	int i;
 	s->name = kstrdup(&ks->name);
@@ -71,10 +71,11 @@ static inline void kseq2bseq(kseq_t *ks, mm_bseq1_t *s, int with_qual)
 		if (s->seq[i] == 'u' || s->seq[i] == 'U')
 			--s->seq[i];
 	s->qual = with_qual && ks->qual.l? kstrdup(&ks->qual) : 0;
+	s->comment = with_comment && ks->comment.l? kstrdup(&ks->comment) : 0;
 	s->l_seq = ks->seq.l;
 }
 
-mm_bseq1_t *mm_bseq_read2(mm_bseq_file_t *fp, int chunk_size, int with_qual, int frag_mode, int *n_)
+mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int with_comment, int frag_mode, int *n_)
 {
 	int64_t size = 0;
 	kvec_t(mm_bseq1_t) a = {0,0,0};
@@ -91,12 +92,12 @@ mm_bseq1_t *mm_bseq_read2(mm_bseq_file_t *fp, int chunk_size, int with_qual, int
 		assert(ks->seq.l <= INT32_MAX);
 		if (a.m == 0) kv_resize(mm_bseq1_t, 0, a, 256);
 		kv_pushp(mm_bseq1_t, 0, a, &s);
-		kseq2bseq(ks, s, with_qual);
+		kseq2bseq(ks, s, with_qual, with_comment);
 		size += s->l_seq;
 		if (size >= chunk_size) {
 			if (frag_mode && a.a[a.n-1].l_seq < CHECK_PAIR_THRES) {
 				while (kseq_read(ks) >= 0) {
-					kseq2bseq(ks, &fp->s, with_qual);
+					kseq2bseq(ks, &fp->s, with_qual, with_comment);
 					if (mm_qname_same(fp->s.name, a.a[a.n-1].name)) {
 						kv_push(mm_bseq1_t, 0, a, fp->s);
 						memset(&fp->s, 0, sizeof(mm_bseq1_t));
@@ -110,12 +111,17 @@ mm_bseq1_t *mm_bseq_read2(mm_bseq_file_t *fp, int chunk_size, int with_qual, int
 	return a.a;
 }
 
+mm_bseq1_t *mm_bseq_read2(mm_bseq_file_t *fp, int chunk_size, int with_qual, int frag_mode, int *n_)
+{
+	return mm_bseq_read3(fp, chunk_size, with_qual, 0, frag_mode, n_);
+}
+
 mm_bseq1_t *mm_bseq_read(mm_bseq_file_t *fp, int chunk_size, int with_qual, int *n_)
 {
 	return mm_bseq_read2(fp, chunk_size, with_qual, 0, n_);
 }
 
-mm_bseq1_t *mm_bseq_read_frag(int n_fp, mm_bseq_file_t **fp, int chunk_size, int with_qual, int *n_)
+mm_bseq1_t *mm_bseq_read_frag2(int n_fp, mm_bseq_file_t **fp, int chunk_size, int with_qual, int with_comment, int *n_)
 {
 	int i;
 	int64_t size = 0;
@@ -136,7 +142,7 @@ mm_bseq1_t *mm_bseq_read_frag(int n_fp, mm_bseq_file_t **fp, int chunk_size, int
 		for (i = 0; i < n_fp; ++i) {
 			mm_bseq1_t *s;
 			kv_pushp(mm_bseq1_t, 0, a, &s);
-			kseq2bseq(fp[i]->ks, s, with_qual);
+			kseq2bseq(fp[i]->ks, s, with_qual, with_comment);
 			size += s->l_seq;
 		}
 		if (size >= chunk_size) break;
@@ -145,6 +151,11 @@ mm_bseq1_t *mm_bseq_read_frag(int n_fp, mm_bseq_file_t **fp, int chunk_size, int
 	return a.a;
 }
 
+mm_bseq1_t *mm_bseq_read_frag(int n_fp, mm_bseq_file_t **fp, int chunk_size, int with_qual, int *n_)
+{
+	return mm_bseq_read_frag2(n_fp, fp, chunk_size, with_qual, 0, n_);
+}
+
 int mm_bseq_eof(mm_bseq_file_t *fp)
 {
 	return (ks_eof(fp->ks->f) && fp->s.seq == 0);

bseq.h

@@ -13,13 +13,15 @@ typedef struct mm_bseq_file_s mm_bseq_file_t;
 
 typedef struct {
 	int l_seq, rid;
-	char *name, *seq, *qual;
+	char *name, *seq, *qual, *comment;
 } mm_bseq1_t;
 
 mm_bseq_file_t *mm_bseq_open(const char *fn);
 void mm_bseq_close(mm_bseq_file_t *fp);
+mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int with_comment, int frag_mode, int *n_);
 mm_bseq1_t *mm_bseq_read2(mm_bseq_file_t *fp, int chunk_size, int with_qual, int frag_mode, int *n_);
 mm_bseq1_t *mm_bseq_read(mm_bseq_file_t *fp, int chunk_size, int with_qual, int *n_);
+mm_bseq1_t *mm_bseq_read_frag2(int n_fp, mm_bseq_file_t **fp, int chunk_size, int with_qual, int with_comment, int *n_);
 mm_bseq1_t *mm_bseq_read_frag(int n_fp, mm_bseq_file_t **fp, int chunk_size, int with_qual, int *n_);
 int mm_bseq_eof(mm_bseq_file_t *fp);
 

cookbook.md

@@ -0,0 +1,243 @@
+## Table of Contents
+
+- [Introduction & Installation](#intro)
+- [Mapping Genomic Reads](#map-reads)
+  * [Mapping long reads](#map-pb)
+  * [Mapping Illumina paired-end reads](#map-sr)
+  * [Evaluating mapping accuracy with simulated reads (for developers)](#mapeval)
+- [Mapping Long RNA-seq Reads](#map-rna)
+  * [Mapping Nanopore 2D cDNA reads](#map-ont-cdna-2d)
+  * [Mapping Nanopore direct-RNA reads](#map-direct-rna)
+  * [Mapping PacBio Iso-seq reads](#map-iso-seq)
+- [Full-Genome Alignment](#genome-aln)
+  * [Intra-species assembly alignment](#asm-to-ref)
+  * [Cross-species full-genome alignment](#x-species)
+  * [Eyeballing alignment](#view-aln)
+  * [Calling variants from assembly-to-reference alignment](#asm-var)
+  * [Constructing self-homology map](#hom-map)
+  * [Lift Over (for developers)](#liftover)
+- [Read Overlap](#read-overlap)
+  * [Long-read overlap](#long-read-overlap)
+  * [Evaluating overlap sensitivity (for developers)](#ov-eval)
+
+## <a name="intro"></a>Introduction & Installation
+
+This cookbook walks you through a variety of applications of minimap2 and its
+companion script `paftools.js`. All data here are freely available from the
+minimap2 release page at version tag [v2.10][v2.10]. Some examples only work
+with v2.10 or later.
+
+To acquire the data used in this cookbook and to install minimap2 and paftools,
+please follow the command lines below:
+```sh
+# install minimap2 executables
+curl -L https://github.com/lh3/minimap2/releases/download/v2.10/minimap2-2.10_x64-linux.tar.bz2 | tar jxf -
+cp minimap2-2.10_x64-linux/{minimap2,k8,paftools.js} .  # copy executables
+export PATH="$PATH:"`pwd`                               # put the current directory on PATH
+# download example datasets
+curl -L https://github.com/lh3/minimap2/releases/download/v2.10/cookbook-data.tgz | tar zxf -
+```
+
+## <a name="map-reads"></a>Mapping Genomic Reads
+
+### <a name="map-pb"></a>Mapping long reads
+```sh
+minimap2 -ax map-pb -t4 ecoli_ref.fa ecoli_p6_25x_canu.fa > mapped.sam
+```
+Alternatively, you can create a minimap2 index first and then map:
+```sh
+minimap2 -x map-pb -d ecoli-pb.mmi ecoli_ref.fa                      # create an index
+minimap2 -ax map-pb ecoli-pb.mmi ecoli_p6_25x_canu.fa > mapped.sam
+```
+This will save you a couple of minutes when you map against the human genome.
+**HOWEVER**, key algorithm parameters such as the k-mer length and window
+size can't be changed after indexing. Minimap2 will give you a warning if
+parameters used in a pre-built index doesn't match parameters on the command
+line. **Please always make sure you are using an intended pre-built index.**
+
+### <a name="map-sr"></a>Mapping Illumina paired-end reads:
+```sh
+minimap2 -ax sr -t4 ecoli_ref.fa ecoli_mason_1.fq ecoli_mason_2.fq > mapped-sr.sam
+```
+
+### <a name="mapeval"></a>Evaluating mapping accuracy with simulated reads (for developers)
+```sh
+minimap2 -ax sr ecoli_ref.fa ecoli_mason_1.fq ecoli_mason_2.fq | paftools.js mapeval -
+```
+The output is:
+```
+Q       60      19712   0       0.000000000     19712
+Q       0       282     219     0.010953286     19994
+U       6
+```
+where a `U`-line gives the number of unmapped reads (for SAM input only); a
+`Q`-line gives:
+
+1. Mapping quality (mapQ) threshold
+2. Number of mapped reads between this threshold and the previous mapQ threshold.
+3. Number of wrong mappings in the same mapQ interval
+4. Accumulative mapping error rate
+5. Accumulative number of mappings
+
+For `paftools.js mapeval` to work, you need to encode the true read positions
+in read names in the right format. For [PBSIM][pbsim] and [mason2][mason2], we
+provide scripts to generate the right format. Simulated reads in this cookbook
+were created with the following command lines:
+```sh
+# in PBSIM source code directory:
+src/pbsim ../ecoli_ref.fa --depth 1 --sample-fastq sample/sample.fastq
+paftools.js pbsim2fq ../ecoli_ref.fa.fai sd_0001.maf > ../ecoli_pbsim.fa
+
+# mason2 simulation
+mason_simulator --illumina-prob-mismatch-scale 2.5 -ir ecoli_ref.fa -n 10000 -o tmp-l.fq -or tmp-r.fq -oa tmp.sam
+paftools.js mason2fq tmp.sam | seqtk seq -1 > ecoli_mason_1.fq
+paftools.js mason2fq tmp.sam | seqtk seq -2 > ecoli_mason_2.fq
+```
+
+
+
+## <a name="map-rna"></a>Mapping Long RNA-seq Reads
+
+### <a name="map-ont-cdna-2d"></a>Mapping Nanopore 2D cDNA reads
+```sh
+minimap2 -ax splice SIRV_E2.fa SIRV_ont-cdna.fa > aln.sam
+```
+You can compare the alignment to the true annotations with:
+```sh
+paftools.js junceval SIRV_E2C.gtf aln.sam
+```
+It gives the percentage of introns found in the annotation. For SIRV data, it
+is possible to achieve higher junction accuracy with
+```sh
+minimap2 -ax splice --splice-flank=no SIRV_E2.fa SIRV_ont-cdna.fa | paftools.js junceval SIRV_E2C.gtf
+```
+This is because minimap2 models one additional evolutionarily conserved base
+around a canonical junction, but SIRV doesn't honor this signal. Option
+`--splice-flank=no` asks minimap2 no to model this additional base.
+
+In the output a tag `ts:A:+` indicates that the read strand is the same as the
+transcript strand; `ts:A:-` indicates the read strand is opposite to the
+transcript strand. This tag is inferred from the GT-AG signal and is thus only
+available to spliced reads.
+
+### <a name="map-direct-rna"></a>Mapping Nanopore direct-RNA reads
+```sh
+minimap2 -ax splice -k14 -uf SIRV_E2.fa SIRV_ont-drna.fa > aln.sam
+```
+Direct-RNA reads are noisier, so we use a shorter k-mer for improved
+sensitivity. Here, option `-uf` forces minimap2 to map reads to the forward
+transcript strand only because direct-RNA reads are stranded. Again, applying
+`--splice-flank=no` helps junction accuracy for SIRV data.
+
+### <a name="map-iso-seq"></a>Mapping PacBio Iso-seq reads
+```sh
+minimap2 -ax splice -uf -C5 SIRV_E2.fa SIRV_iso-seq.fq > aln.sam
+```
+Option `-C5` reduces the penalty on non-canonical splicing sites. It helps
+to align such sites correctly for data with low error rate such as Iso-seq
+reads and traditional cDNAs. On this example, minimap2 makes one junction
+error. Applying `--splice-flank=no` fixes this alignment error.
+
+Note that the command line above is optimized for the final Iso-seq reads.
+PacBio's Iso-seq pipeline produces intermediate sequences at varying quality.
+For example, some intermediate reads are not stranded. For these reads, option
+`-uf` will lead to more errors. Please revise the minimap2 command line
+accordingly.
+
+
+
+## <a name="genome-aln"></a>Full-Genome Alignment
+
+### <a name="asm-to-ref"></a>Intra-species assembly alignment
+```sh
+# option "--cs" is recommended as paftools.js may need it
+minimap2 -cx asm5 --cs ecoli_ref.fa ecoli_canu.fa > ecoli_canu.paf
+```
+Here `ecoli_canu.fa` is the Canu assembly of `ecoli_p6_25x_canu.fa`. This
+command line outputs alignments in the [PAF format][paf]. Use `-a` instead of
+`-c` to get output in the SAM format.
+
+### <a name="x-species"></a>Cross-species full-genome alignment
+```sh
+minimap2 -cx asm20 --cs ecoli_ref.fa ecoli_O104:H4.fa > ecoli_O104:H4.paf
+sort -k6,6 -k8,8n ecoli_O104:H4.paf | paftools.js call -f ecoli_ref.fa -L10000 -l1000 - > out.vcf
+```
+Minimap2 has three presets for full-genome alignment: "asm5" for sequence
+divergence below 1%, "asm10" for divergence around a couple of percent and
+"asm20" for divergence not more than 10%. In theory, with the right setting,
+minimap2 should work for sequence pairs with sequence divergence up to ~15%,
+but this has not been carefully evaluated.
+
+### <a name="view-aln"></a>Eyeballing alignment
+```sh
+# option "--cs" required; minimap2-r741 or higher required for the "asm20" preset
+minimap2 -cx asm20 --cs ecoli_ref.fa ecoli_O104:H4.fa | paftools.js view - | less -S
+```
+This prints the alignment in a BLAST-like format.
+
+### <a name="asm-var"></a>Calling variants from assembly-to-reference alignment
+```sh
+# don't forget the "--cs" option; otherwise it doesn't work
+minimap2 -cx asm5 --cs ecoli_ref.fa ecoli_canu.fa \
+  | sort -k6,6 -k8,8n \
+  | paftools.js call -f ecoli_ref.fa - > out.vcf
+```
+Without option `-f`, `paftools.js call` outputs in a custom format. In this
+format, lines starting with `R` give the regions covered by one contig only.
+This information is not available in the VCF output.
+
+### <a name="hom-map"></a>Constructing self-homology map
+```sh
+minimap2 -DP -k19 -w19 -m200 ecoli_ref.fa ecoli_ref.fa > out.paf
+```
+Option `-D` asks minimap2 to ignore anchors from perfect self match and `-P`
+outputs all chains. For large nomes, we don't recommend to perform base-level
+alignment (with `-c`, `-a` or `--cs`) when `-P` is applied. This is because
+base-alignment is slow and occasionally gives wrong alignments close to the
+diagonal of a dotter plot. For E. coli, though, base-alignment is still fast.
+
+### <a name="liftover"></a>Lift over (for developers)
+```sh
+minimap2 -cx asm5 --cs ecoli_ref.fa ecoli_canu.fa > ecoli_canu.paf
+echo -e 'tig00000001\t200000\t300000' | paftools.js liftover ecoli_canu.paf -
+```
+This lifts over a region on query sequences to one or multiple regions on
+reference sequences. Note that this paftools.js command may not be efficient
+enough to lift millions of regions.
+
+
+
+## <a name="read-overlap"></a>Read Overlap
+
+### <a name="long-read-overlap"></a>Long read overlap
+```sh
+# For pacbio reads:
+minimap2 -x ava-pb ecoli_p6_25x_canu.fa ecoli_p6_25x_canu.fa > overlap.paf
+# For Nanopore reads (ava-ont also works with PacBio but not as good):
+minimap2 -x ava-ont -r 10000 ecoli_p6_25x_canu.fa ecoli_p6_25x_canu.fa > overlap.paf
+# If you have miniasm installed:
+miniasm -f ecoli_p6_25x_canu.fa overlap.paf > asm.gfa
+```
+Here we explicitly applied `-r 10000`. We are considering to set this as the
+default for the `ava-ont` mode as this seems to improve the contiguity for
+nanopore read assembly (Loman, personal communication).
+
+*Minimap2 doesn't work well with short-read overlap.*
+
+### <a name="ov-eval"></a>Evaluating overlap sensitivity (for developers)
+
+```sh
+# read to reference mapping
+minimap2 -cx map-pb ecoli_ref.fa ecoli_p6_25x_canu.fa > to-ref.paf
+# evaluate overlap sensitivity
+sort -k6,6 -k8,8n to-ref.paf | paftools.js ov-eval - overlap.paf
+```
+You can see that for PacBio reads, minimap2 achieves higher overlap sensitivity
+with `-x ava-pb` (99% vs 93% with `-x ava-ont`).
+
+
+
+[pbsim]: https://github.com/pfaucon/PBSIM-PacBio-Simulator
+[mason2]: https://github.com/seqan/seqan/tree/master/apps/mason2
+[paf]: https://github.com/lh3/miniasm/blob/master/PAF.md
+[v2.10]: https://github.com/lh3/minimap2/releases/tag/v2.10

format.c

@@ -118,7 +118,7 @@ void mm_write_sam_hdr(const mm_idx_t *idx, const char *rg, const char *ver, int
 	if (idx) {
 		uint32_t i;
 		for (i = 0; i < idx->n_seq; ++i)
-			printf("@SQ\tSN:%s\tLN:%d\n", idx->seq[i].name, idx->seq[i].len);
+			mm_sprintf_lite(&str, "@SQ\tSN:%s\tLN:%d\n", idx->seq[i].name, idx->seq[i].len);
 	}
 	if (rg) sam_write_rg_line(&str, rg);
 	mm_sprintf_lite(&str, "@PG\tID:minimap2\tPN:minimap2");
@@ -129,36 +129,18 @@ void mm_write_sam_hdr(const mm_idx_t *idx, const char *rg, const char *ver, int
 		for (i = 1; i < argc; ++i)
 			mm_sprintf_lite(&str, " %s", argv[i]);
 	}
-	mm_sprintf_lite(&str, "\n");
-	fputs(str.s, stdout);
+	mm_err_puts(str.s);
 	free(str.s);
 }
 
-static void write_cs(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int no_iden)
+static void write_cs_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq, const mm_reg1_t *r, char *tmp, int no_iden)
 {
-	extern unsigned char seq_nt4_table[256];
 	int i, q_off, t_off;
-	uint8_t *qseq, *tseq;
-	char *tmp;
-	if (r->p == 0) return;
 	mm_sprintf_lite(s, "\tcs:Z:");
-	qseq = (uint8_t*)kmalloc(km, r->qe - r->qs);
-	tseq = (uint8_t*)kmalloc(km, r->re - r->rs);
-	tmp = (char*)kmalloc(km, r->re - r->rs > r->qe - r->qs? r->re - r->rs + 1 : r->qe - r->qs + 1);
-	mm_idx_getseq(mi, r->rid, r->rs, r->re, tseq);
-	if (!r->rev) {
-		for (i = r->qs; i < r->qe; ++i)
-			qseq[i - r->qs] = seq_nt4_table[(uint8_t)t->seq[i]];
-	} else {
-		for (i = r->qs; i < r->qe; ++i) {
-			uint8_t c = seq_nt4_table[(uint8_t)t->seq[i]];
-			qseq[r->qe - i - 1] = c >= 4? 4 : 3 - c;
-		}
-	}
 	for (i = q_off = t_off = 0; i < r->p->n_cigar; ++i) {
 		int j, op = r->p->cigar[i]&0xf, len = r->p->cigar[i]>>4;
 		assert(op >= 0 && op <= 3);
-		if (op == 0) {
+		if (op == 0) { // match
 			int l_tmp = 0;
 			for (j = 0; j < len; ++j) {
 				if (qseq[q_off + j] != tseq[t_off + j]) {
@@ -179,17 +161,17 @@ static void write_cs(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_
 				} else mm_sprintf_lite(s, ":%d", l_tmp);
 			}
 			q_off += len, t_off += len;
-		} else if (op == 1) {
+		} else if (op == 1) { // insertion to ref
 			for (j = 0, tmp[len] = 0; j < len; ++j)
 				tmp[j] = "acgtn"[qseq[q_off + j]];
 			mm_sprintf_lite(s, "+%s", tmp);
 			q_off += len;
-		} else if (op == 2) {
+		} else if (op == 2) { // deletion from ref
 			for (j = 0, tmp[len] = 0; j < len; ++j)
 				tmp[j] = "acgtn"[tseq[t_off + j]];
 			mm_sprintf_lite(s, "-%s", tmp);
 			t_off += len;
-		} else {
+		} else { // intron
 			assert(len >= 2);
 			mm_sprintf_lite(s, "~%c%c%d%c%c", "acgtn"[tseq[t_off]], "acgtn"[tseq[t_off+1]],
 				len, "acgtn"[tseq[t_off+len-2]], "acgtn"[tseq[t_off+len-1]]);
@@ -197,6 +179,59 @@ static void write_cs(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_
 		}
 	}
 	assert(t_off == r->re - r->rs && q_off == r->qe - r->qs);
+}
+
+static void write_MD_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq, const mm_reg1_t *r, char *tmp)
+{
+	int i, q_off, t_off, l_MD = 0;
+	mm_sprintf_lite(s, "\tMD:Z:");
+	for (i = q_off = t_off = 0; i < r->p->n_cigar; ++i) {
+		int j, op = r->p->cigar[i]&0xf, len = r->p->cigar[i]>>4;
+		assert(op >= 0 && op <= 2); // introns (aka reference skips) are not supported
+		if (op == 0) { // match
+			for (j = 0; j < len; ++j) {
+				if (qseq[q_off + j] != tseq[t_off + j]) {
+					mm_sprintf_lite(s, "%d%c", l_MD, "ACGTN"[tseq[t_off + j]]);
+					l_MD = 0;
+				} else ++l_MD;
+			}
+			q_off += len, t_off += len;
+		} else if (op == 1) { // insertion to ref
+			q_off += len;
+		} else if (op == 2) { // deletion from ref
+			for (j = 0, tmp[len] = 0; j < len; ++j)
+				tmp[j] = "ACGTN"[tseq[t_off + j]];
+			mm_sprintf_lite(s, "%d^%s", l_MD, tmp);
+			l_MD = 0;
+			t_off += len;
+		}
+	}
+	if (l_MD > 0) mm_sprintf_lite(s, "%d", l_MD);
+	assert(t_off == r->re - r->rs && q_off == r->qe - r->qs);
+}
+
+static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int no_iden, int is_MD)
+{
+	extern unsigned char seq_nt4_table[256];
+	int i;
+	uint8_t *qseq, *tseq;
+	char *tmp;
+	if (r->p == 0) return;
+	qseq = (uint8_t*)kmalloc(km, r->qe - r->qs);
+	tseq = (uint8_t*)kmalloc(km, r->re - r->rs);
+	tmp = (char*)kmalloc(km, r->re - r->rs > r->qe - r->qs? r->re - r->rs + 1 : r->qe - r->qs + 1);
+	mm_idx_getseq(mi, r->rid, r->rs, r->re, tseq);
+	if (!r->rev) {
+		for (i = r->qs; i < r->qe; ++i)
+			qseq[i - r->qs] = seq_nt4_table[(uint8_t)t->seq[i]];
+	} else {
+		for (i = r->qs; i < r->qe; ++i) {
+			uint8_t c = seq_nt4_table[(uint8_t)t->seq[i]];
+			qseq[r->qe - i - 1] = c >= 4? 4 : 3 - c;
+		}
+	}
+	if (is_MD) write_MD_core(s, tseq, qseq, r, tmp);
+	else write_cs_core(s, tseq, qseq, r, tmp, no_iden);
 	kfree(km, qseq); kfree(km, tseq); kfree(km, tmp);
 }
 
@@ -237,8 +272,10 @@ void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const m
 		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]);
 	}
-	if (r->p && (opt_flag & MM_F_OUT_CS))
-		write_cs(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG));
+	if (r->p && (opt_flag & (MM_F_OUT_CS|MM_F_OUT_MD)))
+		write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD);
+	if ((opt_flag & MM_F_COPY_COMMENT) && t->comment)
+		mm_sprintf_lite(s, "\t%s", t->comment);
 }
 
 static void sam_write_sq(kstring_t *s, char *seq, int l, int rev, int comp)
@@ -434,12 +471,15 @@ 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 (r->p && (opt_flag & (MM_F_OUT_CS|MM_F_OUT_MD)))
+			write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD);
 		if (cigar_in_tag)
 			write_sam_cigar(s, flag, 1, t->l_seq, r, opt_flag);
 	}
 
+	if ((opt_flag & MM_F_COPY_COMMENT) && t->comment)
+		mm_sprintf_lite(s, "\t%s", t->comment);
+
 	s->s[s->l] = 0; // we always have room for an extra byte (see str_enlarge)
 }
 

main.c

@@ -4,9 +4,13 @@
 #include "bseq.h"
 #include "minimap.h"
 #include "mmpriv.h"
+#ifdef HAVE_GETOPT
+#include <getopt.h>
+#else
 #include "getopt.h"
+#endif
 
-#define MM_VERSION "2.9-r720"
+#define MM_VERSION "2.10-r761"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -51,6 +55,8 @@ static struct option long_options[] = {
 	{ "all-chain",      no_argument,       0, 'P' },
 	{ "dual",           required_argument, 0, 0 },   // 26
 	{ "max-clip-ratio", required_argument, 0, 0 },   // 27
+	{ "min-occ-floor",  required_argument, 0, 0 },   // 28
+	{ "MD",             no_argument,       0, 0 },   // 29
 	{ "help",           no_argument,       0, 'h' },
 	{ "max-intron-len", required_argument, 0, 'G' },
 	{ "version",        no_argument,       0, 'V' },
@@ -88,7 +94,7 @@ static inline void yes_or_no(mm_mapopt_t *opt, int flag, int long_idx, const cha
 
 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:O:E:m:N:Qu:R:hF:LC:";
+	const char *opt_str = "2aSDw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:O:E:m:N:Qu:R:hF:LC:y";
 	mm_mapopt_t opt;
 	mm_idxopt_t ipt;
 	int i, c, n_threads = 3, long_idx;
@@ -110,7 +116,7 @@ int main(int argc, char *argv[])
 			}
 			break;
 		}
-	optreset = 1;
+	optind = 0; // for musl getopt, optind=0 has the same effect as optreset=1; older libc doesn't have optreset
 
 	while ((c = getopt_long(argc, argv, opt_str, long_options, &long_idx)) >= 0) {
 		if (c == 'w') ipt.w = atoi(optarg);
@@ -134,6 +140,7 @@ int main(int argc, char *argv[])
 		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 == 'y') opt.flag |= MM_F_COPY_COMMENT;
 		else if (c == 'T') opt.sdust_thres = atoi(optarg);
 		else if (c == 'n') opt.min_cnt = atoi(optarg);
 		else if (c == 'm') opt.min_chain_score = atoi(optarg);
@@ -164,6 +171,8 @@ int main(int argc, char *argv[])
 		else if (c == 0 && long_idx ==22) opt.flag |= MM_F_FOR_ONLY; // --for-only
 		else if (c == 0 && long_idx ==23) opt.flag |= MM_F_REV_ONLY; // --rev-only
 		else if (c == 0 && long_idx ==27) opt.max_clip_ratio = atof(optarg); // --max-clip-ratio
+		else if (c == 0 && long_idx ==28) opt.min_mid_occ = atoi(optarg); // --min-occ-floor
+		else if (c == 0 && long_idx ==29) opt.flag |= MM_F_OUT_MD; // --MD
 		else if (c == 0 && long_idx == 14) { // --frag
 			yes_or_no(&opt, MM_F_FRAG_MODE, long_idx, optarg, 1);
 		} else if (c == 0 && long_idx == 15) { // --secondary
@@ -264,6 +273,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, "    --MD         output the MD tag\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");
@@ -276,7 +286,7 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "                 asm5: -k19 -w19 -A1 -B19 -O39,81 -E3,1 -s200 -z200 (asm to ref mapping; break at 5%% div.)\n");
 		fprintf(fp_help, "                 asm10: -k19 -w19 -A1 -B9 -O16,41 -E2,1 -s200 -z200 (asm to ref mapping; break at 10%% div.)\n");
 		fprintf(fp_help, "                 ava-pb: -Hk19 -Xw5 -m100 -g10000 --max-chain-skip 25 (PacBio read overlap)\n");
-		fprintf(fp_help, "                 ava-ont: -k15 -Xw5 -m100 -g10000 --max-chain-skip 25 (ONT read overlap)\n");
+		fprintf(fp_help, "                 ava-ont: -k15 -Xw5 -m100 -g10000 -r2000 --max-chain-skip 25 (ONT read overlap)\n");
 		fprintf(fp_help, "                 splice: long-read spliced alignment (see minimap2.1 for details)\n");
 		fprintf(fp_help, "                 sr: short single-end reads without splicing (see minimap2.1 for details)\n");
 		fprintf(fp_help, "\nSee `man ./minimap2.1' for detailed description of command-line options.\n");
@@ -330,10 +340,17 @@ int main(int argc, char *argv[])
 	}
 	mm_idx_reader_close(idx_rdr);
 
-	fprintf(stderr, "[M::%s] Version: %s\n", __func__, MM_VERSION);
-	fprintf(stderr, "[M::%s] CMD:", __func__);
-	for (i = 0; i < argc; ++i)
-		fprintf(stderr, " %s", argv[i]);
-	fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec\n", __func__, realtime() - mm_realtime0, cputime());
+	if (fflush(stdout) == EOF) {
+		fprintf(stderr, "[ERROR] failed to write the results\n");
+		exit(EXIT_FAILURE);
+	}
+
+	if (mm_verbose >= 3) {
+		fprintf(stderr, "[M::%s] Version: %s\n", __func__, MM_VERSION);
+		fprintf(stderr, "[M::%s] CMD:", __func__);
+		for (i = 0; i < argc; ++i)
+			fprintf(stderr, " %s", argv[i]);
+		fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec\n", __func__, realtime() - mm_realtime0, cputime());
+	}
 	return 0;
 }

map.c

@@ -442,11 +442,12 @@ static void *worker_pipeline(void *shared, int step, void *in)
     pipeline_t *p = (pipeline_t*)shared;
     if (step == 0) { // step 0: read sequences
 		int with_qual = (!!(p->opt->flag & MM_F_OUT_SAM) && !(p->opt->flag & MM_F_NO_QUAL));
+		int with_comment = !!(p->opt->flag & MM_F_COPY_COMMENT);
 		int frag_mode = (p->n_fp > 1 || !!(p->opt->flag & MM_F_FRAG_MODE));
         step_t *s;
         s = (step_t*)calloc(1, sizeof(step_t));
-		if (p->n_fp > 1) s->seq = mm_bseq_read_frag(p->n_fp, p->fp, p->mini_batch_size, with_qual, &s->n_seq);
-		else s->seq = mm_bseq_read2(p->fp[0], p->mini_batch_size, with_qual, frag_mode, &s->n_seq);
+		if (p->n_fp > 1) s->seq = mm_bseq_read_frag2(p->n_fp, p->fp, p->mini_batch_size, with_qual, with_comment, &s->n_seq);
+		else s->seq = mm_bseq_read3(p->fp[0], p->mini_batch_size, with_qual, with_comment, frag_mode, &s->n_seq);
 		if (s->seq) {
 			s->p = p;
 			for (i = 0; i < s->n_seq; ++i)
@@ -489,11 +490,11 @@ static void *worker_pipeline(void *shared, int step, void *in)
 						mm_write_sam2(&p->str, mi, t, i - seg_st, j, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag);
 					else
 						mm_write_paf(&p->str, mi, t, r, km, p->opt->flag);
-					puts(p->str.s);
+					mm_err_puts(p->str.s);
 				}
-				if (s->n_reg[i] == 0 && (p->opt->flag & MM_F_OUT_SAM)) {
+				if (s->n_reg[i] == 0 && (p->opt->flag & MM_F_OUT_SAM)) { // write an unmapped record
 					mm_write_sam2(&p->str, mi, t, i - seg_st, -1, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag);
-					puts(p->str.s);
+					mm_err_puts(p->str.s);
 				}
 			}
 			for (i = seg_st; i < seg_en; ++i) {

minimap.h

@@ -29,6 +29,8 @@
 #define MM_F_REV_ONLY      0x200000
 #define MM_F_HEAP_SORT     0x400000
 #define MM_F_ALL_CHAINS    0x800000
+#define MM_F_OUT_MD        0x1000000
+#define MM_F_COPY_COMMENT  0x2000000
 
 #define MM_I_HPC          0x1
 #define MM_I_NO_SEQ       0x2
@@ -126,6 +128,7 @@ typedef struct {
 	int pe_ori, pe_bonus;
 
 	float mid_occ_frac;  // only used by mm_mapopt_update(); see below
+	int32_t min_mid_occ;
 	int32_t mid_occ;     // ignore seeds with occurrences above this threshold
 	int32_t max_occ;
 	int mini_batch_size; // size of a batch of query bases to process in parallel

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "24 February 2018" "minimap2-2.9 (r720)" "Bioinformatics tools"
+.TH minimap2 1 "27 March 2018" "minimap2-2.10 (r761)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -123,10 +123,27 @@ provided as the target sequences, options
 will be effectively overridden by the options stored in the index file.
 .SS Mapping options
 .TP 10
-.BI -f \ FLOAT
-Ignore top
+.BI -f \ FLOAT | INT1 [, INT2 ]
+If fraction, ignore top
 .I FLOAT
-fraction of most frequent minimizers [0.0002]
+fraction of most frequent minimizers [0.0002]. If integer,
+ignore minimizers occuring more than
+.I INT1
+times.
+.I INT2
+is only effective in the
+.B --sr
+or
+.B -xsr
+mode, which sets the threshold for a second round of seeding.
+.TP
+.BI --min-occ-floor \ INT
+Force minimap2 to always use k-mers occurring
+.I INT
+times or less [0]. In effect, the max occurrence threshold is set to
+the
+.RI max{ INT ,
+.BR -f }.
 .TP
 .BI -g \ INT
 Stop chain enlongation if there are no minimizers within
@@ -353,6 +370,9 @@ SAM read group line in a format like
 .B @RG\\\\tID:foo\\\\tSM:bar
 [].
 .TP
+.B -y
+Copy input FASTA/Q comments to output.
+.TP
 .B -c
 Generate CIGAR. In PAF, the CIGAR is written to the `cg' custom tag.
 .TP
@@ -371,6 +391,9 @@ is given,
 .I short
 is assumed. [none]
 .TP
+.B --MD
+Output the MD tag (see the SAM spec).
+.TP
 .B -Y
 In SAM output, use soft clipping for supplementary alignments.
 .TP
@@ -433,18 +456,25 @@ is determined by the sequencing error mode.
 .B asm5
 Long assembly to reference mapping
 .RB ( -k19
-.B -w19 -A1 -B19 -O39,81 -E3,1 -s200
-.BR -z200 ).
+.B -w19 -A1 -B19 -O39,81 -E3,1 -s200 -z200
+.BR --min-occ-floor=100 ).
 Typically, the alignment will not extend to regions with 5% or higher sequence
 divergence. Only use this preset if the average divergence is far below 5%.
 .TP
 .B asm10
 Long assembly to reference mapping
 .RB ( -k19
-.B -w19 -A1 -B9 -O16,41 -E2,1 -s200
-.BR -z200 ).
+.B -w19 -A1 -B9 -O16,41 -E2,1 -s200 -z200
+.BR --min-occ-floor=100 ).
 Up to 10% sequence divergence.
 .TP
+.B asm20
+Long assembly to reference mapping
+.RB ( -k19
+.B -w10 -A1 -B6 -O6,26 -E2,1 -s200 -z200
+.BR --min-occ-floor=100 ).
+Up to 20% sequence divergence.
+.TP
 .B ava-pb
 PacBio all-vs-all overlap mapping
 .RB ( -Hk19
@@ -454,7 +484,7 @@ PacBio all-vs-all overlap mapping
 .B ava-ont
 Oxford Nanopore all-vs-all overlap mapping
 .RB ( -k15
-.B -Xw5 -m100 -g10000 --max-chain-skip
+.B -Xw5 -m100 -g10000 -r2000 --max-chain-skip
 .BR 25 ).
 Similarly, the major difference from
 .B ava-pb
@@ -535,6 +565,7 @@ cm	i	Number of minimizers on the chain
 s1	i	Chaining score
 s2	i	Chaining score of the best secondary chain
 NM	i	Total number of mismatches and gaps in the alignment
+MD	Z	To generate the ref sequence in the alignment
 AS	i	DP alignment score
 ms	i	DP score of the max scoring segment in the alignment
 nn	i	Number of ambiguous bases in the alignment

misc.c

@@ -1,3 +1,4 @@
+#include <stdlib.h>
 #include "mmpriv.h"
 
 int mm_verbose = 1;
@@ -120,6 +121,16 @@ double realtime(void)
 	return tp.tv_sec + tp.tv_usec * 1e-6;
 }
 
+void mm_err_puts(const char *str)
+{
+	int ret;
+	ret = puts(str);
+	if (ret == EOF) {
+		fprintf(stderr, "[ERROR] failed to write the results\n");
+		exit(EXIT_FAILURE);
+	}
+}
+
 #include "ksort.h"
 
 #define sort_key_128x(a) ((a).x)

misc/paftools.js

@@ -1,6 +1,6 @@
 #!/usr/bin/env k8
 
-var paftools_version = 'r713';
+var paftools_version = 'r755';
 
 /*****************************
  ***** Library functions *****
@@ -131,6 +131,40 @@ Interval.find_ovlp = function(a, st, en)
  * Reverse and reverse complement *
  **********************************/
 
+function fasta_read(fn)
+{
+	var h = {}, gt = '>'.charCodeAt(0);
+	var file = fn == '-'? new File() : new File(fn);
+	var buf = new Bytes(), seq = null, name = null, seqlen = [];
+	while (file.readline(buf) >= 0) {
+		if (buf[0] == gt) {
+			if (seq != null && name != null) {
+				seqlen.push([name, seq.length]);
+				h[name] = seq;
+				name = seq = null;
+			}
+			var m, line = buf.toString();
+			if ((m = /^>(\S+)/.exec(line)) != null) {
+				name = m[1];
+				seq = new Bytes();
+			}
+		} else seq.set(buf);
+	}
+	if (seq != null && name != null) {
+		seqlen.push([name, seq.length]);
+		h[name] = seq;
+	}
+	buf.destroy();
+	file.close();
+	return [h, seqlen];
+}
+
+function fasta_free(fa)
+{
+	for (var name in fa)
+		fa[name].destroy();
+}
+
 Bytes.prototype.reverse = function()
 {
 	for (var i = 0; i < this.length>>1; ++i) {
@@ -305,14 +339,17 @@ function paf_liftover(args)
 // variant calling
 function paf_call(args)
 {
-	var re_cs = /([:=*+-])(\d+|[A-Za-z]+)/g;
+	var re_cs = /([:=*+-])(\d+|[A-Za-z]+)/g, re_tag = /\t(\S\S:[AZif]):(\S+)/g;
 	var c, min_cov_len = 10000, min_var_len = 50000, gap_thres = 50, min_mapq = 5;
-	while ((c = getopt(args, "l:L:g:q:B:")) != null) {
+	var fa_tmp = null, fa, fa_lens, is_vcf = false;
+	while ((c = getopt(args, "l:L:g:q:B:f:")) != null) {
 		if (c == 'l') min_cov_len = parseInt(getopt.arg);
 		else if (c == 'L') min_var_len = parseInt(getopt.arg);
 		else if (c == 'g') gap_thres = parseInt(getopt.arg);
 		else if (c == 'q') min_mapq = parseInt(getopt.arg);
+		else if (c == 'f') fa_tmp = fasta_read(getopt.arg, fa_lens);
 	}
+	if (fa_tmp != null) fa = fa_tmp[0], fa_lens = fa_tmp[1], is_vcf = true;
 
 	if (args.length == getopt.ind) {
 		print("Usage: sort -k6,6 -k8,8n <with-cs.paf> | paftools.js call [options] -");
@@ -321,6 +358,7 @@ function paf_call(args)
 		print("  -L INT    min alignment length to call variants ["+min_var_len+"]");
 		print("  -q INT    min mapping quality ["+min_mapq+"]");
 		print("  -g INT    short/long gap threshold (for statistics only) ["+gap_thres+"]");
+		print("  -f FILE   reference sequences (enabling VCF output) [null]");
 		exit(1);
 	}
 
@@ -328,6 +366,27 @@ function paf_call(args)
 	var buf = new Bytes();
 	var tot_len = 0, n_sub = [0, 0, 0], n_ins = [0, 0, 0, 0], n_del = [0, 0, 0, 0];
 
+	function print_vcf(o, fa)
+	{
+		var v = null;
+		if (o[3] != 1) return; // coverage is one; skip
+		if (o[5] == '-' && o[6] == '-') return;
+		if (o[5] != '-' && o[6] != '-') { // snp
+			v = [o[0], o[1] + 1, '.', o[5].toUpperCase(), o[6].toUpperCase()];
+		} else if (o[1] > 0) { // shouldn't happen in theory
+			if (fa[o[0]] == null) throw Error('sequence "' + o[0] + '" is absent from the reference FASTA');
+			if (o[1] >= fa[o[0]].length) throw Error('position ' + o[1] + ' exceeds the length of sequence "' + o[0] + '"');
+			var ref = String.fromCharCode(fa[o[0]][o[1]-1]).toUpperCase();
+			if (o[5] == '-') // insertion
+				v = [o[0], o[1], '.', ref, ref + o[6].toUpperCase()];
+			else // deletion
+				v = [o[0], o[1], '.', ref + o[5].toUpperCase(), ref];
+		}
+		v.push(o[4], '.', 'QNAME=' + o[7] + ';QSTART=' + (o[8]+1) + ';QSTRAND=' + (rev? '-' : '+'), 'GT', '1/1');
+		if (v == null) throw Error("unexpected variant: [" + o.join(",") + "]");
+		print(v.join("\t"));
+	}
+
 	function count_var(o)
 	{
 		if (o[3] > 1) return;
@@ -346,46 +405,66 @@ function paf_call(args)
 			else ++n_del[3];
 		} else {
 			++n_sub[0];
-			var s = o[5] + o[6];
+			var s = (o[5] + o[6]).toLowerCase();
 			if (s == 'ag' || s == 'ga' || s == 'ct' || s == 'tc')
 				++n_sub[1];
 			else ++n_sub[2];
 		}
 	}
 
+	if (is_vcf) {
+		print('##fileformat=VCFv4.1');
+		for (var i = 0; i < fa_lens.length; ++i)
+			print('##contig=<ID=' + fa_lens[i][0] + ',length=' + fa_lens[i][1] + '>');
+		print('##INFO=<ID=QNAME,Number=1,Type=String,Description="Query name">');
+		print('##INFO=<ID=QSTART,Number=1,Type=Integer,Description="Query start">');
+		print('##INFO=<ID=QSTRAND,Number=1,Type=String,Description="Query strand">');
+		print('##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">');
+		print('#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	sample');
+	}
+
 	var a = [], out = [];
 	var c1_ctg = null, c1_start = 0, c1_end = 0, c1_counted = false, c1_len = 0;
 	while (file.readline(buf) >= 0) {
 		var line = buf.toString();
-		if (!/\ts2:i:/.test(line)) continue; // skip secondary alignments
 		var m, t = line.split("\t", 12);
 		for (var i = 6; i <= 11; ++i)
 			t[i] = parseInt(t[i]);
 		if (t[10] < min_cov_len || t[11] < min_mapq) continue;
-		print(t[0], t[7], t[8], c1_start, c1_end);
+		//print(t[0], t[7], t[8], c1_start, c1_end);
 		for (var i = 1; i <= 3; ++i)
 			t[i] = parseInt(t[i]);
 		var ctg = t[5], x = t[7], end = t[8];
 		var query = t[0], rev = (t[4] == '-'), y = rev? t[3] : t[2];
+		// collect tags
+		var cs = null, tp = null, have_s1 = false, have_s2 = false;
+		while ((m = re_tag.exec(line)) != null) {
+			if (m[1] == 'cs:Z') cs = m[2];
+			else if (m[1] == 'tp:A') tp = m[2];
+			else if (m[1] == 's1:i') have_s1 = true;
+			else if (m[1] == 's2:i') have_s2 = true;
+		}
+		if (have_s1 && !have_s2) continue;
+		if (tp != null && (tp == 'S' || tp == 'i')) continue;
 		// compute regions covered by 1 contig
 		if (ctg != c1_ctg || x >= c1_end) {
 			if (c1_counted && c1_end > c1_start) {
 				c1_len += c1_end - c1_start;
-				print('R', c1_ctg, c1_start, c1_end);
+				if (!is_vcf) print('R', c1_ctg, c1_start, c1_end);
 			}
 			c1_ctg = ctg, c1_start = x, c1_end = end;
 			c1_counted = (t[10] >= min_var_len);
 		} else if (end > c1_end) { // overlap
 			if (c1_counted && x > c1_start) {
 				c1_len += x - c1_start;
-				print('R', c1_ctg, c1_start, x);
+				if (!is_vcf) print('R', c1_ctg, c1_start, x);
 			}
 			c1_start = c1_end, c1_end = end;
 			c1_counted = (t[10] >= min_var_len);
 		} else if (end > c1_start) { // contained
 			if (c1_counted && x > c1_start) {
 				c1_len += x - c1_start;
-				print('R', c1_ctg, c1_start, x);
+				if (!is_vcf) print('R', c1_ctg, c1_start, x);
 			}
 			c1_start = end;
 		} // else, the alignment precedes the cov1 region; do nothing
@@ -393,7 +472,8 @@ function paf_call(args)
 		while (out.length) {
 			if (out[0][0] != ctg || out[0][2] <= x) {
 				count_var(out[0]);
-				print('V', out[0].join("\t"));
+				if (is_vcf) print_vcf(out[0], fa);
+				else print('V', out[0].join("\t"));
 				out.shift();
 			} else break;
 		}
@@ -409,8 +489,7 @@ function paf_call(args)
 		a.length = k;
 		// core loop
 		if (t[10] >= min_var_len) {
-			if ((m = /\tcs:Z:(\S+)/.exec(line)) == null) continue; // no cs tag
-			var cs = m[1];
+			if (cs == null) continue; // no cs tag
 			var blen = 0, n_diff = 0;
 			tot_len += t[10];
 			while ((m = re_cs.exec(cs)) != null) {
@@ -450,11 +529,12 @@ function paf_call(args)
 	}
 	if (c1_counted && c1_end > c1_start) {
 		c1_len += c1_end - c1_start;
-		print('R', c1_ctg, c1_start, c1_end);
+		if (!is_vcf) print('R', c1_ctg, c1_start, c1_end);
 	}
 	while (out.length) {
 		count_var(out[0]);
-		print('V', out[0].join("\t"));
+		if (is_vcf) print_vcf(out[0], fa);
+		else print('V', out[0].join("\t"));
 		out.shift();
 	}
 
@@ -472,6 +552,7 @@ function paf_call(args)
 
 	buf.destroy();
 	file.close();
+	if (fa != null) fasta_free(fa);
 }
 
 function paf_stat(args)
@@ -912,14 +993,15 @@ function paf_view(args)
 
 function paf_gff2bed(args)
 {
-	var c, fn_ucsc_fai = null, is_short = false;
-	while ((c = getopt(args, "u:s")) != null) {
+	var c, fn_ucsc_fai = null, is_short = false, keep_gff = false;
+	while ((c = getopt(args, "u:sg")) != null) {
 		if (c == 'u') fn_ucsc_fai = getopt.arg;
 		else if (c == 's') is_short = true;
+		else if (c == 'g') keep_gff = true;
 	}
 
 	if (getopt.ind == args.length) {
-		print("Usage: paftools.js gff2bed [-u ucsc-genome.fa.fai] <in.gff>");
+		print("Usage: paftools.js gff2bed [-g] [-u ucsc-genome.fa.fai] <in.gff>");
 		exit(1);
 	}
 
@@ -980,6 +1062,12 @@ function paf_gff2bed(args)
 	var exons = [], cds_st = 1<<30, cds_en = 0, last_id = null;
 	while (file.readline(buf) >= 0) {
 		var t = buf.toString().split("\t");
+		if (keep_gff) {
+			if (t[0].charAt(0) != '#' && ens2ucsc[t[0]] != null)
+				t[0] = ens2ucsc[t[0]];
+			print(t.join("\t"));
+			continue;
+		}
 		if (t[0].charAt(0) == '#') continue;
 		if (t[2] != "CDS" && t[2] != "exon") continue;
 		t[3] = parseInt(t[3]) - 1;
@@ -1028,15 +1116,19 @@ function paf_gff2bed(args)
 
 function paf_sam2paf(args)
 {
-	var c, pri_only = false;
+	var c, pri_only = false, use_eq = false;
 	while ((c = getopt(args, "p")) != null)
 		if (c == 'p') pri_only = true;
+	if (args.length == getopt.ind) {
+		print("Usage: paftools.js sam2paf [-p] <in.sam>");
+		exit(1);
+	}
 
-	var file = args.length == getopt.ind || args[getopt.ind] == "-"? new File() : new File(args[getopt.ind]);
+	var file = args[getopt.ind] == "-"? new File() : new File(args[getopt.ind]);
 	var buf = new Bytes();
-	var re = /(\d+)([MIDSHNX=])/g;
+	var re = /(\d+)([MIDSHNX=])/g, re_MD = /(\d+)|(\^[A-Za-z]+)|([A-Za-z])/g, re_tag = /\t(\S\S:[AZif]):(\S+)/g;
 
-	var len = {}, lineno = 0;
+	var ctg_len = {}, lineno = 0;
 	while (file.readline(buf) >= 0) {
 		var m, n_cigar = 0, line = buf.toString();
 		++lineno;
@@ -1044,37 +1136,52 @@ function paf_sam2paf(args)
 			if (/^@SQ/.test(line)) {
 				var name = (m = /\tSN:(\S+)/.exec(line)) != null? m[1] : null;
 				var l = (m = /\tLN:(\d+)/.exec(line)) != null? parseInt(m[1]) : null;
-				if (name != null && l != null) len[name] = l;
+				if (name != null && l != null) ctg_len[name] = l;
 			}
 			continue;
 		}
-		var t = line.split("\t");
+		var t = line.split("\t", 11);
 		var flag = parseInt(t[1]);
-		if (t[9] != '*' && t[10] != '*' && t[9].length != t[10].length) throw Error("ERROR at line " + lineno + ": inconsistent SEQ and QUAL lengths - " + t[9].length + " != " + t[10].length);
-		if (t[2] == '*' || (flag&4)) continue;
+		if (t[9] != '*' && t[10] != '*' && t[9].length != t[10].length)
+			throw Error("at line " + lineno + ": inconsistent SEQ and QUAL lengths - " + t[9].length + " != " + t[10].length);
+		if (t[2] == '*' || (flag&4) || t[5] == '*') continue;
 		if (pri_only && (flag&0x100)) continue;
-		var tlen = len[t[2]];
-		if (tlen == null) throw Error("ERROR at line " + lineno + ": can't find the length of contig " + t[2]);
-		var nn = (m = /\tnn:i:(\d+)/.exec(line)) != null? parseInt(m[1]) : 0;
-		var NM = (m = /\tNM:i:(\d+)/.exec(line)) != null? parseInt(m[1]) : null;
-		var have_NM = NM == null? false : true;
-		NM += nn;
-		var clip = [0, 0], I = [0, 0], D = [0, 0], M = 0, N = 0, ql = 0, tl = 0, mm = 0, ext_cigar = false;
-		while ((m = re.exec(t[5])) != null) {
-			var l = parseInt(m[1]);
-			if (m[2] == 'M') M += l, ql += l, tl += l, ext_cigar = false;
-			else if (m[2] == 'I') ++I[0], I[1] += l, ql += l;
-			else if (m[2] == 'D') ++D[0], D[1] += l, tl += l;
-			else if (m[2] == 'N') N += l, tl += l;
-			else if (m[2] == 'S') clip[M == 0? 0 : 1] = l, ql += l;
-			else if (m[2] == 'H') clip[M == 0? 0 : 1] = l;
-			else if (m[2] == '=') M += l, ql += l, tl += l, ext_cigar = true;
-			else if (m[2] == 'X') M += l, ql += l, tl += l, mm += l, ext_cigar = true;
-			++n_cigar;
+		var tlen = ctg_len[t[2]];
+		if (tlen == null) throw Error("at line " + lineno + ": can't find the length of contig " + t[2]);
+		// find tags
+		var nn = 0, NM = null, MD = null, md_list = [];
+		while ((m = re_tag.exec(line)) != null) {
+			if (m[1] == "NM:i") NM = parseInt(m[2]);
+			else if (m[1] == "nn:i") nn = parseInt(m[2]);
+			else if (m[1] == "MD:Z") MD = m[2];
 		}
+		if (t[9] == '*') MD = null;
+		// infer various lengths from CIGAR
+		var clip = [0, 0], soft_clip = 0, I = [0, 0], D = [0, 0], M = 0, N = 0, mm = 0, have_M = false, have_ext = false, cigar = [];
+		while ((m = re.exec(t[5])) != null) {
+			var l = parseInt(m[1]), op = m[2];
+			if (op == 'M') M += l, have_M = true;
+			else if (op == 'I') ++I[0], I[1] += l;
+			else if (op == 'D') ++D[0], D[1] += l;
+			else if (op == 'N') N += l;
+			else if (op == 'S') clip[n_cigar == 0? 0 : 1] = l, soft_clip += l;
+			else if (op == 'H') clip[n_cigar == 0? 0 : 1] = l;
+			else if (op == '=') M += l, have_ext = true, op = 'M';
+			else if (op == 'X') M += l, mm += l, have_ext = true, op = 'M';
+			++n_cigar;
+			if (MD != null && op != 'H') {
+				if (cigar.length > 0 && cigar[cigar.length-1][1] == op)
+					cigar[cigar.length-1][0] += l;
+				else cigar.push([l, op]);
+			}
+		}
+		var ql = M + I[1] + soft_clip;
+		var tl = M + D[1] + N;
+		var ts = parseInt(t[3]) - 1, te = ts + tl;
+		// checking coordinate and length consistencies
 		if (n_cigar > 65535)
 			warn("WARNING at line " + lineno + ": " + n_cigar + " CIGAR operations");
-		if (tl + parseInt(t[3]) - 1 > tlen) {
+		if (te > tlen) {
 			warn("WARNING at line " + lineno + ": alignment end position larger than ref length; skipped");
 			continue;
 		}
@@ -1082,24 +1189,78 @@ function paf_sam2paf(args)
 			warn("WARNING at line " + lineno + ": SEQ length inconsistent with CIGAR (" + t[9].length + " != " + ql + "); skipped");
 			continue;
 		}
-		if (!have_NM || ext_cigar) NM = I[1] + D[1] + mm;
-		if (NM < I[1] + D[1] + mm) {
-			warn("WARNING at line " + lineno + ": NM is less than the total number of gaps (" + NM + " < " + (I[1]+D[1]+mm) + ")");
-			NM = I[1] + D[1] + mm;
+		// parse MD
+		var cs = [];
+		if (MD != null) {
+			var k = 0, cx = 0, cy = 0, mx = 0, my = 0;
+			while ((m = re_MD.exec(MD)) != null) {
+				if (m[2] != null) { // deletion from the reference
+					var len = m[2].length - 1;
+					cs.push('-', m[2].substr(1));
+					mx += len, cx += len, ++k;
+				} else { // copy or mismatch
+					var ml = m[1] != null? parseInt(m[1]) : 1;
+					while (k < cigar.length && cigar[k][1] != 'D') {
+						var cl = cigar[k][0], op = cigar[k][1];
+						if (op == 'M') {
+							if (my + ml < cy + cl) {
+								if (ml > 0) {
+									if (m[3] != null) cs.push('*', m[3], t[9][my]);
+									else cs.push(':', ml);
+								}
+								mx += ml, my += ml, ml = 0;
+								break;
+							} else {
+								var dl = cy + cl - my;
+								cs.push(':', dl);
+								cx += cl, cy += cl, ++k;
+								mx += dl, my += dl, ml -= dl;
+							}
+						} else if (op == 'I') {
+							cs.push('+', t[9].substr(cy, cl));
+							cy += cl, my += cl, ++k;
+						} else if (op == 'S') {
+							cy += cl, my += cl, ++k;
+						} else throw Error("at line " + lineno + ": inconsistent MD tag");
+					}
+					if (ml != 0) throw Error("at line " + lineno + ": inconsistent MD tag");
+				}
+			}
+			if (cx != mx || cy != my) throw Error("at line " + lineno + ": inconsistent MD tag");
 		}
-		var extra = ["mm:i:"+(NM-I[1]-D[1]), "io:i:"+I[0], "in:i:"+I[1], "do:i:"+D[0], "dn:i:"+D[1]];
-		var match = M - (NM - I[1] - D[1]);
+		// compute matching length, block length and calibrate NM
+		if (have_ext && !have_M) { // extended CIGAR
+			if (NM != null && NM != I[1] + D[1] + mm)
+				warn("WARNING at line " + lineno + ": NM is different from sum of gaps and mismatches");
+			NM = I[1] + D[1] + mm;
+		} else if (NM != null) { // standard CIGAR; NM present
+			if (NM < I[1] + D[1]) {
+				warn("WARNING at line " + lineno + ": NM is less than the total number of gaps (" + NM + " < " + (I[1]+D[1]) + ")");
+				NM = I[1] + D[1];
+			}
+			mm = NM - (I[1] + D[1]);
+		} else { // no way to compute mm
+			warn("WARNING at line " + lineno + ": unable to find the number of mismatches; assuming zero");
+			mm = 0;
+		}
+		var mlen = M - mm;
 		var blen = M + I[1] + D[1];
+		// find query name, start and end
 		var qlen = M + I[1] + clip[0] + clip[1];
-		var qs, qe;
-		if (flag&16) qs = clip[1], qe = qlen - clip[0];
-		else qs = clip[0], qe = qlen - clip[1];
-		var ts = parseInt(t[3]) - 1, te = ts + M + D[1] + N;
-		var qname = t[0];
+		var qname = t[0], qs, qe;
 		if ((flag&1) && (flag&0x40)) qname += '/1';
 		if ((flag&1) && (flag&0x80)) qname += '/2';
-		var a = [qname, qlen, qs, qe, flag&16? '-' : '+', t[2], tlen, ts, te, match, blen, t[4]];
-		print(a.join("\t"), extra.join("\t"));
+		if (flag&16) qs = clip[1], qe = qlen - clip[0];
+		else qs = clip[0], qe = qlen - clip[1];
+		// optional tags
+		var type = flag&0x100? 'S' : 'P';
+		var tags = ["tp:A:" + type];
+		if (NM != null) tags.push("mm:i:"+mm);
+		tags.push("gn:i:"+(I[1]+D[1]), "go:i:"+(I[0]+D[0]), "cg:Z:" + t[5].replace(/\d+[SH]/g, ''));
+		if (cs.length > 0) tags.push("cs:Z:" + cs.join(""));
+		// print out
+		var a = [qname, qlen, qs, qe, flag&16? '-' : '+', t[2], tlen, ts, te, mlen, blen, t[4]];
+		print(a.join("\t"), tags.join("\t"));
 	}
 
 	buf.destroy();
@@ -1597,26 +1758,28 @@ function paf_pbsim2fq(args)
 
 function paf_junceval(args)
 {
-	var c, l_fuzzy = 0, print_ovlp = false, print_err_only = false, first_only = false;
-	while ((c = getopt(args, "l:ep")) != null) {
+	var c, l_fuzzy = 0, print_ovlp = false, print_err_only = false, first_only = false, chr_only = false;
+	while ((c = getopt(args, "l:epc")) != null) {
 		if (c == 'l') l_fuzzy = parseInt(getopt.arg);
 		else if (c == 'e') print_err_only = print_ovlp = true;
 		else if (c == 'p') print_ovlp = true;
+		else if (c == 'c') chr_only = true;
 	}
 
-	if (args.length - getopt.ind < 2) {
+	if (args.length - getopt.ind < 1) {
 		print("Usage: paftools.js junceval [options] <gene.gtf> <aln.sam>");
 		print("Options:");
 		print("  -l INT    tolerance of junction positions (0 for exact) [0]");
 		print("  -p        print overlapping introns");
 		print("  -e        print erroreous overlapping introns");
+		print("  -c        only consider alignments to /^(chr)?([0-9]+|X|Y)$/");
 		exit(1);
 	}
 
 	var file, buf = new Bytes();
 
 	var tr = {};
-	file = new File(args[getopt.ind]);
+	file = args[getopt.ind] == '-'? new File() : new File(args[getopt.ind]);
 	while (file.readline(buf) >= 0) {
 		var m, t = buf.toString().split("\t");
 		if (t[0].charAt(0) == '#') continue;
@@ -1661,13 +1824,14 @@ function paf_junceval(args)
 	var n_pri = 0, n_unmapped = 0, n_mapped = 0;
 	var n_sgl = 0, n_splice = 0, n_splice_hit = 0, n_splice_novel = 0;
 
-	file = new File(args[getopt.ind+1]);
+	file = getopt.ind+1 >= args.length || args[getopt.ind+1] == '-'? new File() : new File(args[getopt.ind+1]);
 	var last_qname = null;
 	var re_cigar = /(\d+)([MIDNSHX=])/g;
 	while (file.readline(buf) >= 0) {
 		var m, t = buf.toString().split("\t");
 
 		if (t[0].charAt(0) == '@') continue;
+		if (chr_only && !/^(chr)?([0-9]+|X|Y)$/.test(t[2])) continue;
 		var flag = parseInt(t[1]);
 		if (flag&0x100) continue;
 		if (first_only && last_qname == t[0]) continue;

mmpriv.h

@@ -86,6 +86,8 @@ mm_seg_t *mm_seg_gen(void *km, uint32_t hash, int n_segs, const int *qlens, int
 void mm_seg_free(void *km, int n_segs, mm_seg_t *segs);
 void mm_pair(void *km, int max_gap_ref, int dp_bonus, int sub_diff, int match_sc, const int *qlens, int *n_regs, mm_reg1_t **regs);
 
+void mm_err_puts(const char *str);
+
 #ifdef __cplusplus
 }
 #endif

options.c

@@ -51,6 +51,8 @@ void mm_mapopt_update(mm_mapopt_t *opt, const mm_idx_t *mi)
 		opt->flag |= MM_F_SPLICE;
 	if (opt->mid_occ <= 0)
 		opt->mid_occ = mm_idx_cal_max_occ(mi, opt->mid_occ_frac);
+	if (opt->mid_occ < opt->min_mid_occ)
+		opt->mid_occ = opt->min_mid_occ;
 	if (mm_verbose >= 3)
 		fprintf(stderr, "[M::%s::%.3f*%.2f] mid_occ = %d\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), opt->mid_occ);
 }
@@ -74,6 +76,7 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		io->flag |= MM_I_HPC, io->k = 19, io->w = 5;
 		mo->flag |= MM_F_ALL_CHAINS | MM_F_NO_DIAG | MM_F_NO_DUAL | MM_F_NO_LJOIN;
 		mo->min_chain_score = 100, mo->pri_ratio = 0.0f, mo->max_gap = 10000, mo->max_chain_skip = 25;
+		mo->bw = 2000;
 	} else if (strcmp(preset, "map10k") == 0 || strcmp(preset, "map-pb") == 0) {
 		io->flag |= MM_I_HPC, io->k = 19;
 	} else if (strcmp(preset, "map-ont") == 0) {
@@ -81,11 +84,19 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 	} else if (strcmp(preset, "asm5") == 0) {
 		io->flag = 0, io->k = 19, io->w = 19;
 		mo->a = 1, mo->b = 19, mo->q = 39, mo->q2 = 81, mo->e = 3, mo->e2 = 1, mo->zdrop = mo->zdrop_inv = 200;
+		mo->min_mid_occ = 100;
 		mo->min_dp_max = 200;
 		mo->best_n = 50;
 	} else if (strcmp(preset, "asm10") == 0) {
 		io->flag = 0, io->k = 19, io->w = 19;
 		mo->a = 1, mo->b = 9, mo->q = 16, mo->q2 = 41, mo->e = 2, mo->e2 = 1, mo->zdrop = mo->zdrop_inv = 200;
+		mo->min_mid_occ = 100;
+		mo->min_dp_max = 200;
+		mo->best_n = 50;
+	} else if (strcmp(preset, "asm20") == 0) {
+		io->flag = 0, io->k = 19, io->w = 10;
+		mo->a = 1, mo->b = 4, mo->q = 6, mo->q2 = 26, mo->e = 2, mo->e2 = 1, mo->zdrop = mo->zdrop_inv = 200;
+		mo->min_mid_occ = 100;
 		mo->min_dp_max = 200;
 		mo->best_n = 50;
 	} else if (strcmp(preset, "short") == 0 || strcmp(preset, "sr") == 0) {

python/cmappy.h

@@ -119,11 +119,11 @@ static char *mappy_fetch_seq(const mm_idx_t *mi, const char *name, int st, int e
 	*len = 0;
 	rid = mm_idx_name2id(mi, name);
 	if (rid < 0) return 0;
-	if (st >= mi->seq[i].len || st >= en) return 0;
-	if (en < 0 || en > mi->seq[i].len)
-		en = mi->seq[i].len;
+	if (st >= mi->seq[rid].len || st >= en) return 0;
+	if (en < 0 || en > mi->seq[rid].len)
+		en = mi->seq[rid].len;
 	s = (char*)malloc(en - st + 1);
-	*len = mm_idx_getseq(mi, rid, st, en, s);
+	*len = mm_idx_getseq(mi, rid, st, en, (uint8_t*)s);
 	for (i = 0; i < *len; ++i)
 		s[i] = "ACGTN"[(uint8_t)s[i]];
 	s[*len] = 0;

python/cmappy.pxd

@@ -34,6 +34,7 @@ cdef extern from "minimap.h":
 		float max_clip_ratio
 		int pe_ori, pe_bonus
 		float mid_occ_frac
+		int32_t min_mid_occ
 		int32_t mid_occ
 		int32_t max_occ
 		int mini_batch_size
@@ -58,7 +59,8 @@ cdef extern from "minimap.h":
 		mm_idx_seq_t *seq
 		uint32_t *S
 		mm_idx_bucket_t *B
-		void *km, *h
+		void *km
+		void *h
 
 	ctypedef struct mm_idx_reader_t:
 		pass

setup.py

@@ -23,7 +23,7 @@ def readme():
 
 setup(
 	name = 'mappy',
-	version = '2.9',
+	version = '2.10',
 	url = 'https://github.com/lh3/minimap2',
 	description = 'Minimap2 python binding',
 	long_description = readme(),
@@ -39,7 +39,7 @@ setup(
 		depends = ['minimap.h', 'bseq.h', 'kalloc.h', 'kdq.h', 'khash.h', 'kseq.h', 'ksort.h',
 				   'ksw2.h', 'kthread.h', 'kvec.h', 'mmpriv.h', 'sdust.h',
 				   'python/cmappy.h', 'python/cmappy.pxd'],
-		extra_compile_args = ['-DHAVE_KALLOC', '-msse4'], # WARNING: ancient x86_64 CPUs don't have SSE4
+		extra_compile_args = ['-DHAVE_KALLOC', '-msse4.1'], # WARNING: ancient x86_64 CPUs don't have SSE4
 		include_dirs = ['.'],
 		libraries = ['z', 'm', 'pthread'])],
 	classifiers = [

test/MT-orang.fa

@@ -1,4 +1,4 @@
->MT_orang
+>MT_orang co:Z:comment
 GTTTATGTAGCTTATTCTATCCAAAGCAATGCACTGAAAATGTCTCGACGGGCCCACACG
 CCCCATAAACAAATAGGTTTGGTCCTAGCCTTTCTATTAGCTCTTAGTGAGGTTACACAT
 GCAAGCATCCCCGCCCCAGTGAGTCGCCCTCCAAGTCACTCTGACTAAGAGGAGCAAGCA

tex/minimap2.bib

@@ -61,13 +61,6 @@
 	Volume = {32},
 	Year = {2016}}
 
-@misc{Suzuki:2016,
-	title = {Fast and accurate alignment tool for PacBio and Nanopore long reads},
-	author = {Hajime Suzuki},
-	journal = {Unpublished},
-	howpublished = {\href{https://github.com/ocxtal/minialign}{https://github.com/ocxtal/minialign}},
-	year = {2016}}
-
 @misc{Ruan:2016,
 	title = {Ultra-fast de novo assembler using long noisy reads},
 	author = {Jue Ruan},
@@ -172,14 +165,6 @@
 	Volume = {29},
 	Year = {2011}}
 
-@article {Suzuki130633,
-	author = {Suzuki, Hajime and Kasahara, Masahiro},
-	title = {Acceleration Of Nucleotide Semi-Global Alignment With Adaptive Banded Dynamic Programming},
-	year = {2017},
-	note = {doi:10.1101/130633},
-	publisher = {Cold Spring Harbor Labs Journals},
-	journal = {bioRxiv}}
-
 @article{Gotoh:1982aa,
 	Author = {Gotoh, O},
 	Journal = {J Mol Biol},
@@ -337,3 +322,19 @@
 	Title = {{MUMmer4}: A fast and versatile genome alignment system},
 	Volume = {14},
 	Year = {2018}}
+
+@article{Li:2009ys,
+	Author = {Li, Heng and others},
+	Journal = {Bioinformatics},
+	Pages = {2078-9},
+	Title = {The {Sequence Alignment/Map format and SAMtools}},
+	Volume = {25},
+	Year = {2009}}
+
+@article{Suzuki:2018aa,
+	Author = {Suzuki, Hajime and Kasahara, Masahiro},
+	Journal = {BMC Bioinformatics},
+	Pages = {45},
+	Title = {Introducing difference recurrence relations for faster semi-global alignment of long sequences},
+	Volume = {19},
+	Year = {2018}}

tex/minimap2.tex

@@ -19,7 +19,7 @@
 \begin{document}
 \firstpage{1}
 
-\title[Aligning nucleotide sequences with minimap2]{Minimap2: versatile pairwise alignment for nucleotide sequences}
+\title[Aligning nucleotide sequences with minimap2]{Minimap2: pairwise alignment for nucleotide sequences}
 \author[Li]{Heng Li}
 \address{Broad Institute, 415 Main Street, Cambridge, MA 02142, USA}
 
@@ -64,7 +64,7 @@ the thought that 10kb long sequences should be easier to map than 100bp reads
 because we can more effectively skip repetitive regions, which are often the
 bottleneck of short-read alignment. We confirmed our speculation by achieving
 approximate mapping 50 times faster than BWA-MEM~\citep{Li:2016aa}.
-\citet{Suzuki130633} extended our work with a fast and novel algorithm on
+\citet{Suzuki:2018aa} extended our work with a fast and novel algorithm on
 generating base-level alignment, which in turn inspired us to develop minimap2
 with added functionality.
 
@@ -88,7 +88,9 @@ the versatility of minimap2.
 
 Minimap2 follows a typical seed-chain-align procedure as is used by most
 full-genome aligners. It collects minimizers~\citep{Roberts:2004fv} of the
-reference sequences and indexes them in a hash table. Then for each query
+reference sequences and indexes them in a hash table, with the key being the
+hash of a minimizer and the value being a list of locations of the minimizer
+copies. Then for each query
 sequence, minimap2 takes query minimizers as \emph{seeds}, finds exact matches
 (i.e. \emph{anchors}) to the reference, and identifies sets of colinear anchors as
 \emph{chains}. If base-level alignment is requested, minimap2 applies dynamic
@@ -118,9 +120,12 @@ distance between two anchors is too large); otherwise
 \begin{equation}\label{eq:chain-gap}
 \beta(j,i)=\gamma_c\big((y_i-y_j)-(x_i-x_j)\big)
 \end{equation}
-In implementation, a gap of length $l\not=0$ costs
+In implementation, a gap of length $l$ costs
 \[
-\gamma_c(l)=0.01\cdot \bar{w}\cdot|l|+0.5\log_2|l|
+\gamma_c(l)=\left\{\begin{array}{ll}
+0.01\cdot \bar{w}\cdot|l|+0.5\log_2|l| & (l\not=0) \\
+0 & (l=0)
+\end{array}\right.
 \]
 where $\bar{w}$ is the average seed length. For $N$ anchors, directly computing all $f(\cdot)$ with
 Eq.~(\ref{eq:chain}) takes $O(N^2)$ time. Although theoretically faster
@@ -164,7 +169,7 @@ empirical formula:
 \[
 {\rm mapQ}=40\cdot (1-f_2/f_1)\cdot\min\{1,m/10\}\cdot\log f_1
 \]
-where $m$ is the number of anchors on the primary chain, $f_1$ is the chaining
+where $\log$ denotes natural logarithm, $m$ is the number of anchors on the primary chain, $f_1$ is the chaining
 score, and $f_2\le f_1$ is the score of the best chain that is secondary to the
 primary chain. Intuitively, a chain is assigned to a higher mapping quality if
 it is long and its best secondary chain is weak.
@@ -253,7 +258,7 @@ performance of minimap2. Traditional SSE implementations~\citep{Farrar:2007hs}
 based on Eq.~(\ref{eq:ae86}) can achieve 16-way parallelization for short
 sequences, but only 4-way parallelization when the peak alignment score reaches
 32767. Long sequence alignment may exceed this threshold. Inspired by
-\citet{Wu:1996aa} and the following work, \citet{Suzuki130633} proposed a
+\citet{Wu:1996aa} and the following work, \citet{Suzuki:2018aa} proposed a
 difference-based formulation that lifted this limitation.
 In case of 2-piece gap cost, define
 \[
@@ -320,7 +325,7 @@ y_{rt}&=&\max\{0,y_{r-1,t}+u_{r-1,t}-z_{rt}+q\}-q-e\\
 \end{equation*}
 In this formulation, cells with the same diagonal index $r$ are independent of
 each other. This allows us to fully vectorize the computation of all cells on
-the same anti-diagonal in one inner loop. It also simplifies banded alignment,
+the same anti-diagonal in one inner loop. It also simplifies banded alignment (500bp band width by default),
 which would be difficult with striped vectorization~\citep{Farrar:2007hs}.
 
 On the condition that $q+e<\tilde{q}+\tilde{e}$ and $e>\tilde{e}$, the initial
@@ -355,7 +360,7 @@ times as fast as Parasail's 4-way vectorization~\citep{Daily:2016aa}.  Without
 banding, our implementation is slower than Edlib~\citep{Sosic:2017aa}, but with
 a 1000bp band, it is considerably faster. When performing global alignment
 between anchors, we expect the alignment to stay close to the diagonal of the
-DP matrix. Banding is applicable most of time.
+DP matrix. Banding is applicable most of the time.
 
 \subsubsection{The Z-drop heuristic}
 
@@ -478,7 +483,7 @@ both C and Python. It is distributed under the MIT license, free to both
 commercial and academic uses. Minimap2 uses the same base algorithm for all
 applications, but it has to apply different sets of parameters depending on
 input data types. Similar to BWA-MEM, minimap2 introduces `presets' that
-modify multiple parameters with a simple invokation. Detailed settings
+modify multiple parameters with a simple invocation. Detailed settings
 and command-line options can be found in the minimap2 manpage. In addition to
 the applications evaluated in the following sections, minimap2 also retains
 minimap's functionality to find overlaps between long reads and to search
@@ -570,7 +575,7 @@ Peak RAM (GByte)          & 8.9      & 14.5     & 3.2      & 29.2\vspace{1em}\\
 \% approx. introns        & 91.8\%   & 96.9\%   & 92.5\%   & 82.4\% \\
 \botrule
 \end{tabular}
-}{Mouse reads (AC:SRR5286960; R9.4 chemistry) were mapped to the primary assembly of mouse
+}{Mouse cDNA reads (AC:SRR5286960; R9.4 chemistry) were mapped to the primary assembly of mouse
 genome GRCm38 with the following tools and command options: minimap2 (`-ax
 splice'); GMAP (`-n 0 --min-intronlength 30 --cross-species'); SpAln (`-Q7 -LS
 -S3'); STARlong (according to
@@ -579,7 +584,7 @@ compared to the EnsEMBL gene annotation, release 89. A predicted intron
 is \emph{novel} if it has no overlaps with any annotated introns. An intron
 is \emph{exact} if it is identical to an annotated intron. An intron is
 \emph{approximate} if both its 5'- and 3'-end are within 10bp around the ends
-of an annotated intron.}
+of an annotated intron. Chimeric alignments are defined in the SAM spec~\citep{Li:2009ys}.}
 \end{table}
 
 We next aligned real mouse reads~\citep{Byrne:2017aa} with GMAP~(v2017-06-20;
@@ -647,9 +652,9 @@ across the whole genome and have been \emph{de novo} assembled with SMRT reads
 to high quality. This allowed us to construct an independent truth variant
 dataset~\citep{Li223297} for
 ERR1341796. In this evaluation, minimap2 has higher SNP false negative rate
-(FNR; 2.5\% of minimap2 vs 2.2\% of BWA-MEM), but fewer false positive SNPs per
-million bases (FPPM; 3.0 vs 3.9), lower 2--50bp INDEL FNR (7.3\% vs 7.5\%) and
-similar INDEL FPPM (both 1.0). Minimap2 is broadly similar to BWA-MEM in the
+(FNR; 2.6\% of minimap2 vs 2.3\% of BWA-MEM), but fewer false positive SNPs per
+million bases (FPPM; 7.0 vs 8.8), similar INDEL FNR (11.2\% vs 11.3\%) and
+similar INDEL FPPM (6.4 vs 6.5). Minimap2 is broadly comparable to BWA-MEM in the
 context of small variant calling.
 
 \subsection{Aligning long-read assemblies}
@@ -687,7 +692,7 @@ involving $>$100kb introns, which was impractically slow ten years ago.  The
 minimap2 chaining algorithm is fast and highly accurate by itself.  In fact,
 chaining alone is more accurate than all the other long-read mappers in
 Fig.~\ref{fig:eval}a (data not shown). This accuracy helps to reduce downstream
-base-level alignment of candidate chains, which is still times slower than
+base-level alignment of candidate chains, which is still several times slower than
 chaining even with the Suzuki-Kasahara improvement. In addition, taking a
 general form, minimap2 chaining can be adapted to non-typical data types such as
 spliced reads and multiple reads per fragment. This gives us the opportunity to
@@ -712,6 +717,8 @@ Schatz, P. Rescheneder and F. Sedlazeck for pointing out the limitation of
 BWA-MEM. We are also grateful to minimap2 users who have greatly helped to
 suggest features and to fix various issues.
 
+\paragraph{Funding\textcolon} NHGRI 1R01HG010040-01
+
 \bibliography{minimap2}
 
 \end{document}