Fixed collection of bugs in reads aligning to multiple locations.

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

c/bwa/bwa_gateway.cpp

@@ -105,6 +105,8 @@ void BWA::copy_bases_into_sequence(bwa_seq_t* sequence, const char* bases, const
 }
 
 void BWA::create_alignments(bwa_seq_t* sequence, Alignment*& alignments, unsigned& num_alignments) {
+  bool debug = false;
+
   num_alignments = 0;
   for(unsigned i = 0; i < (unsigned)sequence->n_aln; i++)
     num_alignments += (sequence->aln + i)->l - (sequence->aln + i)->k + 1;
@@ -127,8 +129,22 @@ void BWA::create_alignments(bwa_seq_t* sequence, Alignment*& alignments, unsigne
       sequence->aln = &working_alignment_block;
       sequence->n_aln = 1;
 
+      sequence->sa = sa_idx;
+      sequence->strand = alignment_block->a;
+      sequence->score = alignment_block->score;
+
+      // Each time through bwa_refine_gapped, seq gets reversed.  Revert the reverse.
+      // TODO: Fix the interface to bwa_refine_gapped so its easier to work with.
+      if(alignment_idx > 0)
+	seq_reverse(sequence->len, sequence->seq, 0);
+
       // Calculate the local alignment.
       bwa_cal_pac_pos_core(bwts[0],bwts[1],sequence,options.max_diff,options.fnr);
+      if(debug) {
+	printf("alignment_idx = %d, k = %d, l = %d, sa_idx = %d\n", alignment_idx, alignment_block->k, alignment_block->l, sa_idx);
+	printf("sequence->pos = %d\n",sequence->pos);
+      }
+
       bwa_refine_gapped(bns, 1, sequence, 0, NULL);
 
       // Copy the local alignment data into the alignment object.

java/src/org/broadinstitute/sting/alignment/bwa/AlignmentValidationWalker.java

@@ -25,6 +25,11 @@ public class AlignmentValidationWalker extends ReadWalker<Integer,Integer> {
      */
     private BWACAligner aligner = null;
 
+    /**
+     * Number of reads processed.
+     */
+    private int count = 0;
+
     /**
      * Create an aligner object.  The aligner object will load and hold the BWT until close() is called.
      */
@@ -39,6 +44,15 @@ public class AlignmentValidationWalker extends ReadWalker<Integer,Integer> {
                                   prefix + ".rsa");        
     }
 
+    /** Must return true for reads that need to be processed. Reads, for which this method return false will
+     * be skipped by the engine and never passed to the walker.
+     */
+    @Override
+    public boolean filter(char[] ref, SAMRecord read) {
+        return true;
+        //return read.getReadName().contains("SL-XBC:1:76:604:397#0");
+    }
+
     /**
      * Aligns a read to the given reference.
      * @param ref Reference over the read.  Read will most likely be unmapped, so ref will be null.
@@ -47,6 +61,8 @@ public class AlignmentValidationWalker extends ReadWalker<Integer,Integer> {
      */
     @Override
     public Integer map(char[] ref, SAMRecord read) {
+        //logger.info(String.format("examining read %s", read.getReadName()));
+
         byte[] bases = read.getReadBases();
         if(read.getReadNegativeStrandFlag()) bases = BaseUtils.simpleReverseComplement(bases);
 
@@ -85,6 +101,10 @@ public class AlignmentValidationWalker extends ReadWalker<Integer,Integer> {
             throw new StingException(String.format("Read %s mismatches!", read.getReadName()));
         }
 
+        if(++count % 10000 == 0) logger.info(String.format("Processed %d reads", count));
+
+        //logger.info(String.format("validated read %s", read.getReadName()));        
+
         return 1;
     }