Skip to content

Commit 7e48e3c

Browse files
gperteagpertea
committed
re-fix of samtools aux tag parsing bug using htslib code
1 parent 4047408 commit 7e48e3c

3 files changed

Lines changed: 58 additions & 17 deletions

File tree

‎samtools-0.1.18/bam_aux.c‎

Lines changed: 46 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -25,14 +25,53 @@ uint8_t *bam_aux_get_core(bam1_t *b, const char tag[2])
2525
return bam_aux_get(b, tag);
2626
}
2727

28+
inline int aux_type2size(uint8_t type)
29+
{
30+
switch (type) {
31+
case 'A': case 'c': case 'C':
32+
return 1;
33+
case 's': case 'S':
34+
return 2;
35+
case 'i': case 'I': case 'f':
36+
return 4;
37+
case 'd':
38+
return 8;
39+
case 'Z': case 'H': case 'B':
40+
return type;
41+
default:
42+
return 0;
43+
}
44+
}
45+
46+
inline uint8_t* skip_aux(uint8_t* s) {
47+
int size = aux_type2size(*s); ++s; // skip type
48+
uint32_t n;
49+
switch (size) {
50+
case 'Z':
51+
case 'H':
52+
while (*s) ++s;
53+
return s+1;
54+
case 'B':
55+
size = aux_type2size(*s); ++s;
56+
memcpy(&n, s, 4); s += 4;
57+
return s + size * n;
58+
case 0:
59+
abort();
60+
break;
61+
default:
62+
return s + size;
63+
}
64+
}
65+
66+
/*
2867
#define __skip_tag(s) do { \
2968
int type = toupper(*(s)); \
3069
++(s); \
3170
if (type == 'Z' || type == 'H') { while (*(s)) ++(s); ++(s); } \
3271
else if (type == 'B') (s) += 5 + bam_aux_type2size(*(s)) * (*(int32_t*)((s)+1)); \
3372
else (s) += bam_aux_type2size(*(s)); \
3473
} while(0)
35-
74+
*/
3675
uint8_t *bam_aux_get(const bam1_t *b, const char tag[2])
3776
{
3877
uint8_t *s;
@@ -42,7 +81,8 @@ uint8_t *bam_aux_get(const bam1_t *b, const char tag[2])
4281
int x = (int)s[0]<<8 | s[1];
4382
s += 2;
4483
if (x == y) return s;
45-
__skip_tag(s);
84+
//__skip_tag(s);
85+
s=skip_aux(s);
4686
}
4787
return 0;
4888
}
@@ -52,7 +92,8 @@ int bam_aux_del(bam1_t *b, uint8_t *s)
5292
uint8_t *p, *aux;
5393
aux = bam1_aux(b);
5494
p = s - 2;
55-
__skip_tag(s);
95+
//__skip_tag(s);
96+
s=skip_aux(s);
5697
memmove(p, s, b->l_aux - (s - aux));
5798
b->data_len -= s - p;
5899
b->l_aux -= s - p;
@@ -65,7 +106,8 @@ int bam_aux_drop_other(bam1_t *b, uint8_t *s)
65106
uint8_t *p, *aux;
66107
aux = bam1_aux(b);
67108
p = s - 2;
68-
__skip_tag(s);
109+
//__skip_tag(s);
110+
s=skip_aux(s);
69111
memmove(aux, p, s - p);
70112
b->data_len -= b->l_aux - (s - p);
71113
b->l_aux = s - p;

‎samtools-0.1.18/bgzf.c‎

Lines changed: 8 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -73,24 +73,19 @@ static const int GZIP_WINDOW_BITS = -15; // no zlib header
7373
static const int Z_DEFAULT_MEM_LEVEL = 8;
7474

7575

76-
inline
77-
void
78-
packInt16(uint8_t* buffer, uint16_t value)
76+
77+
void packInt16(uint8_t* buffer, uint16_t value)
7978
{
8079
buffer[0] = value;
8180
buffer[1] = value >> 8;
8281
}
8382

84-
inline
85-
int
86-
unpackInt16(const uint8_t* buffer)
83+
int unpackInt16(const uint8_t* buffer)
8784
{
8885
return (buffer[0] | (buffer[1] << 8));
8986
}
9087

91-
inline
92-
void
93-
packInt32(uint8_t* buffer, uint32_t value)
88+
void packInt32(uint8_t* buffer, uint32_t value)
9489
{
9590
buffer[0] = value;
9691
buffer[1] = value >> 8;
@@ -627,11 +622,12 @@ int bgzf_close(BGZF* fp)
627622
if (fp->open_mode == 'w') {
628623
if (bgzf_flush(fp) != 0) return -1;
629624
{ // add an empty block
630-
int count, block_length = deflate_block(fp, 0);
625+
int block_length = deflate_block(fp, 0);
626+
//int count=
631627
#ifdef _USE_KNETFILE
632-
count = fwrite(fp->compressed_block, 1, block_length, fp->x.fpw);
628+
fwrite(fp->compressed_block, 1, block_length, fp->x.fpw);
633629
#else
634-
count = fwrite(fp->compressed_block, 1, block_length, fp->file);
630+
fwrite(fp->compressed_block, 1, block_length, fp->file);
635631
#endif
636632
}
637633
#ifdef _USE_KNETFILE

‎stringtie.cpp‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -442,6 +442,8 @@ if (tstackSize<DEF_TSTACK_SIZE) defStackSize=DEF_TSTACK_SIZE;
442442
xstrand=brec->spliceStrand();
443443
if (xstrand=='.' && brec->exons.Count()>1) {
444444
no_xs++;
445+
//if (verbose)
446+
// GMessage("XS tag missing for read\t%s\t(%d)\t%s\n",brec->name(), brec->start, brec->cigar());
445447
continue; //skip spliced alignments lacking XS tag (e.g. HISAT alignments)
446448
}
447449
if (refseqName==NULL) GError("Error: cannot retrieve target seq name from BAM record!\n");
@@ -713,7 +715,8 @@ if (tstackSize<DEF_TSTACK_SIZE) defStackSize=DEF_TSTACK_SIZE;
713715
if(!mergeMode) {
714716
if(verbose) {
715717
GMessage("Total count of aligned fragments: %g\n", Num_Fragments);
716-
GMessage("Fragment coverage length: %g\n", Frag_Len/Num_Fragments);
718+
if (Num_Fragments)
719+
GMessage("Fragment coverage length: %g\n", Frag_Len/Num_Fragments);
717720
}
718721

719722
f_out=stdout;

0 commit comments

Comments
 (0)