X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bam_aux.c;h=2fa9ac2e83a7a00e2e83a523fe687448e3c4255e;hb=41586185b43962c1ffb3f82e93066a480254aa17;hp=f9be398f0e820d9b48b2002073039af1ebeb95d5;hpb=7bd63d5c2677364126702d76c56e4d75981c62f4;p=samtools.git

diff --git a/bam_aux.c b/bam_aux.c
index f9be398..2fa9ac2 100644
--- a/bam_aux.c
+++ b/bam_aux.c
@@ -3,6 +3,21 @@
 #include "khash.h"
 KHASH_MAP_INIT_STR(s, int)
 
+void bam_aux_append(bam1_t *b, const char tag[2], char type, int len, uint8_t *data)
+{
+	int ori_len = b->data_len;
+	b->data_len += 3 + len;
+	b->l_aux += 3 + len;
+	if (b->m_data < b->data_len) {
+		b->m_data = b->data_len;
+		kroundup32(b->m_data);
+		b->data = (uint8_t*)realloc(b->data, b->m_data);
+	}
+	b->data[ori_len] = tag[0]; b->data[ori_len + 1] = tag[1];
+	b->data[ori_len + 2] = type;
+	memcpy(b->data + ori_len + 3, data, len);
+}
+
 uint8_t *bam_aux_get(bam1_t *b, const char tag[2])
 {
 	uint8_t *s;
@@ -59,7 +74,7 @@ void bam_parse_region(bam_header_t *header, const char *str, int *ref_id, int *b
 
 	bam_init_header_hash(header);
 	h = (khash_t(s)*)header->hash;
-	
+
 	l = strlen(str);
 	p = s = (char*)malloc(l+1);
 	/* squeeze out "," */
@@ -137,3 +152,105 @@ char *bam_aux2Z(const uint8_t *s)
 	if (type == 'Z' || type == 'H') return (char*)s;
 	else return 0;
 }
+
+char bam_aux_getCSi(bam1_t *b, int i)
+{
+	uint8_t *c = bam_aux_get(b, "CS");
+	char *cs = NULL;
+
+	// return the base if the tag was not found
+	if(0 == c) return 0;
+
+	cs = bam_aux2Z(c);
+	// adjust for strandedness and leading adaptor
+	if(bam1_strand(b)) i = strlen(cs) - 1 - i;
+	else i++;
+	return cs[i];
+}
+
+char bam_aux_getCQi(bam1_t *b, int i)
+{
+	uint8_t *c = bam_aux_get(b, "CQ");
+	char *cq = NULL;
+	
+	// return the base if the tag was not found
+	if(0 == c) return 0;
+
+	cq = bam_aux2Z(c);
+	// adjust for strandedness
+	if(bam1_strand(b)) i = strlen(cq) - 1 - i;
+	return cq[i];
+}
+
+char bam_aux_nt2int(char a)
+{
+	switch(toupper(a)) {
+		case 'A':
+			return 0;
+			break;
+		case 'C':
+			return 1;
+			break;
+		case 'G':
+			return 2;
+			break;
+		case 'T':
+			return 3;
+			break;
+		default:
+			return 4;
+			break;
+	}
+}
+
+char bam_aux_ntnt2cs(char a, char b)
+{
+	a = bam_aux_nt2int(a);
+	b = bam_aux_nt2int(b);
+	if(4 == a || 4 == b) return '4';
+	return "0123"[(int)(a ^ b)];
+}
+
+char bam_aux_getCEi(bam1_t *b, int i)
+{
+	int cs_i;
+	uint8_t *c = bam_aux_get(b, "CS");
+	char *cs = NULL;
+	char prev_b, cur_b;
+	char cur_color, cor_color;
+
+	// return the base if the tag was not found
+	if(0 == c) return 0;
+	
+	cs = bam_aux2Z(c);
+
+	// adjust for strandedness and leading adaptor
+	if(bam1_strand(b)) { //reverse strand
+		cs_i = strlen(cs) - 1 - i;
+		// get current color
+		cur_color = cs[cs_i];
+		// get previous base
+		prev_b = (0 == cs_i) ? cs[0] : bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i+1)];
+		// get current base
+		cur_b = bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i)]; 
+	}
+	else {
+		cs_i=i+1;
+		// get current color
+		cur_color = cs[cs_i];
+		// get previous base
+		prev_b = (0 == i) ? cs[0] : bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i-1)];
+		// get current base
+		cur_b = bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i)];
+	}
+
+	// corrected color
+	cor_color = bam_aux_ntnt2cs(prev_b, cur_b);
+
+	if(cur_color == cor_color) { 
+		return '-';
+	}
+	else {
+		return cur_color;
+	}
+}