From 4192b093b88ec31f2d00c63daf71b6fd03bc6936 Mon Sep 17 00:00:00 2001
From: chartl <chartl@348d0f76-0448-11de-a6fe-93d51630548a>
Date: Tue, 27 Oct 2009 17:20:44 +0000
Subject: [PATCH] More robust error handling with parallelization +
 usePreviousBase. Added forceReadBasesToMatchRef to use in conjunction with
 nPreviousReadBases as a less stringent approximation of usePreviousBases
 (requiring previous pileups only had mismatches, and that read mapping
 quality be high was throwing everything away)

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1916 348d0f76-0448-11de-a6fe-93d51630548a
---
 ...seTransitionTableCalculatorJavaWalker.java | 39 ++++++++++++++++---
 1 file changed, 34 insertions(+), 5 deletions(-)

diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/BaseTransitionTableCalculatorJavaWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/BaseTransitionTableCalculatorJavaWalker.java
index 1abbc40e3..746d4cbbd 100644
--- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/BaseTransitionTableCalculatorJavaWalker.java
+++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/BaseTransitionTableCalculatorJavaWalker.java
@@ -47,6 +47,8 @@ public class BaseTransitionTableCalculatorJavaWalker extends LocusWalker<Set<Bas
     boolean useReadGroup = false;
     @Argument(fullName="outputFile", shortName="of", doc="Output to this file rather than standard out. Must be used with -nt.", required = false)
     String outFilePath = null;
+    @Argument(fullName="forcePreviousReadBasesToMatchRef", doc="Forces previous read bases to match the reference", required = false)
+    boolean readBasesMustMatchRef = false;
 
     private UnifiedGenotyper ug;
     // private ReferenceContextWindow refWindow;
@@ -55,7 +57,7 @@ public class BaseTransitionTableCalculatorJavaWalker extends LocusWalker<Set<Bas
     private List<GenomeLoc> previousBaseLoci;
 
     public void initialize() {
-        if ( nPreviousBases > 3 ) {
+        if ( nPreviousBases > 3 || ( nPreviousReadBases > 3 && readBasesMustMatchRef ) ) {
             throw new StingException("You have opted to use a number of previous bases in excess of 3. In order to do this you must change the reference window size in the walker itself.");
         }
         ug = new UnifiedGenotyper();
@@ -116,6 +118,13 @@ public class BaseTransitionTableCalculatorJavaWalker extends LocusWalker<Set<Bas
     }
 
     public Set<BaseTransitionTable> treeReduce( Set<BaseTransitionTable> reduce1, Set<BaseTransitionTable> reduce2 ) {
+        // check to see if this is a truly tree-reducable calculation
+        if ( nPreviousBases >= 1 ) {
+            String errMsg = "Parallelization cannot be used with UsePreviousBases due to the fact that internal walker data specifies whether a previous reference base is usable or not.";
+            String errMsg2 = " This can cause cause concurrency issues and unpredictable behavior when used with parallelization. Either do not specify -nt, or try a the conjunction of ";
+            String errMsg3 = "--usePreviousReadBases and --forcePreviousReadBasesToMatchRef.";
+            throw new StingException(errMsg+errMsg2+errMsg3);
+        }
         return reduce(reduce1,reduce2);
     }
 
@@ -215,11 +224,15 @@ public class BaseTransitionTableCalculatorJavaWalker extends LocusWalker<Set<Bas
             return false;
         } else if ( read.getBaseQualities()[offset] <= minQualityScore ) {
             return false;
-        } else if ( useSecondaryBase && read.getAttribute("SQ") == null )
+        } else if ( useSecondaryBase && read.getAttribute("SQ") == null ) {
             return false;
-        else if ( nPreviousBases >= 1 && previousReadBasesMismatchRef(read, offset, ref) )
+        } else if ( nPreviousBases >= 1 && previousReadBasesMismatchRef(read, offset, ref) ) {
             return false;
-        else {
+        } else if ( nPreviousReadBases >= 1 && readLacksPreviousBases(read,offset,nPreviousReadBases) ) {
+            return false;
+        } else if ( nPreviousReadBases >= 1 && readBasesMustMatchRef && previousReadBasesMismatchRef(read, offset, ref) ) {
+            return false;
+        } else {
             return true;
         }
     }
@@ -241,6 +254,14 @@ public class BaseTransitionTableCalculatorJavaWalker extends LocusWalker<Set<Bas
         return false;
     }
 
+    public boolean readLacksPreviousBases( SAMRecord read, int offset, int prevBases ) {
+        if ( ! read.getReadNegativeStrandFlag() ) {
+            return offset - prevBases < 0;
+        } else {
+            return offset + prevBases + 1 >= read.getReadLength();
+        }
+    }
+
     public List<Comparable> buildConditions( SAMRecord read, int offset, ReferenceContext ref, ReadBackedPileup pileup ) {
         ArrayList<Comparable> conditions = new ArrayList<Comparable>();
 
@@ -254,7 +275,7 @@ public class BaseTransitionTableCalculatorJavaWalker extends LocusWalker<Set<Bas
         }
 
         if ( nPreviousReadBases > 0 ) {
-            conditions.add(read.getReadString().substring(offset-nPreviousReadBases,offset));
+            conditions.add(buildReadString(read, offset, nPreviousReadBases));
         }
 
         if ( usePileupMismatches ) {
@@ -276,6 +297,14 @@ public class BaseTransitionTableCalculatorJavaWalker extends LocusWalker<Set<Bas
         }
     }
 
+    public String buildReadString( SAMRecord read, int offset, int nPreviousReadBases ) {
+        if ( ! read.getReadNegativeStrandFlag() ) {
+            return read.getReadString().substring(offset-nPreviousReadBases,offset);
+        } else {
+            return BaseUtils.simpleReverseComplement( read.getReadString().substring(offset+1,offset+nPreviousReadBases+1) );
+        }
+    }
+
     public String createHeaderFromConditions() {
         String header = "Observed_base\tTrue_base";