Skip to content

Commit 0486a68

Browse files
sirus20x6pd3
authored andcommitted
Fix three bugs in bam2bcf indel calling
- bam2bcf_indel.c: Initialize K=0 at declaration to avoid using an uninitialized variable in the max-deletion loop (line ~783) which runs before K is reset at line ~791. - bam2bcf_iaux.c: Fix memset size in iaux_init_scores() to use n*sizeof(int) instead of n, since read_scores is int*. Without this, only 1/4 of the array (on most platforms) was being zeroed. - bam2bcf_iaux.c: Fix per-sample indel fraction denominator in iaux_init_types() to use ntot (per-sample read count) instead of naux (running cross-sample indel count), which produced incorrect filtering thresholds.
1 parent 944cebd commit 0486a68

2 files changed

Lines changed: 3 additions & 3 deletions

File tree

bam2bcf_iaux.c

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -162,7 +162,7 @@ static int iaux_init_scores(indel_aux_t *iaux, int ismpl)
162162
iaux->mread_scores = n;
163163
iaux->read_scores = tmp;
164164
}
165-
memset(iaux->read_scores,0,n);
165+
memset(iaux->read_scores,0,n*sizeof(int));
166166
return 0;
167167
}
168168

@@ -287,7 +287,7 @@ static int iaux_init_types(indel_aux_t *iaux)
287287
nalt = naux - nalt;
288288
if ( iaux->bca->per_sample_flt )
289289
{
290-
double frac = (double)nalt/naux;
290+
double frac = (double)nalt/ntot;
291291
if ( nalt >= iaux->bca->min_support && frac >= iaux->bca->min_frac ) indel_support_ok = 1;
292292
if ( nalt > iaux->bca->max_support && frac > 0 ) iaux->bca->max_support = nalt, iaux->bca->max_frac = frac;
293293
}

bam2bcf_indel.c

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -702,7 +702,7 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos,
702702

703703
int i, s, j, k, t, n_types, *types, max_rd_len, left, right, max_ins;
704704
int *score, max_ref2;
705-
int N, K, l_run, ref_type, n_alt;
705+
int N, K = 0, l_run, ref_type, n_alt;
706706
char *inscns = 0, *ref2, *query, **ref_sample;
707707

708708
// determine if there is a gap

0 commit comments

Comments
 (0)