Mercurial > repos > lsong10 > psiclass
comparison PsiCLASS-1.0.2/samtools-0.1.19/bcftools/bcf2qcall.c @ 0:903fc43d6227 draft default tip
Uploaded
| author | lsong10 |
|---|---|
| date | Fri, 26 Mar 2021 16:52:45 +0000 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:903fc43d6227 |
|---|---|
| 1 #include <errno.h> | |
| 2 #include <math.h> | |
| 3 #include <string.h> | |
| 4 #include <stdlib.h> | |
| 5 #include "bcf.h" | |
| 6 | |
| 7 static int8_t nt4_table[256] = { | |
| 8 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, | |
| 9 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, | |
| 10 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 /*'-'*/, 4, 4, | |
| 11 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, | |
| 12 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, | |
| 13 4, 4, 4, 4, 3, 4, 4, 4, -1, 4, 4, 4, 4, 4, 4, 4, | |
| 14 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, | |
| 15 4, 4, 4, 4, 3, 4, 4, 4, -1, 4, 4, 4, 4, 4, 4, 4, | |
| 16 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, | |
| 17 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, | |
| 18 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, | |
| 19 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, | |
| 20 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, | |
| 21 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, | |
| 22 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, | |
| 23 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 | |
| 24 }; | |
| 25 | |
| 26 static int read_I16(bcf1_t *b, int anno[16]) | |
| 27 { | |
| 28 char *p; | |
| 29 int i; | |
| 30 if ((p = strstr(b->info, "I16=")) == 0) return -1; | |
| 31 p += 4; | |
| 32 for (i = 0; i < 16; ++i) { | |
| 33 anno[i] = strtol(p, &p, 10); | |
| 34 if (anno[i] == 0 && (errno == EINVAL || errno == ERANGE)) return -2; | |
| 35 ++p; | |
| 36 } | |
| 37 return 0; | |
| 38 } | |
| 39 | |
| 40 int bcf_2qcall(bcf_hdr_t *h, bcf1_t *b) | |
| 41 { | |
| 42 int a[4], k, g[10], l, map[4], k1, j, i, i0, anno[16], dp, mq, d_rest; | |
| 43 char *s; | |
| 44 if (b->ref[1] != 0 || b->n_alleles > 4) return -1; // ref is not a single base | |
| 45 for (i = 0; i < b->n_gi; ++i) | |
| 46 if (b->gi[i].fmt == bcf_str2int("PL", 2)) break; | |
| 47 if (i == b->n_gi) return -1; // no PL | |
| 48 if (read_I16(b, anno) != 0) return -1; // no I16; FIXME: can be improved | |
| 49 d_rest = dp = anno[0] + anno[1] + anno[2] + anno[3]; | |
| 50 if (dp == 0) return -1; // depth is zero | |
| 51 mq = (int)(sqrt((double)(anno[9] + anno[11]) / dp) + .499); | |
| 52 i0 = i; | |
| 53 a[0] = nt4_table[(int)b->ref[0]]; | |
| 54 if (a[0] > 3) return -1; // ref is not A/C/G/T | |
| 55 a[1] = a[2] = a[3] = -2; // -1 has a special meaning | |
| 56 if (b->alt[0] == 0) return -1; // no alternate allele | |
| 57 map[0] = map[1] = map[2] = map[3] = -2; | |
| 58 map[a[0]] = 0; | |
| 59 for (k = 0, s = b->alt, k1 = -1; k < 3 && *s; ++k, s += 2) { | |
| 60 if (s[1] != ',' && s[1] != 0) return -1; // ALT is not single base | |
| 61 a[k+1] = nt4_table[(int)*s]; | |
| 62 if (a[k+1] >= 0) map[a[k+1]] = k+1; | |
| 63 else k1 = k+1; | |
| 64 if (s[1] == 0) break; | |
| 65 } | |
| 66 for (k = 0; k < 4; ++k) | |
| 67 if (map[k] < 0) map[k] = k1; | |
| 68 for (i = 0; i < h->n_smpl; ++i) { | |
| 69 int d; | |
| 70 uint8_t *p = b->gi[i0].data + i * b->gi[i0].len; | |
| 71 for (j = 0; j < b->gi[i0].len; ++j) | |
| 72 if (p[j]) break; | |
| 73 d = (int)((double)d_rest / (h->n_smpl - i) + .499); | |
| 74 if (d == 0) d = 1; | |
| 75 if (j == b->gi[i0].len) d = 0; | |
| 76 d_rest -= d; | |
| 77 for (k = j = 0; k < 4; ++k) { | |
| 78 for (l = k; l < 4; ++l) { | |
| 79 int t, x = map[k], y = map[l]; | |
| 80 if (x > y) t = x, x = y, y = t; // swap | |
| 81 g[j++] = p[y * (y+1) / 2 + x]; | |
| 82 } | |
| 83 } | |
| 84 printf("%s\t%d\t%c", h->ns[b->tid], b->pos+1, *b->ref); | |
| 85 printf("\t%d\t%d\t0", d, mq); | |
| 86 for (j = 0; j < 10; ++j) | |
| 87 printf("\t%d", g[j]); | |
| 88 printf("\t%s\n", h->sns[i]); | |
| 89 } | |
| 90 return 0; | |
| 91 } |