Mercurial > repos > lsong10 > psiclass
view PsiCLASS-1.0.2/samtools-0.1.19/bam_color.c @ 0:903fc43d6227 draft default tip
Uploaded
author | lsong10 |
---|---|
date | Fri, 26 Mar 2021 16:52:45 +0000 |
parents | |
children |
line wrap: on
line source
#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; // adjust for leading hard clip uint32_t cigar = bam1_cigar(b)[0]; if((cigar & BAM_CIGAR_MASK) == BAM_CHARD_CLIP) { i -= cigar >> BAM_CIGAR_SHIFT; } } 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; // adjust for leading hard clip uint32_t cigar = bam1_cigar(b)[0]; if((cigar & BAM_CIGAR_MASK) == BAM_CHARD_CLIP) { i -= (cigar >> BAM_CIGAR_SHIFT); } } 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; // adjust for leading hard clip uint32_t cigar = bam1_cigar(b)[0]; if((cigar & BAM_CIGAR_MASK) == BAM_CHARD_CLIP) { cs_i -= cigar >> BAM_CIGAR_SHIFT; } // 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; } }