X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bam2bcf_indel.c;h=2d0b5c6f8a66a4c5c006dbec560718ab1f7a0377;hb=165f1d8cb88f04a48b7da3f4e5e5fd2a4f074972;hp=de9df3e6c2fb9042b4287ac60e828224aaa651a2;hpb=e6091454768385d722644774042e822658d20713;p=samtools.git

diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c
index de9df3e..2d0b5c6 100644
--- 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;