modified HET transition: equal prob for all clonal clusters

gavinha · Feb 9, 2015 · 977c48a · 977c48a
1 parent 2c76334
commit 977c48a
Show file tree

Hide file tree

Showing 3 changed files with 19 additions and 14 deletions.
diff --git a/R/utils.R b/R/utils.R
@@ -55,7 +55,7 @@ loadDefaultParameters <- function(copyNumber = 5, numberClonalClusters = 1,
         highStates <- c(1,16:length(rt))
         hetState <- c(5, 13, 25, 41)
     }
-    ZS[hetState] <- -1
+    ZS[hetState[1]] <- -1
     rn = rn + skew
     ind <- ct <= copyNumber
     rt <- rt[ind]
@@ -406,10 +406,10 @@ correctReadDepth <- function(tumWig, normWig, gcWig, mapWig,
     	## BUG WITH USE OF FINDOVERLAPS #######
     	## unlist indices are for each space ##
     	#######################################
-        keepInd <- unlist(as.list(findOverlaps(tumour_reads, targetIR, select = "first")))
-        keepInd <- !is.na(keepInd)
-        #hits <- findOverlaps(query = tumour_reads, subject = targetIR)
-        #keepInd <- queryHits(hits)      
+        #keepInd <- unlist(as.list(findOverlaps(tumour_reads, targetIR, select = "first")))
+        #keepInd <- !is.na(keepInd)
+        hits <- findOverlaps(query = tumour_reads, subject = targetIR)
+        keepInd <- unique(queryHits(hits))    
 
         # ind <- tumour_reads$value>10 &
         # normal_reads$value>10 tumThres <-

diff --git a/src/fwd_backC_clonalCN.c b/src/fwd_backC_clonalCN.c
@@ -224,7 +224,7 @@ void preparePositionSpecificMatrix(double * transSlice, unsigned int K, unsigned
         if (OUTLIERSTATE==1){
 					if (i==0){//garbage state 
 						z1 = 0;  unitI = -1;
-						iZS = -1;
+						iZS = -100;
 					}else{
 						z1 = ceil(((double)i)/numUnitStates); /* cluster number */
 						unitI = (int)(i-1)%(numUnitStates); /* unit state */
@@ -239,7 +239,7 @@ void preparePositionSpecificMatrix(double * transSlice, unsigned int K, unsigned
           	if (OUTLIERSTATE==1){
 							if (j==0){//garbage state 
 								z2 = 0;  unitJ = -1;
-								jZS = -1;
+								jZS = -100;
 							}else{
 								z2 = ceil(((double)j)/numUnitStates); /* cluster number */
 								unitJ = (int)(j-1)%(numUnitStates); /* unit state */
@@ -252,19 +252,23 @@ void preparePositionSpecificMatrix(double * transSlice, unsigned int K, unsigned
 						} 
             //printf("i=%d\tj=%d\tz1=%f\tz2=%f\tZS[%d]=%f\tZS[%d]=%f\n",i,j,z1,z2,unitI,iZS,unitJ,jZS);
 						//transitions to same state or same zygosity status
-            //if (i==j || (iZS==jZS)){
+
+            /** GENOTYPE TRANSITION **/
             if ((iZS==jZS)){
                 transSlice[i + j*K] = rhoG; 
             }else{
 	              transSlice[i + j*K] = (1.0-rhoG)/((double)K-1.0); 
             }
 
-	    			//same cluster
-            if(z1 == z2 || jZS==-1){  
-                transSlice[i + j*K] = transSlice[i + j*K] * rhoZ;
-            }else{  //different cluster (except for diploid HET)
-                transSlice[i + j*K] = transSlice[i + j*K] * (1.0-rhoZ);
-            }
+			/** CLONAL CLUSTER TRANSITION **/
+			//if (K > numUnitStates){ //only use if more than 1 cluster
+	    	//same clust or het (-1) 
+				if(z1 == z2 || jZS == -1){ // || CT[unitJ] >= 4){  
+					transSlice[i + j*K] = transSlice[i + j*K] * rhoZ;
+				}else{  //different cluster (except for diploid HET)
+					transSlice[i + j*K] = transSlice[i + j*K] * (1.0-rhoZ);
+				}
+			//}
         }
     }          
 

diff --git a/src/viterbiC_clonalCN.c b/src/viterbiC_clonalCN.c
@@ -30,6 +30,7 @@ void initializeTxnV(double * transSlice, unsigned int K);
 void logMatrixInPlace(double * A, unsigned int K);
 void outputMatrixV(double * A, unsigned int K);
 double distanceTransitionFunctionV(double, double, double);
+void preparePositionSpecificMatrix(double *, unsigned int, unsigned int, double *, double *, double, double, unsigned int, unsigned int);
 
 SEXP viterbiC_clonalCN(SEXP piGiZi, SEXP py, SEXP copyNumKey, SEXP zygosityKey, SEXP numClust, SEXP positions, SEXP zStrength, SEXP txnLen, SEXP useOutlier) {
   double * prior, * obslik, * transSlice;