Merge pull request #601 from szilvajuhos/master

GATK4 Singularity container and fixes for ROI

Author: maxulysse

containers/gatk4/Singularity

@@ -0,0 +1,19 @@
+Bootstrap: docker
+From: ubuntu:18.04
+
+%post
+    apt-get -y update
+    apt-get -y install openjdk-8-jre-headless
+
+%help
+	Custom container for GATK4 to contain only single jar as the result of "./gradlew localJar" in /opt ; you have to change the filename for the created jar below
+
+%files
+   /home/szilva/dev/gatk/build/libs/gatk-package-4.0.4.0-15-g61c6891-SNAPSHOT-local.jar	/opt/gatk4.jar
+
+%environment
+    GATK_HOME=/opt
+    export GATK_HOME
+
+%labels
+   AUTHOR szilveszter.juhos@scilifelab.se

scripts/selectROI.py

@@ -83,6 +83,7 @@ def __init__(self, threads, maxRegions, prefix, width, pad):
         # prefix of the ROI BAM file
         self.prefix = prefix
         # width of region for VCF records
+        self.pad = pad
         # add readlength padding
         self.width = width + pad
 
@@ -118,10 +119,34 @@ def selectROI(self,vcfs,bams,beds):
         # now collapse overlapping intervals
         for chrom in self.callDict:
             self.callDict[chrom] = reduce( self.callDict[chrom] )
+            # we still have to collapse intervals that are too close to each other 
+            # if they are closer than the readlength, there can be cases when a read is picked 
+            # in both interval, so we have to merge them
+            self.mergeCloseIntervals(chrom)
+
         # save intervals to a BAM
         print("Saving reads from BAMs:")
         self.saveCallsToBAM(bams)
 
+    def mergeCloseIntervals(self,chrom):
+        # we are assuming the intervals are ordered (they should to be)
+        mergedChrom = []
+        start = (0,0)
+        for interv in self.callDict[chrom]:
+            # (start[0], start[1])
+            #                          (interv[0], interv[1])
+            # |------------------|     |--------------------|     
+            if interv[0] - start[1] <= self.pad*2:
+                # merge the two intervals
+                start = (start[0],interv[1])
+            elif start[1]!= 0:  # if its not the very first one and have enough space between the intervals
+                # add to the new interval set
+                mergedChrom.append(start)
+                start = interv
+            else: # most of the intervals should fall here
+                start = interv
+        self.callDict[chrom] = mergedChrom
+
     def writeChunk(self,bam,chrom,regions):
         # from the regions list make a list of lists - sublists have 50 entries max
         sublists = [regions[i:i+self.maxRegions] for i  in range(0, len(regions), self.maxRegions)]