samtools/bcftools/mut.c.pysam.c

   1 #include "pysam.h"
   2
   3 #include <stdlib.h>
   4 #include <stdint.h>
   5 #include "bcf.h"
   6
   7 #define MAX_GENO 359
   8
   9 int8_t seq_bitcnt[] = { 4, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4 };
  10 char *seq_nt16rev = "XACMGRSVTWYHKDBN";
  11
  12 uint32_t *bcf_trio_prep(int is_x, int is_son)
  13 {
  14         int i, j, k, n, map[10];
  15         uint32_t *ret;
  16         ret = calloc(MAX_GENO, 4);
  17         for (i = 0, k = 0; i < 4; ++i)
  18                 for (j = i; j < 4; ++j)
  19                         map[k++] = 1<<i|1<<j;
  20         for (i = 0, n = 1; i < 10; ++i) { // father
  21                 if (is_x && seq_bitcnt[map[i]] != 1) continue;
  22                 if (is_x && is_son) {
  23                         for (j = 0; j < 10; ++j) // mother
  24                                 for (k = 0; k < 10; ++k) // child
  25                                         if (seq_bitcnt[map[k]] == 1 && (map[j]&map[k]))
  26                                                 ret[n++] = j<<16 | i<<8 | k;
  27                 } else {
  28                         for (j = 0; j < 10; ++j) // mother
  29                                 for (k = 0; k < 10; ++k) // child
  30                                         if ((map[i]&map[k]) && (map[j]&map[k]) && ((map[i]|map[j])&map[k]) == map[k])
  31                                                 ret[n++] = j<<16 | i<<8 | k;
  32                 }
  33         }
  34         ret[0] = n - 1;
  35         return ret;
  36 }
  37
  38
  39 int bcf_trio_call(const uint32_t *prep, const bcf1_t *b, int *llr, int64_t *gt)
  40 {
  41         int i, j, k;
  42         const bcf_ginfo_t *PL;
  43         uint8_t *gl10;
  44         int map[10];
  45         if (b->n_smpl != 3) return -1; // not a trio
  46         for (i = 0; i < b->n_gi; ++i)
  47                 if (b->gi[i].fmt == bcf_str2int("PL", 2)) break;
  48         if (i == b->n_gi) return -1; // no PL
  49         gl10 = alloca(10 * b->n_smpl);
  50         if (bcf_gl10(b, gl10) < 0) {
  51                 if (bcf_gl10_indel(b, gl10) < 0) return -1;
  52         }
  53         PL = b->gi + i;
  54         for (i = 0, k = 0; i < 4; ++i)
  55                 for (j = i; j < 4; ++j)
  56                         map[k++] = seq_nt16rev[1<<i|1<<j];
  57         for (j = 0; j < 3; ++j) // check if ref hom is the most probable in all members
  58                 if (((uint8_t*)PL->data)[j * PL->len] != 0) break;
  59         if (j < 3) { // we need to go through the complex procedure
  60                 uint8_t *g[3];
  61                 int minc = 1<<30, minc_j = -1, minf = 0, gtf = 0, gtc = 0;
  62                 g[0] = gl10;
  63                 g[1] = gl10 + 10;
  64                 g[2] = gl10 + 20;
  65                 for (j = 1; j <= (int)prep[0]; ++j) { // compute LK with constraint
  66                         int sum = g[0][prep[j]&0xff] + g[1][prep[j]>>8&0xff] + g[2][prep[j]>>16&0xff];
  67                         if (sum < minc) minc = sum, minc_j = j;
  68                 }
  69                 gtc |= map[prep[minc_j]&0xff]; gtc |= map[prep[minc_j]>>8&0xff]<<8; gtc |= map[prep[minc_j]>>16]<<16;
  70                 for (j = 0; j < 3; ++j) { // compute LK without constraint
  71                         int min = 1<<30, min_k = -1;
  72                         for (k = 0; k < 10; ++k)
  73                                 if (g[j][k] < min) min = g[j][k], min_k = k;
  74                         gtf |= map[min_k]<<(j*8);
  75                         minf += min;
  76                 }
  77                 *llr = minc - minf; *gt = (int64_t)gtc<<32 | gtf;
  78         } else *llr = 0, *gt = -1;
  79         return 0;
  80 }
  81
  82 int bcf_pair_call(const bcf1_t *b)
  83 {
  84         int i, j, k;
  85         const bcf_ginfo_t *PL;
  86         if (b->n_smpl != 2) return -1; // not a pair
  87         for (i = 0; i < b->n_gi; ++i)
  88                 if (b->gi[i].fmt == bcf_str2int("PL", 2)) break;
  89         if (i == b->n_gi) return -1; // no PL
  90         PL = b->gi + i;
  91         for (j = 0; j < 2; ++j) // check if ref hom is the most probable in all members
  92                 if (((uint8_t*)PL->data)[j * PL->len] != 0) break;
  93         if (j < 2) { // we need to go through the complex procedure
  94                 uint8_t *g[2];
  95                 int minc = 1<<30, minf = 0;
  96                 g[0] = PL->data;
  97                 g[1] = (uint8_t*)PL->data + PL->len;
  98                 for (j = 0; j < PL->len; ++j) // compute LK with constraint
  99                         minc = minc < g[0][j] + g[1][j]? minc : g[0][j] + g[1][j];
 100                 for (j = 0; j < 2; ++j) { // compute LK without constraint
 101                         int min = 1<<30;
 102                         for (k = 0; k < PL->len; ++k)
 103                                 min = min < g[j][k]? min : g[j][k];
 104                         minf += min;
 105                 }
 106                 return minc - minf;
 107         } else return 0;
 108 }
 109
 110 int bcf_min_diff(const bcf1_t *b)
 111 {
 112         int i, min = 1<<30;
 113         const bcf_ginfo_t *PL;
 114         for (i = 0; i < b->n_gi; ++i)
 115                 if (b->gi[i].fmt == bcf_str2int("PL", 2)) break;
 116         if (i == b->n_gi) return -1; // no PL
 117         PL = b->gi + i;
 118         for (i = 0; i < b->n_smpl; ++i) {
 119                 int m1, m2, j;
 120                 const uint8_t *p = (uint8_t*)PL->data;
 121                 m1 = m2 = 1<<30;
 122                 for (j = 0; j < PL->len; ++j) {
 123                         if ((int)p[j] < m1) m2 = m1, m1 = p[j];
 124                         else if ((int)p[j] < m2) m2 = p[j];
 125                 }
 126                 min = min < m2 - m1? min : m2 - m1;
 127         }
 128         return min;
 129 }