mod README

README.md

@@ -41,16 +41,16 @@ jump-diffusion propagator. It can be used for scale-free smoothing,
 fuzzy data segmentation and data filtering. 
 
 ![cna gof](/images/cna.gof.png "CNA goodness of fit")
-![cna gof](/images/cna.post.png "CNA posterior")
 ![baf gof](/images/baf.gof.png "BAF goodness of fit")
-![cna gof](/images/baf.post.png "BAF posterior")
+![cna post](/images/cna.post.png "CNA posterior")
+![baf post](/images/baf.post.png "BAF posterior")
+![snv gof](/images/snv.gof.png "SNV goodness of fit")
 
 Visualization of the cloneHD output for the simulated data set. From
 top to bottom: (i) the bias corrected read depth data and the cloneHD
-posterior mean emission rate (ii) the total copy number posterior
-distribution for subclone 1 with f1=0.52 and subclone 2 with f2=0.07
-(iii) the BAF and (iv) the minor allele posterior. (Plots created with Wolfram
-[Mathematica](http://www.wolfram.com/mathematica/).)
+posterior mean emission rate (ii) the BAF (iii-iv) the total and minor copy number posterior
+distribution for subclone 1 with f1=0.52 and subclone 2 with f2=0.07 and (v) the SNV goodness of fit, where each SNV is assigned to a genotype according to the its posterior, and the observed SNV freequency is scaled by one half of the local mean total copynumber. 
+(All plots are created with Wolfram [Mathematica](http://www.wolfram.com/mathematica/).)
 
 # filterHD command line arguments
 

run-example.sh

@@ -6,16 +6,17 @@ 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/"
+#data="/Users/af7/Projects/cloneHD/Test/data/"
+#results="/Users/af7/Projects/cloneHD/Test/results/"
 filterHD="./build/filterHD"
 cloneHD="./build/cloneHD"
 
 normalCNA="${data}/normal.cna.txt"
 tumorCNA="${data}/tumor.cna.txt"
 tumorBAF="${data}/tumor.baf.txt"
+tumorSNV="${data}/tumor.snv.txt"
 bias="${results}/normal.cna.posterior-1.txt"
 tumorCNAjumps="${results}/tumor.cna.bias.jumps.txt"
 tumorBAFjumps="${results}/tumor.baf.jumps.txt"
@@ -59,9 +60,14 @@ then
     echo "*** cloneHD ***"
     echo
     cmd="$cloneHD --cna $tumorCNA --baf $tumorBAF --pre ${results}/tumor --bias $bias --seed 123 --trials 2\
- --nmax 3 --force --max-tcn 5 --cna-jumps $tumorCNAjumps --baf-jumps $tumorBAFjumps --min-jump 0.01 --restarts 10"    
+ --nmax 3 --force --max-tcn 4 --cna-jumps $tumorCNAjumps --baf-jumps $tumorBAFjumps --min-jump 0.01 --restarts 10"    
     echo $cmd
     $cmd
     echo
     cat ${results}/tumor.summary.txt
+    cmd="$cloneHD --snv $tumorSNV --pre ${results}/tumor --seed 123 --trials 2\
+ --nmax 3 --force --max-tcn 4 --restarts 10 --mean-tcn ${results}/tumor.mean-tcn.txt --avail-cn ${results}/tumor.avail-cn.txt"    
+    echo $cmd
+    $cmd
+    echo
 fi

src/filterHD.cpp

@@ -179,7 +179,7 @@ int main (int argc, const char * argv[]){
     FILE * total_fp = fopen(buff,"w");
     fprintf(total_fp, "#sample site mean std-dev jump-prob");
     fprintf(total_fp, " posterior %.5e %.5e\n", myJD.myEmit->xmin, myJD.myEmit->xmax);
-    double gof=0, xobs=0,gofNorm=0;
+    double gof=0, xobs=0, gofNorm=0;
     for (int s=0; s < myJD.nSamples; s++){//get posterior distribution with the ML parameters
       myJD.get_posterior(s);
       double mstd = 0.0;
@@ -202,6 +202,21 @@ int main (int argc, const char * argv[]){
 	if (opts.jumps==1){
 	  jumps[s][l] *= exp(myJD.pnojump[s][l]);
 	}
+	//goodness of fit...
+	if (mask==NULL || mask[s][l] == 1){
+	  xobs = double(myEmit.reads[t][s][l]) / double(myEmit.depths[t][s][l]);
+	  if (opts.reflect) xobs = min(xobs,1.0-xobs);	 
+	  double x=0,g=0,dg=0;
+	  double b = (myEmit.bias == NULL) ? 1.0 : myEmit.bias[s][l];
+	  for (int i=0; i<=uidx; i++){
+	    x = myEmit.xgrid[i] * b;
+	    dg = fabs(x - xobs) * post->data[i];
+	    if (i==0||i==uidx) dg *= 0.5;
+	    g += dg;
+	  }
+	  gof += g * myEmit.dx;
+	  gofNorm += 1.0;
+	}
       }
       mstd /= double(myJD.nSites[s]);
       //filter out data points which are not compatible with the emission model
@@ -232,21 +247,6 @@ int main (int argc, const char * argv[]){
 	  }
 	}
 	fprintf(total_fp,"\n");
-	//goodness of fit...
-	if (mask==NULL || mask[s][l] == 1){
-	  xobs = double(myEmit.reads[t][s][l]) / double(myEmit.depths[t][s][l]);
-	  if (opts.reflect) xobs = min(xobs,1.0-xobs);	 
-	  double x=0,g=0,dg=0;
-	  double b = (myEmit.bias == NULL) ? 1.0 : myEmit.bias[s][l];
-	  for (int i=0; i<=uidx; i++){
-	    x = myEmit.xgrid[i] * b;
-	    dg = fabs(x - xobs) * post->data[i];
-	    if (i==0||i==uidx) dg *= 0.5;
-	    g += dg;
-	  }
-	  gof += g * myEmit.dx;
-	  gofNorm += 1.0;
-	}
       }
       gsl_matrix_free(myJD.gamma[s]);
       myJD.gamma[s] = NULL;