Skip to content

Commit

Permalink
bam2bcf_indel tidyups.
Browse files Browse the repository at this point in the history
- Remove STRBZ and BQBZ as unuseful stats.

- Only compute stats once per indel rather than per type.

- Replace some pointless mallocs with fixed arrays as the maximum
  number of types was already hard-coded.

- General code tidyups
  • Loading branch information
jkbonfield committed Mar 9, 2021
1 parent 5aac2fe commit 57fe763
Show file tree
Hide file tree
Showing 6 changed files with 127 additions and 145 deletions.
13 changes: 0 additions & 13 deletions bam2bcf.c
Expand Up @@ -64,8 +64,6 @@ bcf_callaux_t *bcf_call_init(double theta, int min_baseQ)
bca->ialt_mq = (int*) malloc(bca->nqual*sizeof(int));
bca->ref_bq = (int*) malloc(bca->nqual*sizeof(int));
bca->alt_bq = (int*) malloc(bca->nqual*sizeof(int));
bca->iref_bq = (int*) malloc(bca->nqual*sizeof(int));
bca->ialt_bq = (int*) malloc(bca->nqual*sizeof(int));
bca->fwd_mqs = (int*) malloc(bca->nqual*sizeof(int));
bca->rev_mqs = (int*) malloc(bca->nqual*sizeof(int));
return bca;
Expand All @@ -83,7 +81,6 @@ void bcf_call_destroy(bcf_callaux_t *bca)
free(bca->ref_mq); free(bca->alt_mq);
free(bca->iref_mq); free(bca->ialt_mq);
free(bca->ref_bq); free(bca->alt_bq);
free(bca->iref_bq); free(bca->ialt_bq);
free(bca->fwd_mqs); free(bca->rev_mqs);
bca->nqual = 0;
free(bca->bases); free(bca->inscns); free(bca);
Expand Down Expand Up @@ -162,8 +159,6 @@ void bcf_callaux_clean(bcf_callaux_t *bca, bcf_call_t *call)
memset(bca->ialt_mq,0,sizeof(int)*bca->nqual);
memset(bca->ref_bq,0,sizeof(int)*bca->nqual);
memset(bca->alt_bq,0,sizeof(int)*bca->nqual);
memset(bca->iref_bq,0,sizeof(int)*bca->nqual);
memset(bca->ialt_bq,0,sizeof(int)*bca->nqual);
memset(bca->fwd_mqs,0,sizeof(int)*bca->nqual);
memset(bca->rev_mqs,0,sizeof(int)*bca->nqual);
if ( call->ADF ) memset(call->ADF,0,sizeof(int32_t)*(call->n+1)*B2B_MAX_ALLELES);
Expand All @@ -174,8 +169,6 @@ void bcf_callaux_clean(bcf_callaux_t *bca, bcf_call_t *call)
memset(bca->alt_scl, 0, 100*sizeof(int));
memset(bca->iref_scl, 0, 100*sizeof(int));
memset(bca->ialt_scl, 0, 100*sizeof(int));
memset(bca->iref_str, 0, 100*sizeof(int));
memset(bca->ialt_str, 0, 100*sizeof(int));
}

/*
Expand Down Expand Up @@ -885,13 +878,9 @@ int bcf_call_combine(int n, const bcf_callret1_t *calls, bcf_callaux_t *bca, int
bca->npos, 0, 1);
call->mwu_mq = calc_mwu_biasZ(bca->iref_mq, bca->ialt_mq,
bca->nqual,1,1);
call->mwu_bq = calc_mwu_biasZ(bca->iref_bq, bca->ialt_bq,
bca->nqual,0,1);
if ( bca->fmt_flag & B2B_INFO_SCB )
call->mwu_sc = calc_mwu_biasZ(bca->iref_scl, bca->ialt_scl,
100, 0,1);
call->mwu_str = calc_mwu_biasZ(bca->iref_str, bca->ialt_str,
100,0,1);
} else {
if (bca->fmt_flag & B2B_INFO_RPB)
call->mwu_pos = calc_mwu_biasZ(bca->ref_pos, bca->alt_pos,
Expand All @@ -905,7 +894,6 @@ int bcf_call_combine(int n, const bcf_callret1_t *calls, bcf_callaux_t *bca, int
if ( bca->fmt_flag & B2B_INFO_SCB )
call->mwu_sc = calc_mwu_biasZ(bca->ref_scl, bca->alt_scl,
100, 0,1);
call->mwu_str = HUGE_VAL;
}


Expand Down Expand Up @@ -1025,7 +1013,6 @@ int bcf_call2bcf(bcf_call_t *bc, bcf1_t *rec, bcf_callret1_t *bcr, int fmt_flag,
if ( bc->mwu_mqs != HUGE_VAL ) bcf_update_info_float(hdr, rec, "MQSBZ", &bc->mwu_mqs, 1);
if ( bc->mwu_bq != HUGE_VAL ) bcf_update_info_float(hdr, rec, "BQBZ", &bc->mwu_bq, 1);
if ( bc->mwu_sc != HUGE_VAL ) bcf_update_info_float(hdr, rec, "SCBZ", &bc->mwu_sc, 1);
if ( bc->mwu_str != HUGE_VAL ) bcf_update_info_float(hdr, rec, "STRBZ", &bc->mwu_str, 1);
#else
if ( bc->mwu_pos != HUGE_VAL ) bcf_update_info_float(hdr, rec, "RPB", &bc->mwu_pos, 1);
if ( bc->mwu_mq != HUGE_VAL ) bcf_update_info_float(hdr, rec, "MQB", &bc->mwu_mq, 1);
Expand Down
5 changes: 2 additions & 3 deletions bam2bcf.h
Expand Up @@ -91,10 +91,9 @@ typedef struct __bcf_callaux_t {
float max_frac; // for collecting indel candidates
int per_sample_flt; // indel filtering strategy
int *ref_pos, *alt_pos, npos, *ref_mq, *alt_mq, *ref_bq, *alt_bq, *fwd_mqs, *rev_mqs, nqual; // for bias tests
int *iref_pos, *ialt_pos, *iref_mq, *ialt_mq, *iref_bq, *ialt_bq; // for indels
int *iref_pos, *ialt_pos, *iref_mq, *ialt_mq; // for indels
int ref_scl[100], alt_scl[100]; // soft-clip length bias; SNP
int iref_scl[100], ialt_scl[100]; // soft-clip length bias; INDEL
int iref_str[100], ialt_str[100]; // STR score bias; INDEL
// for internal uses
int max_bases;
int indel_types[4]; // indel lengths
Expand Down Expand Up @@ -141,7 +140,7 @@ typedef struct {
int32_t *PL, *DP4, *ADR, *ADF, *SCR, *QS;
uint8_t *fmt_arr;
float vdb; // variant distance bias
float mwu_pos, mwu_mq, mwu_bq, mwu_mqs, mwu_sc, mwu_str;
float mwu_pos, mwu_mq, mwu_bq, mwu_mqs, mwu_sc;
#if CDF_MWU_TESTS
float mwu_pos_cdf, mwu_mq_cdf, mwu_bq_cdf, mwu_mqs_cdf;
#endif
Expand Down

0 comments on commit 57fe763

Please sign in to comment.