X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bam_md.c;h=63c25e961d704b0c555eaea115fe69f38275e17b;hb=addab0d26aa2f0ea4461bfcacf6304556786408a;hp=17b0a4a4fa6943313be1b6d78b18c5255fe44982;hpb=07c8d1180fa6b712eb2bc6c86ecb9e6b3dbed87b;p=samtools.git

diff --git a/bam_md.c b/bam_md.c
index 17b0a4a..63c25e9 100644
--- a/bam_md.c
+++ b/bam_md.c
@@ -5,6 +5,7 @@
 #include "faidx.h"
 #include "sam.h"
 #include "kstring.h"
+#include "kaln.h"
 
 void bam_fillmd1_core(bam1_t *b, char *ref, int is_equal, int max_nm)
 {
@@ -108,19 +109,68 @@ void bam_fillmd1(bam1_t *b, char *ref, int is_equal)
 	bam_fillmd1_core(b, ref, is_equal, 0);
 }
 
+// local realignment
+
+#define MIN_REF_LEN 10
+
+int bam_realn(bam1_t *b, const char *ref)
+{
+	int k, l_ref, score, n_cigar;
+	uint32_t *cigar = bam1_cigar(b);
+	uint8_t *s_ref = 0, *s_read = 0, *seq;
+	ka_param_t par;
+	bam1_core_t *c = &b->core;
+	// set S/W parameters
+	par = ka_param_blast; par.gap_open = 4; par.gap_ext = 1; par.gap_end_open = par.gap_end_ext = 0;
+	// calculate the length of the reference in the alignment
+	for (k = l_ref = 0; k < c->n_cigar; ++k) {
+		if ((cigar[k]&0xf) == BAM_CREF_SKIP) break; // do not do realignment if there is an `N' operation
+		if ((cigar[k]&0xf) == BAM_CMATCH || (cigar[k]&0xf) == BAM_CDEL)
+			l_ref += cigar[k]>>4;
+	}
+	if (k != c->n_cigar || l_ref < MIN_REF_LEN) return -1;
+	for (k = 0; k < l_ref; ++k)
+		if (ref[c->pos + k] == 0) return -1; // the read stands out of the reference
+	// allocate
+	s_ref = calloc(l_ref, 1);
+	s_read = calloc(c->l_qseq, 1);
+	for (k = 0, seq = bam1_seq(b); k < c->l_qseq; ++k)
+		s_read[k] = bam_nt16_nt4_table[bam1_seqi(seq, k)];
+	for (k = 0; k < l_ref; ++k)
+		s_ref[k] = bam_nt16_nt4_table[(int)bam_nt16_table[(int)ref[c->pos + k]]];
+	// do alignment
+	cigar = ka_global_core(s_ref, l_ref, s_read, c->l_qseq, &par, &score, &n_cigar);
+	if (score <= 0) { // realignment failed
+		free(cigar); free(s_ref); free(s_read);
+	}
+	// copy over the alignment
+	if (4 * (n_cigar - (int)c->n_cigar) + b->data_len > b->m_data) { // enlarge b->data
+		b->m_data = 4 * (n_cigar - (int)c->n_cigar) + b->data_len;
+		kroundup32(b->m_data);
+		b->data = realloc(b->data, b->m_data);
+	}
+	if (n_cigar != (int)c->n_cigar) // move data
+		memmove(b->data + c->l_qname + 4 * n_cigar, bam1_seq(b), b->data_len - c->l_qseq - 4 * c->n_cigar);
+	memcpy(bam1_cigar(b), cigar, n_cigar * 4);
+	c->n_cigar = n_cigar;
+	free(s_ref); free(s_read); free(cigar);
+	return 0;
+}
+
 int bam_fillmd(int argc, char *argv[])
 {
-	int c, is_equal = 0, tid = -2, ret, len, is_bam_out, is_sam_in, is_uncompressed, max_nm = 0;
+	int c, is_equal = 0, tid = -2, ret, len, is_bam_out, is_sam_in, is_uncompressed, max_nm = 0, is_realn;
 	samfile_t *fp, *fpout = 0;
 	faidx_t *fai;
 	char *ref = 0, mode_w[8], mode_r[8];
 	bam1_t *b;
 
-	is_bam_out = is_sam_in = is_uncompressed = 0;
+	is_bam_out = is_sam_in = is_uncompressed = is_realn = 0;
 	mode_w[0] = mode_r[0] = 0;
 	strcpy(mode_r, "r"); strcpy(mode_w, "w");
-	while ((c = getopt(argc, argv, "eubSn:")) >= 0) {
+	while ((c = getopt(argc, argv, "reubSn:")) >= 0) {
 		switch (c) {
+		case 'r': is_realn = 1; break;
 		case 'e': is_equal = 1; break;
 		case 'b': is_bam_out = 1; break;
 		case 'u': is_uncompressed = is_bam_out = 1; break;
@@ -135,11 +185,12 @@ int bam_fillmd(int argc, char *argv[])
 	if (is_uncompressed) strcat(mode_w, "u");
 	if (optind + 1 >= argc) {
 		fprintf(stderr, "\n");
-		fprintf(stderr, "Usage:   samtools fillmd [-eubS] <aln.bam> <ref.fasta>\n\n");
+		fprintf(stderr, "Usage:   samtools fillmd [-eubrS] <aln.bam> <ref.fasta>\n\n");
 		fprintf(stderr, "Options: -e       change identical bases to '='\n");
 		fprintf(stderr, "         -u       uncompressed BAM output (for piping)\n");
 		fprintf(stderr, "         -b       compressed BAM output\n");
-		fprintf(stderr, "         -S       the input is SAM with header\n\n");
+		fprintf(stderr, "         -S       the input is SAM with header\n");
+		fprintf(stderr, "         -r       read-independent local realignment\n\n");
 		return 1;
 	}
 	fp = samopen(argv[optind], mode_r, 0);
@@ -162,6 +213,7 @@ int bam_fillmd(int argc, char *argv[])
 					fprintf(stderr, "[bam_fillmd] fail to find sequence '%s' in the reference.\n",
 							fp->header->target_name[tid]);
 			}
+			if (is_realn) bam_realn(b, ref);
 			if (ref) bam_fillmd1_core(b, ref, is_equal, max_nm);
 		}
 		samwrite(fpout, b);