annotate PsiCLASS-1.0.2/samtools-0.1.19/cut_target.c @ 0:903fc43d6227 draft default tip

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