rbh · rbh · Aug 28, 2017 · Mar 31, 2016 · Mar 31, 2016 · Jun 10, 2016
diff --git a/NEWS.md b/NEWS.md
@@ -1,7 +1,24 @@
+Release 0.7.16 (30 July 2017)
+-----------------------------
+
+This release added a couple of minor features and incorporated multiple pull
+requests, including:
+
+ * Added option -5, which is useful to some Hi-C pipelines.
+
+ * Fixed an error with samtools sorting (#129). Updated download link for
+   GRCh38 (#123). Fixed README MarkDown formatting (#70). Addressed multiple
+   issues via a collected pull request #139 by @jmarshall. Avoid malformatted
+   SAM header when -R is used with TAB (#84). Output mate CIGAR (#138).
+
+(0.7.16: 30 July 2017, r1180)
+
+
+
 Release 0.7.15 (31 May 2016)
 ----------------------------
 
-Fixed a long existing bug which potentially leads underestimated insert size
+Fixed a long existing bug which potentially leads to underestimated insert size
 upper bound. This bug should have little effect in practice.
 
 (0.7.15: 31 May 2016, r1140)

diff --git a/README.md b/README.md
@@ -1,14 +1,13 @@
 [![Build Status](https://travis-ci.org/lh3/bwa.svg?branch=dev)](https://travis-ci.org/lh3/bwa)
-[![Build Status](https://drone.io/github.com/lh3/bwa/status.png)](https://drone.io/github.com/lh3/bwa/latest)
-##Getting started
+## Getting started
 
 	git clone https://github.com/lh3/bwa.git
 	cd bwa; make
 	./bwa index ref.fa
 	./bwa mem ref.fa read-se.fq.gz | gzip -3 > aln-se.sam.gz
 	./bwa mem ref.fa read1.fq read2.fq | gzip -3 > aln-pe.sam.gz
 
-##Introduction
+## Introduction
 
 BWA is a software package for mapping DNA sequences against a large reference
 genome, such as the human genome. It consists of three algorithms:
@@ -24,7 +23,7 @@ reference genome (the **index** command). Alignment algorithms are invoked with
 different sub-commands: **aln/samse/sampe** for BWA-backtrack,
 **bwasw** for BWA-SW and **mem** for the BWA-MEM algorithm.
 
-##Availability
+## Availability
 
 BWA is released under [GPLv3][1]. The latest source code is [freely
 available at github][2]. Released packages can [be downloaded][3] at
@@ -37,7 +36,7 @@ In addition to BWA, this self-consistent package also comes with bwa-associated
 and 3rd-party tools for proper BAM-to-FASTQ conversion, mapping to ALT contigs,
 adapter triming, duplicate marking, HLA typing and associated data files.
 
-##Seeking helps
+## Seeking help
 
 The detailed usage is described in the man page available together with the
 source code. You can use `man ./bwa.1` to view the man page in a terminal. The
@@ -46,7 +45,7 @@ have questions about BWA, you may [sign up the mailing list][6] and then send
 the questions to [bio-bwa-help@sourceforge.net][7]. You may also ask questions
 in forums such as [BioStar][8] and [SEQanswers][9].
 
-##Citing BWA
+## Citing BWA
 
 * Li H. and Durbin R. (2009) Fast and accurate short read alignment with
  Burrows-Wheeler transform. *Bioinformatics*, **25**, 1754-1760. [PMID:
@@ -63,7 +62,7 @@ in forums such as [BioStar][8] and [SEQanswers][9].
 Please note that the last reference is a preprint hosted at [arXiv.org][13]. I
 do not have plan to submit it to a peer-reviewed journal in the near future.
 
-##Frequently asked questions (FAQs)
+## Frequently asked questions (FAQs)
 
 1. [What types of data does BWA work with?](#type)
 2. [Why does a read appear multiple times in the output SAM?](#multihit)
@@ -73,7 +72,7 @@ do not have plan to submit it to a peer-reviewed journal in the near future.
 6. [Does BWA work with ALT contigs in the GRCh38 release?](#altctg)
 7. [Can I just run BWA-MEM against GRCh38+ALT without post-processing?](#postalt)
 
-####<a name="type"></a>1. What types of data does BWA work with?
+#### <a name="type"></a>1. What types of data does BWA work with?
 
 BWA works with a variety types of DNA sequence data, though the optimal
 algorithm and setting may vary. The following list gives the recommended
@@ -108,42 +107,42 @@ errors given longer query sequences as the chance of missing all seeds is small.
 As is shown above, with non-default settings, BWA-MEM works with Oxford Nanopore
 reads with a sequencing error rate over 20%.
 
-####<a name="multihit"></a>2. Why does a read appear multiple times in the output SAM?
+#### <a name="multihit"></a>2. Why does a read appear multiple times in the output SAM?
 
 BWA-SW and BWA-MEM perform local alignments. If there is a translocation, a gene
 fusion or a long deletion, a read bridging the break point may have two hits,
 occupying two lines in the SAM output. With the default setting of BWA-MEM, one
 and only one line is primary and is soft clipped; other lines are tagged with
 0x800 SAM flag (supplementary alignment) and are hard clipped.
 
-####<a name="4gb"></a>3. Does BWA work on reference sequences longer than 4GB in total?
+#### <a name="4gb"></a>3. Does BWA work on reference sequences longer than 4GB in total?
 
 Yes. Since 0.6.x, all BWA algorithms work with a genome with total length over
 4GB. However, individual chromosome should not be longer than 2GB.
 
-####<a name="pe0"></a>4. Why can one read in a pair has high mapping quality but the other has zero?
+#### <a name="pe0"></a>4. Why can one read in a pair have a high mapping quality but the other has zero?
 
 This is correct. Mapping quality is assigned for individual read, not for a read
 pair. It is possible that one read can be mapped unambiguously, but its mate
 falls in a tandem repeat and thus its accurate position cannot be determined.
 
-####<a name="endref"></a>5. How can a BWA-backtrack alignment stands out of the end of a chromosome?
+#### <a name="endref"></a>5. How can a BWA-backtrack alignment stand out of the end of a chromosome?
 
 Internally BWA concatenates all reference sequences into one long sequence. A
 read may be mapped to the junction of two adjacent reference sequences. In this
 case, BWA-backtrack will flag the read as unmapped (0x4), but you will see
 position, CIGAR and all the tags. A similar issue may occur to BWA-SW alignment
 as well. BWA-MEM does not have this problem.
 
-####<a name="altctg"></a>6. Does BWA work with ALT contigs in the GRCh38 release?
+#### <a name="altctg"></a>6. Does BWA work with ALT contigs in the GRCh38 release?
 
 Yes, since 0.7.11, BWA-MEM officially supports mapping to GRCh38+ALT.
 BWA-backtrack and BWA-SW don't properly support ALT mapping as of now. Please
 see [README-alt.md][18] for details. Briefly, it is recommended to use
 [bwakit][17], the binary release of BWA, for generating the reference genome
 and for mapping.
 
-####<a name="postalt"></a>7. Can I just run BWA-MEM against GRCh38+ALT without post-processing?
+#### <a name="postalt"></a>7. Can I just run BWA-MEM against GRCh38+ALT without post-processing?
 
 If you are not interested in hits to ALT contigs, it is okay to run BWA-MEM
 without post-processing. The alignments produced this way are very close to

diff --git a/bntseq.c b/bntseq.c
@@ -299,9 +299,9 @@ int64_t bns_fasta2bntseq(gzFile fp_fa, const char *prefix, int for_only)
 	// read sequences
 	while (kseq_read(seq) >= 0) pac = add1(seq, bns, pac, &m_pac, &m_seqs, &m_holes, &q);
 	if (!for_only) { // add the reverse complemented sequence
-		m_pac = (bns->l_pac * 2 + 3) / 4 * 4;
-		pac = realloc(pac, m_pac/4);
-		memset(pac + (bns->l_pac+3)/4, 0, (m_pac - (bns->l_pac+3)/4*4) / 4);
+		int64_t ll_pac = (bns->l_pac * 2 + 3) / 4 * 4;
+		if (ll_pac > m_pac) pac = realloc(pac, ll_pac/4);
+		memset(pac + (bns->l_pac+3)/4, 0, (ll_pac - (bns->l_pac+3)/4*4) / 4);
 		for (l = bns->l_pac - 1; l >= 0; --l, ++bns->l_pac)
 			_set_pac(pac, bns->l_pac, 3-_get_pac(pac, l));
 	}

diff --git a/bwa.1 b/bwa.1
@@ -1,4 +1,4 @@
-.TH bwa 1 "31 May 2016" "bwa-0.7.15-r1140" "Bioinformatics tools"
+.TH bwa 1 "30 July 2017" "bwa-0.7.16-r1180" "Bioinformatics tools"
 .SH NAME
 .PP
 bwa - Burrows-Wheeler Alignment Tool
@@ -107,6 +107,8 @@ appropriate algorithm will be chosen automatically.
 .IR clipPen ]
 .RB [ -U
 .IR unpairPen ]
+.RB [ -x
+.IR readType ]
 .RB [ -R
 .IR RGline ]
 .RB [ -H
@@ -256,7 +258,23 @@ Penalty for an unpaired read pair. BWA-MEM scores an unpaired read pair as
 and scores a paired as scoreRead1+scoreRead2-insertPenalty. It compares these
 two scores to determine whether we should force pairing. A larger value leads to
 more aggressive read pair. [17]
-
+.TP
+.BI -x \ STR
+Read type. Changes multiple parameters unless overriden [null]
+.RS
+.TP 10
+.BR pacbio :
+.B -k17 -W40 -r10 -A1 -B1 -O1 -E1 -L0
+(PacBio reads to ref)
+.TP
+.BR ont2d :
+.B -k14 -W20 -r10 -A1 -B1 -O1 -E1 -L0
+(Oxford Nanopore 2D-reads to ref)
+.TP
+.BR intractg :
+.B -B9 -O16 -L5
+(intra-species contigs to ref)
+.RE
 .TP
 .B INPUT/OUTPUT OPTIONS:
 .TP
@@ -277,6 +295,29 @@ If ARG starts with @, it is interpreted as a string and gets inserted into the
 output SAM header; otherwise, ARG is interpreted as a file with all lines
 starting with @ in the file inserted into the SAM header. [null]
 .TP
+.BI -o \ FILE
+Write the output SAM file to
+.IR FILE .
+For compatibility with other BWA commands, this option may also be given as
+.B -f
+.IR FILE .
+[standard ouptut]
+.TP
+.B -5
+For split alignment, mark the segment with the smallest coordinate as the
+primary. This option may help some Hi-C pipelines. By default, BWA-MEM marks
+highest scoring segment as primary.
+.TP
+.B -K \ INT
+Process
+.I INT
+input bases in each batch regardless of the number of threads in use
+.RI [10000000* nThreads ].
+By default, the batch size is proportional to the number of threads in use.
+Because the inferred insert size distribution slightly depends on the batch
+size, using different number of threads may produce different output.
+Specifying this option helps reproducibility.
+.TP
 .BI -T \ INT
 Don't output alignment with score lower than
 .IR INT .
@@ -302,7 +343,7 @@ Output all found alignments for single-end or unpaired paired-end reads. These
 alignments will be flagged as secondary alignments.
 .TP
 .B -C
-Append append FASTA/Q comment to SAM output. This option can be used to
+Append FASTA/Q comment to SAM output. This option can be used to
 transfer read meta information (e.g. barcode) to the SAM output. Note that the
 FASTA/Q comment (the string after a space in the header line) must conform the SAM
 spec (e.g. BC:Z:CGTAC). Malformated comments lead to incorrect SAM output.
@@ -316,7 +357,7 @@ supplementary alignments.
 Mark shorter split hits as secondary (for Picard compatibility).
 .TP
 .BI -v \ INT
-Control the verbose level of the output. This option has not been fully
+Control the verbosity level of the output. This option has not been fully
 supported throughout BWA. Ideally, a value 0 for disabling all the output to
 stderr; 1 for outputting errors only; 2 for warnings and errors; 3 for
 all normal messages; 4 or higher for debugging. When this option takes value

diff --git a/bwa.c b/bwa.c
@@ -30,13 +30,23 @@ static inline void trim_readno(kstring_t *s)
 		s->l -= 2, s->s[s->l] = 0;
 }
 
+static inline char *dupkstring(const kstring_t *str, int dupempty)
+{
+	char *s = (str->l > 0 || dupempty)? malloc(str->l + 1) : NULL;
+	if (!s) return NULL;
+
+	memcpy(s, str->s, str->l);
+	s[str->l] = '\0';
+	return s;
+}
+
 static inline void kseq2bseq1(const kseq_t *ks, bseq1_t *s)
 { // TODO: it would be better to allocate one chunk of memory, but probably it does not matter in practice
-	s->name = strdup(ks->name.s);
-	s->comment = ks->comment.l? strdup(ks->comment.s) : 0;
-	s->seq = strdup(ks->seq.s);
-	s->qual = ks->qual.l? strdup(ks->qual.s) : 0;
-	s->l_seq = strlen(s->seq);
+	s->name = dupkstring(&ks->name, 1);
+	s->comment = dupkstring(&ks->comment, 0);
+	s->seq = dupkstring(&ks->seq, 1);
+	s->qual = dupkstring(&ks->qual, 0);
+	s->l_seq = ks->seq.l;
 }
 
 bseq1_t *bseq_read(int chunk_size, int *n_, void *ks1_, void *ks2_)
@@ -144,9 +154,9 @@ uint32_t *bwa_gen_cigar2(const int8_t mat[25], int o_del, int e_del, int o_ins,
 		max_del = (int)((double)(((l_query+1)>>1) * mat[0] - o_del) / e_del + 1.);
 		max_gap = max_ins > max_del? max_ins : max_del;
 		max_gap = max_gap > 1? max_gap : 1;
-		w = (max_gap + abs(rlen - l_query) + 1) >> 1;
+		w = (max_gap + abs((int)rlen - l_query) + 1) >> 1;
 		w = w < w_? w : w_;
-		min_w = abs(rlen - l_query) + 3;
+		min_w = abs((int)rlen - l_query) + 3;
 		w = w > min_w? w : min_w;
 		// NW alignment
 		if (bwa_verbose >= 4) {
@@ -414,10 +424,14 @@ char *bwa_set_rg(const char *s)
 		if (bwa_verbose >= 1) fprintf(stderr, "[E::%s] the read group line is not started with @RG\n", __func__);
 		goto err_set_rg;
 	}
+	if (strstr(s, "\t") != NULL) {
+		if (bwa_verbose >= 1) fprintf(stderr, "[E::%s] the read group line contained literal <tab> characters -- replace with escaped tabs: \\t\n", __func__);
+		goto err_set_rg;
+	}
 	rg_line = strdup(s);
 	bwa_escape(rg_line);
 	if ((p = strstr(rg_line, "\tID:")) == 0) {
-		if (bwa_verbose >= 1) fprintf(stderr, "[E::%s] no ID at the read group line\n", __func__);
+		if (bwa_verbose >= 1) fprintf(stderr, "[E::%s] no ID within the read group line\n", __func__);
 		goto err_set_rg;
 	}
 	p += 4;

diff --git a/bwakit/run-bwamem b/bwakit/run-bwamem
@@ -5,7 +5,7 @@ use warnings;
 use Getopt::Std;
 
 my %opts = (t=>1);
-getopts("PSadskHo:R:x:t:", \%opts);
+getopts("MPSadskHo:R:x:t:", \%opts);
 
 die('
 Usage:   run-bwamem [options] <idxbase> <file1> [file2]
@@ -24,6 +24,7 @@ Options: -o STR    prefix for output files                       [inferred from
          -S        for BAM input, don\'t shuffle
          -s        sort the output alignment (via samtools; requring more RAM)
          -k        keep temporary files generated by typeHLA
+         -M        mark shorter split hits as secondary
 
 Examples:
 
@@ -143,7 +144,7 @@ if ($is_bam) {
 	$cmd = "cat $ARGV[1] \\\n";
 }
 
-my $bwa_opts = "-p " . ($opts{t} > 1? "-t$opts{t} " : "") . (defined($opts{x})? "-x $opts{x} " : "") . (defined($opts{R})? "-R'$opts{R}' " : "");
+my $bwa_opts = "-p " . ($opts{t} > 1? "-t$opts{t} " : "") . (defined($opts{x})? "-x $opts{x} " : "") . (defined($opts{R})? "-R'$opts{R}' " : "") . (defined($opts{M})? "-M " : "");
 $bwa_opts .= join(" ", @RG_lines) . " -C " if @RG_lines > 0;
 
 $cmd .= "  | $root/trimadap 2> $prefix.log.trim \\\n" if defined($opts{a});
@@ -163,7 +164,7 @@ if (-f "$ARGV[0].alt" && !defined($opts{P})) {
 }
 
 my $t_sort = $opts{t} < 4? $opts{t} : 4;
-$cmd .= defined($opts{s})? "  | $root/samtools sort -@ $t_sort -m1G - $prefix.aln;\n" : "  | $root/samtools view -1 - > $prefix.aln.bam;\n";
+$cmd .= defined($opts{s})? "  | $root/samtools sort -@ $t_sort -m1G - -o $prefix.aln.bam;\n" : "  | $root/samtools view -1 - > $prefix.aln.bam;\n";
 
 if ($has_hla && defined($opts{H}) && (!defined($opts{x}) || $opts{x} eq 'intractg')) {
 	$cmd .= "$root/run-HLA ". (defined($opts{x}) && $opts{x} eq 'intractg'? "-A " : "") . "$prefix.hla > $prefix.hla.top 2> $prefix.log.hla;\n";

diff --git a/bwakit/run-gen-ref b/bwakit/run-gen-ref
@@ -2,7 +2,7 @@
 
 root=`dirname $0`
 
-url38="ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_full_analysis_set.fna.gz"
+url38="ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_full_analysis_set.fna.gz"
 url37d5="ftp://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz"
 
 if [ $# -eq 0 ]; then

diff --git a/bwamem.c b/bwamem.c
@@ -809,6 +809,19 @@ static inline int get_rlen(int n_cigar, const uint32_t *cigar)
 	return l;
 }
 
+static inline void add_cigar(const mem_opt_t *opt, mem_aln_t *p, kstring_t *str, int which)
+{
+	int i;
+	if (p->n_cigar) { // aligned
+		for (i = 0; i < p->n_cigar; ++i) {
+			int c = p->cigar[i]&0xf;
+			if (!(opt->flag&MEM_F_SOFTCLIP) && !p->is_alt && (c == 3 || c == 4))
+				c = which? 4 : 3; // use hard clipping for supplementary alignments
+			kputw(p->cigar[i]>>4, str); kputc("MIDSH"[c], str);
+		}
+	} else kputc('*', str); // having a coordinate but unaligned (e.g. when copy_mate is true)
+}
+
 void mem_aln2sam(const mem_opt_t *opt, const bntseq_t *bns, kstring_t *str, bseq1_t *s, int n, const mem_aln_t *list, int which, const mem_aln_t *m_)
 {
 	int i, l_name;
@@ -835,14 +848,7 @@ void mem_aln2sam(const mem_opt_t *opt, const bntseq_t *bns, kstring_t *str, bseq
 		kputs(bns->anns[p->rid].name, str); kputc('\t', str); // RNAME
 		kputl(p->pos + 1, str); kputc('\t', str); // POS
 		kputw(p->mapq, str); kputc('\t', str); // MAPQ
-		if (p->n_cigar) { // aligned
-			for (i = 0; i < p->n_cigar; ++i) {
-				int c = p->cigar[i]&0xf;
-				if (!(opt->flag&MEM_F_SOFTCLIP) && !p->is_alt && (c == 3 || c == 4))
-					c = which? 4 : 3; // use hard clipping for supplementary alignments
-				kputw(p->cigar[i]>>4, str); kputc("MIDSH"[c], str);
-			}
-		} else kputc('*', str); // having a coordinate but unaligned (e.g. when copy_mate is true)
+		add_cigar(opt, p, str, which);
 	} else kputsn("*\t0\t0\t*", 7, str); // without coordinte
 	kputc('\t', str);
 
@@ -899,6 +905,7 @@ void mem_aln2sam(const mem_opt_t *opt, const bntseq_t *bns, kstring_t *str, bseq
 		kputsn("\tNM:i:", 6, str); kputw(p->NM, str);
 		kputsn("\tMD:Z:", 6, str); kputs((char*)(p->cigar + p->n_cigar), str);
 	}
+	if (m && m->n_cigar) { kputsn("\tMC:Z:", 6, str); add_cigar(opt, m, str, which); }
 	if (p->score >= 0) { kputsn("\tAS:i:", 6, str); kputw(p->score, str); }
 	if (p->sub >= 0) { kputsn("\tXS:i:", 6, str); kputw(p->sub, str); }
 	if (bwa_rg_id[0]) { kputsn("\tRG:Z:", 6, str); kputs(bwa_rg_id, str); }