bam2bcf.c

   1 #include <math.h>
   2 #include <stdint.h>
   3 #include <assert.h>
   4 #include "bam.h"
   5 #include "kstring.h"
   6 #include "bam2bcf.h"
   7 #include "errmod.h"
   8 #include "bcftools/bcf.h"
   9
  10 extern  void ks_introsort_uint32_t(size_t n, uint32_t a[]);
  11
  12 #define CALL_ETA 0.03f
  13 #define CALL_MAX 256
  14 #define CALL_DEFTHETA 0.83f
  15 #define DEF_MAPQ 20
  16
  17 #define CAP_DIST 25
  18
  19 bcf_callaux_t *bcf_call_init(double theta, int min_baseQ)
  20 {
  21         bcf_callaux_t *bca;
  22         if (theta <= 0.) theta = CALL_DEFTHETA;
  23         bca = calloc(1, sizeof(bcf_callaux_t));
  24         bca->capQ = 60;
  25         bca->openQ = 40; bca->extQ = 20; bca->tandemQ = 100;
  26         bca->min_baseQ = min_baseQ;
  27         bca->e = errmod_init(1. - theta);
  28         bca->min_frac = 0.002;
  29         bca->min_support = 1;
  30     bca->per_sample_flt = 0;
  31         return bca;
  32 }
  33
  34
  35 int dist_from_end(const bam_pileup1_t *p)
  36 {
  37     int icig, n_tot_bases = 0, iread = 0, edist = p->qpos + 1;
  38     for (icig=0; icig<p->b->core.n_cigar; icig++)
  39     {
  40         // Conversion from uint32_t to MIDNSHP
  41         //  0123456
  42         //  MIDNSHP
  43         int cig  = bam1_cigar(p->b)[icig] & BAM_CIGAR_MASK;
  44         int ncig = bam1_cigar(p->b)[icig] >> BAM_CIGAR_SHIFT;
  45         if ( cig==0 )
  46         {
  47             n_tot_bases += ncig;
  48             iread += ncig;
  49         }
  50         else if ( cig==1 )
  51         {
  52             n_tot_bases += ncig;
  53             iread += ncig;
  54         }
  55         else if ( cig==4 )
  56         {
  57             iread += ncig;
  58             if ( iread<=p->qpos ) edist -= ncig;
  59         }
  60     }
  61     // distance from either end
  62     if ( edist > n_tot_bases/2. ) edist = n_tot_bases - edist + 1;
  63     if ( edist<0 ) { fprintf(stderr,"uh: %d %d -> %d (%d)\n", p->qpos,n_tot_bases,edist,p->b->core.pos+1); exit(1); }
  64     return edist;
  65 }
  66
  67 void bcf_call_destroy(bcf_callaux_t *bca)
  68 {
  69         if (bca == 0) return;
  70         errmod_destroy(bca->e);
  71     if (bca->npos) { free(bca->ref_pos); free(bca->alt_pos); bca->npos = 0; }
  72         free(bca->bases); free(bca->inscns); free(bca);
  73 }
  74 /* ref_base is the 4-bit representation of the reference base. It is
  75  * negative if we are looking at an indel. */
  76 int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base, bcf_callaux_t *bca, bcf_callret1_t *r)
  77 {
  78         int i, n, ref4, is_indel, ori_depth = 0;
  79         memset(r, 0, sizeof(bcf_callret1_t));
  80         if (ref_base >= 0) {
  81                 ref4 = bam_nt16_nt4_table[ref_base];
  82                 is_indel = 0;
  83         } else ref4 = 4, is_indel = 1;
  84         if (_n == 0) return -1;
  85         // enlarge the bases array if necessary
  86         if (bca->max_bases < _n) {
  87                 bca->max_bases = _n;
  88                 kroundup32(bca->max_bases);
  89                 bca->bases = (uint16_t*)realloc(bca->bases, 2 * bca->max_bases);
  90         }
  91         // fill the bases array
  92     bca->read_len = 0;
  93         for (i = n = r->n_supp = 0; i < _n; ++i) {
  94                 const bam_pileup1_t *p = pl + i;
  95                 int q, b, mapQ, baseQ, is_diff, min_dist, seqQ;
  96                 // set base
  97                 if (p->is_del || p->is_refskip || (p->b->core.flag&BAM_FUNMAP)) continue;
  98                 ++ori_depth;
  99                 baseQ = q = is_indel? p->aux&0xff : (int)bam1_qual(p->b)[p->qpos]; // base/indel quality
 100                 seqQ = is_indel? (p->aux>>8&0xff) : 99;
 101                 if (q < bca->min_baseQ) continue;
 102                 if (q > seqQ) q = seqQ;
 103                 mapQ = p->b->core.qual < 255? p->b->core.qual : DEF_MAPQ; // special case for mapQ==255
 104                 mapQ = mapQ < bca->capQ? mapQ : bca->capQ;
 105                 if (q > mapQ) q = mapQ;
 106                 if (q > 63) q = 63;
 107                 if (q < 4) q = 4;
 108                 if (!is_indel) {
 109                         b = bam1_seqi(bam1_seq(p->b), p->qpos); // base
 110                         b = bam_nt16_nt4_table[b? b : ref_base]; // b is the 2-bit base
 111                         is_diff = (ref4 < 4 && b == ref4)? 0 : 1;
 112                 } else {
 113                         b = p->aux>>16&0x3f;
 114                         is_diff = (b != 0);
 115                 }
 116                 if (is_diff) ++r->n_supp;
 117                 bca->bases[n++] = q<<5 | (int)bam1_strand(p->b)<<4 | b;
 118                 // collect annotations
 119                 if (b < 4) r->qsum[b] += q;
 120                 ++r->anno[0<<2|is_diff<<1|bam1_strand(p->b)];
 121                 min_dist = p->b->core.l_qseq - 1 - p->qpos;
 122                 if (min_dist > p->qpos) min_dist = p->qpos;
 123                 if (min_dist > CAP_DIST) min_dist = CAP_DIST;
 124                 r->anno[1<<2|is_diff<<1|0] += baseQ;
 125                 r->anno[1<<2|is_diff<<1|1] += baseQ * baseQ;    // FIXME: signed int is not enough for thousands of samples
 126                 r->anno[2<<2|is_diff<<1|0] += mapQ;
 127                 r->anno[2<<2|is_diff<<1|1] += mapQ * mapQ;              // FIXME: signed int is not enough for thousands of samples
 128                 r->anno[3<<2|is_diff<<1|0] += min_dist;
 129                 r->anno[3<<2|is_diff<<1|1] += min_dist * min_dist;
 130
 131         // collect read positions for ReadPosBias
 132         int epos = dist_from_end(p);
 133         int npos = p->b->core.l_qseq;
 134         if ( bca->npos < npos )
 135         {
 136             bca->ref_pos = realloc(bca->ref_pos, sizeof(int)*npos);
 137             bca->alt_pos = realloc(bca->alt_pos, sizeof(int)*npos);
 138             int j;
 139             for (j=bca->npos; j<npos; j++) bca->ref_pos[j] = 0;
 140             for (j=bca->npos; j<npos; j++) bca->alt_pos[j] = 0;
 141             bca->npos = npos;
 142         }
 143         if ( bam1_seqi(bam1_seq(p->b),p->qpos) == ref_base )
 144             bca->ref_pos[epos]++;
 145         else
 146             bca->alt_pos[epos]++;
 147         bca->read_len += npos;
 148         //fprintf(stderr,"qpos=%d  epos=%d  fwd=%d  var=%d  len=%d  %s\n", p->qpos, epos, p->b->core.flag&BAM_FREVERSE ? 0 : 1, bam1_seqi(bam1_seq(p->b),p->qpos) == ref_base ? 1 : 0, p->b->core.l_qseq, bam1_qname(p->b));
 149         }
 150     if ( n ) bca->read_len /= n;
 151         r->depth = n; r->ori_depth = ori_depth;
 152         // glfgen
 153         errmod_cal(bca->e, n, 5, bca->bases, r->p);
 154         return r->depth;
 155 }
 156
 157 double mann_whitney_1947(int n, int m, int U)
 158 {
 159     if (U<0) return 0;
 160     if (n==0||m==0) return U==0 ? 1 : 0;
 161     return (double)n/(n+m)*mann_whitney_1947(n-1,m,U-m) + (double)m/(n+m)*mann_whitney_1947(n,m-1,U);
 162 }
 163
 164 void calc_ReadPosBias(bcf_callaux_t *bca, bcf_call_t *call)
 165 {
 166     int i, nref = 0, nalt = 0;
 167     unsigned long int U = 0;
 168     for (i=0; i<bca->npos; i++)
 169     {
 170         nref += bca->ref_pos[i];
 171         nalt += bca->alt_pos[i];
 172         U += nref*bca->alt_pos[i];
 173         bca->ref_pos[i] = 0;
 174         bca->alt_pos[i] = 0;
 175     }
 176     if ( !nref || !nalt )
 177     {
 178         call->read_pos_bias = -1;
 179         return;
 180     }
 181
 182     if ( nref>=8 || nalt>=8 )
 183     {
 184         // normal approximation
 185         double mean = ((double)nref*nalt+1.0)/2.0;
 186         double var2 = (double)nref*nalt*(nref+nalt+1.0)/12.0;
 187         double z    = (U-mean)/sqrt(var2);
 188         call->read_pos_bias = z;
 189         //fprintf(stderr,"nref=%d  nalt=%d  U=%ld  mean=%e  var=%e  zval=%e\n", nref,nalt,U,mean,sqrt(var2),call->read_pos_bias);
 190     }
 191     else
 192     {
 193         double p = mann_whitney_1947(nalt,nref,U);
 194         // biased form claimed by GATK to behave better empirically
 195         // double var2 = (1.0+1.0/(nref+nalt+1.0))*(double)nref*nalt*(nref+nalt+1.0)/12.0;
 196         double var2 = (double)nref*nalt*(nref+nalt+1.0)/12.0;
 197         double z;
 198         if ( p >= 1./sqrt(var2*2*M_PI) ) z = 0;   // equal to mean
 199         else
 200         {
 201             if ( U >= nref*nalt/2. ) z = sqrt(-2*log(sqrt(var2*2*M_PI)*p));
 202             else z = -sqrt(-2*log(sqrt(var2*2*M_PI)*p));
 203         }
 204         call->read_pos_bias = z;
 205         //fprintf(stderr,"nref=%d  nalt=%d  U=%ld  p=%e var2=%e  zval=%e\n", nref,nalt,U, p,var2,call->read_pos_bias);
 206     }
 207 }
 208
 209 float mean_diff_to_prob(float mdiff, int dp, int readlen)
 210 {
 211     if ( dp==2 )
 212     {
 213         if ( mdiff==0 )
 214             return (2.0*readlen + 4.0*(readlen-1.0))/((float)readlen*readlen);
 215         else
 216             return 8.0*(readlen - 4.0*mdiff)/((float)readlen*readlen);
 217     }
 218
 219     // This is crude empirical approximation and is not very accurate for
 220     // shorter read lengths (<100bp). There certainly is a room for
 221     // improvement.
 222     const float mv[24][2] = { {0,0}, {0,0}, {0,0},
 223         { 9.108, 4.934}, { 9.999, 3.991}, {10.273, 3.485}, {10.579, 3.160},
 224         {10.828, 2.889}, {11.014, 2.703}, {11.028, 2.546}, {11.244, 2.391},
 225         {11.231, 2.320}, {11.323, 2.138}, {11.403, 2.123}, {11.394, 1.994},
 226         {11.451, 1.928}, {11.445, 1.862}, {11.516, 1.815}, {11.560, 1.761},
 227         {11.544, 1.728}, {11.605, 1.674}, {11.592, 1.652}, {11.674, 1.613},
 228         {11.641, 1.570} };
 229
 230     float m, v;
 231     if ( dp>=24 )
 232     {
 233         m = (int)readlen/8.;
 234         if (dp>100) dp = 100;
 235         v = 1.476/(0.182*pow(dp,0.514));
 236         v = v*(readlen/100.);
 237     }
 238     else
 239     {
 240         m = mv[dp][0];
 241         v = mv[dp][1];
 242         m = (int)m*readlen/100.;
 243         v = v*readlen/100.;
 244         v *= 1.2;   // allow more variability
 245     }
 246     return 1.0/(v*sqrt(2*M_PI)) * exp(-0.5*((mdiff-m)/v)*((mdiff-m)/v));
 247 }
 248
 249 void calc_vdb(bcf_callaux_t *bca, bcf_call_t *call)
 250 {
 251     int i, dp = 0;
 252     float mean_pos = 0, mean_diff = 0;
 253     for (i=0; i<bca->npos; i++)
 254     {
 255         if ( !bca->alt_pos[i] ) continue;
 256         dp += bca->alt_pos[i];
 257         mean_pos += bca->alt_pos[i]*i;
 258     }
 259     if ( !dp )
 260     {
 261         call->vdb = -1;
 262         return;
 263     }
 264     mean_pos /= dp;
 265     for (i=0; i<bca->npos; i++)
 266     {
 267         if ( !bca->alt_pos[i] ) continue;
 268         mean_diff += bca->alt_pos[i] * fabs(i - mean_pos);
 269     }
 270     mean_diff /= dp;
 271     call->vdb = mean_diff_to_prob(mean_diff, dp, bca->read_len);
 272 }
 273
 274 /**
 275  *  bcf_call_combine() - sets the PL array and VDB, RPB annotations, finds the top two alleles
 276  *  @n:         number of samples
 277  *  @calls:     each sample's calls
 278  *  @bca:       auxiliary data structure for holding temporary values
 279  *  @ref_base:  the reference base
 280  *  @call:      filled with the annotations
 281  */
 282 int bcf_call_combine(int n, const bcf_callret1_t *calls, bcf_callaux_t *bca, int ref_base /*4-bit*/, bcf_call_t *call)
 283 {
 284         int ref4, i, j, qsum[4];
 285         int64_t tmp;
 286         if (ref_base >= 0) {
 287                 call->ori_ref = ref4 = bam_nt16_nt4_table[ref_base];
 288                 if (ref4 > 4) ref4 = 4;
 289         } else call->ori_ref = -1, ref4 = 0;
 290         // calculate qsum
 291         memset(qsum, 0, 4 * sizeof(int));
 292         for (i = 0; i < n; ++i)
 293                 for (j = 0; j < 4; ++j)
 294                         qsum[j] += calls[i].qsum[j];
 295     int qsum_tot=0;
 296     for (j=0; j<4; j++) { qsum_tot += qsum[j]; call->qsum[j] = 0; }
 297         for (j = 0; j < 4; ++j) qsum[j] = qsum[j] << 2 | j;
 298         // find the top 2 alleles
 299         for (i = 1; i < 4; ++i) // insertion sort
 300                 for (j = i; j > 0 && qsum[j] < qsum[j-1]; --j)
 301                         tmp = qsum[j], qsum[j] = qsum[j-1], qsum[j-1] = tmp;
 302         // set the reference allele and alternative allele(s)
 303         for (i = 0; i < 5; ++i) call->a[i] = -1;
 304         call->unseen = -1;
 305         call->a[0] = ref4;
 306         for (i = 3, j = 1; i >= 0; --i) {
 307                 if ((qsum[i]&3) != ref4) {
 308                         if (qsum[i]>>2 != 0)
 309             {
 310                 if ( j<4 ) call->qsum[j] = (float)(qsum[i]>>2)/qsum_tot; // ref N can make j>=4
 311                 call->a[j++]  = qsum[i]&3;
 312             }
 313                         else break;
 314                 }
 315         else
 316             call->qsum[0] = (float)(qsum[i]>>2)/qsum_tot;
 317         }
 318         if (ref_base >= 0) { // for SNPs, find the "unseen" base
 319                 if (((ref4 < 4 && j < 4) || (ref4 == 4 && j < 5)) && i >= 0)
 320                         call->unseen = j, call->a[j++] = qsum[i]&3;
 321                 call->n_alleles = j;
 322         } else {
 323                 call->n_alleles = j;
 324                 if (call->n_alleles == 1) return -1; // no reliable supporting read. stop doing anything
 325         }
 326         // set the PL array
 327         if (call->n < n) {
 328                 call->n = n;
 329                 call->PL = realloc(call->PL, 15 * n);
 330         }
 331         {
 332                 int x, g[15], z;
 333                 double sum_min = 0.;
 334                 x = call->n_alleles * (call->n_alleles + 1) / 2;
 335                 // get the possible genotypes
 336                 for (i = z = 0; i < call->n_alleles; ++i)
 337                         for (j = 0; j <= i; ++j)
 338                                 g[z++] = call->a[j] * 5 + call->a[i];
 339                 for (i = 0; i < n; ++i) {
 340                         uint8_t *PL = call->PL + x * i;
 341                         const bcf_callret1_t *r = calls + i;
 342                         float min = 1e37;
 343                         for (j = 0; j < x; ++j)
 344                                 if (min > r->p[g[j]]) min = r->p[g[j]];
 345                         sum_min += min;
 346                         for (j = 0; j < x; ++j) {
 347                                 int y;
 348                                 y = (int)(r->p[g[j]] - min + .499);
 349                                 if (y > 255) y = 255;
 350                                 PL[j] = y;
 351                         }
 352                 }
 353 //              if (ref_base < 0) fprintf(stderr, "%d,%d,%f,%d\n", call->n_alleles, x, sum_min, call->unseen);
 354                 call->shift = (int)(sum_min + .499);
 355         }
 356         // combine annotations
 357         memset(call->anno, 0, 16 * sizeof(int));
 358         for (i = call->depth = call->ori_depth = 0, tmp = 0; i < n; ++i) {
 359                 call->depth += calls[i].depth;
 360                 call->ori_depth += calls[i].ori_depth;
 361                 for (j = 0; j < 16; ++j) call->anno[j] += calls[i].anno[j];
 362         }
 363
 364     calc_vdb(bca, call);
 365     calc_ReadPosBias(bca, call);
 366
 367         return 0;
 368 }
 369
 370 int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b, bcf_callret1_t *bcr, int fmt_flag,
 371                                  const bcf_callaux_t *bca, const char *ref)
 372 {
 373         extern double kt_fisher_exact(int n11, int n12, int n21, int n22, double *_left, double *_right, double *two);
 374         kstring_t s;
 375         int i, j;
 376         b->n_smpl = bc->n;
 377         b->tid = tid; b->pos = pos; b->qual = 0;
 378         s.s = b->str; s.m = b->m_str; s.l = 0;
 379         kputc('\0', &s);
 380         if (bc->ori_ref < 0) { // an indel
 381                 // write REF
 382                 kputc(ref[pos], &s);
 383                 for (j = 0; j < bca->indelreg; ++j) kputc(ref[pos+1+j], &s);
 384                 kputc('\0', &s);
 385                 // write ALT
 386                 kputc(ref[pos], &s);
 387                 for (i = 1; i < 4; ++i) {
 388                         if (bc->a[i] < 0) break;
 389                         if (i > 1) {
 390                                 kputc(',', &s); kputc(ref[pos], &s);
 391                         }
 392                         if (bca->indel_types[bc->a[i]] < 0) { // deletion
 393                                 for (j = -bca->indel_types[bc->a[i]]; j < bca->indelreg; ++j)
 394                                         kputc(ref[pos+1+j], &s);
 395                         } else { // insertion; cannot be a reference unless a bug
 396                                 char *inscns = &bca->inscns[bc->a[i] * bca->maxins];
 397                                 for (j = 0; j < bca->indel_types[bc->a[i]]; ++j)
 398                                         kputc("ACGTN"[(int)inscns[j]], &s);
 399                                 for (j = 0; j < bca->indelreg; ++j) kputc(ref[pos+1+j], &s);
 400                         }
 401                 }
 402                 kputc('\0', &s);
 403         } else { // a SNP
 404                 kputc("ACGTN"[bc->ori_ref], &s); kputc('\0', &s);
 405                 for (i = 1; i < 5; ++i) {
 406                         if (bc->a[i] < 0) break;
 407                         if (i > 1) kputc(',', &s);
 408                         kputc(bc->unseen == i? 'X' : "ACGT"[bc->a[i]], &s);
 409                 }
 410                 kputc('\0', &s);
 411         }
 412         kputc('\0', &s);
 413         // INFO
 414         if (bc->ori_ref < 0) ksprintf(&s,"INDEL;IS=%d,%f;", bca->max_support, bca->max_frac);
 415         kputs("DP=", &s); kputw(bc->ori_depth, &s); kputs(";I16=", &s);
 416         for (i = 0; i < 16; ++i) {
 417                 if (i) kputc(',', &s);
 418                 kputw(bc->anno[i], &s);
 419         }
 420     ksprintf(&s,";QS=%f,%f,%f,%f", bc->qsum[0],bc->qsum[1],bc->qsum[2],bc->qsum[3]);
 421     if (bc->vdb != -1)
 422         ksprintf(&s, ";VDB=%.4f", bc->vdb);
 423     if (bc->read_pos_bias != -1 )
 424         ksprintf(&s, ";RPB=%.4f", bc->read_pos_bias);
 425         kputc('\0', &s);
 426         // FMT
 427         kputs("PL", &s);
 428         if (bcr && fmt_flag) {
 429                 if (fmt_flag & B2B_FMT_DP) kputs(":DP", &s);
 430                 if (fmt_flag & B2B_FMT_DV) kputs(":DV", &s);
 431                 if (fmt_flag & B2B_FMT_SP) kputs(":SP", &s);
 432         }
 433         kputc('\0', &s);
 434         b->m_str = s.m; b->str = s.s; b->l_str = s.l;
 435         bcf_sync(b);
 436         memcpy(b->gi[0].data, bc->PL, b->gi[0].len * bc->n);
 437         if (bcr && fmt_flag) {
 438                 uint16_t *dp = (fmt_flag & B2B_FMT_DP)? b->gi[1].data : 0;
 439                 uint16_t *dv = (fmt_flag & B2B_FMT_DV)? b->gi[1 + ((fmt_flag & B2B_FMT_DP) != 0)].data : 0;
 440                 int32_t  *sp = (fmt_flag & B2B_FMT_SP)? b->gi[1 + ((fmt_flag & B2B_FMT_DP) != 0) + ((fmt_flag & B2B_FMT_DV) != 0)].data : 0;
 441                 for (i = 0; i < bc->n; ++i) {
 442                         bcf_callret1_t *p = bcr + i;
 443                         if (dp) dp[i] = p->depth  < 0xffff? p->depth  : 0xffff;
 444                         if (dv) dv[i] = p->n_supp < 0xffff? p->n_supp : 0xffff;
 445                         if (sp) {
 446                                 if (p->anno[0] + p->anno[1] < 2 || p->anno[2] + p->anno[3] < 2
 447                                         || p->anno[0] + p->anno[2] < 2 || p->anno[1] + p->anno[3] < 2)
 448                                 {
 449                                         sp[i] = 0;
 450                                 } else {
 451                                         double left, right, two;
 452                                         int x;
 453                                         kt_fisher_exact(p->anno[0], p->anno[1], p->anno[2], p->anno[3], &left, &right, &two);
 454                                         x = (int)(-4.343 * log(two) + .499);
 455                                         if (x > 255) x = 255;
 456                                         sp[i] = x;
 457                                 }
 458                         }
 459                 }
 460         }
 461         return 0;
 462 }