X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=pysam.git;a=blobdiff_plain;f=samtools%2Fbam_color.c.pysam.c;fp=samtools%2Fbam_color.c.pysam.c;h=f71c65f214b0ec0587807482cc59c8a82e2a7a90;hp=0000000000000000000000000000000000000000;hb=bd0c3067c187d1f718004fb38acc093af8810a02;hpb=1b740fc70684c92a5e2293013217d5a2fd661d8a

diff --git a/samtools/bam_color.c.pysam.c b/samtools/bam_color.c.pysam.c
new file mode 100644
index 0000000..f71c65f
--- /dev/null
+++ b/samtools/bam_color.c.pysam.c
@@ -0,0 +1,129 @@
+#include "pysam.h"
+
+#include <ctype.h>
+#include "bam.h"
+
+/*!
+ @abstract     Get the color encoding the previous and current base
+ @param b      pointer to an alignment
+ @param i      The i-th position, 0-based
+ @return       color
+
+ @discussion   Returns 0 no color information is found.
+ */
+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];
+}
+
+/*!
+ @abstract     Get the color quality of the color encoding the previous and current base
+ @param b      pointer to an alignment
+ @param i      The i-th position, 0-based
+ @return       color quality
+
+ @discussion   Returns 0 no color information is found.
+ */
+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)];
+}
+
+/*!
+ @abstract     Get the color error profile at the give position    
+ @param b      pointer to an alignment
+ @return       the original color if the color was an error, '-' (dash) otherwise
+
+ @discussion   Returns 0 no color information is found.
+ */
+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.  Note: must rc adaptor
+		prev_b = (cs_i == 1) ? "TGCAN"[(int)bam_aux_nt2int(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;
+	}
+}