#include <assert.h>
#include <ctype.h>
+#include <string.h>
#include "bam.h"
#include "bam2bcf.h"
#include "ksort.h"
#include "kaln.h"
+#include "khash.h"
+KHASH_SET_INIT_STR(rg)
#define MINUS_CONST 0x10000000
#define INDEL_WINDOW_SIZE 50
+void *bcf_call_add_rg(void *_hash, const char *hdtext, const char *list)
+{
+ const char *s, *p, *q, *r, *t;
+ khash_t(rg) *hash;
+ if (list == 0 || hdtext == 0) return _hash;
+ if (_hash == 0) _hash = kh_init(rg);
+ hash = (khash_t(rg)*)_hash;
+ if ((s = strstr(hdtext, "@RG\t")) == 0) return hash;
+ do {
+ t = strstr(s + 4, "@RG\t"); // the next @RG
+ if ((p = strstr(s, "\tID:")) != 0) p += 4;
+ if ((q = strstr(s, "\tPL:")) != 0) q += 4;
+ if (p && q && (t == 0 || (p < t && q < t))) { // ID and PL are both present
+ int lp, lq;
+ char *x;
+ for (r = p; *r && *r != '\t' && *r != '\n'; ++r); lp = r - p;
+ for (r = q; *r && *r != '\t' && *r != '\n'; ++r); lq = r - q;
+ x = calloc((lp > lq? lp : lq) + 1, 1);
+ for (r = q; *r && *r != '\t' && *r != '\n'; ++r) x[r-q] = *r;
+ if (strstr(list, x)) { // insert ID to the hash table
+ khint_t k;
+ int ret;
+ for (r = p; *r && *r != '\t' && *r != '\n'; ++r) x[r-p] = *r;
+ x[r-p] = 0;
+ k = kh_get(rg, hash, x);
+ if (k == kh_end(hash)) k = kh_put(rg, hash, x, &ret);
+ else free(x);
+ } else free(x);
+ }
+ s = t;
+ } while (s);
+ return hash;
+}
+
+void bcf_call_del_rghash(void *_hash)
+{
+ khint_t k;
+ khash_t(rg) *hash = (khash_t(rg)*)_hash;
+ if (hash == 0) return;
+ for (k = kh_begin(hash); k < kh_end(hash); ++k)
+ if (kh_exist(hash, k))
+ free((char*)kh_key(hash, k));
+ kh_destroy(rg, hash);
+}
+
static int tpos2qpos(const bam1_core_t *c, const uint32_t *cigar, int32_t tpos, int is_left, int32_t *_tpos)
{
int k, x = c->pos, y = 0, last_y = 0;
return max_i - pos;
}
-int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref)
+int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref,
+ const void *rghash)
{
extern void ks_introsort_uint32_t(int, uint32_t*);
- int i, s, j, k, t, n_types, *types, max_rd_len, left, right, max_ins, *score, N, K, l_run, ref_type;
+ int i, s, j, k, t, n_types, *types, max_rd_len, left, right, max_ins, *score, N, K, l_run, ref_type, n_alt;
char *inscns = 0, *ref2, *query;
+ khash_t(rg) *hash = (khash_t(rg)*)rghash;
if (ref == 0 || bca == 0) return -1;
+ // mark filtered reads
+ if (rghash) {
+ N = 0;
+ for (s = N = 0; s < n; ++s) {
+ for (i = 0; i < n_plp[s]; ++i) {
+ bam_pileup1_t *p = plp[s] + i;
+ const uint8_t *rg = bam_aux_get(p->b, "RG");
+ p->aux = 1; // filtered by default
+ if (rg) {
+ khint_t k = kh_get(rg, hash, (const char*)(rg + 1));
+ if (k != kh_end(hash)) p->aux = 0, ++N; // not filtered
+ }
+ }
+ }
+ if (N == 0) return -1; // no reads left
+ }
// determine if there is a gap
for (s = N = 0; s < n; ++s) {
- N += n_plp[s]; // N is the total number of reads
for (i = 0; i < n_plp[s]; ++i)
if (plp[s][i].indel != 0) break;
if (i < n_plp[s]) break;
}
if (s == n) return -1; // there is no indel at this position.
+ for (s = N = 0; s < n; ++s) N += n_plp[s]; // N is the total number of reads
{ // find out how many types of indels are present
int m;
uint32_t *aux;
for (s = 0; s < n; ++s) {
for (i = 0; i < n_plp[s]; ++i) {
const bam_pileup1_t *p = plp[s] + i;
- if (p->indel != 0)
+ if (p->indel != 0 && (rghash == 0 || p->aux == 0))
aux[m++] = MINUS_CONST + p->indel;
j = bam_cigar2qlen(&p->b->core, bam1_cigar(p->b));
if (j > max_rd_len) max_rd_len = j;
// squeeze out identical types
for (i = 1, n_types = 1; i < m; ++i)
if (aux[i] != aux[i-1]) ++n_types;
- assert(n_types > 1); // there must at least one type of non-reference indel
+ if (n_types == 1) { // no indels
+ free(aux); return -1;
+ }
types = (int*)calloc(n_types, sizeof(int));
t = 0;
types[t++] = aux[0] - MINUS_CONST;
query[l - qbeg] = bam_nt16_nt4_table[bam1_seqi(seq, l)];
// do alignment; this takes most of computing time for indel calling
sc = ka_global_score((uint8_t*)ref2 + tbeg - left, tend - tbeg + abs(types[t]),
- (uint8_t*)query + qbeg, qend - qbeg, &ap);
+ (uint8_t*)query, qend - qbeg, &ap);
score[K*n_types + t] = -sc;
/*
for (l = 0; l < tend - tbeg + abs(types[t]); ++l)
fputc("ACGTN"[(int)ref2[tbeg-left+l]], stderr);
fputc('\n', stderr);
- for (l = 0; l < qend - qbeg; ++l) fputc("ACGTN"[(int)query[qbeg + l]], stderr);
+ for (l = 0; l < qend - qbeg; ++l) fputc("ACGTN"[(int)query[l]], stderr);
fputc('\n', stderr);
- fprintf(stderr, "pos=%d type=%d read=%d:%d name=%s score=%d\n", pos, types[t], s, i, bam1_qname(p->b), sc);
+ fprintf(stderr, "pos=%d type=%d read=%d:%d name=%s qbeg=%d tbeg=%d score=%d\n", pos, types[t], s, i, bam1_qname(p->b), qbeg, tbeg, sc);
*/
}
}
for (s = K = 0; s < n; ++s) {
for (i = 0; i < n_plp[s]; ++i, ++K) {
bam_pileup1_t *p = plp[s] + i;
- int *sct = &score[K*n_types], indelQ;
+ int *sct = &score[K*n_types], indelQ, seqQ;
for (t = 0; t < n_types; ++t) sc[t] = sct[t]<<6 | t;
for (t = 1; t < n_types; ++t) // insertion sort
for (j = t; j > 0 && sc[j] < sc[j-1]; --j)
*/
if ((sc[0]&0x3f) == ref_type) {
indelQ = (sc[1]>>6) - (sc[0]>>6);
- tmp = est_seqQ(bca, types[sc[1]&0x3f], l_run);
+ seqQ = est_seqQ(bca, types[sc[1]&0x3f], l_run);
} else {
for (t = 0; t < n_types; ++t) // look for the reference type
if ((sc[t]&0x3f) == ref_type) break;
indelQ = (sc[t]>>6) - (sc[0]>>6);
- tmp = est_seqQ(bca, types[sc[0]&0x3f], l_run);
+ seqQ = est_seqQ(bca, types[sc[0]&0x3f], l_run);
}
- if (indelQ > tmp) indelQ = tmp;
- p->aux = (sc[0]&0x3f)<<8 | indelQ;
- sumq[sc[0]&0x3f] += indelQ;
+ p->aux = (sc[0]&0x3f)<<16 | seqQ<<8 | indelQ;
+ sumq[sc[0]&0x3f] += indelQ < seqQ? indelQ : seqQ;
// fprintf(stderr, "pos=%d read=%d:%d name=%s call=%d q=%d\n", pos, s, i, bam1_qname(p->b), types[sc[0]&0x3f], indelQ);
}
}
memcpy(&bca->inscns[t * bca->maxins], &inscns[(sumq[t]&0x3f) * max_ins], bca->maxins);
}
// update p->aux
- for (s = K = 0; s < n; ++s) {
- for (i = 0; i < n_plp[s]; ++i, ++K) {
+ for (s = n_alt = 0; s < n; ++s) {
+ for (i = 0; i < n_plp[s]; ++i) {
bam_pileup1_t *p = plp[s] + i;
- int x = types[p->aux>>8&0x3f];
+ int x = types[p->aux>>16&0x3f];
for (j = 0; j < 4; ++j)
if (x == bca->indel_types[j]) break;
- p->aux = j<<8 | (j == 4? 0 : (p->aux&0xff));
-// fprintf(stderr, "pos=%d read=%d:%d name=%s call=%d q=%d\n", pos, s, i, bam1_qname(p->b), p->aux>>8&63, p->aux&0xff);
+ p->aux = j<<16 | (j == 4? 0 : (p->aux&0xffff));
+ if ((p->aux>>16&0x3f) > 0) ++n_alt;
+// fprintf(stderr, "pos=%d read=%d:%d name=%s call=%d q=%d\n", pos, s, i, bam1_qname(p->b), p->aux>>16&63, p->aux&0xff);
}
}
}
free(score);
// free
free(types); free(inscns);
- return 0;
+ return n_alt > 0? 0 : -1;
}
projects / samtools.git / blobdiff