trim.c

#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <stdint.h>
#include <ctype.h>
#include <stdio.h>
#include "kvec.h"
#include "kstring.h"

/********************
 * Global variables *
 ********************/

const char *lt_adapter1 = "AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC"; // Illumina 3'-end adapter
const char *lt_adapter2 = "AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT";

#define AD_MAX_TYPE 40

enum ad_type_e {
	AD_UNKNOWN = 0,
	AD_NO_BARCODE = 11,
	AD_SHORT_SE = 12,
	AD_SHORT_PE = 13,
	AD_COMPLETE_MERGE = 21,
	AD_PARTIAL_MERGE = 22,
	AD_AMBI_MERGE = 31,
	AD_NO_MERGE = 32
};

typedef struct {
	int n_threads;
	int chunk_size;
	int min_seq_len;
	int max_qual;
	int max_ovlp_pen, min_ovlp_len;
	int max_adap_pen, min_adap_len;
	int bc_len[2];
	int trim_len;
	int tab_out, retain_pe;
} lt_opt_t;

static void lt_opt_init(lt_opt_t *opt)
{
	memset(opt, 0, sizeof(lt_opt_t));
	opt->n_threads = 2;
	opt->chunk_size = 10000000;
	opt->max_qual = 50;
	opt->min_seq_len = 30;
	opt->max_ovlp_pen = 3;
	opt->min_ovlp_len = 10;
	opt->max_adap_pen = 2;
	opt->min_adap_len = 5;
	opt->bc_len[0] = opt->bc_len[1] = 0;
	opt->trim_len = 0;
}

/******************
 * K-mer matching *
 ******************/

#include "khash.h"
KHASH_SET_INIT_INT64(s64)
typedef khash_t(s64) lt_seqcloud1_t;

typedef struct {
	int l;
	uint64_t s;
	lt_seqcloud1_t *mm;
} lt_seqcloud_t;

#define KMER_LETTER_NUM_VALUES 5
#define KMER_LETTER_BITWIDTH 3

/* A, a -> 0
 * C, c -> 1
 * G, g -> 2
 * T, t -> 3
 * N, n -> 4
 * all others -> 5
 */
unsigned char seq_nt4_table[256] = {
	5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5, 
	5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5, 
	5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5 /*'-'*/, 5, 5,
	5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5, 
	5, 0, 5, 1,  5, 5, 5, 2,  5, 5, 5, 5,  5, 5, 4, 5, 
	5, 5, 5, 5,  3, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5, 
	5, 0, 5, 1,  5, 5, 5, 2,  5, 5, 5, 5,  5, 5, 4, 5, 
	5, 5, 5, 5,  3, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5, 
	5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5, 
	5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5, 
	5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5, 
	5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5, 
	5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5, 
	5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5, 
	5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5, 
	5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5
};

static lt_seqcloud_t *lt_sc_init(void)
{
	lt_seqcloud_t *sc;
	sc = (lt_seqcloud_t*)calloc(1, sizeof(lt_seqcloud_t));
	sc->mm = kh_init(s64);
	return sc;
}

void lt_sc_destroy(lt_seqcloud_t *sc)
{
	kh_destroy(s64, sc->mm);
	free(sc);
}

void lt_sc_add_core(lt_seqcloud_t *sc, uint64_t s)
{
	int i, absent;
	sc->s = s = s & ((1ULL<<sc->l * KMER_LETTER_BITWIDTH) - 1);
	for (i = 0; i < sc->l; ++i) {
		int i2 = i * KMER_LETTER_BITWIDTH;
		int c = (s >> i2) & 7;
		int a;
		for (a = 1; a < KMER_LETTER_NUM_VALUES; ++a) {
			uint64_t x = (s & ~(7ULL << i2)) | (uint64_t)((a+c) % KMER_LETTER_NUM_VALUES) << i2;
			kh_put(s64, sc->mm, x, &absent);
		}
	}
}

lt_seqcloud_t *lt_sc_gen(const char *s)
{
	lt_seqcloud_t *sc;
	uint64_t x = 0;
	int i;
	sc = lt_sc_init();
	sc->l = strlen(s);
	for (i = 0; s[i] && i < sc->l; ++i) {
		int c = seq_nt4_table[(uint8_t)s[i]];
		if (c >= KMER_LETTER_NUM_VALUES) {
			lt_sc_destroy(sc);
			return 0;
		}
		x = x << KMER_LETTER_BITWIDTH | c;
	}
	lt_sc_add_core(sc, x);
	return sc;
}

typedef struct {
	uint32_t pos:30, type:2;
} lt_sc_hit_t;

int lt_sc_test(const lt_seqcloud_t *sc, const char *seq, int max_hits, lt_sc_hit_t *hits)
{
	int i, l, n = 0;
	uint64_t x = 0, mask = (1ULL << sc->l * KMER_LETTER_BITWIDTH) - 1;
	for (i = l = 0; seq[i]; ++i) {
		int c = seq_nt4_table[(uint8_t)seq[i]];
		if (c < KMER_LETTER_NUM_VALUES) {
			x = (x << KMER_LETTER_BITWIDTH | c) & mask;
			if (++l >= sc->l) {
				if (x == sc->s || kh_get(s64, sc->mm, x) != kh_end(sc->mm)) {
					hits[n].pos = i - (sc->l - 1);
					hits[n++].type = x == sc->s? 0 : 1;
					if (n == max_hits) return n;
				}
			}
		} else l = 0, x = 0;
	}
	return n;
}

/**********************
 * Reverse complement *
 **********************/

char comp_tab[] = {
	  0,   1,	2,	 3,	  4,   5,	6,	 7,	  8,   9,  10,	11,	 12,  13,  14,	15,
	 16,  17,  18,	19,	 20,  21,  22,	23,	 24,  25,  26,	27,	 28,  29,  30,	31,
	 32,  33,  34,	35,	 36,  37,  38,	39,	 40,  41,  42,	43,	 44,  45,  46,	47,
	 48,  49,  50,	51,	 52,  53,  54,	55,	 56,  57,  58,	59,	 60,  61,  62,	63,
	 64, 'T', 'V', 'G', 'H', 'E', 'F', 'C', 'D', 'I', 'J', 'M', 'L', 'K', 'N', 'O',
	'P', 'Q', 'Y', 'S', 'A', 'A', 'B', 'W', 'X', 'R', 'Z',	91,	 92,  93,  94,	95,
	 64, 't', 'v', 'g', 'h', 'e', 'f', 'c', 'd', 'i', 'j', 'm', 'l', 'k', 'n', 'o',
	'p', 'q', 'y', 's', 'a', 'a', 'b', 'w', 'x', 'r', 'z', 123, 124, 125, 126, 127
};

void lt_seq_rev(int l, const char *f, char *r)
{
	int i;
	for (i = 0; i < l; ++i)
		r[l - i - 1] = f[i];
	r[l] = 0;
}

void lt_seq_revcomp(int l, const char *f, char *r)
{
	int i;
	for (i = 0; i < l; ++i)
		r[l - i - 1] = (uint8_t)f[i] >= 128? 'N' : comp_tab[(uint8_t)f[i]];
	r[l] = 0;
}

/**********************
 * Ungapped extension *
 **********************/

#define LT_QUAL_THRES 53 // =33+20
#define LT_HIGH_PEN 3
#define LT_LOW_PEN  1

int lt_ue_for1(int l1, const char *s1, const char *q1, int l2, const char *s2, const char *q2, int min_len, int max_pen)
{
	int i, pen = 0;
	for (i = 0; i < l1 && i < l2; ++i) {
		if (s1[i] != s2[i]) {
			pen += q1[i] >= LT_QUAL_THRES && (q2 == 0 || q2[i] >= LT_QUAL_THRES)? LT_HIGH_PEN : LT_LOW_PEN;
			if (i <= min_len && pen > max_pen) break;
			if (i > min_len && pen * min_len > i * max_pen) break; // in effect: pen > max_pen * ((double)i / min_len)
		}
	}
	return i;
}

int lt_ue_rev1(int l1, const char *s1, const char *q1, int l2, const char *s2, const char *q2, int min_len, int max_pen)
{
	int i, pen = 0;
	for (i = 0; i < l1 && i < l2; ++i) {
		if (s1[l1-1-i] != s2[l2-1-i]) {
			pen += q1[l1-1-i] >= LT_QUAL_THRES && (q2 == 0 || q2[l2-1-i] >= LT_QUAL_THRES)? LT_HIGH_PEN : LT_LOW_PEN;
			if (i <= min_len && pen > max_pen) break;
			if (i > min_len && pen * min_len > i * max_pen) break;
		}
	}
	return i;
}

int lt_ue_for(int l1, const char *s1, const char *q1, int l2, const char *s2, const char *q2, int max_pen, int min_len, int max_pos, uint64_t *pos)
{
	int i, n = 0;
	for (i = min_len; i <= l1; ++i) {
		int l;
		l = lt_ue_for1(i, s1 + l1 - i, q1 + l1 - i, l2, s2, q2, min_len, max_pen);
		if (l >= min_len && (l == i || l == l2)) {
			pos[n++] = (uint64_t)(l1 - i) << 32 | l;
			if (n == max_pos) return n;
		}
	}
	return n;
}

int lt_ue_rev(int l1, const char *s1, const char *q1, int l2, const char *s2, const char *q2, int max_pen, int min_len, int max_pos, uint64_t *pos)
{
	int i, n = 0;
	for (i = min_len; i <= l1; ++i) {
		int l;
		l = lt_ue_rev1(i, s1, q1, l2, s2, q2, min_len, max_pen);
		if (l >= min_len && (l == i || l == l2)) {
			pos[n++] = (uint64_t)(l1 - i) << 32 | l;
			if (n == max_pos) return n;
		}
	}
	return n;
}

int lt_ue_contained(int l1, const char *s1, const char *q1, int l2, const char *s2, const char *q2, int max_pen, int max_pos, uint64_t *pos)
{
	int i, n = 0;
	for (i = 1; i < l2 - l1; ++i) {
		int l;
		l = lt_ue_for1(l1, s1, q1, l2 - i, s2 + i, q2 + i, l1, max_pen);
		if (l == l1) {
			pos[n++] = (uint64_t)i << 32 | l;
			if (n == max_pos) return n;
		}
	}
	return n;
}

/**********************
 * Batch FASTQ reader *
 **********************/

#include <zlib.h>
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

typedef struct {
	uint32_t l_seq:31, dbl_bind:1;
	enum ad_type_e type;
	char *name, *seq, *qual, *bc, *trim;
} bseq1_t;

bseq1_t *bseq_read(kseq_t *ks, int chunk_size, int *n_)
{
	int size = 0, m, n;
	bseq1_t *seqs;
	m = n = 0; seqs = 0;
	while (kseq_read(ks) >= 0) {
		bseq1_t *s;
		if (n >= m) {
			m = m? m<<1 : 256;
			seqs = realloc(seqs, m * sizeof(bseq1_t));
		}
		s = &seqs[n];
		s->name = strdup(ks->name.s);
		s->seq = strdup(ks->seq.s);
		s->qual = ks->qual.l? strdup(ks->qual.s) : 0;
		s->bc = 0;
		s->l_seq = ks->seq.l;
		s->dbl_bind = 0;
		s->type = AD_UNKNOWN;
		size += seqs[n++].l_seq;
		if (size >= chunk_size && (n&1) == 0) break;
	}
	*n_ = n;
	return seqs;
}

/****************
 * Read barcode *
 ****************/

typedef struct {
	int n[2], len[2];
	lt_seqcloud_t **bc[2];
	char **seq[2];
} ad_barcode_t;

typedef char *cstr_t;

ad_barcode_t *ad_bc_read(const char *fn)
{
	gzFile fp;
	kstream_t *ks;
	kstring_t s = {0,0,0};
	ad_barcode_t *bc = 0;
	int dret, i, len[2];
	kvec_t(cstr_t) x[2];

	fp = gzopen(fn, "r");
	if (fp == 0) return 0;
	ks = ks_init(fp);
	kv_init(x[0]); kv_init(x[1]);
	len[0] = len[1] = 0;
	for (i = 0; i < 2; ++i) {
		int j, j0;
		if (ks_getuntil(ks, KS_SEP_LINE, &s, &dret) < 0) break;
		for (j = j0 = 0; j <= s.l; ++j) {
			if (j == s.l || ispunct(s.s[j])) {
				s.s[j] = 0;
				if (len[i] == 0) len[i] = j - j0;
				else if (len[i] != j - j0) {
					j0 = j + 1;
					continue;
				}
				if (len[i] > 0) kv_push(cstr_t, x[i], strdup(&s.s[j0]));
				j0 = j + 1;
			}
		}
	}
	free(s.s);
	ks_destroy(ks);
	gzclose(fp);

	fprintf(stderr, "[M::%s] read %ld read1 barcodes and %ld read2 barcodes\n", __func__, x[0].n, x[1].n);

	bc = (ad_barcode_t*)calloc(1, sizeof(ad_barcode_t));
	for (i = 0; i < 2; ++i) {
		int j;
		bc->n[i] = x[i].n, bc->len[i] = len[i];
		bc->bc[i] = (lt_seqcloud_t**)calloc(bc->n[i], sizeof(lt_seqcloud_t*));
		for (j = 0; j < x[i].n; ++j)
			bc->bc[i][j] = lt_sc_gen(x[i].a[j]);
		bc->seq[i] = x[i].a;
	}
	return bc;
}

/*********************************
 * Core trimming/merging routine *
 *********************************/

typedef struct {
	lt_opt_t opt;
	kseq_t *ks;
	gzFile fp_pe[2];
	uint64_t types[AD_MAX_TYPE+1];
	ad_barcode_t *bc;
} lt_global_t;

void lt_global_init(lt_global_t *g)
{
	memset(g, 0, sizeof(lt_global_t));
	lt_opt_init(&g->opt);
}

#define MAX_BINDING_HITS 3

static inline void trim_bseq_5(bseq1_t *s, int l)
{
	memmove(s->seq, s->seq + l, s->l_seq - l);
	memmove(s->qual, s->qual + l, s->l_seq - l);
	s->l_seq -= l;
	s->seq[s->l_seq] = s->qual[s->l_seq] = 0;
}

static inline int merge_base(int max_qual, char fc, char fq, char rc, char rq)
{
	int y;
	if (fc == rc) {
		int q = fq > rq? fq - 33 : rq - 33;
		y = toupper(fc) | (33 + (q < max_qual? q : max_qual)) << 8;
	} else {
		if (fq > rq) y = toupper(fc) | (33 + (fq - rq)) << 8;
		else y = toupper(rc) | (33 + (rq - fq)) << 8;
	}
	return y;
}

static inline void trim_adap(bseq1_t *s, const char *adap, int is_5, int min_len, int max_pen, int allow_contained)
{
	int n_hits, l_adap;
	uint64_t hits[4];
	l_adap = strlen(adap);
	if (is_5) n_hits = lt_ue_rev(s->l_seq, s->seq, s->qual, l_adap, adap, 0, max_pen, min_len, 4, hits);
	else n_hits = lt_ue_for(s->l_seq, s->seq, s->qual, l_adap, adap, 0, max_pen, min_len, 4, hits);
	if (n_hits > 0 && (allow_contained || (hits[0]>>32) + (uint32_t)hits[0] == s->l_seq || (hits[n_hits-1]>>32) + (uint32_t)hits[n_hits-1] == s->l_seq)) {
		int len = s->l_seq - (hits[n_hits-1]>>32); // trim the longest hit
		if (is_5) {
			if (len > min_len) trim_bseq_5(s, len); // trim
			else memset(s->qual, 33+1, len); // reduce baseQ
		} else {
			if (len > min_len) s->l_seq -= len, s->seq[s->l_seq] = s->qual[s->l_seq] = 0; // trim
			else memset(s->qual + (s->l_seq - len), 33+1, len); // reduce baseQ
		}
	}
}

void lt_process(const lt_global_t *g, bseq1_t s[2])
{
	int i, k, mlen, bc_id[2];
	char *rseq, *rqual, *xseq, *xqual;

	mlen = s[0].l_seq > s[1].l_seq? s[0].l_seq : s[1].l_seq;
	rseq = (char*)alloca(mlen + 1);
	rqual = (char*)alloca(mlen + 1);
	xseq = (char*)alloca(s[0].l_seq + s[1].l_seq + 1);
	xqual = (char*)alloca(s[0].l_seq + s[1].l_seq + 1);

	// trim trailing N
	for (k = 0; k < 2; ++k) {
		bseq1_t *sk = &s[k];
		for (i = sk->l_seq - 1; i >= 0; --i) // trim trailing "N"
			if (sk->seq[i] != 'N') break;
		sk->l_seq = i + 1;
		sk->seq[sk->l_seq] = sk->qual[sk->l_seq] = 0;
	}
	// test barcode
	if (g->bc && g->bc->n[0] + g->bc->n[1] > 0) {
		char *bc;
		bc = (char*)alloca(g->bc->len[0] + g->bc->len[1] + 1);
		for (k = 0; k < 2; ++k) {
			lt_sc_hit_t hits[3];
			strncpy(bc, s[k].seq, g->bc->len[k]);
			bc[g->bc->len[k]] = 0;
			for (i = 0; i < g->bc->n[k]; ++i)
				if (lt_sc_test(g->bc->bc[k][i], bc, 3, hits))
					break;
			if (i == g->bc->n[k]) break;
			bc_id[k] = i;
		}
		if (k < 2) {
			s[0].type = s[1].type = AD_NO_BARCODE;
			return;
		}
		// write barcode
		strncpy(bc, g->bc->seq[0][bc_id[0]], g->bc->len[0]);
		strncpy(bc + g->bc->len[0], g->bc->seq[1][bc_id[1]], g->bc->len[1]);
		bc[g->bc->len[0] + g->bc->len[1]] = 0;
		s->bc = strdup(bc);
		if (g->bc->len[0]) trim_bseq_5(&s[0], g->bc->len[0]);
		if (g->bc->len[1]) trim_bseq_5(&s[1], g->bc->len[1]);
	} else if (g->opt.bc_len[0] + g->opt.bc_len[1] > 0) {
		char *bc;
		bc = (char*)alloca(g->opt.bc_len[0] + g->opt.bc_len[1] + 1);
		strncpy(bc, s[0].seq, g->opt.bc_len[0]);
		strncpy(bc + g->opt.bc_len[0], s[1].seq, g->opt.bc_len[1]);
		bc[g->opt.bc_len[0] + g->opt.bc_len[1]] = 0;
		s->bc = strdup(bc);
		if (g->opt.bc_len[0]) trim_bseq_5(&s[0], g->opt.bc_len[0]);
		if (g->opt.bc_len[1]) trim_bseq_5(&s[1], g->opt.bc_len[1]);
	}
	// find end overlaps
	{
		int n_fh, n_rh, n_ch;
		uint64_t fh[2], rh[2], ch[2];
		// reverse the other read
		lt_seq_revcomp(s[1].l_seq, s[1].seq, rseq);
		lt_seq_rev(s[1].l_seq, s[1].qual, rqual);
		// find overlaps
		n_fh = lt_ue_for(s[0].l_seq, &s[0].seq[0], &s[0].qual[0], s[1].l_seq, rseq, rqual, g->opt.max_ovlp_pen, g->opt.min_ovlp_len, 2, fh);
		n_rh = lt_ue_rev(s[0].l_seq, &s[0].seq[0], &s[0].qual[0], s[1].l_seq, rseq, rqual, g->opt.max_ovlp_pen, g->opt.min_ovlp_len, 2, rh);
		n_ch = lt_ue_contained(s[0].l_seq, &s[0].seq[0], &s[0].qual[0], s[1].l_seq, rseq, rqual, g->opt.max_ovlp_pen, 2, ch);
		if (n_fh + n_rh + n_ch > 1) {
			s[0].type = s[1].type = AD_AMBI_MERGE;
		} else if (n_fh + n_rh + n_ch == 0) {
			s[0].type = s[1].type = AD_NO_MERGE;
		} else {
			int x = 0;
			if (n_fh == 1) {
				int l = (uint32_t)fh[0], st = fh[0]>>32;
				s[0].type = s[1].type = AD_PARTIAL_MERGE;
				for (i = 0; i < st; ++i)
					xseq[x] = s[0].seq[i], xqual[x++] = s[0].qual[i];
				for (i = 0; i < l; ++i) {
					int j = st + i, y;
					y = merge_base(g->opt.max_qual, s[0].seq[j], s[0].qual[j], rseq[i], rqual[i]);
					xseq[x] = (uint8_t)y, xqual[x++] = y>>8;
				}
				if (l < s[1].l_seq) {
					for (i = l; i < s[1].l_seq; ++i)
						xseq[x] = rseq[i], xqual[x++] = rqual[i];
				} else {
					for (i = l; i < s[0].l_seq; ++i)
						xseq[x] = s[0].seq[i], xqual[x++] = s[0].qual[i];
				}
			} else if (n_rh == 1) {
				int l = (uint32_t)rh[0], st = rh[0]>>32;
				s[0].type = s[1].type = AD_COMPLETE_MERGE;
				for (i = 0; i < s[0].l_seq - st - l; ++i)
					xseq[x] = s[0].seq[i], xqual[x++] = s[0].qual[i];
				for (i = s[1].l_seq - l; i < s[1].l_seq; ++i) {
					int j = i - s[1].l_seq + (s[0].l_seq - st), y;
					y = merge_base(g->opt.max_qual, s[0].seq[j], s[0].qual[j], rseq[i], rqual[i]);
					xseq[x] = (uint8_t)y, xqual[x++] = y>>8;
				}
			} else {
				int j, st = ch[0]>>32;
				s[0].type = s[1].type = AD_COMPLETE_MERGE;
				for (j = 0; j < s[0].l_seq; ++j) {
					int i = j + st, y;
					y = merge_base(g->opt.max_qual, s[0].seq[j], s[0].qual[j], rseq[i], rqual[i]);
					xseq[x] = (uint8_t)y, xqual[x++] = y>>8;
				}
			}
			xseq[x] = xqual[x] = 0;
			if (x < g->opt.min_seq_len) s[0].type = s[1].type = AD_SHORT_SE;
			free(s[0].seq); free(s[0].qual);
			s[0].seq = strdup(xseq);
			s[0].qual = strdup(xqual);
			s[0].l_seq = x;
			s[1].l_seq = 0;
		}
	}
	if (s->type == AD_NO_MERGE || s->type == AD_AMBI_MERGE) {
		if (s->bc) s[1].bc = strdup(s->bc);
		// trim Illumina PE adapters
		trim_adap(&s[0], lt_adapter1, 0, g->opt.min_adap_len, g->opt.max_adap_pen, 1);
		trim_adap(&s[1], lt_adapter2, 0, g->opt.min_adap_len, g->opt.max_adap_pen, 1);
		if (s[0].l_seq < g->opt.min_seq_len || s[1].l_seq < g->opt.min_seq_len)
			s[0].type = s[1].type = AD_SHORT_PE;
	}
	if (g->opt.trim_len) {
		if (s->type == AD_PARTIAL_MERGE || s->type == AD_COMPLETE_MERGE) {
			s->l_seq -= g->opt.trim_len;
			s->seq[s->l_seq] = s->qual[s->l_seq] = 0;
			trim_bseq_5(s, g->opt.trim_len);
		} else if (s->type == AD_NO_MERGE || s->type == AD_AMBI_MERGE) {
			trim_bseq_5(&s[0], g->opt.trim_len);
			trim_bseq_5(&s[1], g->opt.trim_len);
		}
	}
}

/**********************
 * Callback functions *
 **********************/

void kt_for(int n_threads, void (*func)(void*,long,int), void *data, long n);
void kt_pipeline(int n_threads, void *(*func)(void*, int, void*), void *shared_data, int n_steps);

typedef struct {
	int n_seqs;
	bseq1_t *seqs;
	lt_global_t *g;
} data_for_t;

static void worker_for(void *_data, long i, int tid)
{
	data_for_t *data = (data_for_t*)_data;
	lt_process(data->g, &data->seqs[i<<1]);
}

static void *worker_pipeline(void *shared, int step, void *_data)
{
	int i;
	lt_global_t *g = (lt_global_t*)shared;
	if (step == 0) {
		data_for_t *ret;
		ret = calloc(1, sizeof(data_for_t));
		ret->seqs = bseq_read(g->ks, g->opt.chunk_size, &ret->n_seqs);
		assert((ret->n_seqs&1) == 0);
		ret->g = g;
		if (ret->seqs) return ret;
		else free(ret);
	} else if (step == 1) {
		data_for_t *data = (data_for_t*)_data;
		kt_for(g->opt.n_threads, worker_for, data, data->n_seqs>>1);
		return data;
	} else if (step == 2) {
		data_for_t *data = (data_for_t*)_data;
		if (g->opt.tab_out) { // tabular output
			for (i = 0; i < data->n_seqs; i += 2) {
				bseq1_t *s = &data->seqs[i];
				printf("%s\t%d\n", s->name, s->type);
			}
		} else { // FASTQ output (FASTA not supported yet)
			kstring_t pe[2] = {{0,0,0},{0,0,0}};
			for (i = 0; i < data->n_seqs; i += 2)
				++g->types[data->seqs[i].type];
			for (i = 0; i < data->n_seqs; ++i) {
				bseq1_t *s = &data->seqs[i];
				if (s->l_seq <= 0) continue;
				if (g->fp_pe[0] && g->fp_pe[1]) {
					if (s->type == AD_PARTIAL_MERGE || s->type == AD_COMPLETE_MERGE) {
						printf("%c%s_%d:%s\n%s\n", s->qual? '@' : '>', s->name, s->type, s->bc? s->bc : "*", s->seq);
						if (s->qual) { puts("+"); puts(s->qual); }
					} else if (s->type == AD_NO_MERGE || s->type == AD_AMBI_MERGE) {
						kstring_t *str = &pe[i&1];
						ksprintf(str, "%c%s_%d:%s\n%s\n", s->qual? '@' : '>', s->name, s->type, s->bc? s->bc : "*", s->seq);
						if (s->qual) { kputs("+\n", str); kputs(s->qual, str); kputc('\n', str); }
					}
				} else {
					if (s->type == AD_PARTIAL_MERGE || s->type == AD_COMPLETE_MERGE || s->type == AD_NO_MERGE || s->type == AD_AMBI_MERGE) {
						if (!g->opt.retain_pe && (s->type == AD_NO_MERGE || s->type == AD_AMBI_MERGE)) continue;
						printf("%c%s_%d:%s", s->qual? '@' : '>', s->name, s->type, s->bc? s->bc : "*");
						if (s->type != AD_PARTIAL_MERGE && s->type != AD_COMPLETE_MERGE) {
							putchar('/'); putchar("12"[i&1]);
						}
						putchar('\n');
						puts(s->seq);
						if (s->qual) { puts("+"); puts(s->qual); }
					}
				}
			}
			for (i = 0; i < 2; ++i) {
				gzwrite(g->fp_pe[i], pe[i].s, pe[i].l);
				free(pe[i].s);
			}
		}
		for (i = 0; i < data->n_seqs; ++i) { // deallocate
			bseq1_t *s = &data->seqs[i];
			free(s->bc); free(s->seq); free(s->qual); free(s->name);
		}
		free(data->seqs); free(data);
	}
	return 0;
}

#include <unistd.h>

int main(int argc, char *argv[])
{
	int c;
	lt_global_t g;
	gzFile fp;
	char *p, *pe_prefix = 0, *fn_bc = 0;

	lt_global_init(&g);
	while ((c = getopt(argc, argv, "dt:b:l:o:p:B:T:P")) >= 0) {
		if (c == 't') g.opt.n_threads = atoi(optarg);
		else if (c == 'd') g.opt.tab_out = 1;
		else if (c == 'P') g.opt.retain_pe = 1;
		else if (c == 'l') g.opt.min_seq_len = atoi(optarg);
		else if (c == 'o') g.opt.min_ovlp_len = atoi(optarg);
		else if (c == 'p') pe_prefix = optarg;
		else if (c == 'b') fn_bc = optarg;
		else if (c == 'T') g.opt.trim_len = atoi(optarg);
		else if (c == 'B') {
			g.opt.bc_len[0] = g.opt.bc_len[1] = strtol(optarg, &p, 10);
			if (*p == ',') g.opt.bc_len[1] = strtol(p+1, &p, 10);
		}
	}
	if (argc - optind < 1) {
		fprintf(stderr, "Usage: seqtk mergepe <read1.fq> <read2.fq> | adna-trim [options] -\n");
		fprintf(stderr, "Options:\n");
		fprintf(stderr, "  -b FILE      barcode file (overriding -B) []\n");
		fprintf(stderr, "  -B INT[,INT] barcode length [%d,%d]\n", g.opt.bc_len[0], g.opt.bc_len[1]);
		fprintf(stderr, "  -t INT       number of threads [%d]\n", g.opt.n_threads);
		fprintf(stderr, "  -l INT       min read/fragment length to output [%d]\n", g.opt.min_seq_len);
		fprintf(stderr, "  -o INT       min overlap length [%d]\n", g.opt.min_ovlp_len);
		fprintf(stderr, "  -p STR       output PE reads to STR.R[12].fq.gz (overriding -P) [discard pe]\n");
		fprintf(stderr, "  -T INT       trim INT-bp [%d]\n", g.opt.trim_len);
		fprintf(stderr, "  -P           output paired reads also to stdout [discard pe]\n");
		fprintf(stderr, "  -d           tabular output for debugging\n");
		return 1;
	}

	fp = strcmp(argv[optind], "-")? gzopen(argv[optind], "r") : gzdopen(fileno(stdin), "r");
	g.ks = kseq_init(fp);
	if (fn_bc) g.bc = ad_bc_read(fn_bc), g.opt.bc_len[0] = g.opt.bc_len[1] = 0;
	if (pe_prefix) {
		kstring_t str = {0,0,0};
		for (c = 0; c < 2; ++c) {
			str.l = 0;
			ksprintf(&str, "%s.R%d.fq.gz", pe_prefix, c+1);
			g.fp_pe[c] = gzopen(str.s, "w1");
		}
		g.opt.retain_pe = 0;
	}

	kt_pipeline(2, worker_pipeline, &g, 3);
	for (c = 0; c <= AD_MAX_TYPE; ++c)
		if (g.types[c] > 0)
			fprintf(stderr, "[M::%s] %2d  %ld\n", __func__, c, (long)g.types[c]);

	if (g.fp_pe[0]) gzclose(g.fp_pe[0]);
	if (g.fp_pe[1]) gzclose(g.fp_pe[1]);
	kseq_destroy(g.ks);
	gzclose(fp);
	return 0;
}