From c63af32fc78d4c1ed9a9658f1fcb680aa7d2dc2e Mon Sep 17 00:00:00 2001
From: hanna <hanna@348d0f76-0448-11de-a6fe-93d51630548a>
Date: Tue, 3 Nov 2009 00:01:55 +0000
Subject: [PATCH] The BWA/C bindings were triggering the local aligner to
 repeatedly reload the ref genome.  Make sure the reference genome is cached.

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1961 348d0f76-0448-11de-a6fe-93d51630548a
---
 c/bwa/bwa_gateway.cpp                         | 19 ++++++++++---------
 c/bwa/bwa_gateway.h                           |  1 +
 .../bwa/AlignmentValidationWalker.java        |  8 ++++++++
 3 files changed, 19 insertions(+), 9 deletions(-)

diff --git a/c/bwa/bwa_gateway.cpp b/c/bwa/bwa_gateway.cpp
index 538b8b16a..a5fb045ec 100644
--- a/c/bwa/bwa_gateway.cpp
+++ b/c/bwa/bwa_gateway.cpp
@@ -12,7 +12,13 @@ BWA::BWA(const char* ann_filename,
 	 const char* reverse_bwt_filename, 
 	 const char* reverse_sa_filename) 
 {
+  // Load the bns (?) and reference
   bns = bns_restore_core(ann_filename,amb_filename,pac_filename);
+  reference = new ubyte_t[bns->l_pac/4+1];
+  rewind(bns->fp_pac);
+  fread(reference, 1, bns->l_pac/4+1, bns->fp_pac);
+
+  // Load the BWTs (both directions) and suffix arrays (both directions)
   bwts[0] = bwt_restore_bwt(forward_bwt_filename);
   bwt_restore_sa(forward_sa_filename, bwts[0]);
   bwts[1] = bwt_restore_bwt(reverse_bwt_filename);
@@ -21,6 +27,8 @@ BWA::BWA(const char* ann_filename,
 }
 
 BWA::~BWA() {
+  delete[] reference;
+  bns_destroy(bns);
   bwt_destroy(bwts[0]);
   bwt_destroy(bwts[1]);
 }
@@ -105,8 +113,6 @@ 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;
@@ -138,14 +144,9 @@ void BWA::create_alignments(bwa_seq_t* sequence, Alignment*& alignments, unsigne
       if(alignment_idx > 0)
 	seq_reverse(sequence->len, sequence->seq, 0);
 
-      // Calculate the local alignment.
+      // Calculate the local coordinate and 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);
+      bwa_refine_gapped(bns, 1, sequence, reference, NULL);
 
       // Copy the local alignment data into the alignment object.
       Alignment& alignment = *(alignments + alignment_idx);
diff --git a/c/bwa/bwa_gateway.h b/c/bwa/bwa_gateway.h
index d517acb6a..bb8f72d68 100644
--- a/c/bwa/bwa_gateway.h
+++ b/c/bwa/bwa_gateway.h
@@ -12,6 +12,7 @@ class Alignment;
 class BWA {
  private:
   bntseq_t *bns;
+  ubyte_t* reference;
   bwt_t* bwts[2];
   gap_opt_t options;
 
diff --git a/java/src/org/broadinstitute/sting/alignment/bwa/AlignmentValidationWalker.java b/java/src/org/broadinstitute/sting/alignment/bwa/AlignmentValidationWalker.java
index 72c566f6a..50cba3151 100644
--- a/java/src/org/broadinstitute/sting/alignment/bwa/AlignmentValidationWalker.java
+++ b/java/src/org/broadinstitute/sting/alignment/bwa/AlignmentValidationWalker.java
@@ -30,6 +30,11 @@ public class AlignmentValidationWalker extends ReadWalker<Integer,Integer> {
      */
     private int count = 0;
 
+    /**
+     * Max delta in mapping quality.
+     */
+    private int maxMapQDelta = 0;
+
     /**
      * Create an aligner object.  The aligner object will load and hold the BWT until close() is called.
      */
@@ -77,6 +82,8 @@ public class AlignmentValidationWalker extends ReadWalker<Integer,Integer> {
                 matches &= (alignment.getAlignmentStart() == read.getAlignmentStart());
                 matches &= (alignment.isNegativeStrand() == read.getReadNegativeStrandFlag());
                 matches &= (alignment.getCigar().equals(read.getCigar()));
+                int mapQDelta = Math.abs(alignment.getMappingQuality()-read.getMappingQuality());
+                maxMapQDelta = Math.max(mapQDelta,maxMapQDelta);
                 //matches &= (alignment.getMappingQuality() == read.getMappingQuality());
                 if(matches) break;
             }
@@ -133,6 +140,7 @@ public class AlignmentValidationWalker extends ReadWalker<Integer,Integer> {
     @Override
     public void onTraversalDone(Integer result) {
         aligner.close();
+        out.printf("Max mapping quality delta: %s%n", maxMapQDelta);
         super.onTraversalDone(result);
     }