Added bedtools coverage to counts script

biopet-framework/examples/sage.json

@@ -9,5 +9,15 @@
         "seqtk": {"exe":"/data/DIV5/SASC/common/programs/seqtk/seqtk/seqtk"},
         "cutadapt": {"exe":"/home/pjvan_thof/.local/bin/cutadapt"},
         "sickle": {"exe":"/data/DIV5/SASC/pjvan_thof/bin/sickle"}
-    }
+    },
+    "samples": {
+        "test": {
+            "libraries": {
+                "1": {
+                    "R1": "test.fastq"
+                }
+            }
+        }
+    },
+    "bowtie": {"exe": "test"}
 }

biopet-framework/nbactions.xml

@@ -10,7 +10,7 @@
                 <goal>org.codehaus.mojo:exec-maven-plugin:1.2.1:exec</goal>
             </goals>
             <properties>
-                <exec.args>-classpath %classpath nl.lumc.sasc.biopet.pipelines.sage.CreateTagCounts -I raw.count -taglib tag.lib -sense sense.counts -allsense all.sense.counts -antisense antisense.counts -allantisense all.antisense.counts</exec.args>
+                <exec.args>-classpath %classpath nl.lumc.sasc.biopet.pipelines.sage.Sage -outDir /home/pjvan_thof/pipelines/test --transcriptome bla.fa --countbed bla.bed -config /home/pjvan_thof/pipelines/biopet/biopet-framework/examples/sage.json</exec.args>
                 <exec.executable>java</exec.executable>
                 <exec.workingdir>/home/pjvan_thof/pipelines/test</exec.workingdir>
             </properties>
@@ -25,7 +25,7 @@
                 <goal>org.codehaus.mojo:exec-maven-plugin:1.2.1:exec</goal>
             </goals>
             <properties>
-                <exec.args>-Xdebug -Xrunjdwp:transport=dt_socket,server=n,address=${jpda.address} -classpath %classpath nl.lumc.sasc.biopet.pipelines.sage.CreateTagCounts -I raw.count -taglib tag.lib -sense sense.counts -allsense all.sense.counts -antisense antisense.counts -allantisense all.antisense.counts</exec.args>
+                <exec.args>-Xdebug -Xrunjdwp:transport=dt_socket,server=n,address=${jpda.address} -classpath %classpath nl.lumc.sasc.biopet.pipelines.sage.Sage -outDir /home/pjvan_thof/pipelines/test --transcriptome bla.fa --countbed bla.bed -config /home/pjvan_thof/pipelines/biopet/biopet-framework/examples/sage.json</exec.args>
                 <exec.executable>java</exec.executable>
                 <jpda.listen>true</jpda.listen>
                 <exec.workingdir>/home/pjvan_thof/pipelines/test</exec.workingdir>
@@ -41,7 +41,7 @@
                 <goal>org.codehaus.mojo:exec-maven-plugin:1.2.1:exec</goal>
             </goals>
             <properties>
-                <exec.args>-classpath %classpath nl.lumc.sasc.biopet.pipelines.sage.CreateTagCounts -I raw.count -taglib tag.lib -sense sense.counts -allsense all.sense.counts -antisense antisense.counts -allantisense all.antisense.counts</exec.args>
+                <exec.args>-classpath %classpath nl.lumc.sasc.biopet.pipelines.sage.Sage -outDir /home/pjvan_thof/pipelines/test --transcriptome bla.fa --countbed bla.bed -config /home/pjvan_thof/pipelines/biopet/biopet-framework/examples/sage.json</exec.args>
                 <exec.executable>java</exec.executable>
                 <exec.workingdir>/home/pjvan_thof/pipelines/test</exec.workingdir>
             </properties>

biopet-framework/src/main/scala/nl/lumc/sasc/biopet/pipelines/sage/BedtoolsCoverageToCounts.scala

+package nl.lumc.sasc.biopet.pipelines.sage
+
+import java.io.File
+import java.io.PrintWriter
+import nl.lumc.sasc.biopet.core.BiopetJavaCommandLineFunction
+import nl.lumc.sasc.biopet.core.config.Configurable
+import org.broadinstitute.gatk.utils.commandline.{ Input, Output }
+import scala.collection.JavaConversions._
+import scala.collection.SortedMap
+import scala.collection.mutable.Map
+import scala.io.Source
+
+class BedtoolsCoverageToCounts(val root: Configurable) extends BiopetJavaCommandLineFunction {
+  javaMainClass = getClass.getName
+
+  @Input(doc = "Input fasta", shortName = "input", required = true)
+  var input: File = _
+  
+  @Output(doc = "Output tag library", shortName = "output", required = true)
+  var output: File = _
+  
+  override val defaultVmem = "8G"
+  memoryLimit = Option(4.0)
+    
+  override def commandLine = super.commandLine + 
+    required("-I", input) + 
+    required("-o", output)
+}
+
+object BedtoolsCoverageToCounts {
+  var input: File = _
+  var output: File = _
+  
+  /**
+   * @param args the command line arguments
+   */
+  def main(args: Array[String]): Unit = {
+    for (t <- 0 until args.size) {
+      args(t) match {
+        case "-I" => input = new File(args(t+1))
+        case "-o" => output = new File(args(t+1))
+        case _ =>
+      }
+    }
+    if (input == null || !input.exists) throw new IllegalStateException("Input file not found, use -I")
+    if (output == null) throw new IllegalStateException("Output file not found, use -o")
+    
+    val counts:Map[String, Long] = Map()
+    for (line <- Source.fromFile(input).getLines) {
+      val values = line.split("\t")
+      val gene = values(3)
+      val count = values(6).toLong
+      if (counts.contains(gene)) counts(gene) += count
+      else counts += gene -> count
+    }
+  
+    val sortedCounts:SortedMap[String, Long] = SortedMap(counts.toArray:_*)
+    
+    val writer = new PrintWriter(output)
+    for ((seq,count) <- sortedCounts) {
+      if (count > 0) writer.println(seq + "\t" + count)
+    }
+    writer.close
+  }
+}
\ No newline at end of file

biopet-framework/src/main/scala/nl/lumc/sasc/biopet/pipelines/sage/Sage.scala

@@ -142,12 +142,25 @@ class Sage(val root: Configurable) extends QScript with MultiSampleQScript {
   }
   
   def addBedtoolsCounts(bamFile:File, outputPrefix: String, outputDir: String) {
-    add(BedtoolsCoverage(this, bamFile, squishedCountBed, outputDir + outputPrefix + ".genome.sense.counts", 
-                         depth = false, sameStrand = true, diffStrand = false))
-    add(BedtoolsCoverage(this, bamFile, squishedCountBed, outputDir + outputPrefix + ".genome.antisense.counts", 
-                         depth = false, sameStrand = false, diffStrand = true))
-    add(BedtoolsCoverage(this, bamFile, squishedCountBed, outputDir + outputPrefix + ".genome.counts", 
-                         depth = false, sameStrand = false, diffStrand = false))
+    val bedtoolsSense = BedtoolsCoverage(this, bamFile, squishedCountBed, outputDir + outputPrefix + ".genome.sense.coverage", 
+                         depth = false, sameStrand = true, diffStrand = false)
+    val countSense = new BedtoolsCoverageToCounts(this)
+    countSense.input = bedtoolsSense.output
+    countSense.output = outputDir + outputPrefix + ".genome.sense.counts"
+    
+    val bedtoolsAntisense = BedtoolsCoverage(this, bamFile, squishedCountBed, outputDir + outputPrefix + ".genome.antisense.coverage", 
+                         depth = false, sameStrand = false, diffStrand = true)
+    val countAntisense = new BedtoolsCoverageToCounts(this)
+    countAntisense.input = bedtoolsAntisense.output
+    countAntisense.output = outputDir + outputPrefix + ".genome.antisense.counts"
+    
+    val bedtools = BedtoolsCoverage(this, bamFile, squishedCountBed, outputDir + outputPrefix + ".genome.coverage", 
+                         depth = false, sameStrand = false, diffStrand = false)
+    val count = new BedtoolsCoverageToCounts(this)
+    count.input = bedtools.output
+    count.output = outputDir + outputPrefix + ".genome.counts"
+    
+    add(bedtoolsSense, countSense, bedtoolsAntisense, countAntisense, bedtools, count)
   }
   
   def addTablibCounts(fastq:File, outputPrefix: String, outputDir: String) {