Change in Methods development pipeline: dbsnp130 can't be used for anything, changed it to dbsnp129. Optimization for HaplotypeScore and the to-be-committed ReadRosRankSumTest in AlignmentUtils

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@5301 348d0f76-0448-11de-a6fe-93d51630548a

java/src/org/broadinstitute/sting/gatk/walkers/recalibration/TableRecalibrationWalker.java

@@ -63,7 +63,7 @@ import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
  *
  * @author rpoplin
  * @since Nov 3, 2009
- * @help.summary Second pass of the recalibration. Uses the table generated by CountCovariates to update the base quality scores of the input bam file using a sequential table calculation making the base quality scores more accurately reflect the actual quality of the bases as mesaured by reference mismatch rate.
+ * @help.summary Second pass of the recalibration. Uses the table generated by CountCovariates to update the base quality scores of the input bam file using a sequential table calculation making the base quality scores more accurately reflect the actual quality of the bases as measured by reference mismatch rate.
  */
 
 @BAQMode(QualityMode = BAQ.QualityMode.ADD_TAG, ApplicationTime = BAQ.ApplicationTime.ON_OUTPUT)

java/src/org/broadinstitute/sting/utils/sam/AlignmentUtils.java

@@ -31,6 +31,7 @@ import net.sf.samtools.Cigar;
 import net.sf.samtools.CigarElement;
 import net.sf.samtools.util.StringUtil;
 import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
+import org.broadinstitute.sting.utils.BaseUtils;
 import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
 import org.broadinstitute.sting.utils.pileup.*;
 import org.broadinstitute.sting.utils.Utils;
@@ -368,6 +369,18 @@ public class AlignmentUtils {
     	return n;
     }
 
+    public static int getNumAlignedBases(final SAMRecord r) {
+    	int n = 0;
+        final Cigar cigar = r.getCigar();
+        if (cigar == null) return 0;
+
+        for (final CigarElement e : cigar.getCigarElements()) {
+        	if (e.getOperator() == CigarOperator.M ) { n += e.getLength(); }
+        }
+
+    	return n;
+    }
+
     @Deprecated
     public static char[] alignmentToCharArray( final Cigar cigar, final char[] read, final char[] ref ) {
 
@@ -468,37 +481,30 @@ public class AlignmentUtils {
             switch( ce.getOperator() ) {
             case I:
             case S:
-                for ( int jjj = 0; jjj < elementLength; jjj++ ) {
+                pos += elementLength;
-                    if( pos == pileupOffset ) {
+                if( pos >= pileupOffset ) {
-                        return alignmentPos;
+                    return alignmentPos;
-                    }
-                    pos++;
                 }
                 break;
             case D:
             case N:
                 if(!atDeletion) {
-                    for ( int jjj = 0; jjj < elementLength; jjj++ ) {
+                    alignmentPos += elementLength;
-                        alignmentPos++;
-                    }
                 } else {
-                    for ( int jjj = 0; jjj < elementLength; jjj++ ) {
+                    if( pos + elementLength >= pileupOffset ) {
-                        if( pos == pileupOffset ) {
+                        return alignmentPos + (pileupOffset - pos);
-                            return alignmentPos;
+                    } else {
-                        }
+                        pos += elementLength;
-                        pos++;
+                        alignmentPos += elementLength;
-                        alignmentPos++;
                     }
-
                 }
                 break;
             case M:
-                for ( int jjj = 0; jjj < elementLength; jjj++ ) {
+                if( pos + elementLength >= pileupOffset ) {
-                    if( pos == pileupOffset ) {
+                    return alignmentPos + (pileupOffset - pos);
-                        return alignmentPos;
+                } else {
-                    }
+                    pos += elementLength;
-                    pos++;
+                    alignmentPos += elementLength;
-                    alignmentPos++;
                 }
                 break;
             case H:
@@ -548,6 +554,14 @@ public class AlignmentUtils {
 
             switch( ce.getOperator() ) {
             case I:
+                /*
+                if( alignPos > 0 ) {
+                    if( alignment[alignPos-1] == BaseUtils.A ) { alignment[alignPos-1] = PileupElement.A_FOLLOWED_BY_INSERTION_BASE; }
+                    else if( alignment[alignPos-1] == BaseUtils.C ) { alignment[alignPos-1] = PileupElement.C_FOLLOWED_BY_INSERTION_BASE; }
+                    else if( alignment[alignPos-1] == BaseUtils.T ) { alignment[alignPos-1] = PileupElement.T_FOLLOWED_BY_INSERTION_BASE; }
+                    else if( alignment[alignPos-1] == BaseUtils.G ) { alignment[alignPos-1] = PileupElement.G_FOLLOWED_BY_INSERTION_BASE; }
+                }
+                */
             case S:
                 for ( int jjj = 0; jjj < elementLength; jjj++ ) {
                     readPos++;

scala/qscript/core/MethodsDevelopmentCallingPipeline.scala

@@ -73,9 +73,8 @@ class MethodsDevelopmentCallingPipeline extends QScript {
   val hg18 = new File("/seq/references/Homo_sapiens_assembly18/v0/Homo_sapiens_assembly18.fasta")
   val b36 = new File("/humgen/1kg/reference/human_b36_both.fasta")
   val b37 = new File("/humgen/1kg/reference/human_g1k_v37.fasta")
-  val dbSNP_hg18 = "/humgen/gsa-hpprojects/GATK/data/Comparisons/Validated/dbSNP/dbsnp_130_hg18.rod"                // Special case for NA12878 collections that can't use 132 because they are part of it.
+  val dbSNP_hg18_129 = "/humgen/gsa-hpprojects/GATK/data/Comparisons/Validated/dbSNP/dbsnp_129_hg18.rod"            // Special case for NA12878 collections that can't use 132 because they are part of it.
-  val dbSNP_hg18_129 = "/humgen/gsa-hpprojects/GATK/data/Comparisons/Validated/dbSNP/dbsnp_129_hg18.rod"
+  val dbSNP_b36 = "/humgen/gsa-hpprojects/GATK/data/Comparisons/Validated/dbSNP/dbsnp_129_b36.rod"
-  val dbSNP_b36 = "/humgen/gsa-hpprojects/GATK/data/Comparisons/Validated/dbSNP/dbsnp_130_b36.rod"
   val dbSNP_b37 = "/humgen/gsa-hpprojects/GATK/data/Comparisons/Validated/dbSNP/dbsnp_132_b37.leftAligned.vcf"
   val dbSNP_b37_129 = "/humgen/gsa-hpprojects/GATK/data/Comparisons/Validated/dbSNP/dbsnp_129_b37.rod"              // Special case for NA12878 collections that can't use 132 because they are part of it.
   val hapmap_hg18 = "/humgen/gsa-hpprojects/GATK/data/Comparisons/Validated/HapMap/3.3/sites_r27_nr.hg18_fwd.vcf"