init subclone fraction 0.5

bioinform · Jun 4, 2020 · 3a7a185 · 3a7a185
1 parent da638a6
commit 3a7a185
Show file tree

Hide file tree

Showing 3 changed files with 11 additions and 9 deletions.
diff --git a/R/output.R b/R/output.R
@@ -128,6 +128,7 @@ outputParametersToFile <- function(hmmResults, file){
       write.table(paste0("Gaussian variance: ", paste0(signif(x$vars[,s], digits = 2), collapse = ", ")), file = fc, col.names = FALSE, row.names = FALSE, quote = FALSE, sep = "\t")
     }
     write.table(paste0("GC-Map correction MAD:\t", format(mad(diff(2^as.numeric(hmmResults$cna[[s]][,"logR"])), na.rm=T), digits=4)), file = fc, col.names = FALSE, row.names = FALSE, quote = FALSE, sep = "\t")
+    write.table(paste0("Log-likelihood:\t", format(as.numeric(hmmResults.cor$results$loglik[i, "loglik"]), digits=4)), file = fc, col.names = FALSE, row.names = FALSE, quote = FALSE, sep = "\t")
     write.table("\n", file = fc, col.names = FALSE, row.names = FALSE, quote = FALSE, sep = "\t")
   }
   write.table(x$loglik, file = fc, col.names = TRUE, row.names = FALSE, quote = FALSE, sep = "\t")

diff --git a/R/segmentation.R b/R/segmentation.R
@@ -140,7 +140,7 @@ getDefaultParameters <- function(x, maxCN = 5, ct.sc = c(1,3), n_0 = 0.5, ploidy
 		ct.sc.status = c(rep(FALSE, length(ct)), rep(TRUE, length(ct.sc))),
 		phi_0 = ploidy_0, alphaPhi = 4, betaPhi = 0.75,
 		n_0 = n_0, alphaN = 1, betaN = 1,
-		sp_0 = 0.01, alphaSp = 2, betaSp = 2,
+		sp_0 = 0.5, alphaSp = 2, betaSp = 2,
 		lambda = as.matrix(rep(100, length(ct)+length(ct.sc)), ncol=1),
 		nu = 2.1,
 		kappa = rep(75, length(ct)), 

diff --git a/scripts/runIchorCNA.R b/scripts/runIchorCNA.R
@@ -414,15 +414,16 @@ save.image(outImage)
 ### SORT SOLUTIONS BY DECREASING LIKELIHOOD ###
 fracAltSub[is.nan(fracAltSub$V1), ] <- 0
 if (estimateScPrevalence){ ## sort but excluding solutions with too large % subclonal 
-    if (numSamples > 1){
-        fracInd <- which(rowSums(fracAltSub <= maxFracCNASubclone) == S &
-                        rowSums(fracGenomeSub <= maxFracGenomeSubclone) == S)
-    }else{
-        fracInd <- which(fracAltSub <= maxFracCNASubclone &
-                        fracGenomeSub <= maxFracGenomeSubclone)
-    }
+  if (numSamples > 1){
+      fracInd <- which(rowSums(fracAltSub <= maxFracCNASubclone) == S &
+                      rowSums(fracGenomeSub <= maxFracGenomeSubclone) == S)
+  }else{
+      fracInd <- which(fracAltSub <= maxFracCNASubclone &
+                      fracGenomeSub <= maxFracGenomeSubclone)
+  }
 	if (length(fracInd) > 0){ ## if there is a solution satisfying % subclonal
-		ind <- fracInd[order(loglik[fracInd, "loglik"], decreasing=TRUE)]
+    # sort by highest loglik and then by lowest fracCNAsubclonal (if ties)
+		ind <- fracInd[order(-loglik[fracInd, "loglik"], loglik[fracInd, "Frac_CNA_subclonal"])]
 	}else{ # otherwise just take largest likelihood
 		ind <- order(as.numeric(loglik[, "loglik"]), decreasing=TRUE) 
 	}