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