r388: cleanup mem_process_seqs() interface

Print output outside the function and allow to feed insert size distribution.

bwamem.c

@@ -951,7 +951,7 @@ static void *worker2(void *data)
 	return 0;
 }
 
-void mem_process_seqs(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t *bns, const uint8_t *pac, int n, bseq1_t *seqs)
+void mem_process_seqs(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t *bns, const uint8_t *pac, int n, bseq1_t *seqs, const mem_pestat_t *pes0)
 {
 	int i;
 	worker_t *w;
@@ -967,29 +967,30 @@ void mem_process_seqs(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t *bn
 		p->seqs = seqs; p->regs = regs;
 		p->pes = &pes[0];
 	}
+
 #ifdef HAVE_PTHREAD
 	if (opt->n_threads == 1) {
+#endif
 		worker1(w);
-		if (opt->flag&MEM_F_PE) mem_pestat(opt, bns->l_pac, n, regs, pes);
+		if (opt->flag&MEM_F_PE) { // paired-end mode
+			if (pes0) memcpy(pes, pes0, 4 * sizeof(mem_pestat_t)); // if pes0 != NULL, set the insert-size distribution as pes0
+			else mem_pestat(opt, bns->l_pac, n, regs, pes); // otherwise, infer the insert size distribution from data
+		}
 		worker2(w);
+#ifdef HAVE_PTHREAD
 	} else {
 		pthread_t *tid;
 		tid = (pthread_t*)calloc(opt->n_threads, sizeof(pthread_t));
 		for (i = 0; i < opt->n_threads; ++i) pthread_create(&tid[i], 0, worker1, &w[i]);
 		for (i = 0; i < opt->n_threads; ++i) pthread_join(tid[i], 0);
-		if (opt->flag&MEM_F_PE) mem_pestat(opt, bns->l_pac, n, regs, pes);
+		if (opt->flag&MEM_F_PE) {
+			if (pes0) memcpy(pes, pes0, 4 * sizeof(mem_pestat_t));
+			else mem_pestat(opt, bns->l_pac, n, regs, pes);
+		}
 		for (i = 0; i < opt->n_threads; ++i) pthread_create(&tid[i], 0, worker2, &w[i]);
 		for (i = 0; i < opt->n_threads; ++i) pthread_join(tid[i], 0);
 		free(tid);
 	}
-#else
-	worker1(w);
-	if (opt->flag&MEM_F_PE) mem_pestat(opt, bns->l_pac, n, regs, pes);
-	worker2(w);
 #endif
-	for (i = 0; i < n; ++i) {
-		fputs(seqs[i].sam, stdout);
-		free(seqs[i].name); free(seqs[i].comment); free(seqs[i].seq); free(seqs[i].qual); free(seqs[i].sam);
-	}
 	free(regs); free(w);
 }

bwamem.h

@@ -57,8 +57,9 @@ typedef struct {
 typedef struct { size_t n, m; mem_alnreg_t *a; } mem_alnreg_v;
 
 typedef struct {
-	int low, high, failed;
-	double avg, std;
+	int low, high;   // lower and upper bounds within which a read pair is considered to be properly paired
+	int failed;      // non-zero if the orientation is not supported by sufficient data
+	double avg, std; // mean and stddev of the insert size distribution
 } mem_pestat_t;
 
 typedef struct { // This struct is only used for the convenience of API.
@@ -103,8 +104,10 @@ extern "C" {
 	 * @param pac    2-bit encoded reference
 	 * @param n      number of query sequences
 	 * @param seqs   query sequences; $seqs[i].seq/sam to be modified after the call
+	 * @param pes0   insert-size info; if NULL, infer from data; if not NULL, it should be an array with 4 elements,
+	 *               corresponding to each FF, FR, RF and RR orientation. See mem_pestat() for more info.
 	 */
-	void mem_process_seqs(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t *bns, const uint8_t *pac, int n, bseq1_t *seqs);
+	void mem_process_seqs(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t *bns, const uint8_t *pac, int n, bseq1_t *seqs, const mem_pestat_t *pes0);
 
 	/**
 	 * Find the aligned regions for one query sequence

fastmap.c

@@ -120,7 +120,11 @@ int main_mem(int argc, char *argv[])
 		for (i = 0; i < n; ++i) size += seqs[i].l_seq;
 		if (bwa_verbose >= 3)
 			fprintf(stderr, "[M::%s] read %d sequences (%ld bp)...\n", __func__, n, (long)size);
-		mem_process_seqs(opt, idx->bwt, idx->bns, idx->pac, n, seqs);
+		mem_process_seqs(opt, idx->bwt, idx->bns, idx->pac, n, seqs, 0);
+		for (i = 0; i < n; ++i) {
+			fputs(seqs[i].sam, stdout);
+			free(seqs[i].name); free(seqs[i].comment); free(seqs[i].seq); free(seqs[i].qual); free(seqs[i].sam);
+		}
 		free(seqs);
 	}
 

main.c

@@ -3,7 +3,7 @@
 #include "utils.h"
 
 #ifndef PACKAGE_VERSION
-#define PACKAGE_VERSION "0.7.4-r386-beta"
+#define PACKAGE_VERSION "0.7.4-r388-beta"
 #endif
 
 int bwa_fa2pac(int argc, char *argv[]);