trinity-plugins/bwa-0.5.7-patched_multi_map/bwtsw2_aux.c
static void bsw2_aln_core(int tid, bsw2seq_t *_seq, const bsw2opt_t *_opt, const bntseq_t *bns, uint8_t *pac, bwt_t * const target[2])
{
	int x;
	bsw2opt_t opt = *_opt;
	bsw2global_t *pool = bsw2_global_init();
	for (x = 0; x < _seq->n; ++x) {
		bsw2seq1_t *p = _seq->seq + x;
		uint8_t *seq[2], *rseq[2];
		int i, l, k;
		bwtsw2_t *b[2];
		l = p->l;

#ifdef HAVE_PTHREAD
		if (_opt->n_threads > 1) {
			pthread_mutex_lock(&g_dbwtsw_lock);
			if (p->tid < 0) p->tid = tid;
			else if (p->tid != tid) {
				pthread_mutex_unlock(&g_dbwtsw_lock);
				continue;
			} // in pinciple else should not happen
			pthread_mutex_unlock(&g_dbwtsw_lock);
		}
#endif

		// set opt->t
		opt.t = _opt->t;
		if (opt.t < log(l) * opt.coef) opt.t = (int)(log(l) * opt.coef + .499);
		if (pool->max_l < l) { // then enlarge working space for aln_extend_core()
			int tmp = ((l + 1) / 2 * opt.a + opt.r) / opt.r + l;
			pool->max_l = l;
			pool->aln_mem = realloc(pool->aln_mem, (tmp + 2) * 24);
		}
		// set opt->bw
		opt.bw = _opt->bw;
		k = (l * opt.a - 2 * opt.q) / (2 * opt.r + opt.a);
		i = (l * opt.a - opt.a - opt.t) / opt.r;
		if (k > i) k = i;
		if (k < 1) k = 1; // I do not know if k==0 causes troubles
		opt.bw = _opt->bw < k? _opt->bw : k;
		// set seq[2] and rseq[2]
		seq[0] = calloc(l * 4, 1);
		seq[1] = seq[0] + l;
		rseq[0] = seq[1] + l; rseq[1] = rseq[0] + l;
		// convert sequences to 2-bit representation
		for (i = k = 0; i < l; ++i) {
			int c = nst_nt4_table[(int)p->seq[i]];
			if (c >= 4) { c = (int)(drand48() * 4); ++k; } // FIXME: ambiguous bases are not properly handled
			seq[0][i] = c;
			seq[1][l-1-i] = 3 - c;
			rseq[0][l-1-i] = c;
			rseq[1][i] = 3 - c;
		}
		if (l - k < opt.t) {
			free(seq[0]); continue;
		}
		// alignment
		b[0] = bsw2_aln1_core(&opt, bns, pac, target[0], l, seq, 0, pool);
		for (k = 0; k < b[0]->n; ++k)
			if (b[0]->hits[k].n_seeds < opt.t_seeds) break;
		if (k < b[0]->n) {
			b[1] = bsw2_aln1_core(&opt, bns, pac, target[1], l, rseq, 1, pool);
			for (i = 0; i < b[1]->n; ++i) {
				bsw2hit_t *p = b[1]->hits + i;
				int x = p->beg;
				p->beg = l - p->end;
				p->end = l - x;
				if (p->l == 0) p->k = bns->l_pac - (p->k + p->len);
			}
			flag_fr(b);
			merge_hits(b, l, 0);
			bsw2_resolve_duphits(0, b[0], 0);
			bsw2_resolve_query_overlaps(b[0], opt.mask_level);
		} else b[1] = 0;
		// generate CIGAR and print SAM
		gen_cigar(&opt, l, seq, pac, b[0]);
		print_hits(bns, &opt, p, b[0]);
		// free
		free(seq[0]);
		bsw2_destroy(b[0]);
	}
	bsw2_global_destroy(pool);
}
/* Core routine to align reads in _seq. It is separated from
 * process_seqs() to realize multi-threading */
Core routine to align reads in _seq. It is separated from
process_seqs() to realize multi-threading
Core routine to align reads in _seq . It is separated from process_seqs () to realize multi - threading
static void bsw2_aln_core(int tid, bsw2seq_t *_seq, const bsw2opt_t *_opt, const bntseq_t *bns, uint8_t *pac, bwt_t * const target[2])
{
	int x;
	bsw2opt_t opt = *_opt;
	bsw2global_t *pool = bsw2_global_init();
	for (x = 0; x < _seq->n; ++x) {
		bsw2seq1_t *p = _seq->seq + x;
		uint8_t *seq[2], *rseq[2];
		int i, l, k;
		bwtsw2_t *b[2];
		l = p->l;
#ifdef HAVE_PTHREAD
		if (_opt->n_threads > 1) {
			pthread_mutex_lock(&g_dbwtsw_lock);
			if (p->tid < 0) p->tid = tid;
			else if (p->tid != tid) {
				pthread_mutex_unlock(&g_dbwtsw_lock);
				continue;
			} 
			pthread_mutex_unlock(&g_dbwtsw_lock);
		}
#endif
		opt.t = _opt->t;
		if (opt.t < log(l) * opt.coef) opt.t = (int)(log(l) * opt.coef + .499);
		if (pool->max_l < l) { 
			int tmp = ((l + 1) / 2 * opt.a + opt.r) / opt.r + l;
			pool->max_l = l;
			pool->aln_mem = realloc(pool->aln_mem, (tmp + 2) * 24);
		}
		opt.bw = _opt->bw;
		k = (l * opt.a - 2 * opt.q) / (2 * opt.r + opt.a);
		i = (l * opt.a - opt.a - opt.t) / opt.r;
		if (k > i) k = i;
		if (k < 1) k = 1; 
		opt.bw = _opt->bw < k? _opt->bw : k;
		seq[0] = calloc(l * 4, 1);
		seq[1] = seq[0] + l;
		rseq[0] = seq[1] + l; rseq[1] = rseq[0] + l;
		for (i = k = 0; i < l; ++i) {
			int c = nst_nt4_table[(int)p->seq[i]];
			if (c >= 4) { c = (int)(drand48() * 4); ++k; } 
			seq[0][i] = c;
			seq[1][l-1-i] = 3 - c;
			rseq[0][l-1-i] = c;
			rseq[1][i] = 3 - c;
		}
		if (l - k < opt.t) {
			free(seq[0]); continue;
		}
		b[0] = bsw2_aln1_core(&opt, bns, pac, target[0], l, seq, 0, pool);
		for (k = 0; k < b[0]->n; ++k)
			if (b[0]->hits[k].n_seeds < opt.t_seeds) break;
		if (k < b[0]->n) {
			b[1] = bsw2_aln1_core(&opt, bns, pac, target[1], l, rseq, 1, pool);
			for (i = 0; i < b[1]->n; ++i) {
				bsw2hit_t *p = b[1]->hits + i;
				int x = p->beg;
				p->beg = l - p->end;
				p->end = l - x;
				if (p->l == 0) p->k = bns->l_pac - (p->k + p->len);
			}
			flag_fr(b);
			merge_hits(b, l, 0);
			bsw2_resolve_duphits(0, b[0], 0);
			bsw2_resolve_query_overlaps(b[0], opt.mask_level);
		} else b[1] = 0;
		gen_cigar(&opt, l, seq, pac, b[0]);
		print_hits(bns, &opt, p, b[0]);
		free(seq[0]);
		bsw2_destroy(b[0]);
	}
	bsw2_global_destroy(pool);
}
static void bsw2_aln_core ( int tid , bsw2seq_t * _seq , const bsw2opt_t * _opt , const bntseq_t * bns , uint8_t * pac , bwt_t * const target [ 2 ] ) { int x ; bsw2opt_t opt = * _opt ; bsw2global_t * pool = bsw2_global_init ( ) ; for ( x = 0 ; x < _seq -> n ; ++ x ) { bsw2seq1_t * p = _seq -> seq + x ; uint8_t * seq [ 2 ] , * rseq [ 2 ] ; int i , l , k ; bwtsw2_t * b [ 2 ] ; l = p -> l ; #ifdef HAVE_PTHREAD if ( _opt -> n_threads > 1 ) { pthread_mutex_lock ( & g_dbwtsw_lock ) ; if ( p -> tid < 0 ) p -> tid = tid ; else if ( p -> tid != tid ) { pthread_mutex_unlock ( & g_dbwtsw_lock ) ; continue ; } pthread_mutex_unlock ( & g_dbwtsw_lock ) ; } #endif opt . t = _opt -> t ; if ( opt . t < log ( l ) * opt . coef ) opt . t = ( int ) ( log ( l ) * opt . coef + .499 ) ; if ( pool -> max_l < l ) { int tmp = ( ( l + 1 ) / 2 * opt . a + opt . r ) / opt . r + l ; pool -> max_l = l ; pool -> aln_mem = realloc ( pool -> aln_mem , ( tmp + 2 ) * 24 ) ; } opt . bw = _opt -> bw ; k = ( l * opt . a - 2 * opt . q ) / ( 2 * opt . r + opt . a ) ; i = ( l * opt . a - opt . a - opt . t ) / opt . r ; if ( k > i ) k = i ; if ( k < 1 ) k = 1 ; opt . bw = _opt -> bw < k ? _opt -> bw : k ; seq [ 0 ] = calloc ( l * 4 , 1 ) ; seq [ 1 ] = seq [ 0 ] + l ; rseq [ 0 ] = seq [ 1 ] + l ; rseq [ 1 ] = rseq [ 0 ] + l ; for ( i = k = 0 ; i < l ; ++ i ) { int c = nst_nt4_table [ ( int ) p -> seq [ i ] ] ; if ( c >= 4 ) { c = ( int ) ( drand48 ( ) * 4 ) ; ++ k ; } seq [ 0 ] [ i ] = c ; seq [ 1 ] [ l - 1 - i ] = 3 - c ; rseq [ 0 ] [ l - 1 - i ] = c ; rseq [ 1 ] [ i ] = 3 - c ; } if ( l - k < opt . t ) { free ( seq [ 0 ] ) ; continue ; } b [ 0 ] = bsw2_aln1_core ( & opt , bns , pac , target [ 0 ] , l , seq , 0 , pool ) ; for ( k = 0 ; k < b [ 0 ] -> n ; ++ k ) if ( b [ 0 ] -> hits [ k ] . n_seeds < opt . t_seeds ) break ; if ( k < b [ 0 ] -> n ) { b [ 1 ] = bsw2_aln1_core ( & opt , bns , pac , target [ 1 ] , l , rseq , 1 , pool ) ; for ( i = 0 ; i < b [ 1 ] -> n ; ++ i ) { bsw2hit_t * p = b [ 1 ] -> hits + i ; int x = p -> beg ; p -> beg = l - p -> end ; p -> end = l - x ; if ( p -> l == 0 ) p -> k = bns -> l_pac - ( p -> k + p -> len ) ; } flag_fr ( b ) ; merge_hits ( b , l , 0 ) ; bsw2_resolve_duphits ( 0 , b [ 0 ] , 0 ) ; bsw2_resolve_query_overlaps ( b [ 0 ] , opt . mask_level ) ; } else b [ 1 ] = 0 ; gen_cigar ( & opt , l , seq , pac , b [ 0 ] ) ; print_hits ( bns , & opt , p , b [ 0 ] ) ; free ( seq [ 0 ] ) ; bsw2_destroy ( b [ 0 ] ) ; } bsw2_global_destroy ( pool ) ; }
Core routine to align reads in _seq.
Core routine to align reads in _seq .
// in pinciple else should not happen // set opt->t // then enlarge working space for aln_extend_core() // set opt->bw // I do not know if k==0 causes troubles // set seq[2] and rseq[2] // convert sequences to 2-bit representation // FIXME: ambiguous bases are not properly handled // alignment // generate CIGAR and print SAM // free