* samtools-0.1.9-8 (r807)

[samtools.git] / bam2bcf_indel.c

diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c
index de9df3e6c2fb9042b4287ac60e828224aaa651a2..2d0b5c6f8a66a4c5c006dbec560718ab1f7a0377 100644 (file)
--- a/bam2bcf_indel.c
+++ b/bam2bcf_indel.c
@@ -1,13 +1,61 @@
 #include <assert.h>
 #include <ctype.h>
+#include <string.h>
 #include "bam.h"
 #include "bam2bcf.h"
 #include "ksort.h"
 #include "kaln.h"
+#include "khash.h"
+KHASH_SET_INIT_STR(rg)
 
 #define MINUS_CONST 0x10000000
 #define INDEL_WINDOW_SIZE 50
 
+void *bcf_call_add_rg(void *_hash, const char *hdtext, const char *list)
+{
+       const char *s, *p, *q, *r, *t;
+       khash_t(rg) *hash;
+       if (list == 0 || hdtext == 0) return _hash;
+       if (_hash == 0) _hash = kh_init(rg);
+       hash = (khash_t(rg)*)_hash;
+       if ((s = strstr(hdtext, "@RG\t")) == 0) return hash;
+       do {
+               t = strstr(s + 4, "@RG\t"); // the next @RG
+               if ((p = strstr(s, "\tID:")) != 0) p += 4;
+               if ((q = strstr(s, "\tPL:")) != 0) q += 4;
+               if (p && q && (t == 0 || (p < t && q < t))) { // ID and PL are both present
+                       int lp, lq;
+                       char *x;
+                       for (r = p; *r && *r != '\t' && *r != '\n'; ++r); lp = r - p;
+                       for (r = q; *r && *r != '\t' && *r != '\n'; ++r); lq = r - q;
+                       x = calloc((lp > lq? lp : lq) + 1, 1);
+                       for (r = q; *r && *r != '\t' && *r != '\n'; ++r) x[r-q] = *r;
+                       if (strstr(list, x)) { // insert ID to the hash table
+                               khint_t k;
+                               int ret;
+                               for (r = p; *r && *r != '\t' && *r != '\n'; ++r) x[r-p] = *r;
+                               x[r-p] = 0;
+                               k = kh_get(rg, hash, x);
+                               if (k == kh_end(hash)) k = kh_put(rg, hash, x, &ret);
+                               else free(x);
+                       } else free(x);
+               }
+               s = t;
+       } while (s);
+       return hash;
+}
+
+void bcf_call_del_rghash(void *_hash)
+{
+       khint_t k;
+       khash_t(rg) *hash = (khash_t(rg)*)_hash;
+       if (hash == 0) return;
+       for (k = kh_begin(hash); k < kh_end(hash); ++k)
+               if (kh_exist(hash, k))
+                       free((char*)kh_key(hash, k));
+       kh_destroy(rg, hash);
+}
+
 static int tpos2qpos(const bam1_core_t *c, const uint32_t *cigar, int32_t tpos, int is_left, int32_t *_tpos)
 {
        int k, x = c->pos, y = 0, last_y = 0;
@@ -58,12 +106,30 @@ static inline int est_indelreg(int pos, const char *ref, int l, char *ins4)
        return max_i - pos;
 }
 
-int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref)
+int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref,
+                                         const void *rghash)
 {
        extern void ks_introsort_uint32_t(int, uint32_t*);
        int i, s, j, k, t, n_types, *types, max_rd_len, left, right, max_ins, *score, N, K, l_run, ref_type;
        char *inscns = 0, *ref2, *query;
+       khash_t(rg) *hash = (khash_t(rg)*)rghash;
        if (ref == 0 || bca == 0) return -1;
+       // mark filtered reads
+       if (rghash) {
+               N = 0;
+               for (s = N = 0; s < n; ++s) {
+                       for (i = 0; i < n_plp[s]; ++i) {
+                               bam_pileup1_t *p = plp[s] + i;
+                               const uint8_t *rg = bam_aux_get(p->b, "RG");
+                               p->aux = 1; // filtered by default
+                               if (rg) {
+                                       khint_t k = kh_get(rg, hash, (const char*)(rg + 1));
+                                       if (k != kh_end(hash)) p->aux = 0, ++N; // not filtered
+                               }
+                       }
+               }
+               if (N == 0) return -1; // no reads left
+       }
        // determine if there is a gap
        for (s = N = 0; s < n; ++s) {
                for (i = 0; i < n_plp[s]; ++i)
@@ -81,7 +147,7 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
                for (s = 0; s < n; ++s) {
                        for (i = 0; i < n_plp[s]; ++i) {
                                const bam_pileup1_t *p = plp[s] + i;
-                               if (p->indel != 0)
+                               if (p->indel != 0 && (rghash == 0 || p->aux == 0))
                                        aux[m++] = MINUS_CONST + p->indel;
                                j = bam_cigar2qlen(&p->b->core, bam1_cigar(p->b));
                                if (j > max_rd_len) max_rd_len = j;
@@ -91,7 +157,9 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
                // squeeze out identical types
                for (i = 1, n_types = 1; i < m; ++i)
                        if (aux[i] != aux[i-1]) ++n_types;
-               assert(n_types > 1); // there must at least one type of non-reference indel
+               if (n_types == 1) { // no indels
+                       free(aux); return -1;
+               }
                types = (int*)calloc(n_types, sizeof(int));
                t = 0;
                types[t++] = aux[0] - MINUS_CONST;