* aggressive gapped aligner is implemented in calmd.

[samtools.git] / bam_md.c

diff --git a/bam_md.c b/bam_md.c
index 17b0a4a4fa6943313be1b6d78b18c5255fe44982..63c25e961d704b0c555eaea115fe69f38275e17b 100644 (file)
--- 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);