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

@@ -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) {
+      bwa_refine_gapped(bns, 1, sequence, reference, NULL);
-	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.
       Alignment& alignment = *(alignments + alignment_idx);

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;
 

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);
     }