diff --git a/align.cpp b/align.cpp index 3b4fa2b..4502f94 100644 --- a/align.cpp +++ b/align.cpp @@ -166,8 +166,8 @@ void SingleAlign::ConvertBinaySeq() { void SingleAlign::SnpAlign(RefSeq &ref, bit32_t mode) { - bit32_t i,j,jj1,jj,m, modeindex, cmodeindex, mc, h, read_chain_index_mask, *_refloc0; - Hit prefetch_hit, _hHit; + bit32_t i,j,jj,m, modeindex, cmodeindex, mc, h, read_chain_index_mask, *_refloc0; + Hit _hHit; //cout<<_pread->seq<seq.size(); for(i=0; iseq.size()-4;pos++){ @@ -416,7 +416,7 @@ int SingleAlign::RunAlign(RefSeq &ref) { } void SingleAlign::ReorderSeed(RefSeq &ref) { - bit32_t i,ii,s, total, tt; + bit32_t i,ii, total, tt; //for(i=0;i>offset)¶m.XC64(*s)^(*s))&((q[FIXELEMENT])>>offset)))>snp_thres) return; + if((tmp_snp=N_count+param.XM64(((((*q)>>offset)¶m.XC64(*s))^(*s))&((q[FIXELEMENT])>>offset)))>snp_thres) return; for(bit32_t i=1;i>offset)¶m.XC64(s[i])^s[i])&(((q[FIXELEMENT-1+i]<<1)<<(63-offset))|q[i+FIXELEMENT]>>offset)))>snp_thres) return; + if((tmp_snp+=param.XM64((((((q[i-1]<<1)<<(63-offset))|q[i]>>offset)¶m.XC64(s[i]))^s[i])&(((q[FIXELEMENT-1+i]<<1)<<(63-offset))|q[i+FIXELEMENT]>>offset)))>snp_thres) return; } /* tmp_snp=N_count; @@ -134,7 +134,7 @@ inline void SingleAlign::CountMismatch(bit64_t *q, bit32_t offset, bit64_t *s) { inline bit32_t SingleAlign::MismatchPattern0(bit64_t *q, bit64_t *s, bit32_t offset) { register bit64_t tmp; - bit32_t *mm_array; int i, j,jj, ss=0; + bit32_t *mm_array,i;int j,jj, ss=0; //cout<<"in MMpattern() gap_index:"< "< "<seq.size()-param.index_interval+1)/param.seed_size),_sa.read_max_snp_num+1); _sb.seedseg_num=min((bit32_t)((_sb._pread->seq.size()-param.index_interval+1)/param.seed_size),_sb.read_max_snp_num+1); - seedseg_num=max(_sa.seedseg_num,_sb.seedseg_num); _sa.ConvertBinaySeq(); _sb.ConvertBinaySeq(); _sa.snp_thres=_sa.read_max_snp_num; diff --git a/param.cpp b/param.cpp index f7796be..ff12d9c 100644 --- a/param.cpp +++ b/param.cpp @@ -1,49 +1,49 @@ -#include "param.h" -#include - -using namespace std; - -Param::Param() -{ +#include "param.h" +#include + +using namespace std; + +Param::Param() +{ num_procs=sysconf(_SC_NPROCESSORS_ONLN); - if(num_procs>8) num_procs=8; - - max_dbseq_size=0x1000000; //16Mb - append_dbseq_size=0x1000000; //16Mb - - max_ns = 5; - trim_lowQ=0; - - zero_qual= '!'; - qual_threshold= 0; - default_qual=40; - - min_insert= 28; - max_insert= 1000; - - SetSeedSize(16); - //seed_size= 16; - //seed_bits=(1<<(seed_size*2))-1; - - max_snp_num = 108; - max_num_hits = MAXHITS>100?100:MAXHITS; - max_kmer_ratio = 5e-7; - - min_read_size=seed_size; - input_format=gz_input=gz_ref=-1; - n_adapter=0; - - report_repeat_hits = 1; - - useful_nt="ACGTacgt"; - nx_nt="NXnx"; - - for(bit32_t i=0; i8) num_procs=8; + + max_dbseq_size=0x1000000; //16Mb + append_dbseq_size=0x1000000; //16Mb + + max_ns = 5; + trim_lowQ=0; + + zero_qual= '!'; + qual_threshold= 0; + default_qual=40; + + min_insert= 28; + max_insert= 1000; + + SetSeedSize(16); + //seed_size= 16; + //seed_bits=(1<<(seed_size*2))-1; + + max_snp_num = 108; + max_num_hits = MAXHITS>100?100:MAXHITS; + max_kmer_ratio = 5e-7; + + min_read_size=seed_size; + input_format=gz_input=gz_ref=-1; + n_adapter=0; + + report_repeat_hits = 1; + + useful_nt="ACGTacgt"; + nx_nt="NXnx"; + + for(bit32_t i=0; i0); + int IUPAC_count[256], i, j, ii[256], l; + char IUPAC[256][4]; + for(i=0;i<256;++i) {IUPAC_count[i]=0; for(j=0;j<4;++j) IUPAC[i][j]=(char) 0;} + IUPAC['A'][0]='A'; IUPAC['C'][0]='C'; IUPAC['G'][0]='G'; IUPAC['T'][0]='T'; + IUPAC['N'][0]='A'; IUPAC['N'][1]='C'; IUPAC['N'][2]='G'; IUPAC['N'][3]='T'; + strcpy(IUPAC['R'], "AG"); strcpy(IUPAC['Y'], "CT"); strcpy(IUPAC['S'], "CG"); + strcpy(IUPAC['W'], "AT"); strcpy(IUPAC['K'], "GT"); strcpy(IUPAC['M'], "AC"); + strcpy(IUPAC['B'], "CGT"); strcpy(IUPAC['D'], "AGT"); strcpy(IUPAC['H'], "ACT"); strcpy(IUPAC['V'], "ACG"); + for(i=0;i<256;++i) for(j=0;j<4;++j) IUPAC_count[i]+=(IUPAC[i][j]>0); //cout<=IUPAC_count[ds[j]]&&j=IUPAC_count[(unsigned char)ds[j]]&&j16||n<10) {cerr<<"seed size must be between 10 and 16\n"; exit(1);} - seed_size=n; - seed_bits_lz=(SEGLEN-seed_size)*2; - min_read_size=seed_size+index_interval-1; - seed_bits=0; - for(bit32_t i=0; i16||n<10) {cerr<<"seed size must be between 10 and 16\n"; exit(1);} + seed_size=n; + seed_bits_lz=(SEGLEN-seed_size)*2; + min_read_size=seed_size+index_interval-1; + seed_bits=0; + for(bit32_t i=0; i 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 diff --git a/samtools/bam_pileup.c b/samtools/bam_pileup.c index 57434e0..26d4c8c 100644 --- a/samtools/bam_pileup.c +++ b/samtools/bam_pileup.c @@ -72,11 +72,10 @@ static inline int resolve_cigar2(bam_pileup1_t *p, uint32_t pos, cstate_t *s) bam1_t *b = p->b; bam1_core_t *c = &b->core; uint32_t *cigar = bam1_cigar(b); - int k, is_head = 0; + int k; // determine the current CIGAR operation // fprintf(stderr, "%s\tpos=%d\tend=%d\t(%d,%d,%d)\n", bam1_qname(b), pos, s->end, s->k, s->x, s->y); if (s->k == -1) { // never processed - is_head = 1; if (c->n_cigar == 1) { // just one operation, save a loop if (_cop(cigar[0]) == BAM_CMATCH || _cop(cigar[0]) == BAM_CEQUAL || _cop(cigar[0]) == BAM_CDIFF) s->k = 0, s->x = c->pos, s->y = 0; } else { // find the first match or deletion diff --git a/samtools/bam_reheader.c b/samtools/bam_reheader.c index 0b52267..9b1d6ec 100644 --- a/samtools/bam_reheader.c +++ b/samtools/bam_reheader.c @@ -8,12 +8,11 @@ int bam_reheader(BGZF *in, const bam_header_t *h, int fd) { BGZF *fp; - bam_header_t *old; int len; uint8_t *buf; if (in->open_mode != 'r') return -1; buf = malloc(BUF_SIZE); - old = bam_header_read(in); + bam_header_read(in); fp = bgzf_fdopen(fd, "w"); bam_header_write(fp, h); if (in->block_offset < in->block_length) { diff --git a/samtools/bam_tview.c b/samtools/bam_tview.c index 4eea955..58f1b9d 100644 --- a/samtools/bam_tview.c +++ b/samtools/bam_tview.c @@ -431,7 +431,6 @@ int bam_tview_main(int argc, char *argv[]) } #else // #ifdef _HAVE_CURSES #include -#warning "No curses library is available; tview is disabled." int bam_tview_main(int argc, char *argv[]) { fprintf(stderr, "[bam_tview_main] The ncurses library is unavailable; tview is not compiled.\n"); diff --git a/samtools/bcftools/call1.c b/samtools/bcftools/call1.c index 3cc4649..47aa2af 100644 --- a/samtools/bcftools/call1.c +++ b/samtools/bcftools/call1.c @@ -281,13 +281,13 @@ static void write_header(bcf_hdr_t *h) kputs("##FORMAT=\n", &str); if (!strstr(str.s, "##FORMAT=\n", &str); - if (!strstr(str.s, "##FORMAT=\n", &str); - if (!strstr(str.s, "##FORMAT=\n", &str); - if (!strstr(str.s, "##FORMAT=\n", &str); - h->l_txt = str.l + 1; h->txt = str.s; + if (!strstr(str.s, "##FORMAT=\n", &str); + if (!strstr(str.s, "##FORMAT=\n", &str); + if (!strstr(str.s, "##FORMAT=\n", &str); + h->l_txt = str.l + 1; h->txt = str.s; } double bcf_pair_freq(const bcf1_t *b0, const bcf1_t *b1, double f[4]); diff --git a/samtools/bcftools/em.c b/samtools/bcftools/em.c index b7dfe1a..dd024e3 100644 --- a/samtools/bcftools/em.c +++ b/samtools/bcftools/em.c @@ -171,13 +171,13 @@ static double lk_ratio_test(int n, int n1, const double *pdg, double f3[3][3]) int bcf_em1(const bcf1_t *b, int n1, int flag, double x[10]) { double *pdg; - int i, n, n2; + int i, n; if (b->n_alleles < 2) return -1; // one allele only // initialization if (n1 < 0 || n1 > b->n_smpl) n1 = 0; if (flag & 1<<7) flag |= 7<<5; // compute group freq if LRT is required if (flag & 0xf<<1) flag |= 0xf<<1; - n = b->n_smpl; n2 = n - n1; + n = b->n_smpl; pdg = get_pdg3(b); if (pdg == 0) return -1; for (i = 0; i < 10; ++i) x[i] = -1.; // set to negative diff --git a/samtools/bgzf.c b/samtools/bgzf.c index 216cd04..e847da3 100644 --- a/samtools/bgzf.c +++ b/samtools/bgzf.c @@ -627,11 +627,11 @@ int bgzf_close(BGZF* fp) if (fp->open_mode == 'w') { if (bgzf_flush(fp) != 0) return -1; { // add an empty block - int count, block_length = deflate_block(fp, 0); + int block_length = deflate_block(fp, 0); #ifdef _USE_KNETFILE - count = fwrite(fp->compressed_block, 1, block_length, fp->x.fpw); + fwrite(fp->compressed_block, 1, block_length, fp->x.fpw); #else - count = fwrite(fp->compressed_block, 1, block_length, fp->file); + fwrite(fp->compressed_block, 1, block_length, fp->file); #endif } #ifdef _USE_KNETFILE diff --git a/samtools/cut_target.c b/samtools/cut_target.c index 26f434f..ec5cfff 100644 --- a/samtools/cut_target.c +++ b/samtools/cut_target.c @@ -93,7 +93,7 @@ static void process_cns(bam_header_t *h, int tid, int l, uint16_t *cns) s = b[i]>>s&1; } // print - for (i = 0, s = -1; i <= l; ++i) { + for (i = 0, s = -1; ; ++i) { if (i == l || ((b[i]>>2&3) == 0 && s >= 0)) { if (s >= 0) { int j; @@ -111,6 +111,8 @@ static void process_cns(bam_header_t *h, int tid, int l, uint16_t *cns) //if (s >= 0) printf("%s\t%d\t%d\t%d\n", h->target_name[tid], s, i, i - s); s = -1; } else if ((b[i]>>2&3) && s < 0) s = i; + if (i==l) + break; } free(b); } @@ -134,7 +136,7 @@ static int read_aln(void *data, bam1_t *b) int main_cut_target(int argc, char *argv[]) { - int c, tid, pos, n, lasttid = -1, lastpos = -1, l, max_l; + int c, tid, pos, n, lasttid = -1, l, max_l; const bam_pileup1_t *p; bam_plp_t plp; uint16_t *cns; @@ -177,7 +179,6 @@ int main_cut_target(int argc, char *argv[]) lasttid = tid; } cns[pos] = gencns(&g, n, p); - lastpos = pos; } process_cns(g.h, lasttid, l, cns); free(cns); diff --git a/samtools/errmod.c b/samtools/errmod.c index fba9a8d..abb7656 100644 --- a/samtools/errmod.c +++ b/samtools/errmod.c @@ -76,7 +76,6 @@ int errmod_cal(const errmod_t *em, int n, int m, uint16_t *bases, float *q) call_aux_t aux; int i, j, k, w[32]; - if (m > m) return -1; memset(q, 0, m * m * sizeof(float)); if (n == 0) return 0; // calculate aux.esum and aux.fsum diff --git a/samtools/faidx.c b/samtools/faidx.c index f0798fc..d293b08 100644 --- a/samtools/faidx.c +++ b/samtools/faidx.c @@ -275,8 +275,12 @@ faidx_t *fai_load(const char *fn) } } else -#endif + { + fp = fopen(str, "rb"); + } +#else fp = fopen(str, "rb"); +#endif if (fp == 0) { fprintf(stderr, "[fai_load] build FASTA index.\n"); fai_build(fn); diff --git a/samtools/kprobaln.c b/samtools/kprobaln.c index 894a2ae..1c390d1 100644 --- a/samtools/kprobaln.c +++ b/samtools/kprobaln.c @@ -72,10 +72,11 @@ kpa_par_t kpa_par_alt = { 0.0001, 0.01, 10 }; int kpa_glocal(const uint8_t *_ref, int l_ref, const uint8_t *_query, int l_query, const uint8_t *iqual, const kpa_par_t *c, int *state, uint8_t *q) { - double **f, **b = 0, *s, m[9], sI, sM, bI, bM, pb; + double **f, **b = 0, *s, m[9], sI, sM, bI, bM; + double pb __attribute__((unused));; float *qual, *_qual; const uint8_t *ref, *query; - int bw, bw2, i, k, is_diff = 0, is_backward = 1, Pr; + int bw, bw2, i, k, is_backward = 1, Pr; /*** initialization ***/ is_backward = state && q? 1 : 0; @@ -214,7 +215,6 @@ int kpa_glocal(const uint8_t *_ref, int l_ref, const uint8_t *_query, int l_quer set_u(k, bw, 0, 0); pb = b[0][k] = sum / s[0]; // if everything works as is expected, pb == 1.0 } - is_diff = fabs(pb - 1.) > 1e-7? 1 : 0; /*** MAP ***/ for (i = 1; i <= l_query; ++i) { double sum = 0., *fi = f[i], *bi = b[i], max = 0.; diff --git a/samtools/ksort.h b/samtools/ksort.h index fa850ab..65c7f6c 100644 --- a/samtools/ksort.h +++ b/samtools/ksort.h @@ -141,7 +141,7 @@ typedef struct { tmp = *l; *l = l[i]; l[i] = tmp; ks_heapadjust_##name(0, i, l); \ } \ } \ - inline void __ks_insertsort_##name(type_t *s, type_t *t) \ + void __ks_insertsort_##name(type_t *s, type_t *t) \ { \ type_t *i, *j, swap_tmp; \ for (i = s + 1; i < t; ++i) \ diff --git a/samtools/misc/md5.c b/samtools/misc/md5.c index 55ae181..108eba9 100644 --- a/samtools/misc/md5.c +++ b/samtools/misc/md5.c @@ -152,7 +152,7 @@ void MD5Final(digest, ctx) MD5Transform(ctx->buf, (uint32_t *) ctx->in); byteReverse((unsigned char *) ctx->buf, 4); memcpy(digest, ctx->buf, 16); - memset(ctx, 0, sizeof(ctx)); /* In case it's sensitive */ + memset(ctx, 0, sizeof(*ctx)); /* In case it's sensitive */ } diff --git a/samtools/misc/seqtk.c b/samtools/misc/seqtk.c index 591ddff..0061d40 100644 --- a/samtools/misc/seqtk.c +++ b/samtools/misc/seqtk.c @@ -140,9 +140,9 @@ int stk_comp(int argc, char *argv[]) } for (k = 0; p && k < p->n; ++k) { int beg = p->a[k]>>32, end = p->a[k]&0xffffffff; - int la, lb, lc, na, nb, nc, cnt[11]; - if (beg > 0) la = seq->seq.s[beg-1], lb = seq_nt16_table[la], lc = bitcnt_table[lb]; - else la = 'a', lb = -1, lc = 0; + int la, lb, na, nb, nc, cnt[11]; + if (beg > 0) la = seq->seq.s[beg-1], lb = seq_nt16_table[la]; + else la = 'a', lb = -1; na = seq->seq.s[beg]; nb = seq_nt16_table[na]; nc = bitcnt_table[nb]; memset(cnt, 0, 11 * sizeof(int)); for (i = beg; i < end; ++i) { @@ -166,7 +166,7 @@ int stk_comp(int argc, char *argv[]) if (b == 10 || b == 5) ++cnt[10]; } } - la = a; lb = b; lc = c; + la = a; lb = b; } if (h) printf("%s\t%d\t%d", seq->name.s, beg, end); else printf("%s\t%d", seq->name.s, l); diff --git a/utilities.cpp b/utilities.cpp index 6da73ad..6cf04f9 100644 --- a/utilities.cpp +++ b/utilities.cpp @@ -49,11 +49,11 @@ bit32_t myrand(int i, bit32_t* rseed) { bool HitComp(gHit a, gHit b) { - return ((a.chr