Skip to content
Snippets Groups Projects
Commit 162f02ed authored by Peter van 't Hof's avatar Peter van 't Hof
Browse files

Added test on main method

parent c7d28c03
Branches
Tags
No related merge requests found
...@@ -56,7 +56,7 @@ object BaseCounter extends ToolCommand { ...@@ -56,7 +56,7 @@ object BaseCounter extends ToolCommand {
val overlapGenes = groupGenesOnOverlap(geneReader.getAll) val overlapGenes = groupGenesOnOverlap(geneReader.getAll)
logger.info("Start reading bamFile") logger.info("Start reading bamFile")
val counts = (for (genes <- overlapGenes.values.flatten.par) yield bamToGeneCount(cmdArgs.bamFile, genes)).toList val counts = (for (genes <- overlapGenes.values.flatten.par) yield runThread(cmdArgs.bamFile, genes)).toList
logger.info("Done reading bamFile") logger.info("Done reading bamFile")
writeGeneCounts(counts.flatMap(_.geneCounts), cmdArgs.outputDir, cmdArgs.prefix) writeGeneCounts(counts.flatMap(_.geneCounts), cmdArgs.outputDir, cmdArgs.prefix)
...@@ -281,7 +281,7 @@ object BaseCounter extends ToolCommand { ...@@ -281,7 +281,7 @@ object BaseCounter extends ToolCommand {
nonStrandedMetaExonCounts: List[(String, RegionCount)], nonStrandedMetaExonCounts: List[(String, RegionCount)],
strandedMetaExonCounts: List[(String, RegionCount)]) strandedMetaExonCounts: List[(String, RegionCount)])
def bamToGeneCount(bamFile: File, genes: List[Gene]): ThreadOutput = { private[tools] def runThread(bamFile: File, genes: List[Gene]): ThreadOutput = {
val counts = genes.map(gene => gene -> new GeneCount(gene)).toMap val counts = genes.map(gene => gene -> new GeneCount(gene)).toMap
val bamReader = SamReaderFactory.makeDefault().open(bamFile) val bamReader = SamReaderFactory.makeDefault().open(bamFile)
...@@ -304,25 +304,27 @@ object BaseCounter extends ToolCommand { ...@@ -304,25 +304,27 @@ object BaseCounter extends ToolCommand {
} }
def createMetaExonCounts(genes: List[Gene]): List[(String, RegionCount)] = { def createMetaExonCounts(genes: List[Gene]): List[(String, RegionCount)] = {
val regions = genes.map(gene => gene.getName -> generateMergedExonRegions(gene).sorted) if (genes.nonEmpty) {
val chr = genes.head.getContig val regions = genes.map(gene => gene.getName -> generateMergedExonRegions(gene).sorted)
val begin = regions.map(_._2.allRecords.head.start).min val chr = genes.head.getContig
val end = regions.map(_._2.allRecords.last.end).max val begin = regions.map(_._2.allRecords.head.start).min
val end = regions.map(_._2.allRecords.last.end).max
val posibleEnds = (regions.flatMap(_._2.allRecords.map(_.end)) ++ regions.flatMap(_._2.allRecords.map(_.start))).distinct.sorted
val posibleEnds = (regions.flatMap(_._2.allRecords.map(_.end)) ++ regions.flatMap(_._2.allRecords.map(_.start))).distinct.sorted
def mergeRegions(newBegin: Int, output: List[(String, RegionCount)] = Nil): List[(String, RegionCount)] = {
val newEnds = posibleEnds.filter(_ > newBegin) def mergeRegions(newBegin: Int, output: List[(String, RegionCount)] = Nil): List[(String, RegionCount)] = {
if (newBegin > end || newEnds.isEmpty) output val newEnds = posibleEnds.filter(_ > newBegin)
else { if (newBegin > end || newEnds.isEmpty) output
val newEnd = newEnds.min else {
val record = BedRecord(chr, newBegin, newEnd) val newEnd = newEnds.min
val names = regions.filter(_._2.overlapWith(record).nonEmpty).map(_._1) val record = BedRecord(chr, newBegin, newEnd)
if (names.nonEmpty) mergeRegions(newEnd, (names.mkString(","), new RegionCount(record.start + 1, record.end)) :: output) val names = regions.filter(_._2.overlapWith(record).nonEmpty).map(_._1)
else mergeRegions(newEnd, output) if (names.nonEmpty) mergeRegions(newEnd, (names.mkString(","), new RegionCount(record.start + 1, record.end)) :: output)
else mergeRegions(newEnd, output)
}
} }
} mergeRegions(begin)
mergeRegions(begin) } else Nil
} }
def bamRecordBasesOverlap(samRecord: SAMRecord, start: Int, end: Int): Int = { def bamRecordBasesOverlap(samRecord: SAMRecord, start: Int, end: Int): Int = {
......
geneA t1 chrQ + 200 500 225 475 3 200,300,400 250,350,500
File added
File suppressed by a .gitattributes entry or the file's encoding is unsupported.
package nl.lumc.sasc.biopet.tools package nl.lumc.sasc.biopet.tools
import java.io.File
import java.nio.file.Paths
import com.google.common.io.Files
import htsjdk.samtools.{SAMReadGroupRecord, SAMSequenceRecord, SAMLineParser, SAMFileHeader} import htsjdk.samtools.{SAMReadGroupRecord, SAMSequenceRecord, SAMLineParser, SAMFileHeader}
import org.scalatest.Matchers import org.scalatest.Matchers
import org.scalatest.testng.TestNGSuite import org.scalatest.testng.TestNGSuite
...@@ -164,6 +168,22 @@ class BaseCounterTest extends TestNGSuite with Matchers { ...@@ -164,6 +168,22 @@ class BaseCounterTest extends TestNGSuite with Matchers {
assert(ab.exists(x => x._1 == "geneB" && x._2.start == 201 && x._2.end == 210)) assert(ab.exists(x => x._1 == "geneB" && x._2.start == 201 && x._2.end == 210))
assert(ab.exists(x => x._1 == "geneB" && x._2.start == 221 && x._2.end == 250)) assert(ab.exists(x => x._1 == "geneB" && x._2.start == 221 && x._2.end == 250))
} }
private def resourcePath(p: String): String = {
Paths.get(getClass.getResource(p).toURI).toString
}
@Test
def testMain: Unit = {
val outputDir = Files.createTempDir()
outputDir.deleteOnExit()
val prefix = "test"
val bamFile = new File(resourcePath("/empty.bam"))
val refflat = new File(resourcePath("/chrQ.refflat"))
main(Array("-o", outputDir.getAbsolutePath, "-p", prefix,
"-b", bamFile.getAbsolutePath, "-r", refflat.getAbsolutePath))
outputDir.list().size shouldBe 34
}
} }
object BaseCounterTest { object BaseCounterTest {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment