Include log-likelihood-ratio filter and normalization

Snake.config.json

@@ -19,10 +19,9 @@
     "R_reference"   : "BSgenome.Hsapiens.UCSC.hg38",
 
     "bp_density"    : {
-        "few"       : 0.1,
-        "medium"    : 0.25,
-        "more"      : 0.33,
-        "many"      : 0.5
+        "few"       : 0.2,
+        "medium"    : 0.4,
+        "many"      : 0.6
     },
 
     "paired_end"    : true,

Snakefile

@@ -18,7 +18,7 @@ import os.path
 # * calculate a segmentation into potential SVs using Mosaicatcher
 
 
-METHODS = ["simpleCalls", "biAllelic"]
+METHODS = ["simpleCalls_llr1", "simpleCalls_llr4", "biAllelic_llr1", "biAllelic_llr4"]
 
 rule all:
     input:
@@ -27,7 +27,7 @@ rule all:
                sample = SAMPLE,
                chrom = config["chromosomes"],
                window = [50000, 100000],
-               bpdens = ["few","medium","more","many"],
+               bpdens = ["few","medium","many"],
                method = METHODS)
 
 
@@ -349,9 +349,9 @@ rule mosaiClassifier_make_call:
     input:
         probs = "sv_probabilities/{sample}/{windows}.{bpdens}/probabilities.Rdata"
     output:
-        "sv_calls/{sample}/{windows}.{bpdens}/simpleCalls.txt"
+        "sv_calls/{sample}/{windows}.{bpdens}/simpleCalls_llr{llr}.txt"
     log:
-        "log/{sample}/mosaiClassifier_make_call.{windows}.{bpdens}.txt"
+        "log/{sample}/mosaiClassifier_make_call.{windows}.{bpdens}.{llr}.txt"
     script:
         "utils/mosaiClassifier_call.snakemake.R"
 
@@ -379,9 +379,9 @@ rule mosaiClassifier_make_call_biallelic:
     input:
         probs = "sv_probabilities/{sample}/{windows}.{bpdens}/probabilities.Rdata"
     output:
-        "sv_calls/{sample}/{windows}.{bpdens}/biAllelic.txt"
+        "sv_calls/{sample}/{windows}.{bpdens}/biAllelic_llr{llr}.txt"
     log:
-        "log/{sample}/mosaiClassifier_make_call_biallelic.{windows}.{bpdens}.txt"
+        "log/{sample}/mosaiClassifier_make_call_biallelic.{windows}.{bpdens}.{llr}.txt"
     script:
         "utils/mosaiClassifier_call_biallelic.snakemake.R"
 

utils/mosaiClassifier.snakemake.R

@@ -12,7 +12,8 @@ segs   = fread(snakemake@input[["bp"]])
 
 
 # is there a normalization file given?
-if ("norm" %in% snakemake@input) {
+if ("norm" %in% names(snakemake@input)) {
+  message("[MosaiClassifier] Read normalization from ", snakemake@input[["norm"]])
   normalization = fread(snakemake@input[["norm"]])
 } else {
   normalization = NULL

utils/mosaiClassifier/mosaiClassifier.R

@@ -127,6 +127,7 @@ mosaiClassifierPrepare <- function(counts, info, strand, segs, normVector = NULL
 
   # Add normalization factors to the expected counts ("scalar")
   if (!is.null(normVector)) {
+    message("[MosaiClassifier] Normalize coverage expectation")
     addNormalizationScalar(probs, counts, normVector)
   } else {
     probs[, scalar := 1.0]

utils/mosaiClassifier_call.snakemake.R

@@ -3,9 +3,9 @@ library(data.table)
 source("utils/mosaiClassifier/makeSVcalls.R")
 
 probs = readRDS(snakemake@input[["probs"]])
-llr   = as.numeric(snakemake@params[["llr"]])
+llr   = as.numeric(snakemake@wildcards[["llr"]])
 
 probs <- mosaiClassifierPostProcessing(probs)
-tab <- makeSVCallSimple(probs)
+tab <- makeSVCallSimple(probs, llr_thr = llr)
 
 write.table(tab, file = snakemake@output[[1]], sep = "\t", quote=F, row.names = F, col.names = T)

utils/mosaiClassifier_call_biallelic.snakemake.R

@@ -3,10 +3,10 @@ library(data.table)
 source("utils/mosaiClassifier/makeSVcalls.R")
 
 probs = readRDS(snakemake@input[["probs"]])
-llr   = as.numeric(snakemake@params[["llr"]])
+llr   = as.numeric(snakemake@wildcards[["llr"]])
 
 probs <- mosaiClassifierPostProcessing(probs)
 probs <- forceBiallelic(probs)
-tab <- makeSVCallSimple(probs)
+tab <- makeSVCallSimple(probs, llr_thr = llr)
 
 write.table(tab, file = snakemake@output[[1]], sep = "\t", quote=F, row.names = F, col.names = T)