diff --git a/run-example.sh b/run-example.sh index 57f0577..0e176e5 100644 --- a/run-example.sh +++ b/run-example.sh @@ -1,64 +1,65 @@ # RUN filterHD & cloneHD FOR A SIMULATED EXAMPLE DATA SET # fix the number of threads -export OMP_NUM_THREADS=1; +export OMP_NUM_THREADS=4; + +part=$1 # input data -#data="./test/data/" -#results="./test/results/" -data="/Users/af7/Projects/cloneHD/Test/data/" -results="/Users/af7/Projects/cloneHD/Test/results/" +data="./test/data/" +results="./test/results/" filterHD="./build/filterHD" cloneHD="./build/cloneHD" -### filterHD ### -echo "*** filterHD ***" -echo +normalCNA="${data}/normal.cna.txt" +tumorCNA="${data}/tumor.cna.txt" +tumorBAF="${data}/tumor.baf.txt" +bias="${results}/normal.cna.posterior.1.txt" +tumorCNAjumps="${results}/tumor.cna.bias.jumps.txt" +tumorBAFjumps="${results}/tumor.baf.jumps.txt" +### filterHD ### +if [ -z $part ] || [ $part -eq 1 ] +then + echo "*** filterHD ***" + echo #emission modes: # 1: Binomial # 2: Beta-Binomial # 3: Poisson # 4: Negative Binomial - +# # normal read depth -normalCNA="${data}/normal.cna.txt" -cmd="$filterHD --data $normalCNA --mode 3 --pre ${results}/normal.cna --rnd 0" -echo $cmd -$cmd -echo - -# tumour read depth without bias -tumorCNA="${data}/tumor.cna.txt" -cmd="$filterHD --data $tumorCNA --mode 3 --pre ${results}/tumor.cna --rnd 0" -echo $cmd -$cmd -echo - -# tumour read depth with bias from normal -bias="${results}/normal.cna.posterior.1.txt" -cmd="$filterHD --data $tumorCNA --mode 3 --pre ${results}/tumor.cna.bias --bias $bias --sigma 0 --rnd 0 --jumps 1" -echo $cmd -$cmd -echo -tumorCNAjumps="${results}/tumor.cna.bias.jumps.txt" - -# tumour BAF -tumorBAF="${data}/tumor.baf.txt" -cmd="$filterHD --data $tumorBAF --mode 1 --pre ${results}/tumor.baf --sigma 0 --jumps 1 --reflect 1 --dist 1" -echo $cmd -$cmd -echo -tumorBAFjumps="${results}/tumor.baf.jumps.txt" - + cmd="$filterHD --data $normalCNA --mode 3 --pre ${results}/normal.cna --rnd 0" + echo $cmd + $cmd + echo +# tumor read depth without bias + cmd="$filterHD --data $tumorCNA --mode 3 --pre ${results}/tumor.cna --rnd 0" + echo $cmd + $cmd + echo +# tumor read depth with bias from normal + cmd="$filterHD --data $tumorCNA --mode 3 --pre ${results}/tumor.cna.bias --bias $bias --sigma 0 --rnd 0 --jumps 1" + echo $cmd + $cmd + echo +# tumor BAF + cmd="$filterHD --data $tumorBAF --mode 1 --pre ${results}/tumor.baf --sigma 0 --jumps 1 --reflect 1 --dist 1 --rnd 0" + echo $cmd + $cmd + echo +fi +if [ -z $part ] || [ $part -eq 2 ] +then ### cloneHD ### -echo "*** cloneHD ***" -echo - -cmd="$cloneHD --cna $tumorCNA --baf $tumorBAF --pre ${results}/tumor --bias $bias --seed 123 --trials 1 --nmax 3\ - --force --max-tcn 4 --cna-jumps $tumorCNAjumps --min-jump 0.01 --restarts 20" -echo $cmd -$cmd -echo -cat ${results}/tumor.clonal.txt \ No newline at end of file + echo "*** cloneHD ***" + echo + cmd="$cloneHD --cna $tumorCNA --baf $tumorBAF --pre ${results}/tumor --bias $bias --seed 123 --trials 1\ + --nmax 3 --force --max-tcn 4 --cna-jumps $tumorCNAjumps --baf-jumps $tumorBAFjumps --min-jump 0.01 --restarts 20" + echo $cmd + $cmd + echo + cat ${results}/tumor.clonal.txt +fi \ No newline at end of file diff --git a/src/clone-prior.cpp b/src/clone-prior.cpp index 5fae22d..31651f7 100644 --- a/src/clone-prior.cpp +++ b/src/clone-prior.cpp @@ -144,7 +144,7 @@ void Clone::set_snv_prior_map(){//either via BAF or else via CNA snv_prior_from_cna_map = gsl_matrix_alloc( maxtcn+1, maxtcn+1); } gsl_matrix_set_zero( snv_prior_from_cna_map); - p = (snvEmit->connect) ? 1.0 : snv_pen;// penalty for high genotypes + double p = (snvEmit->connect) ? 1.0 : snv_pen;// penalty for high genotypes for (int cn=0; cn <= maxtcn; cn++){ for (int i=0; i <= cn; i++){ gsl_matrix_set( snv_prior_from_cna_map, i, cn, pow(p,i));