Lines Matching refs:args
110 static args_t args; variable
118 static void init_hmm_trio(args_t *args);
119 static void init_hmm_unrelated(args_t *args);
150 memset(&args,0,sizeof(args_t)); in init()
151 args.prev_rid = -1; in init()
152 args.hdr = in; in init()
153 args.pij = 2e-8; in init()
154 args.pgt_err = 1e-9; in init()
168 case 'p': args.prefix = optarg; break; in init()
179 if ( !args.prefix ) error("Expected the -p option\n"); in init()
181 int ret = bcf_hdr_set_samples(args.hdr, trio_samples ? trio_samples : unrelated_samples, 0); in init()
190 args.imother = bcf_hdr_id2int(args.hdr, BCF_DT_SAMPLE, list[0]); in init()
191 args.ifather = bcf_hdr_id2int(args.hdr, BCF_DT_SAMPLE, list[1]); in init()
192 args.ichild = bcf_hdr_id2int(args.hdr, BCF_DT_SAMPLE, list[2]); in init()
195 args.set_observed_prob = set_observed_prob_trio; in init()
196 args.mode = C_TRIO; in init()
197 init_hmm_trio(&args); in init()
204 args.isample = bcf_hdr_id2int(args.hdr, BCF_DT_SAMPLE, list[0]); in init()
205 args.jsample = bcf_hdr_id2int(args.hdr, BCF_DT_SAMPLE, list[1]); in init()
208 args.set_observed_prob = set_observed_prob_unrelated; in init()
209 args.mode = C_UNRL; in init()
210 init_hmm_unrelated(&args); in init()
215 static void init_hmm_trio(args_t *args) in init_hmm_trio() argument
218 args->nstates = 8; in init_hmm_trio()
219 args->tprob = (double*) malloc(sizeof(double)*args->nstates*args->nstates); in init_hmm_trio()
221 for (i=0; i<args->nstates; i++) in init_hmm_trio()
222 for (j=0; j<args->nstates; j++) hap_switch[i][j] = 0; in init_hmm_trio()
250 for (i=0; i<args->nstates; i++) in init_hmm_trio()
252 for (j=0; j<args->nstates; j++) in init_hmm_trio()
254 if ( !hap_switch[i][j] ) MAT(args->tprob,args->nstates,i,j) = 0; in init_hmm_trio()
257 MAT(args->tprob,args->nstates,i,j) = 1; in init_hmm_trio()
258 if ( hap_switch[i][j] & SW_MOTHER ) MAT(args->tprob,args->nstates,i,j) *= args->pij; in init_hmm_trio()
259 if ( hap_switch[i][j] & SW_FATHER ) MAT(args->tprob,args->nstates,i,j) *= args->pij; in init_hmm_trio()
263 for (i=0; i<args->nstates; i++) in init_hmm_trio()
266 for (j=0; j<args->nstates; j++) in init_hmm_trio()
268 if ( i!=j ) sum += MAT(args->tprob,args->nstates,i,j); in init_hmm_trio()
270 MAT(args->tprob,args->nstates,i,i) = 1 - sum; in init_hmm_trio()
274 for (i=0; i<args->nstates; i++) in init_hmm_trio()
276 for (j=0; j<args->nstates; j++) in init_hmm_trio()
280 for (i=0; i<args->nstates; i++) in init_hmm_trio()
282 for (j=0; j<args->nstates; j++) in init_hmm_trio()
283 fprintf(stderr,"\t%e",MAT(args->tprob,args->nstates,j,i)); in init_hmm_trio()
288 args->hmm = hmm_init(args->nstates, args->tprob, 10000); in init_hmm_trio()
290 static void init_hmm_unrelated(args_t *args) in init_hmm_unrelated() argument
293 args->nstates = 7; in init_hmm_unrelated()
294 args->tprob = (double*) malloc(sizeof(double)*args->nstates*args->nstates); in init_hmm_unrelated()
296 for (i=0; i<args->nstates; i++) in init_hmm_unrelated()
298 for (j=0; j<args->nstates; j++) in init_hmm_unrelated()
299 MAT(args->tprob,args->nstates,i,j) = args->pij; in init_hmm_unrelated()
301 MAT(args->tprob,args->nstates,UNRL_0101,UNRL_xxxx) = args->pij*args->pij; in init_hmm_unrelated()
302 MAT(args->tprob,args->nstates,UNRL_0110,UNRL_xxxx) = args->pij*args->pij; in init_hmm_unrelated()
303 MAT(args->tprob,args->nstates,UNRL_x0x0,UNRL_0x0x) = args->pij*args->pij; in init_hmm_unrelated()
304 MAT(args->tprob,args->nstates,UNRL_0110,UNRL_0x0x) = args->pij*args->pij; in init_hmm_unrelated()
305 MAT(args->tprob,args->nstates,UNRL_x00x,UNRL_0xx0) = args->pij*args->pij; in init_hmm_unrelated()
306 MAT(args->tprob,args->nstates,UNRL_0101,UNRL_0xx0) = args->pij*args->pij; in init_hmm_unrelated()
307 MAT(args->tprob,args->nstates,UNRL_0101,UNRL_x00x) = args->pij*args->pij; in init_hmm_unrelated()
308 MAT(args->tprob,args->nstates,UNRL_0110,UNRL_x0x0) = args->pij*args->pij; in init_hmm_unrelated()
309 MAT(args->tprob,args->nstates,UNRL_0110,UNRL_0101) = args->pij*args->pij; in init_hmm_unrelated()
311 for (i=0; i<args->nstates; i++) in init_hmm_unrelated()
313 for (j=i+1; j<args->nstates; j++) in init_hmm_unrelated()
314 MAT(args->tprob,args->nstates,i,j) = MAT(args->tprob,args->nstates,j,i); in init_hmm_unrelated()
316 for (i=0; i<args->nstates; i++) in init_hmm_unrelated()
319 for (j=0; j<args->nstates; j++) in init_hmm_unrelated()
320 if ( i!=j ) sum += MAT(args->tprob,args->nstates,i,j); in init_hmm_unrelated()
321 MAT(args->tprob,args->nstates,i,i) = 1 - sum; in init_hmm_unrelated()
325 for (i=0; i<args->nstates; i++) in init_hmm_unrelated()
327 for (j=0; j<args->nstates; j++) in init_hmm_unrelated()
328 fprintf(stderr,"\t%e",MAT(args->tprob,args->nstates,j,i)); in init_hmm_unrelated()
333 args->hmm = hmm_init(args->nstates, args->tprob, 10000); in init_hmm_unrelated()
337 return a==b ? 1-args.pgt_err : args.pgt_err; in prob_shared()
349 int ngt = bcf_get_genotypes(args.hdr, rec, &args.gt_arr, &args.ngt_arr); in set_observed_prob_unrelated()
354 a = args.gt_arr[2*args.isample]; in set_observed_prob_unrelated()
355 b = args.gt_arr[2*args.isample+1]; in set_observed_prob_unrelated()
356 c = args.gt_arr[2*args.jsample]; in set_observed_prob_unrelated()
357 d = args.gt_arr[2*args.jsample+1]; in set_observed_prob_unrelated()
367 int m = args.msites; in set_observed_prob_unrelated()
368 args.nsites++; in set_observed_prob_unrelated()
369 hts_expand(uint32_t,args.nsites,args.msites,args.sites); in set_observed_prob_unrelated()
370 …if ( m!=args.msites ) args.eprob = (double*) realloc(args.eprob, sizeof(double)*args.msites*args.n… in set_observed_prob_unrelated()
372 args.sites[args.nsites-1] = rec->pos; in set_observed_prob_unrelated()
373 double *prob = args.eprob + args.nstates*(args.nsites-1); in set_observed_prob_unrelated()
397 int ngt = bcf_get_genotypes(args.hdr, rec, &args.gt_arr, &args.ngt_arr); in set_observed_prob_trio()
402 a = args.gt_arr[2*args.imother]; in set_observed_prob_trio()
403 b = args.gt_arr[2*args.imother+1]; in set_observed_prob_trio()
404 c = args.gt_arr[2*args.ifather]; in set_observed_prob_trio()
405 d = args.gt_arr[2*args.ifather+1]; in set_observed_prob_trio()
406 e = args.gt_arr[2*args.ichild]; in set_observed_prob_trio()
407 f = args.gt_arr[2*args.ichild+1]; in set_observed_prob_trio()
426 if ( a!=b ) args.nhet_mother++; in set_observed_prob_trio()
427 if ( c!=d ) args.nhet_father++; in set_observed_prob_trio()
429 int m = args.msites; in set_observed_prob_trio()
430 args.nsites++; in set_observed_prob_trio()
431 hts_expand(uint32_t,args.nsites,args.msites,args.sites); in set_observed_prob_trio()
432 …if ( m!=args.msites ) args.eprob = (double*) realloc(args.eprob, sizeof(double)*args.msites*args.n… in set_observed_prob_trio()
434 args.sites[args.nsites-1] = rec->pos; in set_observed_prob_trio()
435 double *prob = args.eprob + args.nstates*(args.nsites-1); in set_observed_prob_trio()
446 void flush_viterbi(args_t *args) in flush_viterbi() argument
449 if ( args->mode==C_UNRL ) in flush_viterbi()
451 s1 = bcf_hdr_int2id(args->hdr,BCF_DT_SAMPLE,args->isample); in flush_viterbi()
452 s2 = bcf_hdr_int2id(args->hdr,BCF_DT_SAMPLE,args->jsample); in flush_viterbi()
454 else if ( args->mode==C_TRIO ) in flush_viterbi()
456 s1 = bcf_hdr_int2id(args->hdr,BCF_DT_SAMPLE,args->imother); in flush_viterbi()
457 s3 = bcf_hdr_int2id(args->hdr,BCF_DT_SAMPLE,args->ifather); in flush_viterbi()
458 s2 = bcf_hdr_int2id(args->hdr,BCF_DT_SAMPLE,args->ichild); in flush_viterbi()
462 if ( !args->fp ) in flush_viterbi()
465 kputs(args->prefix, &str); in flush_viterbi()
467 args->fp = fopen(str.s,"w"); in flush_viterbi()
468 if ( !args->fp ) error("%s: %s\n", str.s,strerror(errno)); in flush_viterbi()
470 …fprintf(args->fp,"# SG, shared segment\t[2]Chromosome\t[3]Start\t[4]End\t[5]%s:1\t[6]%s:2\n",s2,s2… in flush_viterbi()
471 …fprintf(args->fp,"# SW, number of switches\t[3]Sample\t[4]Chromosome\t[5]nHets\t[5]nSwitches\t[6]s… in flush_viterbi()
474 hmm_run_viterbi(args->hmm,args->nsites,args->eprob,args->sites); in flush_viterbi()
475 uint8_t *vpath = hmm_get_viterbi_path(args->hmm); in flush_viterbi()
476 int i, iprev = -1, prev_state = -1, nstates = hmm_get_nstates(args->hmm); in flush_viterbi()
478 for (i=0; i<args->nsites; i++) in flush_viterbi()
481 if ( state!=prev_state || i+1==args->nsites ) in flush_viterbi()
483 uint32_t start = iprev>=0 ? args->sites[iprev]+1 : 1, end = i>0 ? args->sites[i-1] : 1; in flush_viterbi()
484 const char *chr = bcf_hdr_id2name(args->hdr,args->prev_rid); in flush_viterbi()
485 if ( args->mode==C_UNRL ) in flush_viterbi()
490 fprintf(args->fp,"SG\t%s\t%d\t%d\t%s:1\t-\n", chr,start,end,s1); break; in flush_viterbi()
492 fprintf(args->fp,"SG\t%s\t%d\t%d\t-\t%s:1\n", chr,start,end,s1); break; in flush_viterbi()
494 fprintf(args->fp,"SG\t%s\t%d\t%d\t%s:2\t-\n", chr,start,end,s1); break; in flush_viterbi()
496 fprintf(args->fp,"SG\t%s\t%d\t%d\t-\t%s:2\n", chr,start,end,s1); break; in flush_viterbi()
498 … fprintf(args->fp,"SG\t%s\t%d\t%d\t%s:1\t%s:2\n", chr,start,end,s1,s1); break; in flush_viterbi()
500 … fprintf(args->fp,"SG\t%s\t%d\t%d\t%s:2\t%s:1\n", chr,start,end,s1,s1); break; in flush_viterbi()
503 else if ( args->mode==C_TRIO ) in flush_viterbi()
508 … fprintf(args->fp,"SG\t%s\t%d\t%d\t%s:1\t%s:1\n", chr,start,end,s1,s3); break; in flush_viterbi()
510 … fprintf(args->fp,"SG\t%s\t%d\t%d\t%s:1\t%s:2\n", chr,start,end,s1,s3); break; in flush_viterbi()
512 … fprintf(args->fp,"SG\t%s\t%d\t%d\t%s:2\t%s:1\n", chr,start,end,s1,s3); break; in flush_viterbi()
514 … fprintf(args->fp,"SG\t%s\t%d\t%d\t%s:2\t%s:2\n", chr,start,end,s1,s3); break; in flush_viterbi()
516 … fprintf(args->fp,"SG\t%s\t%d\t%d\t%s:1\t%s:1\n", chr,start,end,s3,s1); break; in flush_viterbi()
518 … fprintf(args->fp,"SG\t%s\t%d\t%d\t%s:2\t%s:1\n", chr,start,end,s3,s1); break; in flush_viterbi()
520 … fprintf(args->fp,"SG\t%s\t%d\t%d\t%s:1\t%s:2\n", chr,start,end,s3,s1); break; in flush_viterbi()
522 … fprintf(args->fp,"SG\t%s\t%d\t%d\t%s:2\t%s:2\n", chr,start,end,s3,s1); break; in flush_viterbi()
531 float mrate = args->nhet_mother>1 ? (float)nswitch_mother/(args->nhet_mother-1) : 0; in flush_viterbi()
532 float frate = args->nhet_father>1 ? (float)nswitch_father/(args->nhet_father-1) : 0; in flush_viterbi()
533 …fprintf(args->fp,"SW\t%s\t%s\t%d\t%d\t%f\n", s1,bcf_hdr_id2name(args->hdr,args->prev_rid),args->nh… in flush_viterbi()
534 …fprintf(args->fp,"SW\t%s\t%s\t%d\t%d\t%f\n", s3,bcf_hdr_id2name(args->hdr,args->prev_rid),args->nh… in flush_viterbi()
535 args->nsites = 0; in flush_viterbi()
536 args->nhet_father = args->nhet_mother = 0; in flush_viterbi()
541 if ( args.prev_rid==-1 ) args.prev_rid = rec->rid; in process()
542 if ( args.prev_rid!=rec->rid ) flush_viterbi(&args); in process()
543 args.prev_rid = rec->rid; in process()
544 args.set_observed_prob(rec); in process()
550 flush_viterbi(&args); in destroy()
551 fclose(args.fp); in destroy()
553 free(args.gt_arr); in destroy()
554 free(args.tprob); in destroy()
555 free(args.sites); in destroy()
556 free(args.eprob); in destroy()
557 hmm_destroy(args.hmm); in destroy()