diff --git a/PosteriorProb.h b/PosteriorProb.h deleted file mode 100644 index 7870335..0000000 --- a/PosteriorProb.h +++ /dev/null @@ -1,8 +0,0 @@ -#ifndef POSTERPROB_H -#define POSTERPROB_H -double NoPMDGivenAnc_b(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv,double Tol); -double NoPMDGivenAnc_nb(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv,double Tol); -double AncProb(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv, uint8_t seqError[], int model, double eps, double anc_mu, double anc_si, double mod_mu, double mod_si, int isrecal, int len_limit, int len_min,double Tol); -double PMDProb(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv, uint8_t seqError[], int model,double Tol); -bam_hdr_t* calc_pp_pmd_prob(char *refName,char *fname, char* ofname, int mapped_only,int se_only, int mapq, faidx_t *seq_ref, int len_limit, int len_min, int model, double eps, double lambda, double delta, double delta_s, double nv, double anc_mu, double anc_si, double mod_mu, double mod_si, int isrecal, kstring_t *str_cli,double **deamRateCT,double **deamRateGA,double Tol); -#endif diff --git a/likelihood.cpp b/likelihood.cpp index 695c2d7..68c474a 100644 --- a/likelihood.cpp +++ b/likelihood.cpp @@ -1,4 +1,4 @@ -#include "likelihood.h" + #include #include #include @@ -20,11 +20,11 @@ #include "profile.h" #include "bfgs.h" #include "htslib/bgzf.h" - +#include "likelihood.h" #include "misc.h" #include "recalibration.h" #include "ngsBriggs.h" -#include "PosteriorProb.h" +//#include "PosteriorProb.h" double PhredError[255]; double PhredErrorAThird[255]; @@ -1475,7 +1475,7 @@ void like_hess_master(const double *xs,double **y){ } } -double ErrorLik(char reffrag[], char frag[], int L, uint8_t seqError[]){ +double ErrorLik(char reffrag[], char frag[], int L, uint8_t seqError[],int l_check){ double l1 = 0; for(int i=0;0&&i multiplicity model + speed up -double PMDLik_b(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv, uint8_t seqError[],double Tol){ +double PMDLik_b(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv, uint8_t seqError[],double Tol,int l_check){ #if 0 for(int i=0;i mulplicative -double PMDLik_nb(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv, uint8_t seqError[],double Tol){ +double PMDLik_nb(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv, uint8_t seqError[],double Tol,int l_check){ double l_pmd=0; // Likelihood double p = 0; // Accumulated prob of (l,r) double exterm_s = 0; diff --git a/likelihood.h b/likelihood.h index 9c6ae55..ac72fd2 100644 --- a/likelihood.h +++ b/likelihood.h @@ -77,10 +77,10 @@ void like_grad_master(const double *xs,double *y,const void *); void like_hess_master(const double *xs,double **y); -double ErrorLik(char reffrag[], char frag[], int L, uint8_t seqError[]); +double ErrorLik(char reffrag[], char frag[], int L, uint8_t seqError[],int l_check); -double PMDLik_b(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv, uint8_t seqError[],double Tol); +double PMDLik_b(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv, uint8_t seqError[],double Tol,int l_check); -double PMDLik_nb(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv, uint8_t seqError[],double Tol); +double PMDLik_nb(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv, uint8_t seqError[],double Tol,int l_check); #endif diff --git a/ngsBriggs.cpp b/ngsBriggs.cpp index 487928d..a6cfa59 100644 --- a/ngsBriggs.cpp +++ b/ngsBriggs.cpp @@ -16,7 +16,7 @@ #include "misc.h" #include "recalibration.h" #include "likelihood.h" -#include "PosteriorProb.h" +#include "pp.h" #include "ngsBriggs_cli.h" #include "ngsBriggs.h" #include "bdamagereader.h" @@ -24,7 +24,7 @@ // definining all global variables used across multiple scripts -double l_check = 15; + tsk_struct *my_tsk_struct = NULL; int tsk_nthreads = -1; @@ -32,7 +32,7 @@ int tsk_nthreads = -1; // defining our main ngsBriggs function int main(int argc, char **argv){ - + double l_check = 15; double** mm5p, **mm3p; double Tol = 1.0E-8; // Tolerance double **deamRateCT; @@ -475,7 +475,7 @@ int main(int argc, char **argv){ } bam_hdr_t *hdr = sam_hdr_read(in); int ndim = 0; - double **mat = read_all_reads(in,hdr,seq_ref,len_limit,invec2[0],invec2[1],invec2[2],invec2[3],Tol,ndim,model); + double **mat = read_all_reads(in,hdr,seq_ref,len_limit,invec2[0],invec2[1],invec2[2],invec2[3],Tol,ndim,model,l_check); //tsk stop //double lbd2[4] = {30,1e-8,30,1e-8}; @@ -533,7 +533,7 @@ int main(int argc, char **argv){ kstr3->l = kstr3->m = 0; if (mypars->ihts!=NULL){ - bam_hdr_t* hdr = calc_pp_pmd_prob(refName,mypars->hts,mypars->ohts,mapped_only,se_only,mapq,seq_ref,len_limit,len_min,model,Contam_eps,invec2[0],invec2[1],invec2[2],invec2[3],distparam[0],distparam[1],distparam[2], distparam[3], mypars->dorecal,&str_cli,deamRateCT,deamRateGA,Tol); + bam_hdr_t* hdr = calc_pp_pmd_prob(refName,mypars->hts,mypars->ohts,mapped_only,se_only,mapq,seq_ref,len_limit,len_min,model,Contam_eps,invec2[0],invec2[1],invec2[2],invec2[3],distparam[0],distparam[1],distparam[2], distparam[3], mypars->dorecal,&str_cli,deamRateCT,deamRateGA,Tol,l_check); sam_hdr_destroy(hdr); } diff --git a/ngsBriggs.h b/ngsBriggs.h index f37f540..51a3730 100644 --- a/ngsBriggs.h +++ b/ngsBriggs.h @@ -1,10 +1,6 @@ #ifndef NGSBRIGGS_H #define NGSBRIGGS_H extern int MAXLENGTH; -extern double l_check; - extern char refToChar[256]; extern char com[5]; - #endif - diff --git a/PosteriorProb.cpp b/pp.cpp similarity index 96% rename from PosteriorProb.cpp rename to pp.cpp index 353b4d9..defa8ab 100644 --- a/PosteriorProb.cpp +++ b/pp.cpp @@ -12,14 +12,14 @@ #include #include "ngsBriggs.h" #include "likelihood.h" -#include "PosteriorProb.h" +#include "pp.h" #include "misc.h" #include "profile.h" typedef unsigned char uchar; -double NoPMDGivenAnc_b(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv,double Tol){ +double NoPMDGivenAnc_b(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv,double Tol,int l_check){ double np_pmd_anc=0; // probability of no pmd given ancient double p = 0; //Investigate each possible left and right 5' overhang pair with (l,r) <= l_check (15) @@ -402,7 +402,7 @@ double NoPMDGivenAnc_b(char reffrag[], char frag[], int L, double lambda, double } -double NoPMDGivenAnc_nb(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv,double Tol){ +double NoPMDGivenAnc_nb(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv,double Tol,int l_check){ double np_pmd_anc=0; // Likelihood double p = 0; // Accumulated prob of (l,r) //Investigate each possible left and right 5' overhang pair with (l,r) <= l_check (15) @@ -814,8 +814,8 @@ double NoPMDGivenAnc_nb(char reffrag[], char frag[], int L, double lambda, doubl } // The function below is for calculating the posterior prob of being ancient -double AncProb(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv, uint8_t seqError[], int model, double eps, double anc_mu, double anc_si, double mod_mu, double mod_si, int isrecal, int len_limit, int len_min,double Tol){ - double l_err = ErrorLik(reffrag, frag, L, seqError); // Modern Likelihood/Only Sequencing-error Likelihood +double AncProb(char reffrag[], char frag[], int L, double lambda, double delta, double delta_s, double nv, uint8_t seqError[], int model, double eps, double anc_mu, double anc_si, double mod_mu, double mod_si, int isrecal, int len_limit, int len_min,double Tol,int l_check){ + double l_err = ErrorLik(reffrag, frag, L, seqError,l_check); // Modern Likelihood/Only Sequencing-error Likelihood double prior_anc = 1-eps; double post_anc; double l_anc_l = 1; @@ -833,12 +833,12 @@ double AncProb(char reffrag[], char frag[], int L, double lambda, double delta, l_err_l = NormalINC(y_max2, y_min2, x_max2, x_min2); } if (model==0){ - double l_anc_b = PMDLik_b(reffrag, frag, L, lambda, delta, delta_s, nv, seqError,Tol); // Ancient Likelihood based on biotin model + double l_anc_b = PMDLik_b(reffrag, frag, L, lambda, delta, delta_s, nv, seqError,Tol,l_check); // Ancient Likelihood based on biotin model post_anc = prior_anc * l_anc_b * l_anc_l/(prior_anc * l_anc_b * l_anc_l+ (1-prior_anc) * l_err * l_err_l); //double test = exp(-10); //cout<core.l_qseq, lambda, delta, delta_s, nv, yourqual, model ,eps, anc_mu, anc_si, mod_mu, mod_si, 0, len_limit, len_min,Tol); + double PostAncProb1 = AncProb(yourrefe, yourread, b->core.l_qseq, lambda, delta, delta_s, nv, yourqual, model ,eps, anc_mu, anc_si, mod_mu, mod_si, 0, len_limit, len_min,Tol,l_check); double PostAncProb2 = 0; if (isrecal==1) { - PostAncProb2 = AncProb(yourrefe, yourread, b->core.l_qseq, lambda, delta, delta_s, nv, yourqual, model ,eps, anc_mu, anc_si, mod_mu, mod_si, 1, len_limit, len_min,Tol); + PostAncProb2 = AncProb(yourrefe, yourread, b->core.l_qseq, lambda, delta, delta_s, nv, yourqual, model ,eps, anc_mu, anc_si, mod_mu, mod_si, 1, len_limit, len_min,Tol,l_check); } - double PostPMDProb = PMDProb(yourrefe, yourread, b->core.l_qseq, lambda, delta, delta_s, nv, yourqual, model,Tol); + double PostPMDProb = PMDProb(yourrefe, yourread, b->core.l_qseq, lambda, delta, delta_s, nv, yourqual, model,Tol,l_check); //cout << "length 0" << " " << b->core.l_qseq << "\n"; -#if 0 - // Read Name - std::cout << bam_get_qname(b) << "\t"; - cout << bam_aux_get(b,"FL") << "\t"; - for (int cycle=0;cyclecore.l_qseq;cycle++) - std::cout<core.l_qseq << "\n"; - for (int cycle=0;cyclecore.l_qseq;cycle++) - std::cout<core.l_qseq;cycle++) - std::cout<<(char)(bam_get_qual(b)[cycle]+33); - if (isrecal==1) - std::cout<<"\t"<<"AO:f:"<core.l_qseq, yourqual); - double l_anc = PMDLik_b(yourrefe, yourread, b->core.l_qseq, lambda, delta, delta_s, nv, yourqual,Tol); + double l_err = ErrorLik(yourrefe, yourread, b->core.l_qseq, yourqual,l_check); + double l_anc = PMDLik_b(yourrefe, yourread, b->core.l_qseq, lambda, delta, delta_s, nv, yourqual,Tol,l_check); // fprintf(stderr,"pmdlik_b: %f nb: %f model: %d\n",PMDLik_b(yourrefe, yourread, b->core.l_qseq, lambda, delta, delta_s, nv, yourqual,Tol),PMDLik_nb(yourrefe, yourread, b->core.l_qseq, lambda, delta, delta_s, nv, yourqual,Tol),model); if(model==1) - l_anc = 0.5*l_anc + 0.5*PMDLik_nb(yourrefe, yourread, b->core.l_qseq, lambda, delta, delta_s, nv, yourqual,Tol); + l_anc = 0.5*l_anc + 0.5*PMDLik_nb(yourrefe, yourread, b->core.l_qseq, lambda, delta, delta_s, nv, yourqual,Tol,l_check); // fprintf(stderr,"lanc: %f\n",l_anc);exit(0); vec.push_back(l_err); vec.push_back(l_anc); diff --git a/read_all_reads.h b/read_all_reads.h index a42f951..b253e30 100644 --- a/read_all_reads.h +++ b/read_all_reads.h @@ -8,5 +8,5 @@ #include #include -double **read_all_reads(samFile *in,sam_hdr_t *hdr,faidx_t *seq_ref,int len_limit,double lambda,double delta,double delta_s,double nv,double Tol,int &ndim,int model); +double **read_all_reads(samFile *in,sam_hdr_t *hdr,faidx_t *seq_ref,int len_limit,double lambda,double delta,double delta_s,double nv,double Tol,int &ndim,int model,int l_check); #endif