Mauricio pointed out to me that dynamic merging the unmapped regions …

…of multiple BAMs ('-L unmapped' with a BAM list) was completely broken. Sorry about this! Fixed.
eroller · Sep 13, 2011 · e63d9d8 · e63d9d8
1 parent 367bbee
commit e63d9d8
Show file tree

Hide file tree

Showing 3 changed files with 45 additions and 18 deletions.
diff --git a/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/BAMScheduler.java b/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/BAMScheduler.java
@@ -26,6 +26,7 @@
 
 import net.sf.picard.util.PeekableIterator;
 import net.sf.samtools.GATKBAMFileSpan;
+import net.sf.samtools.GATKChunk;
 import org.broadinstitute.sting.utils.GenomeLoc;
 import org.broadinstitute.sting.utils.GenomeLocSortedSet;
 
@@ -84,7 +85,7 @@ private void advance() {
             if(currentLocus == GenomeLoc.UNMAPPED) {
                 nextFilePointer = new FilePointer(GenomeLoc.UNMAPPED);
                 for(SAMReaderID id: dataSource.getReaderIDs())
-                    nextFilePointer.addFileSpans(id,new GATKBAMFileSpan());
+                    nextFilePointer.addFileSpans(id,new GATKBAMFileSpan(new GATKChunk(indexFiles.get(id).getStartOfLastLinearBin(),Long.MAX_VALUE)));
                 currentLocus = null;
                 continue;
             }

diff --git a/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/GATKBAMIndex.java b/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/GATKBAMIndex.java
@@ -215,6 +215,45 @@ public int getLastLocusInBin(final Bin bin) {
         return (new GATKBin(bin).getBinNumber()-levelStart+1)*(BIN_GENOMIC_SPAN /levelSize);
     }
 
+    /**
+     * Use to get close to the unmapped reads at the end of a BAM file.
+     * @return The file offset of the first record in the last linear bin, or -1
+     * if there are no elements in linear bins (i.e. no mapped reads).
+     */
+    public long getStartOfLastLinearBin() {
+        openIndexFile();
+
+        seek(4);
+
+        final int sequenceCount = readInteger();
+        // Because no reads may align to the last sequence in the sequence dictionary,
+        // grab the last element of the linear index for each sequence, and return
+        // the last one from the last sequence that has one.
+        long lastLinearIndexPointer = -1;
+        for (int i = 0; i < sequenceCount; i++) {
+            // System.out.println("# Sequence TID: " + i);
+            final int nBins = readInteger();
+            // System.out.println("# nBins: " + nBins);
+            for (int j1 = 0; j1 < nBins; j1++) {
+                // Skip bin #
+                skipBytes(4);
+                final int nChunks = readInteger();
+                // Skip chunks
+                skipBytes(16 * nChunks);
+            }
+            final int nLinearBins = readInteger();
+            if (nLinearBins > 0) {
+                // Skip to last element of list of linear bins
+                skipBytes(8 * (nLinearBins - 1));
+                lastLinearIndexPointer = readLongs(1)[0];
+            }
+        }
+
+        closeIndexFile();
+
+        return lastLinearIndexPointer;
+    }
+
     /**
      * Gets the possible number of bins for a given reference sequence.
      * @return How many bins could possibly be used according to this indexing scheme to index a single contig.

diff --git a/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/ReadShardStrategy.java b/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/ReadShardStrategy.java
@@ -134,24 +134,11 @@ public void advance() {
             Map<SAMReaderID,SAMFileSpan> selectedReaders = new HashMap<SAMReaderID,SAMFileSpan>();
             while(selectedReaders.size() == 0 && currentFilePointer != null) {
                 shardPosition = currentFilePointer.fileSpans;
+
                 for(SAMReaderID id: shardPosition.keySet()) {
-                    // If the region contains location information (in other words, it is not at
-                    // the start of the unmapped region), add the region.
-                    if(currentFilePointer.isRegionUnmapped) {
-                        // If the region is unmapped and no location data exists, add a null as an indicator to
-                        // start at the next unmapped region.
-                        if(!isIntoUnmappedRegion) {
-                            selectedReaders.put(id,null);
-                            isIntoUnmappedRegion = true;
-                        }
-                        else
-                            selectedReaders.put(id,position.get(id));
-                    }
-                    else {
-                        SAMFileSpan fileSpan = shardPosition.get(id).removeContentsBefore(position.get(id));
-                        if(!fileSpan.isEmpty())
-                            selectedReaders.put(id,fileSpan);
-                    }
+                    SAMFileSpan fileSpan = shardPosition.get(id).removeContentsBefore(position.get(id));
+                    if(!fileSpan.isEmpty())
+                        selectedReaders.put(id,fileSpan);
                 }
 
                 if(selectedReaders.size() > 0) {