adding reading reports and methylation and coverage stats generation

.test/comparison/config.yaml

@@ -7,17 +7,12 @@ resources:
     genome: ".test/resources/genome/Malus_domestica_cultivar_Golden_Delicious-chr4.fasta"
   annot: ".test/resources/annot/WT_vs_Meth13-CG.bed" #optional must be "" empty string if not being used
 params:
-  method: bins
-  binSize: 100
-  kernelFunction: triangular
-  test: score
-  pseudocountM: 1
-  pseudocountN: 2
-  pValueThreshold: 0.01
-  minCytosinesCount: 4
-  minProportionDifference: 0.4
-  minGap: 0
-  minSize: 50
-  minReadsPerCytosine: 4
-  cores: 1
+  bins: TRUE
+  windowSize: 1000
+  stepSize: 1000
+  test: F # default F, Chisq for overdispersion
+  overdispersion: none # default none, NM for overdispersion
+  qValue: 0.01
+  minCov: 4
+  minDiff: 0.4
   context: ["CG", "CHG", "CHH"]

workflow/scripts/methylkit.R

@@ -26,3 +26,43 @@ qValue <- 0.01 #snakemake@params[["qValue"]]
 
 ### inferred params  ###
 method <- if(bins) "bins" else "bases"
+
+### create output directory if it doesn't already exist
+dir.create(file.path(outdir), showWarnings = FALSE)
+
+### reading bismark's CX reports, by default using the tabix database for minimal memory usage
+readingReports=methRead(files,
+           sample.id=c(rep("control",length(control)),rep("treatment",length(treatment))),
+           assembly="bismark",
+           treatment=c(rep(0,length(control)),rep(1,length(treatment))),
+           context=context,
+           mincov = minCov,
+           pipeline="bismarkCytosineReport",
+           dbtype = "tabix",
+           dbdir = "methylDB"
+           )
+
+
+### meth stats and plot ###
+for (i in 1:length(files)) {
+	fileTxt<-file(snakemake@output[["methylationStatsTxt"]])
+    sink(fileTxt)
+    getMethylationStats(readingReports[[i]],plot=FALSE,both.strands=FALSE)
+    sink()
+    close(fileTxt)
+	fileTxt<-file(snakemake@output[["methylationStatsPdf"]])
+    getMethylationStats(readingReports[[i]],plot=TRUE,both.strands=FALSE)
+    dev.off()
+    }
+
+### coverage stats and plot ###
+for (i in 1:length(files)) {
+	fileTxt<-file(snakemake@output[["coverageStatsTxt"]])
+    sink(fileTxt)
+    getCoverageStats(readingReports[[i]],plot=FALSE,both.strands=FALSE)
+    sink()
+    close(fileTxt)
+	fileTxt<-file(snakemake@output[["coverageStatsPdf"]])
+    getCoverageStats(readingReports[[i]],plot=TRUE,both.strands=FALSE)
+    dev.off()
+    }