From c9a9c2742eecde8e097f70efdb6d599374ca325e Mon Sep 17 00:00:00 2001 From: Heng Li Date: Wed, 27 Oct 2010 15:01:11 +0000 Subject: [PATCH] * samtools-0.1.8-21 (r780) * minor speedup to pileup --- bam.h | 2 +- bam_plcmd.c | 33 ++++++++++++++++++++------------- bamtk.c | 2 +- 3 files changed, 22 insertions(+), 15 deletions(-) diff --git a/bam.h b/bam.h index 5943123..4c8536f 100644 --- a/bam.h +++ b/bam.h @@ -230,7 +230,7 @@ typedef struct __bam_iter_t *bam_iter_t; @param b pointer to an alignment @return pointer to quality string */ -#define bam1_qual(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname + ((b)->core.l_qseq + 1)/2) +#define bam1_qual(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname + (((b)->core.l_qseq + 1)>>1)) /*! @function @abstract Get a base on read diff --git a/bam_plcmd.c b/bam_plcmd.c index f0e63ce..6d42471 100644 --- a/bam_plcmd.c +++ b/bam_plcmd.c @@ -183,10 +183,15 @@ static inline void pileup_seq(const bam_pileup1_t *p, int pos, int ref_len, cons putchar(p->b->core.qual > 93? 126 : p->b->core.qual + 33); } if (!p->is_del) { - int rb, c = bam_nt16_rev_table[bam1_seqi(bam1_seq(p->b), p->qpos)]; - rb = (ref && pos < ref_len)? ref[pos] : 'N'; - if (c == '=' || toupper(c) == toupper(rb)) c = bam1_strand(p->b)? ',' : '.'; - else c = bam1_strand(p->b)? tolower(c) : toupper(c); + int c = bam_nt16_rev_table[bam1_seqi(bam1_seq(p->b), p->qpos)]; + if (ref) { + int rb = pos < ref_len? ref[pos] : 'N'; + if (c == '=' || bam_nt16_table[c] == bam_nt16_table[rb]) c = bam1_strand(p->b)? ',' : '.'; + else c = bam1_strand(p->b)? tolower(c) : toupper(c); + } else { + if (c == '=') c = bam1_strand(p->b)? ',' : '.'; + else c = bam1_strand(p->b)? tolower(c) : toupper(c); + } putchar(c); } else putchar(p->is_refskip? (bam1_strand(p->b)? '<' : '>') : '*'); if (p->indel > 0) { @@ -286,16 +291,18 @@ static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *p } // print pileup sequences printw(n, stdout); putchar('\t'); - rms_aux = 0; // we need to recalculate rms_mapq when -c is not flagged on the command line - for (i = 0; i < n; ++i) { - const bam_pileup1_t *p = pu + i; - int tmp = p->b->core.qual < d->c->cap_mapQ? p->b->core.qual : d->c->cap_mapQ; - rms_aux += tmp * tmp; - pileup_seq(p, pos, d->len, d->ref); - } + for (i = 0; i < n; ++i) + pileup_seq(pu + i, pos, d->len, d->ref); // finalize rms_mapq - rms_aux = (uint64_t)(sqrt((double)rms_aux / n) + .499); - if (rms_mapq < 0) rms_mapq = rms_aux; + if (d->format & BAM_PLF_CNS) { + for (i = rms_aux = 0; i < n; ++i) { + const bam_pileup1_t *p = pu + i; + int tmp = p->b->core.qual < d->c->cap_mapQ? p->b->core.qual : d->c->cap_mapQ; + rms_aux += tmp * tmp; + } + rms_aux = (uint64_t)(sqrt((double)rms_aux / n) + .499); + if (rms_mapq < 0) rms_mapq = rms_aux; + } putchar('\t'); // print quality for (i = 0; i < n; ++i) { diff --git a/bamtk.c b/bamtk.c index fb0e388..7bc8756 100644 --- a/bamtk.c +++ b/bamtk.c @@ -9,7 +9,7 @@ #endif #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.1.8-20 (r778)" +#define PACKAGE_VERSION "0.1.8-21 (r780)" #endif int bam_taf2baf(int argc, char *argv[]); -- 2.39.2