generate test and update workflow

workflows/epigenetics/peaks-consensus/generate_fake_reads.R

@@ -0,0 +1,73 @@
+set.seed(1)
+peak.width <- 100
+frag.length <- 300
+read.length <- 40
+gauss.centers <- c(74670000, 74668000, 74666000)
+reads.per.peaks <- c(100, 50, 150)
+inputs.df <- data.frame(gauss.center = gauss.centers, reads.per.peak = reads.per.peaks)
+
+rand.center <- do.call(c, apply(inputs.df, 1, function(v){
+  rnorm(n = v["reads.per.peak"], mean = v["gauss.center"], sd = peak.width)
+}))
+rand.frag.length <- rnorm(n = sum(reads.per.peaks), mean = frag.length, sd = frag.length / 4)
+
+frag <- data.frame(start = as.integer(rand.center - rand.frag.length / 2),
+                   end = as.integer(rand.center + rand.frag.length / 2))
+frag$start1 <- frag$start
+frag$end1 <- frag$start + read.length
+frag$start2 <- frag$end - read.length
+frag$end2 <- frag$end
+
+frag$end1[frag$end1 > frag$frag$end] <- frag$end[frag$end1 > frag$frag$end] - 1
+frag$start2[frag$start2 < frag$frag$start] <- frag$start[frag$start2 < frag$frag$start] + 1
+
+frag$chr <- "chr2"
+frag$name <- paste0("read", 1:nrow(frag))
+frag$plus <- "+"
+frag$minus <- "-"
+frag$score <- 0
+write.table(frag[, c("chr", "start1", "end1", "name", "score", "plus")],
+            file = "/tmp/sample1_R1.bed",
+            sep = "\t", col.names = F, row.names = F, quote = F)
+write.table(frag[, c("chr", "start2", "end2", "name", "score", "minus")],
+            file = "/tmp/sample1_R2.bed",
+            sep = "\t", col.names = F, row.names = F, quote = F)
+set.seed(2)
+peak.width <- 100
+frag.length <- 350
+read.length <- 40
+gauss.centers <- c(74670050, 74668300, 74665000)
+reads.per.peaks <- c(1000, 3000, 1000)
+inputs.df <- data.frame(gauss.center = gauss.centers, reads.per.peak = reads.per.peaks)
+
+rand.center <- do.call(c, apply(inputs.df, 1, function(v){
+  rnorm(n = v["reads.per.peak"], mean = v["gauss.center"], sd = peak.width)
+}))
+rand.frag.length <- rnorm(n = sum(reads.per.peaks), mean = frag.length, sd = frag.length / 4)
+
+frag <- data.frame(start = as.integer(rand.center - rand.frag.length / 2),
+                   end = as.integer(rand.center + rand.frag.length / 2))
+frag$start1 <- frag$start
+frag$end1 <- frag$start + read.length
+frag$start2 <- frag$end - read.length
+frag$end2 <- frag$end
+
+frag$end1[frag$end1 > frag$frag$end] <- frag$end[frag$end1 > frag$frag$end] - 1
+frag$start2[frag$start2 < frag$frag$start] <- frag$start[frag$start2 < frag$frag$start] + 1
+
+frag$chr <- "chr2"
+frag$name <- paste0("read", 1:nrow(frag))
+frag$plus <- "+"
+frag$minus <- "-"
+frag$score <- 0
+write.table(frag[, c("chr", "start1", "end1", "name", "score", "plus")],
+            file = "/tmp/sample2_R1.bed",
+            sep = "\t", col.names = F, row.names = F, quote = F)
+write.table(frag[, c("chr", "start2", "end2", "name", "score", "minus")],
+            file = "/tmp/sample2_R2.bed",
+            sep = "\t", col.names = F, row.names = F, quote = F)
+
+# Then in bash
+# for sample in 1 2; do for i in 1 2; do bedtools getfasta -nameOnly -s -fi ~/igv/genomes/seq/mm10.fa -bed sample${sample}_R${i}.bed | seqtk seq -F 'A' - > sample${sample}_R${i}.fastq; done; done
+# Then in galaxy bowtie2 on mm10
+

workflows/epigenetics/peaks-consensus/peaks-consensus-chip-pe-tests.yml

@@ -0,0 +1,48 @@
+- doc: Test for peaks-consensus-chip-pe.ga
+  job:
+    2 rmDup BAMPE:
+      class: Collection
+      collection_type: list
+      elements:
+      - class: File
+        identifier: rep1
+        path: test-data/rep1.bam
+        dbkey: mm10
+      - class: File
+        identifier: rep2
+        path: test-data/rep2.bam
+        dbkey: mm10
+    effective_genome_size: 1870000000
+    bin_size: 50
+  outputs:
+    individual_macs2_narrowPeaks:
+      element_tests:
+        rep1:
+          asserts:
+            has_n_lines:
+              n: 3
+        rep2:
+          asserts:
+            has_n_lines:
+              n: 3
+    strict_intersect_peaks:
+      asserts:
+        has_n_lines:
+          n: 2
+    average_bigwig:
+      asserts:
+        has_size:
+          value: 1290
+          delta: 100
+    merged_macs2_narrowPeaks:
+      asserts:
+        has_n_lines:
+          n: 3
+    multiqc_output:
+      asserts:
+        has_text:
+          text: "<span class=\"mqc_table_tooltip\" title=\"MACS2: Fragment Length\">Fragment Length</span>"
+    shared_narrowPeak:
+      asserts:
+        has_n_lines:
+          n: 1