Now one can specify the expected insert size distribution as supposed… (

#19)

* Now one can specify the expected insert size distribution as supposed to allow
BWA to infer it on the fly. This might have some value when re-aligning a small number of reads (inference may fail)
or when one knows that there is a heterozygous variant and we actually now the expected insert size distribution; the
input data would have a mixture of gaussians rather than the expected single gaussian.

src/main/c/Makefile

@@ -19,7 +19,7 @@ JNI_BASE_NAME=org_broadinstitute_hellbender_utils_bwa_BwaMemIndex
 
 all: libbwa.$(LIB_EXT)
 
-libbwa.$(LIB_EXT): $(JNI_BASE_NAME).o jnibwa.o bwa/libbwa.a
+libbwa.$(LIB_EXT): $(JNI_BASE_NAME).o init.o jnibwa.o bwa/libbwa.a
 	$(CC) -ggdb -dynamiclib -shared -o $@ $^ -lm -lz -lpthread
 
 bwa:
@@ -31,7 +31,9 @@ bwa/libbwa.a: bwa
 
 $(JNI_BASE_NAME).o: $(JNI_BASE_NAME).c jnibwa.h bwa
 
-jnibwa.o: jnibwa.c jnibwa.h bwa
+jnibwa.o: jnibwa.c jnibwa.h init.h bwa
+
+init.o: init.c init.h
 
 clean:
 	rm -rf bwa *.o *.$(LIB_EXT)

src/main/c/init.c

@@ -0,0 +1,30 @@
+#include "init.h"
+
+jfieldID peStatClass_failedID;
+jfieldID peStatClass_lowID;
+jfieldID peStatClass_highID;
+jfieldID peStatClass_averageID;
+jfieldID peStatClass_stdID;
+
+
+#define FAIL_IF_NULL(x) if (!(x)) return JNI_ERR;
+
+jint JNI_OnLoad(JavaVM* vm, void* reserved) {
+
+    JNIEnv* env;
+    jclass peStatClass;
+    if ((*vm)->GetEnv(vm, &env, JNI_VERSION_1_8) != JNI_OK) {
+         return JNI_ERR;
+    }
+
+    FAIL_IF_NULL(peStatClass = (*env)->FindClass(env, "org/broadinstitute/hellbender/utils/bwa/BwaMemPairEndStats"));
+    FAIL_IF_NULL(peStatClass_failedID = (*env)->GetFieldID(env, peStatClass, "failed", "Z"));
+    FAIL_IF_NULL(peStatClass_lowID = (*env)->GetFieldID(env, peStatClass, "low", "I"));
+    FAIL_IF_NULL(peStatClass_highID = (*env)->GetFieldID(env, peStatClass, "high", "I"));
+    FAIL_IF_NULL(peStatClass_averageID = (*env)->GetFieldID(env, peStatClass, "average", "D"));
+    FAIL_IF_NULL(peStatClass_stdID = (*env)->GetFieldID(env, peStatClass, "std", "D"));
+    (*env)->DeleteLocalRef(env, peStatClass);
+
+    return JNI_VERSION_1_8;
+}
+

src/main/c/init.h

@@ -0,0 +1,15 @@
+#include <jni.h>
+
+#ifndef INIT_H
+#define INIT_H
+
+extern jfieldID peStatClass_failedID;
+extern jfieldID peStatClass_lowID;
+extern jfieldID peStatClass_highID;
+extern jfieldID peStatClass_averageID;
+extern jfieldID peStatClass_stdID;
+
+
+jint JNI_OnLoad(JavaVM* vm, void* reserved);
+
+#endif

src/main/c/jnibwa.c

@@ -194,7 +194,7 @@ void* jnibwa_getRefContigNames( bwaidx_t* pIdx, size_t* pBufSize ) {
 	return bufMem;
 }
 
-void* jnibwa_createAlignments( bwaidx_t* pIdx, mem_opt_t* pOpts, char* pSeq, size_t* pBufSize) {
+void* jnibwa_createAlignments( bwaidx_t* pIdx, mem_opt_t* pOpts, mem_pestat_t* pPestat, char* pSeq, size_t* pBufSize) {
 	char c = 0;
 	char* emptyString = &c;
 	uint32_t nSeqs = *(uint32_t*)pSeq;
@@ -211,7 +211,7 @@ void* jnibwa_createAlignments( bwaidx_t* pIdx, mem_opt_t* pOpts, char* pSeq, siz
 		pSeq += seqLen + 1;
 	}
 
-	mem_process_seqs(pOpts, pIdx->bwt, pIdx->bns, pIdx->pac, 0, nSeqs, pSeq1Beg, 0);
+	mem_process_seqs(pOpts, pIdx->bwt, pIdx->bns, pIdx->pac, 0, nSeqs, pSeq1Beg, pPestat);
 
 	size_t nInts = 0;
 	for ( pSeq1 = pSeq1Beg; pSeq1 != pSeq1End; ++pSeq1 ) {

src/main/c/jnibwa.h

@@ -13,6 +13,6 @@ int jnibwa_createIndexFile( char const* refName, char const* imgSuffix );
 bwaidx_t* jnibwa_openIndex( int fd );
 int jnibwa_destroyIndex( bwaidx_t* pIdx );
 void* jnibwa_getRefContigNames( bwaidx_t* pIdx, size_t* pBufSize );
-void* jnibwa_createAlignments( bwaidx_t* pIdx, mem_opt_t* pOpts, char* pSeq, size_t* pBufSize);
+void* jnibwa_createAlignments( bwaidx_t* pIdx, mem_opt_t* pOpts, mem_pestat_t* peStats, char* pSeq, size_t* pBufSize);
 
 #endif /* JNIBWA_H_ */

src/main/c/org_broadinstitute_hellbender_utils_bwa_BwaMemIndex.c

@@ -2,6 +2,7 @@
 #include <fcntl.h>
 #include <stdlib.h>
 #include "jnibwa.h"
+#include "init.h"
 #include "bwa/bwa_commit.h"
 
 
@@ -17,6 +18,27 @@ jint throwIllegalArgumentException(JNIEnv* env, char* message) {
    return (*env)->ThrowNew(env, iaeClass, message);
 }
 
+int jobject_to_mem_pestat_t(JNIEnv* env, jobject in, mem_pestat_t *out) {
+   if (in == NULL) {
+     return 0;
+   }
+   memset(out, 0, sizeof(mem_pestat_t) * 4);
+   for (int i = 0; i < 4; i++, out++) {
+      if (i == 1) {
+         out->failed = (int) (*env)->GetBooleanField(env, in, peStatClass_failedID);
+         if (!out->failed) {
+            out->low = (int) (*env)->GetIntField(env, in, peStatClass_lowID);
+            out->high = (int) (*env)->GetIntField(env, in, peStatClass_highID);
+            out->avg = (double) (*env)->GetDoubleField(env, in, peStatClass_averageID);
+            out->std = (double) (*env)->GetDoubleField(env, in, peStatClass_stdID);
+         }
+      } else {
+         out->failed = 1;
+      }
+   }
+   return 1;
+}
+
 JNIEXPORT jboolean JNICALL
 Java_org_broadinstitute_hellbender_utils_bwa_BwaMemIndex_createReferenceIndex( JNIEnv* env, jclass cls, jstring jReferenceFileName, jstring jIndexPrefix, jstring jAlgoName) {
 
@@ -119,12 +141,14 @@ typedef struct {
 */
 JNIEXPORT jobject JNICALL
 Java_org_broadinstitute_hellbender_utils_bwa_BwaMemIndex_createAlignments(
-				JNIEnv* env, jclass cls, jobject seqsBuf, jlong idxAddr, jobject optsBuf ) {
+				JNIEnv* env, jclass cls, jobject seqsBuf, jlong idxAddr, jobject optsBuf, jobject frPEStats ) {
 	bwaidx_t* pIdx = (bwaidx_t*)idxAddr;
 	mem_opt_t* pOpts = (*env)->GetDirectBufferAddress(env, optsBuf);
+	mem_pestat_t peStats[4];
+	int pestatProvided = jobject_to_mem_pestat_t(env, frPEStats, peStats);
 	char* pSeq = (*env)->GetDirectBufferAddress(env, seqsBuf);
 	size_t bufSize = 0;
-	void* bufMem = jnibwa_createAlignments(pIdx, pOpts, pSeq, &bufSize);
+	void* bufMem = jnibwa_createAlignments(pIdx, pOpts, pestatProvided ? peStats : 0, pSeq, &bufSize);
 	jobject alnBuf = (*env)->NewDirectByteBuffer(env, bufMem, bufSize);
 	if ( !alnBuf ) free(bufMem);
 	return alnBuf;

src/main/java/org/broadinstitute/hellbender/utils/bwa/BwaMemAligner.java

@@ -3,6 +3,7 @@
 import java.nio.ByteBuffer;
 import java.nio.ByteOrder;
 import java.util.ArrayList;
+import java.util.Arrays;
 import java.util.List;
 import java.util.function.Function;
 
@@ -20,13 +21,16 @@ public final class BwaMemAligner implements AutoCloseable {
     private final BwaMemIndex index;
     private ByteBuffer opts;
 
+    private BwaMemPairEndStats pairEndStats;
+
     public BwaMemAligner( final BwaMemIndex index ) {
         this.index = index;
         if ( !index.isOpen() ) {
             throw new IllegalStateException("Can't create aligner: bwa-mem index has been closed");
         }
         opts = BwaMemIndex.createDefaultOptions();
         opts.order(ByteOrder.nativeOrder()).position(0).limit(opts.capacity());
+        pairEndStats = null;
     }
 
     public boolean isOpen() { return opts != null; }
@@ -141,6 +145,30 @@ public void setIntraCtgOptions() {
         setClip3PenaltyOption(5);
     }
 
+    /**
+     * Indicate that we want bwa-mem to infer pair-end inter-size stat from the sequences given.
+     * Only has effect in paired alignment.
+     */
+    public void inferPairEndStats() {
+        pairEndStats = null;
+    }
+
+    /**
+     * Tells the aligner to avoid trying to infer inter-size stats and just go with bwamem defaults when
+     * that information is not-available.
+     */
+    public void dontInferPairEndStats() {
+        pairEndStats = BwaMemPairEndStats.DO_NOT_INFER;
+    }
+
+    /**
+     * Indicate the pair-end inter size stats for "properly" oriented read-pairs.
+     * @param stats
+     */
+    public void setProperPairEndStats(final BwaMemPairEndStats stats) {
+        pairEndStats = stats;
+    }
+
     public BwaMemIndex getIndex() {
         return index;
     }
@@ -179,7 +207,7 @@ public <T> List<List<BwaMemAlignment>> alignSeqs( final Iterable<T> iterable, fi
                 contigBuf.put(func.apply(ele)).put((byte) 0);
             }
             contigBuf.flip();
-            alignsBuf = index.doAlignment(contigBuf, tmpOpts);
+            alignsBuf = index.doAlignment(contigBuf, tmpOpts, pairEndStats);
         }
         finally {
             index.deRefIndex();

src/main/java/org/broadinstitute/hellbender/utils/bwa/BwaMemIndex.java

@@ -407,8 +407,8 @@ public static String getBWAVersion() {
         return getVersion();
     }
 
-    ByteBuffer doAlignment( final ByteBuffer seqs, final ByteBuffer opts ) {
-        final ByteBuffer alignments = createAlignments(seqs, indexAddress, opts);
+    ByteBuffer doAlignment( final ByteBuffer seqs, final ByteBuffer opts, final BwaMemPairEndStats peStats) {
+        final ByteBuffer alignments = createAlignments(seqs, indexAddress, opts, peStats);
         if ( alignments == null ) {
             throw new IllegalStateException("Unable to get alignments from bwa-mem index "+indexImageFile+": We don't know why.");
         }
@@ -482,7 +482,7 @@ private static void loadNativeLibrary() {
     private static native int destroyIndex( long indexAddress );
     static native ByteBuffer createDefaultOptions();
     private static native ByteBuffer getRefContigNames( long indexAddress );
-    private static native ByteBuffer createAlignments( ByteBuffer seqs, long indexAddress, ByteBuffer opts );
+    private static native ByteBuffer createAlignments( ByteBuffer seqs, long indexAddress, ByteBuffer opts, BwaMemPairEndStats peStats);
     static native void destroyByteBuffer( ByteBuffer alignments );
     private static native String getVersion();
 }