From 9cae33902ab50f91d82bf0ffa7a7b0aad0af7d49 Mon Sep 17 00:00:00 2001
From: Peter van 't Hof <p.j.van_t_hof@lumc.nl>
Date: Wed, 27 Jan 2016 19:44:47 +0100
Subject: [PATCH] Fixed createMetaExonCounts method
---
.../lumc/sasc/biopet/tools/BaseCounter.scala | 29 ++++++++++++-------
1 file changed, 18 insertions(+), 11 deletions(-)
diff --git a/public/biopet-tools/src/main/scala/nl/lumc/sasc/biopet/tools/BaseCounter.scala b/public/biopet-tools/src/main/scala/nl/lumc/sasc/biopet/tools/BaseCounter.scala
index 9cfb85103..5f251506c 100644
--- a/public/biopet-tools/src/main/scala/nl/lumc/sasc/biopet/tools/BaseCounter.scala
+++ b/public/biopet-tools/src/main/scala/nl/lumc/sasc/biopet/tools/BaseCounter.scala
@@ -248,18 +248,25 @@ object BaseCounter extends ToolCommand {
counts.values.toList
}
- //FIXME: This is not yet working
- def createMetaExonCounts(genes: List[Gene]): List[String, RegionCount] = {
- val regions = (for (gene <- genes; region <- generateMergedExonRegions(gene).allRecords)
- yield gene.getName -> region).sortBy(_._2.start)
-
- def mergeRegions(todo: List[(String, BedRecord)],
- inOverlap: List[(String, BedRecord)] = Nil,
- output: List[(String, RegionCount)] = Nil): List[(String, RegionCount)] = {
- if (todo.head._2.overlapWith(todo.tail.head._2)) Nil
- else Nil
+ def createMetaExonCounts(genes: List[Gene]): List[(String, RegionCount)] = {
+ val regions = genes.map(gene => gene.getName -> generateMergedExonRegions(gene).sorted)
+ val chr = genes.head.getContig
+ 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
+
+ def mergeRegions(newBegin: Int, output: List[(String, RegionCount)] = Nil): List[(String, RegionCount)] = {
+ if (newBegin > end) output
+ else {
+ val newEnd = posibleEnds.filter(_ > begin).min
+ val record = BedRecord(chr, newBegin, newEnd)
+ val names = regions.filter(_._2.overlapWith(record).size > 0).map(_._1)
+ if (names.nonEmpty) mergeRegions(newEnd, (names.mkString(","), new RegionCount(record.start + 1, record.end)) :: output)
+ else mergeRegions(newEnd, output)
+ }
}
- Nil
+ mergeRegions(begin)
}
def bamRecordBasesOverlap(samRecord: SAMRecord, start: Int, end: Int): Int = {
--
GitLab