A bug fix, "samtools view" is now working again.

[samtools.git] / bam_rmdupse.c

diff --git a/bam_rmdupse.c b/bam_rmdupse.c
index ed0ef47e2a69209cfb4900eb7d7d5f6697a1cd75..cf1b7bd7b71461938bf43033b8c452eb4a914eaf 100644 (file)
--- a/bam_rmdupse.c
+++ b/bam_rmdupse.c
@@ -9,8 +9,7 @@ typedef struct {
 
 KHASH_MAP_INIT_INT(32, listelem_t)
 
-//#define BLOCK_SIZE 65536
-#define BLOCK_SIZE 100
+#define BLOCK_SIZE 65536
 
 typedef struct {
        bam1_t *b;
@@ -64,12 +63,13 @@ static int fill_buf(samfile_t *in, buffer_t *buf)
                for (i = 0; i < c->l_qseq; ++i) e->score += qual[i] + 1;
                e->score = (double)e->score / sqrt(c->l_qseq + 1);
                is_mapped = (c->tid < 0 || c->tid >= in->header->n_targets || (c->flag&BAM_FUNMAP))? 0 : 1;
+               if (!is_mapped) e->score = -1;
                if (is_mapped && (c->flag & BAM_FREVERSE)) {
                        e->rpos = b->core.pos + bam_calend(&b->core, bam1_cigar(b));
                        if (min_rpos > e->rpos) min_rpos = e->rpos;
                }
                if (buf->n >= capacity) {
-                       if (c->pos <= min_rpos) capacity += BLOCK_SIZE;
+                       if (is_mapped && c->pos <= min_rpos) capacity += BLOCK_SIZE;
                        else break;
                }
        }
@@ -161,7 +161,8 @@ int bam_rmdupse(int argc, char *argv[])
        samfile_t *in, *out;
        buffer_t *buf;
        if (argc < 3) {
-               fprintf(stderr, "Usage: samtools rmdupse <in.bam> <out.bam>\n");
+               fprintf(stderr, "Usage: samtools rmdupse <in.bam> <out.bam>\n\n");
+               fprintf(stderr, "Note: Picard is recommended for this task.\n");
                return 1;
        }
        buf = calloc(1, sizeof(buffer_t));