samtools/cut_target.c.pysam.c

   1 #include "pysam.h"
   2
   3 #include <unistd.h>
   4 #include <stdlib.h>
   5 #include <string.h>
   6 #include "bam.h"
   7 #include "errmod.h"
   8 #include "faidx.h"
   9
  10 #define ERR_DEP 0.83f
  11
  12 typedef struct {
  13         int e[2][3], p[2][2];
  14 } score_param_t;
  15
  16 /* Note that although the two matrics have 10 parameters in total, only 4
  17  * (probably 3) are free.  Changing the scoring matrices in a sort of symmetric
  18  * way will not change the result. */
  19 static score_param_t g_param = { {{0,0,0},{-4,1,6}}, {{0,-14000}, {0,0}} };
  20
  21 typedef struct {
  22         int min_baseQ, tid, max_bases;
  23         uint16_t *bases;
  24         bamFile fp;
  25         bam_header_t *h;
  26         char *ref;
  27         faidx_t *fai;
  28         errmod_t *em;
  29 } ct_t;
  30
  31 static uint16_t gencns(ct_t *g, int n, const bam_pileup1_t *plp)
  32 {
  33         int i, j, ret, tmp, k, sum[4], qual;
  34         float q[16];
  35         if (n > g->max_bases) { // enlarge g->bases
  36                 g->max_bases = n;
  37                 kroundup32(g->max_bases);
  38                 g->bases = realloc(g->bases, g->max_bases * 2);
  39         }
  40         for (i = k = 0; i < n; ++i) {
  41                 const bam_pileup1_t *p = plp + i;
  42                 uint8_t *seq;
  43                 int q, baseQ, b;
  44                 if (p->is_refskip || p->is_del) continue;
  45                 baseQ = bam1_qual(p->b)[p->qpos];
  46                 if (baseQ < g->min_baseQ) continue;
  47                 seq = bam1_seq(p->b);
  48                 b = bam_nt16_nt4_table[bam1_seqi(seq, p->qpos)];
  49                 if (b > 3) continue;
  50                 q = baseQ < p->b->core.qual? baseQ : p->b->core.qual;
  51                 if (q < 4) q = 4;
  52                 if (q > 63) q = 63;
  53                 g->bases[k++] = q<<5 | bam1_strand(p->b)<<4 | b;
  54         }
  55         if (k == 0) return 0;
  56         errmod_cal(g->em, k, 4, g->bases, q);
  57         for (i = 0; i < 4; ++i) sum[i] = (int)(q[i<<2|i] + .499) << 2 | i;
  58         for (i = 1; i < 4; ++i) // insertion sort
  59                 for (j = i; j > 0 && sum[j] < sum[j-1]; --j)
  60                         tmp = sum[j], sum[j] = sum[j-1], sum[j-1] = tmp;
  61         qual = (sum[1]>>2) - (sum[0]>>2);
  62         k = k < 256? k : 255;
  63         ret = (qual < 63? qual : 63) << 2 | (sum[0]&3);
  64         return ret<<8|k;
  65 }
  66
  67 static void process_cns(bam_header_t *h, int tid, int l, uint16_t *cns)
  68 {
  69         int i, f[2][2], *prev, *curr, *swap_tmp, s;
  70         uint8_t *b; // backtrack array
  71         b = calloc(l, 1);
  72         f[0][0] = f[0][1] = 0;
  73         prev = f[0]; curr = f[1];
  74         // fill the backtrack matrix
  75         for (i = 0; i < l; ++i) {
  76                 int c = (cns[i] == 0)? 0 : (cns[i]>>8 == 0)? 1 : 2;
  77                 int tmp0, tmp1;
  78                 // compute f[0]
  79                 tmp0 = prev[0] + g_param.e[0][c] + g_param.p[0][0]; // (s[i+1],s[i])=(0,0)
  80                 tmp1 = prev[1] + g_param.e[0][c] + g_param.p[1][0]; // (0,1)
  81                 if (tmp0 > tmp1) curr[0] = tmp0, b[i] = 0;
  82                 else curr[0] = tmp1, b[i] = 1;
  83                 // compute f[1]
  84                 tmp0 = prev[0] + g_param.e[1][c] + g_param.p[0][1]; // (s[i+1],s[i])=(1,0)
  85                 tmp1 = prev[1] + g_param.e[1][c] + g_param.p[1][1]; // (1,1)
  86                 if (tmp0 > tmp1) curr[1] = tmp0, b[i] |= 0<<1;
  87                 else curr[1] = tmp1, b[i] |= 1<<1;
  88                 // swap
  89                 swap_tmp = prev; prev = curr; curr = swap_tmp;
  90         }
  91         // backtrack
  92         s = prev[0] > prev[1]? 0 : 1;
  93         for (i = l - 1; i > 0; --i) {
  94                 b[i] |= s<<2;
  95                 s = b[i]>>s&1;
  96         }
  97         // print
  98         for (i = 0, s = -1; i <= l; ++i) {
  99                 if (i == l || ((b[i]>>2&3) == 0 && s >= 0)) {
 100                         if (s >= 0) {
 101                                 int j;
 102                                 printf("%s:%d-%d\t0\t%s\t%d\t60\t%dM\t*\t0\t0\t", h->target_name[tid], s+1, i, h->target_name[tid], s+1, i-s);
 103                                 for (j = s; j < i; ++j) {
 104                                         int c = cns[j]>>8;
 105                                         if (c == 0) putchar('N');
 106                                         else putchar("ACGT"[c&3]);
 107                                 }
 108                                 putchar('\t');
 109                                 for (j = s; j < i; ++j)
 110                                         putchar(33 + (cns[j]>>8>>2));
 111                                 putchar('\n');
 112                         }
 113                         //if (s >= 0) printf("%s\t%d\t%d\t%d\n", h->target_name[tid], s, i, i - s);
 114                         s = -1;
 115                 } else if ((b[i]>>2&3) && s < 0) s = i;
 116         }
 117         free(b);
 118 }
 119
 120 static int read_aln(void *data, bam1_t *b)
 121 {
 122         extern int bam_prob_realn_core(bam1_t *b, const char *ref, int flag);
 123         ct_t *g = (ct_t*)data;
 124         int ret, len;
 125         ret = bam_read1(g->fp, b);
 126         if (ret >= 0 && g->fai && b->core.tid >= 0 && (b->core.flag&4) == 0) {
 127                 if (b->core.tid != g->tid) { // then load the sequence
 128                         free(g->ref);
 129                         g->ref = fai_fetch(g->fai, g->h->target_name[b->core.tid], &len);
 130                         g->tid = b->core.tid;
 131                 }
 132                 bam_prob_realn_core(b, g->ref, 1<<1|1);
 133         }
 134         return ret;
 135 }
 136
 137 int main_cut_target(int argc, char *argv[])
 138 {
 139         int c, tid, pos, n, lasttid = -1, lastpos = -1, l, max_l;
 140         const bam_pileup1_t *p;
 141         bam_plp_t plp;
 142         uint16_t *cns;
 143         ct_t g;
 144
 145         memset(&g, 0, sizeof(ct_t));
 146         g.min_baseQ = 13; g.tid = -1;
 147         while ((c = getopt(argc, argv, "f:Q:i:o:0:1:2:")) >= 0) {
 148                 switch (c) {
 149                         case 'Q': g.min_baseQ = atoi(optarg); break; // quality cutoff
 150                         case 'i': g_param.p[0][1] = -atoi(optarg); break; // 0->1 transition (in) PENALTY
 151                         case '0': g_param.e[1][0] = atoi(optarg); break; // emission SCORE
 152                         case '1': g_param.e[1][1] = atoi(optarg); break;
 153                         case '2': g_param.e[1][2] = atoi(optarg); break;
 154                         case 'f': g.fai = fai_load(optarg);
 155                                 if (g.fai == 0) fprintf(pysamerr, "[%s] fail to load the fasta index.\n", __func__);
 156                                 break;
 157                 }
 158         }
 159         if (argc == optind) {
 160                 fprintf(pysamerr, "Usage: samtools targetcut [-Q minQ] [-i inPen] [-0 em0] [-1 em1] [-2 em2] [-f ref] <in.bam>\n");
 161                 return 1;
 162         }
 163         l = max_l = 0; cns = 0;
 164         g.fp = strcmp(argv[optind], "-")? bam_open(argv[optind], "r") : bam_dopen(fileno(stdin), "r");
 165         g.h = bam_header_read(g.fp);
 166         g.em = errmod_init(1 - ERR_DEP);
 167         plp = bam_plp_init(read_aln, &g);
 168         while ((p = bam_plp_auto(plp, &tid, &pos, &n)) != 0) {
 169                 if (tid < 0) break;
 170                 if (tid != lasttid) { // change of chromosome
 171                         if (cns) process_cns(g.h, lasttid, l, cns);
 172                         if (max_l < g.h->target_len[tid]) {
 173                                 max_l = g.h->target_len[tid];
 174                                 kroundup32(max_l);
 175                                 cns = realloc(cns, max_l * 2);
 176                         }
 177                         l = g.h->target_len[tid];
 178                         memset(cns, 0, max_l * 2);
 179                         lasttid = tid;
 180                 }
 181                 cns[pos] = gencns(&g, n, p);
 182                 lastpos = pos;
 183         }
 184         process_cns(g.h, lasttid, l, cns);
 185         free(cns);
 186         bam_header_destroy(g.h);
 187         bam_plp_destroy(plp);
 188         bam_close(g.fp);
 189         if (g.fai) {
 190                 fai_destroy(g.fai); free(g.ref);
 191         }
 192         errmod_destroy(g.em);
 193         free(g.bases);
 194         return 0;
 195 }