zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Temporary commit before moving to new native library

This commit is contained in:
Karthik Gururaj 2014-01-24 16:29:35 -08:00
parent 936e9e175e
commit 81bdfbd00d
15 changed files with 902 additions and 529 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.gitignore vendored
View File

@ -36,4 +36,5 @@ org/
package-list
resources/
velocity.log
perf
verify

7
PairHMM_JNI/LoadTimeInitializer.cc
View File

@ -13,6 +13,10 @@ LoadTimeInitializer::LoadTimeInitializer() //will be called when library is loa
m_sumSquareNumHaplotypes = 0;
m_sumNumTestcases = 0;
m_sumSquareNumTestcases = 0;
m_sumReadLengths = 0;
m_sumHaplotypeLengths = 0;
m_sumProductReadLengthHaplotypeLength = 0;
m_sumSquareProductReadLengthHaplotypeLength = 0;
m_maxNumTestcases = 0;
m_num_invocations = 0;
m_compute_time = 0;
@ -40,6 +44,9 @@ void LoadTimeInitializer::print_profiling()
mean_val = m_sumNumTestcases/m_num_invocations;
cout << "numtestcases\t"<<m_sumNumTestcases<<"\t"<<m_sumSquareNumTestcases<<"\t"<<mean_val<<"\t"<<
(m_sumSquareNumTestcases/m_num_invocations)-mean_val*mean_val<<"\t"<<m_maxNumTestcases<<"\n";
mean_val = m_sumProductReadLengthHaplotypeLength/m_sumNumTestcases;
cout <<"productReadLengthHaplotypeLength\t"<<m_sumProductReadLengthHaplotypeLength<<"\t"<<m_sumSquareProductReadLengthHaplotypeLength<<"\t"
<<mean_val<<"\t"<<(m_sumSquareProductReadLengthHaplotypeLength/m_sumNumTestcases)-mean_val*mean_val<<"\n";
cout.flush();
}

4
PairHMM_JNI/LoadTimeInitializer.h
View File

@ -23,6 +23,10 @@ class LoadTimeInitializer
double m_sumSquareNumHaplotypes;
double m_sumNumTestcases;
double m_sumSquareNumTestcases;
double m_sumReadLengths;
double m_sumHaplotypeLengths;
double m_sumProductReadLengthHaplotypeLength;
double m_sumSquareProductReadLengthHaplotypeLength;
unsigned m_maxNumTestcases;
unsigned m_num_invocations;
//timing

71
PairHMM_JNI/org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM.h 100644
View File

@ -0,0 +1,71 @@
/* DO NOT EDIT THIS FILE - it is machine generated */
#include <jni.h>
/* Header for class org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM */
#ifndef _Included_org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM
#define _Included_org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM
#ifdef __cplusplus
extern "C" {
#endif
#undef org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_TRISTATE_CORRECTION
#define org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_TRISTATE_CORRECTION 3.0
#undef org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_matchToMatch
#define org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_matchToMatch 0L
#undef org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_indelToMatch
#define org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_indelToMatch 1L
#undef org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_matchToInsertion
#define org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_matchToInsertion 2L
#undef org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_insertionToInsertion
#define org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_insertionToInsertion 3L
#undef org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_matchToDeletion
#define org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_matchToDeletion 4L
#undef org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_deletionToDeletion
#define org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_deletionToDeletion 5L
#undef org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_debug
#define org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_debug 0L
#undef org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_verify
#define org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_verify 0L
#undef org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_debug0_1
#define org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_debug0_1 0L
#undef org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_debug1
#define org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_debug1 0L
#undef org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_debug2
#define org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_debug2 0L
#undef org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_debug3
#define org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_debug3 0L
/*
* Class: org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM
* Method: jniInitialize
* Signature: (II)V
*/
JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_jniInitialize
(JNIEnv *, jobject, jint, jint);
/*
* Class: org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM
* Method: jniInitializeProbabilities
* Signature: ([[D[B[B[B)V
*/
JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_jniInitializeProbabilities
(JNIEnv *, jclass, jobjectArray, jbyteArray, jbyteArray, jbyteArray);
/*
* Class: org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM
* Method: jniInitializePriorsAndUpdateCells
* Signature: (ZII[B[B[BI)D
*/
JNIEXPORT jdouble JNICALL Java_org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_jniInitializePriorsAndUpdateCells
(JNIEnv *, jobject, jboolean, jint, jint, jbyteArray, jbyteArray, jbyteArray, jint);
/*
* Class: org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM
* Method: jniSubComputeReadLikelihoodGivenHaplotypeLog10
* Signature: (II[B[B[B[B[B[BI)D
*/
JNIEXPORT jdouble JNICALL Java_org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_jniSubComputeReadLikelihoodGivenHaplotypeLog10
(JNIEnv *, jobject, jint, jint, jbyteArray, jbyteArray, jbyteArray, jbyteArray, jbyteArray, jbyteArray, jint);
#ifdef __cplusplus
}
#endif
#endif

10
PairHMM_JNI/org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.cc
View File

@ -176,6 +176,9 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPai
#ifdef DEBUG3
for(unsigned k=0;k<env->GetArrayLength(haplotypeBases);++k)
g_load_time_initializer.debug_dump("haplotype_bases_jni.txt",to_string((int)haplotypeBasesArray[k]),true);
#endif
#ifdef DO_PROFILING
g_load_time_initializer.m_sumHaplotypeLengths += env->GetArrayLength(haplotypeBasesGlobalRef);
#endif
}
}
@ -272,6 +275,10 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPai
tc_array[tc_idx].i = (char*)insertionGOPArray;
tc_array[tc_idx].d = (char*)deletionGOPArray;
tc_array[tc_idx].c = (char*)overallGCPArray;
#ifdef DO_PROFILING
g_load_time_initializer.m_sumProductReadLengthHaplotypeLength += (readLength*haplotypeLength);
g_load_time_initializer.m_sumSquareProductReadLengthHaplotypeLength += ((readLength*haplotypeLength)*(readLength*haplotypeLength));
#endif
++tc_idx;
}
//Release read arrays at end because they are used by compute_full_prob
@ -283,6 +290,9 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPai
readBasesArrayVector[i][2] = make_pair(insertionGOP, insertionGOPArray);
readBasesArrayVector[i][3] = make_pair(deletionGOP, deletionGOPArray);
readBasesArrayVector[i][4] = make_pair(overallGCP, overallGCPArray);
#ifdef DO_PROFILING
g_load_time_initializer.m_sumReadLengths += readLength;
#endif
}
#ifdef DO_PROFILING
g_load_time_initializer.m_data_transfer_time += get_time();

79
PairHMM_JNI/org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM.h 100644
View File

@ -0,0 +1,79 @@
/* DO NOT EDIT THIS FILE - it is machine generated */
#include <jni.h>
/* Header for class org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM */
#ifndef _Included_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM
#define _Included_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM
#ifdef __cplusplus
extern "C" {
#endif
#undef org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_TRISTATE_CORRECTION
#define org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_TRISTATE_CORRECTION 3.0
#undef org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_matchToMatch
#define org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_matchToMatch 0L
#undef org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_indelToMatch
#define org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_indelToMatch 1L
#undef org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_matchToInsertion
#define org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_matchToInsertion 2L
#undef org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_insertionToInsertion
#define org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_insertionToInsertion 3L
#undef org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_matchToDeletion
#define org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_matchToDeletion 4L
#undef org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_deletionToDeletion
#define org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_deletionToDeletion 5L
#undef org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_sse42Mask
#define org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_sse42Mask 1LL
#undef org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_avxMask
#define org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_avxMask 2LL
/*
* Class: org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM
* Method: jniGetMachineType
* Signature: ()J
*/
JNIEXPORT jlong JNICALL Java_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_jniGetMachineType
(JNIEnv *, jobject);
/*
* Class: org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM
* Method: jniClose
* Signature: ()V
*/
JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_jniClose
(JNIEnv *, jobject);
/*
* Class: org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM
* Method: jniGlobalInit
* Signature: (Ljava/lang/Class;Ljava/lang/Class;J)V
*/
JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_jniGlobalInit
(JNIEnv *, jobject, jclass, jclass, jlong);
/*
* Class: org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM
* Method: jniInitializeHaplotypes
* Signature: (I[Lorg/broadinstitute/sting/utils/pairhmm/VectorLoglessPairHMM/JNIHaplotypeDataHolderClass;)V
*/
JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_jniInitializeHaplotypes
(JNIEnv *, jobject, jint, jobjectArray);
/*
* Class: org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM
* Method: jniFinalizeRegion
* Signature: ()V
*/
JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_jniFinalizeRegion
(JNIEnv *, jobject);
/*
* Class: org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM
* Method: jniComputeLikelihoods
* Signature: (II[Lorg/broadinstitute/sting/utils/pairhmm/VectorLoglessPairHMM/JNIReadDataHolderClass;[Lorg/broadinstitute/sting/utils/pairhmm/VectorLoglessPairHMM/JNIHaplotypeDataHolderClass;[DI)V
*/
JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_jniComputeLikelihoods
(JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jdoubleArray, jint);
#ifdef __cplusplus
}
#endif
#endif

6
PairHMM_JNI/run.sh
View File

@ -2,17 +2,17 @@
rm -f *.txt *.log
GSA_ROOT_DIR=/home/karthikg/broad/gsa-unstable
#-Djava.library.path is needed if you are using JNI_LOGLESS_CACHING, else not needed
java -Djava.library.path=${GSA_ROOT_DIR}/PairHMM_JNI -jar $GSA_ROOT_DIR/dist/GenomeAnalysisTK.jar -T HaplotypeCaller \
java -Djava.library.path=${GSA_ROOT_DIR}/PairHMM_JNI/ -jar $GSA_ROOT_DIR/dist/GenomeAnalysisTK.jar -T HaplotypeCaller \
-R /data/broad/samples/joint_variant_calling/broad_reference/human_g1k_v37_decoy.fasta \
-I /data/broad/samples/joint_variant_calling/NA12878_low_coverage_alignment/NA12878.chrom11.ILLUMINA.bwa.CEU.low_coverage.20121211.bam \
--dbsnp /data/broad/samples/joint_variant_calling/dbSNP/00-All.vcf \
-stand_call_conf 50.0 \
-stand_emit_conf 10.0 \
--pair_hmm_implementation JNI_LOGLESS_CACHING \
-o output.raw.snps.indels.vcf
-o output.raw.snps.indels.vcf
#-XL unmapped \
#--pair_hmm_implementation JNI_LOGLESS_CACHING \
#-XL unmapped \
#-I /data/simulated/sim1M_pairs_final.bam \
#-I /data/broad/samples/joint_variant_calling/NA12878_low_coverage_alignment/NA12878.chrom11.ILLUMINA.bwa.CEU.low_coverage.20121211.bam \
#-I /data/broad/samples/joint_variant_calling/NA12878_high_coverage_alignment/NA12878.mapped.ILLUMINA.bwa.CEU.high_coverage_pcr_free.20130906.bam \

16
PairHMM_JNI/utils.cc
View File

@ -33,17 +33,27 @@ bool is_sse42_supported()
return ((ecx >> 20)&1) == 1;
}
void initialize_function_pointers()
uint64_t get_machine_capabilities()
{
uint64_t machine_mask = 0ull;
if(is_avx_supported())
machine_mask |= (1 << AVX_CUSTOM_IDX);
if(is_sse42_supported())
machine_mask |= (1 << SSE42_CUSTOM_IDX);
return machine_mask;
}
void initialize_function_pointers(uint64_t mask)
{
//if(false)
if(is_avx_supported())
if(is_avx_supported() && (mask & (1<< AVX_CUSTOM_IDX)))
{
cout << "Using AVX accelerated implementation of PairHMM\n";
g_compute_full_prob_float = compute_full_prob_avxs<float>;
g_compute_full_prob_double = compute_full_prob_avxd<double>;
}
else
if(is_sse42_supported())
if(is_sse42_supported() && (mask & (1<< SSE42_CUSTOM_IDX)))
{
cout << "Using SSE4.2 accelerated implementation of PairHMM\n";
g_compute_full_prob_float = compute_full_prob_sses<float>;

9
PairHMM_JNI/utils.h
View File

@ -25,7 +25,14 @@ extern double (*g_compute_full_prob_double)(testcase *tc, double* before_last_lo
void debug_dump(std::string filename, std::string s, bool to_append, bool add_newline);
template<class NUMBER>
NUMBER compute_full_prob(testcase *tc, NUMBER *before_last_log=0);
void initialize_function_pointers();
void initialize_function_pointers(uint64_t mask=0xFFFFFFFFFFFFFFFFull);
double getCurrClk();
uint64_t get_time(struct timespec* x=0);
enum ProcessorCapabilitiesEnum
{
SSE42_CUSTOM_IDX=0,
AVX_CUSTOM_IDX
};
uint64_t get_machine_capabilities();
#endif

6
protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/PairHMMLikelihoodCalculationEngine.java
View File

@ -91,8 +91,10 @@ public class PairHMMLikelihoodCalculationEngine implements LikelihoodCalculation
return new LoglessPairHMM();
else
return new CnyPairHMM();
case JNI_LOGLESS_CACHING:
return new JNILoglessPairHMM();
case VECTOR_LOGLESS_CACHING:
return new VectorLoglessPairHMM();
case DEBUG_JNI_LOGLESS_CACHING:
return new DebugJNILoglessPairHMM(hmmType);
case ARRAY_LOGLESS:
if (noFpga || !CnyPairHMM.isAvailable())
return new ArrayLoglessPairHMM();

496
protected/java/src/org/broadinstitute/sting/utils/pairhmm/DebugJNILoglessPairHMM.java 100644
View File

@ -0,0 +1,496 @@
/*
* By downloading the PROGRAM you agree to the following terms of use:
*
* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
*
* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
*
* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
*
* 1. DEFINITIONS
* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
*
* 2. LICENSE
* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
*
* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
* Copyright 2012 Broad Institute, Inc.
* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
*
* 4. INDEMNIFICATION
* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
*
* 5. NO REPRESENTATIONS OR WARRANTIES
* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
*
* 6. ASSIGNMENT
* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
*
* 7. MISCELLANEOUS
* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
*/
package org.broadinstitute.sting.utils.pairhmm;
import com.google.java.contract.Ensures;
import com.google.java.contract.Requires;
import org.broadinstitute.sting.utils.QualityUtils;
import org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap;
import org.broadinstitute.sting.utils.haplotype.Haplotype;
import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
import org.broadinstitute.variant.variantcontext.Allele;
import org.broadinstitute.sting.utils.exceptions.UserException;
import java.util.List;
import java.util.Map;
import java.util.HashMap;
import java.io.File;
import java.io.FileWriter;
import java.io.BufferedWriter;
import java.util.Map;
import java.util.HashMap;
import java.io.IOException;
/**
* Created with IntelliJ IDEA.
* User: rpoplin, carneiro
* Date: 10/16/12
*/
public class DebugJNILoglessPairHMM extends LoglessPairHMM {
private static final boolean debug = false; //simulates ifdef
private static final boolean verify = debug || false; //simulates ifdef
private static final boolean debug0_1 = false; //simulates ifdef
private static final boolean debug1 = false; //simulates ifdef
private static final boolean debug2 = false;
private static final boolean debug3 = false;
//Debugging stats
private int numCalls = 0;
private int numComputeLikelihoodCalls = 0;
protected HashMap<String, BufferedWriter> filenameToWriter = new HashMap<String, BufferedWriter>();
private JNILoglessPairHMM jniPairHMM = null;
public DebugJNILoglessPairHMM(final PairHMM.HMM_IMPLEMENTATION hmmType) {
super();
switch(hmmType) {
case VECTOR_LOGLESS_CACHING:
jniPairHMM = new VectorLoglessPairHMM();
break;
default:
throw new UserException.BadArgumentValue("pairHMM","Specified JNIPairHMM implementation is unrecognized or incompatible with the HaplotypeCaller. Acceptable options are VECTOR_LOGLESS_CACHING");
}
}
@Override
public void close()
{
jniPairHMM.close();
debugClose();
}
//Used only when testing parts of the compute kernel
/**
* {@inheritDoc}
*/
@Override
public void initialize( final int readMaxLength, final int haplotypeMaxLength ) {
if(verify)
super.initialize(readMaxLength, haplotypeMaxLength);
if(debug3)
{
System.out.println("Java: alloc initialized readMaxLength : "+readMaxLength+" haplotypeMaxLength : "+haplotypeMaxLength);
debugDump("lengths_java.txt", String.format("%d %d\n",readMaxLength, haplotypeMaxLength),
true);
}
if(debug2)
jniInitialize(readMaxLength, haplotypeMaxLength);
}
private HashMap<Haplotype,Integer> haplotypeToHaplotypeListIdxMap = null;
//Used to transfer data to JNI
//Since the haplotypes are the same for all calls to computeLikelihoods within a region, transfer the haplotypes only once to the JNI per region
/**
* {@inheritDoc}
*/
@Override
public void initialize( final List<Haplotype> haplotypes, final Map<String, List<GATKSAMRecord>> perSampleReadList,
final int readMaxLength, final int haplotypeMaxLength ) {
if(verify)
super.initialize(haplotypes, perSampleReadList, readMaxLength, haplotypeMaxLength);
jniPairHMM.initialize(haplotypes, perSampleReadList, readMaxLength, haplotypeMaxLength);
haplotypeToHaplotypeListIdxMap = jniPairHMM.getHaplotypeToHaplotypeListIdxMap();
}
/**
* {@inheritDoc}
*/
@Override
public void finalizeRegion()
{
jniPairHMM.finalizeRegion();
}
/**
* {@inheritDoc}
*/
@Override
public PerReadAlleleLikelihoodMap computeLikelihoods( final List<GATKSAMRecord> reads, final Map<Allele, Haplotype> alleleHaplotypeMap, final Map<GATKSAMRecord, byte[]> GCPArrayMap ) {
// (re)initialize the pairHMM only if necessary
final int readMaxLength = verify ? findMaxReadLength(reads) : 0;
final int haplotypeMaxLength = verify ? findMaxHaplotypeLength(alleleHaplotypeMap) : 0;
if(verify)
{
if (!initialized || readMaxLength > maxReadLength || haplotypeMaxLength > maxHaplotypeLength)
{ initialize(readMaxLength, haplotypeMaxLength); }
if ( ! initialized )
throw new IllegalStateException("Must call initialize before calling jniComputeLikelihoods in debug/verify mode");
}
int readListSize = reads.size();
int numHaplotypes = alleleHaplotypeMap.size();
int numTestcases = readListSize*numHaplotypes;
if(debug0_1)
System.out.println("Java numReads "+readListSize+" numHaplotypes "+numHaplotypes);
int idx = 0;
for(GATKSAMRecord read : reads)
{
byte [] overallGCP = GCPArrayMap.get(read);
if(debug0_1)
System.out.println("Java read length "+read.getReadBases().length);
if(debug3)
{
for(int i=0;i<read.getReadBases().length;++i)
{
debugDump("reads_java.txt",String.format("%d\n",(int)read.getReadBases()[i]),true);
debugDump("reads_java.txt",String.format("%d\n",(int)read.getBaseQualities()[i]),true);
debugDump("reads_java.txt",String.format("%d\n",(int)read.getBaseInsertionQualities()[i]),true);
debugDump("reads_java.txt",String.format("%d\n",(int)read.getBaseDeletionQualities()[i]),true);
debugDump("reads_java.txt",String.format("%d\n",(int)overallGCP[i]),true);
}
}
++idx;
}
if(verify)
{
idx = 0;
for (Map.Entry<Allele, Haplotype> currEntry : alleleHaplotypeMap.entrySet()) //order is important - access in same order always
{
byte[] haplotypeBases = currEntry.getValue().getBases();
if(debug0_1)
System.out.println("Java haplotype length "+haplotypeBases.length);
if(debug3)
{
for(int i=0;i<haplotypeBases.length;++i)
debugDump("haplotype_bases_java.txt",String.format("%d\n",(int)haplotypeBases[i]),true);
}
++idx;
}
}
jniPairHMM.computeLikelihoods(reads, alleleHaplotypeMap, GCPArrayMap);
double[] likelihoodArray = jniPairHMM.getLikelihoodArray();
//to compare values
final PerReadAlleleLikelihoodMap likelihoodMap = super.computeLikelihoods(reads, alleleHaplotypeMap, GCPArrayMap);
if(verify)
{
//re-order values in likelihoodArray
double[] tmpArray = new double[numHaplotypes];
idx = 0;
int idxInsideHaplotypeList = 0;
int readIdx = 0;
for(GATKSAMRecord read : reads)
{
for(int j=0;j<numHaplotypes;++j)
tmpArray[j] = likelihoodArray[readIdx+j];
for (Map.Entry<Allele, Haplotype> currEntry : alleleHaplotypeMap.entrySet())//order is important - access in same order always
{
idxInsideHaplotypeList = haplotypeToHaplotypeListIdxMap.get(currEntry.getValue());
likelihoodArray[idx] = tmpArray[idxInsideHaplotypeList];
++idx;
}
readIdx += numHaplotypes;
}
//for floating point values, no exact equality
//check whether numbers are close in terms of abs_error or relative_error
//For very large values, relative_error is relevant
//For very small values, abs_error is relevant
boolean toDump = false;
for(int i=0;i<likelihoodArray.length;++i)
{
double abs_error = Math.abs(likelihoodArray[i] - mLikelihoodArray[i]);
double relative_error = 0;
if(mLikelihoodArray[i] == 0)
relative_error = 0;
else
relative_error = Math.abs(abs_error/mLikelihoodArray[i]);
if(abs_error > 1e-5 && relative_error > 1e-5)
{
toDump = true;
break;
}
}
//if numbers are not close, then dump out the data that produced the inconsistency
if(toDump)
{
idx = 0;
System.out.println("Dump : Java numReads "+readListSize+" numHaplotypes "+numHaplotypes);
for(GATKSAMRecord read : reads)
{
byte [] overallGCP = GCPArrayMap.get(read);
for (Map.Entry<Allele, Haplotype> currEntry : alleleHaplotypeMap.entrySet()) //order is important - access in same order always
{
byte[] haplotypeBases = currEntry.getValue().getBases();
debugDump("debug_dump.txt",new String(haplotypeBases)+" ",true);
debugDump("debug_dump.txt",new String(read.getReadBases())+" ",true);
for(int k=0;k<read.getReadBases().length;++k)
debugDump("debug_dump.txt",String.format("%d",(int)(read.getBaseQualities()[k]))+" ",true);
for(int k=0;k<read.getReadBases().length;++k)
debugDump("debug_dump.txt",String.format("%d",(int)(read.getBaseInsertionQualities()[k]))+" ",true);
for(int k=0;k<read.getReadBases().length;++k)
debugDump("debug_dump.txt",String.format("%d",(int)(read.getBaseDeletionQualities()[k]))+" ",true);
for(int k=0;k<read.getReadBases().length;++k)
debugDump("debug_dump.txt",String.format("%d",(int)(overallGCP[k]))+" ",true);
debugDump("debug_dump.txt","\n",true);
debugDump("debug_results.txt",String.format("%e %e\n",mLikelihoodArray[idx],likelihoodArray[idx]),true);
++idx;
}
}
}
debugClose();
}
++numComputeLikelihoodCalls;
//if(numComputeLikelihoodCalls == 5)
//jniPairHMM.jniClose();
//System.exit(0);
return likelihoodMap;
}
//Used to test parts of the compute kernel separately
private native void jniInitialize( final int readMaxLength, final int haplotypeMaxLength);
private native static void jniInitializeProbabilities( final double[][] transition, final byte[] insertionGOP,
final byte[] deletionGOP, final byte[] overallGCP);
private native double jniInitializePriorsAndUpdateCells( boolean doInitialization, final int paddedReadLength,
final int paddedHaplotypeLength, final byte[] readBases, final byte[] haplotypeBases, final byte[] readQuals,
final int hapStartIndex);
private native double jniSubComputeReadLikelihoodGivenHaplotypeLog10( final int readLength, final int haplotypeLength,
final byte[] readBases, final byte[] haplotypeBases, final byte[] readQuals, final byte[] insertionGOP,
final byte[] deletionGOP, final byte[] overallGCP, final int hapStartIndex);
/**
* {@inheritDoc}
*/
@Override
public double subComputeReadLikelihoodGivenHaplotypeLog10( final byte[] haplotypeBases, final byte[] readBases,
final byte[] readQuals, final byte[] insertionGOP, final byte[] deletionGOP, final byte[] overallGCP,
final int hapStartIndex, final boolean recacheReadValues, final int nextHapStartIndex) {
//System.out.println("#### START STACK TRACE ####");
//for (StackTraceElement ste : Thread.currentThread().getStackTrace()) {
//System.out.println(ste);
//}
//System.out.println("#### END STACK TRACE ####");
//
if(debug1)
jniSubComputeReadLikelihoodGivenHaplotypeLog10(readBases.length, haplotypeBases.length,
readBases, haplotypeBases, readQuals,
insertionGOP, deletionGOP, overallGCP,
hapStartIndex);
boolean doInitialization = (previousHaplotypeBases == null || previousHaplotypeBases.length != haplotypeBases.length);
if (doInitialization) {
final double initialValue = INITIAL_CONDITION / haplotypeBases.length;
// set the initial value (free deletions in the beginning) for the first row in the deletion matrix
for( int j = 0; j < paddedHaplotypeLength; j++ ) {
deletionMatrix[0][j] = initialValue;
}
}
if ( ! constantsAreInitialized || recacheReadValues ) {
initializeProbabilities(transition, insertionGOP, deletionGOP, overallGCP);
if(debug3)
{
System.out.println("Java: initializeProbabilities lengths : "+insertionGOP.length+" padded "+paddedReadLength+" "+paddedHaplotypeLength);
for(int i=0;i<insertionGOP.length;++i)
for(int j=0;j<6;++j)
debugDump("transitions_java.txt",String.format("%e\n",transition[i+1][j]),true);
}
if(debug2)
jniInitializeProbabilities(transition, insertionGOP, deletionGOP, overallGCP);
// note that we initialized the constants
constantsAreInitialized = true;
}
if(debug3)
System.out.println("Java: initializePriors : lengths "+readBases.length+" "+haplotypeBases.length+" padded "+paddedReadLength+" "+paddedHaplotypeLength + " doNotUseTristateCorrection "+doNotUseTristateCorrection);
initializePriors(haplotypeBases, readBases, readQuals, hapStartIndex);
for (int i = 1; i < paddedReadLength; i++) {
// +1 here is because hapStartIndex is 0-based, but our matrices are 1 based
for (int j = hapStartIndex+1; j < paddedHaplotypeLength; j++) {
updateCell(i, j, prior[i][j], transition[i]);
}
}
// final probability is the log10 sum of the last element in the Match and Insertion state arrays
// this way we ignore all paths that ended in deletions! (huge)
// but we have to sum all the paths ending in the M and I matrices, because they're no longer extended.
final int endI = paddedReadLength - 1;
double finalSumProbabilities = 0.0;
for (int j = 1; j < paddedHaplotypeLength; j++) {
finalSumProbabilities += matchMatrix[endI][j] + insertionMatrix[endI][j];
}
if(debug2)
jniInitializePriorsAndUpdateCells(doInitialization, paddedReadLength, paddedHaplotypeLength,
readBases, haplotypeBases, readQuals,
hapStartIndex);
if(debug)
debugDump("return_values_java.txt",String.format("%e\n",Math.log10(finalSumProbabilities) - INITIAL_CONDITION_LOG10),true);
++numCalls;
//if(numCalls == 100)
//{
//debugClose();
//System.exit(0);
//}
return Math.log10(finalSumProbabilities) - INITIAL_CONDITION_LOG10;
}
/**
* Initializes the matrix that holds all the constants related to the editing
* distance between the read and the haplotype.
*
* @param haplotypeBases the bases of the haplotype
* @param readBases the bases of the read
* @param readQuals the base quality scores of the read
* @param startIndex where to start updating the distanceMatrix (in case this read is similar to the previous read)
*/
protected void initializePriors(final byte[] haplotypeBases, final byte[] readBases, final byte[] readQuals, final int startIndex) {
// initialize the pBaseReadLog10 matrix for all combinations of read x haplotype bases
// Abusing the fact that java initializes arrays with 0.0, so no need to fill in rows and columns below 2.
if(debug3)
System.out.println("hapStartIndex "+startIndex);
for (int i = 0; i < readBases.length; i++) {
final byte x = readBases[i];
final byte qual = readQuals[i];
for (int j = startIndex; j < haplotypeBases.length; j++) {
final byte y = haplotypeBases[j];
prior[i+1][j+1] = ( x == y || x == (byte) 'N' || y == (byte) 'N' ?
QualityUtils.qualToProb(qual) : (QualityUtils.qualToErrorProb(qual) / (doNotUseTristateCorrection ? 1.0 : TRISTATE_CORRECTION)) );
if(debug3)
debugDump("priors_java.txt",String.format("%e\n",prior[i+1][j+1]),true);
}
}
}
/**
* Initializes the matrix that holds all the constants related to quality scores.
*
* @param insertionGOP insertion quality scores of the read
* @param deletionGOP deletion quality scores of the read
* @param overallGCP overall gap continuation penalty
*/
@Requires({
"insertionGOP != null",
"deletionGOP != null",
"overallGCP != null"
})
@Ensures("constantsAreInitialized")
protected static void initializeProbabilities(final double[][] transition, final byte[] insertionGOP, final byte[] deletionGOP, final byte[] overallGCP) {
for (int i = 0; i < insertionGOP.length; i++) {
final int qualIndexGOP = Math.min(insertionGOP[i] + deletionGOP[i], Byte.MAX_VALUE);
transition[i+1][matchToMatch] = QualityUtils.qualToProb((byte) qualIndexGOP);
transition[i+1][indelToMatch] = QualityUtils.qualToProb(overallGCP[i]);
transition[i+1][matchToInsertion] = QualityUtils.qualToErrorProb(insertionGOP[i]);
transition[i+1][insertionToInsertion] = QualityUtils.qualToErrorProb(overallGCP[i]);
transition[i+1][matchToDeletion] = QualityUtils.qualToErrorProb(deletionGOP[i]);
transition[i+1][deletionToDeletion] = QualityUtils.qualToErrorProb(overallGCP[i]);
//TODO it seems that it is not always the case that matchToMatch + matchToDeletion + matchToInsertion == 1.
//TODO We have detected cases of 1.00002 which can cause problems downstream. This are typically masked
//TODO by the fact that we always add a indelToMatch penalty to all PairHMM likelihoods (~ -0.1)
//TODO This is in fact not well justified and although it does not have any effect (since is equally added to all
//TODO haplotypes likelihoods) perhaps we should just remove it eventually and fix this != 1.0 issue here.
}
}
/**
* Updates a cell in the HMM matrix
*
* The read and haplotype indices are offset by one because the state arrays have an extra column to hold the
* initial conditions
* @param indI row index in the matrices to update
* @param indJ column index in the matrices to update
* @param prior the likelihood editing distance matrix for the read x haplotype
* @param transition an array with the six transition relevant to this location
*/
protected void updateCell( final int indI, final int indJ, final double prior, final double[] transition) {
matchMatrix[indI][indJ] = prior * ( matchMatrix[indI - 1][indJ - 1] * transition[matchToMatch] +
insertionMatrix[indI - 1][indJ - 1] * transition[indelToMatch] +
deletionMatrix[indI - 1][indJ - 1] * transition[indelToMatch] );
insertionMatrix[indI][indJ] = matchMatrix[indI - 1][indJ] * transition[matchToInsertion] + insertionMatrix[indI - 1][indJ] * transition[insertionToInsertion];
deletionMatrix[indI][indJ] = matchMatrix[indI][indJ - 1] * transition[matchToDeletion] + deletionMatrix[indI][indJ - 1] * transition[deletionToDeletion];
if(debug3)
{
debugDump("matrices_java.txt",String.format("%e\n",matchMatrix[indI][indJ]),true);
debugDump("matrices_java.txt",String.format("%e\n",insertionMatrix[indI][indJ]),true);
debugDump("matrices_java.txt",String.format("%e\n",deletionMatrix[indI][indJ]),true);
}
}
protected void debugDump( String filename, String s, boolean toAppend ) {
try {
File file = new File(filename);
if (!file.exists())
file.createNewFile();
BufferedWriter currWriter = filenameToWriter.get(filename);
if(currWriter == null)
{
FileWriter fw = new FileWriter(file, toAppend);
currWriter = new BufferedWriter(fw);
filenameToWriter.put(filename, currWriter);
}
currWriter.write(s);
}
catch(IOException e)
{
e.printStackTrace();
}
}
protected void debugClose() {
for(Map.Entry<String, BufferedWriter> currEntry : filenameToWriter.entrySet()) {
BufferedWriter currWriter = currEntry.getValue();
try
{
currWriter.flush();
currWriter.close();
}
catch(IOException e)
{
e.printStackTrace();
}
}
filenameToWriter.clear();
}
}

471
protected/java/src/org/broadinstitute/sting/utils/pairhmm/JNILoglessPairHMM.java
View File

@ -46,17 +46,8 @@
package org.broadinstitute.sting.utils.pairhmm;
import com.google.java.contract.Ensures;
import com.google.java.contract.Requires;
import org.broadinstitute.sting.utils.QualityUtils;
import org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap;
import org.broadinstitute.sting.utils.haplotype.Haplotype;
import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
import org.broadinstitute.variant.variantcontext.Allele;
import java.util.List;
import java.util.Map;
import java.util.HashMap;
@ -65,464 +56,6 @@ import java.util.HashMap;
* User: rpoplin, carneiro
* Date: 10/16/12
*/
public class JNILoglessPairHMM extends LoglessPairHMM {
private static final boolean debug = false; //simulates ifdef
private static final boolean verify = debug || false; //simulates ifdef
private static final boolean debug0_1 = false; //simulates ifdef
private static final boolean debug1 = false; //simulates ifdef
private static final boolean debug2 = false;
private static final boolean debug3 = false;
private int numComputeLikelihoodCalls = 0;
//Used to copy references to byteArrays to JNI from reads
protected class JNIReadDataHolderClass
{
public byte[] readBases = null;
public byte[] readQuals = null;
public byte[] insertionGOP = null;
public byte[] deletionGOP = null;
public byte[] overallGCP = null;
}
//Used to copy references to byteArrays to JNI from haplotypes
protected class JNIHaplotypeDataHolderClass
{
public byte[] haplotypeBases = null;
}
private native void jniGlobalInit(Class<?> readDataHolderClass, Class<?> haplotypeDataHolderClass);
private native void jniClose();
private static boolean isJNILoglessPairHMMLibraryLoaded = false;
public JNILoglessPairHMM()
{
super();
if(!isJNILoglessPairHMMLibraryLoaded)
{
System.loadLibrary("JNILoglessPairHMM");
isJNILoglessPairHMMLibraryLoaded = true;
jniGlobalInit(JNIReadDataHolderClass.class, JNIHaplotypeDataHolderClass.class); //need to do this only once
}
}
@Override
public void close()
{
jniClose();
debugClose();
}
//Used to test parts of the compute kernel separately
private native void jniInitialize(final int readMaxLength, final int haplotypeMaxLength);
private native static void jniInitializeProbabilities(final double[][] transition, final byte[] insertionGOP, final byte[] deletionGOP, final byte[] overallGCP);
private native double jniInitializePriorsAndUpdateCells(
boolean doInitialization,
final int paddedReadLength, final int paddedHaplotypeLength,
final byte[] readBases, final byte[] haplotypeBases, final byte[] readQuals,
final int hapStartIndex);
private native double jniSubComputeReadLikelihoodGivenHaplotypeLog10(
final int readLength, final int haplotypeLength,
final byte[] readBases, final byte[] haplotypeBases, final byte[] readQuals,
final byte[] insertionGOP, final byte[] deletionGOP, final byte[] overallGCP,
final int hapStartIndex);
//Used only when testing parts of the compute kernel
/**
* {@inheritDoc}
*/
@Override
public void initialize(final int readMaxLength, final int haplotypeMaxLength)
{
if(verify)
super.initialize(readMaxLength, haplotypeMaxLength);
if(debug3)
{
System.out.println("Java: alloc initialized readMaxLength : "+readMaxLength+" haplotypeMaxLength : "+haplotypeMaxLength);
debugDump("lengths_java.txt", String.format("%d %d\n",readMaxLength, haplotypeMaxLength),
true);
}
if(debug2)
jniInitialize(readMaxLength, haplotypeMaxLength);
}
private HashMap<Haplotype,Integer> haplotypeToHaplotypeListIdxMap = new HashMap<Haplotype,Integer>();
private native void jniInitializeHaplotypes(final int numHaplotypes, JNIHaplotypeDataHolderClass[] haplotypeDataArray);
//Used to transfer data to JNI
//Since the haplotypes are the same for all calls to computeLikelihoods within a region, transfer the haplotypes only once to the JNI per region
/**
* {@inheritDoc}
*/
@Override
public void initialize( final List<Haplotype> haplotypes, final Map<String, List<GATKSAMRecord>> perSampleReadList, final int readMaxLength, final int haplotypeMaxLength ) {
if(verify)
super.initialize(haplotypes, perSampleReadList, readMaxLength, haplotypeMaxLength);
int numHaplotypes = haplotypes.size();
JNIHaplotypeDataHolderClass[] haplotypeDataArray = new JNIHaplotypeDataHolderClass[numHaplotypes];
int idx = 0;
for(final Haplotype currHaplotype : haplotypes)
{
haplotypeDataArray[idx] = new JNIHaplotypeDataHolderClass();
haplotypeDataArray[idx].haplotypeBases = currHaplotype.getBases();
haplotypeToHaplotypeListIdxMap.put(currHaplotype, idx);
++idx;
}
jniInitializeHaplotypes(numHaplotypes, haplotypeDataArray);
}
private native void jniFinalizeRegion();
//Tell JNI to release arrays - really important if native code is directly accessing Java memory, if not
//accessing Java memory directly, still important to release memory from C++
/**
* {@inheritDoc}
*/
@Override
public void finalizeRegion()
{
jniFinalizeRegion();
}
//Real compute kernel
private native void jniComputeLikelihoods(int numReads, int numHaplotypes,
JNIReadDataHolderClass[] readDataArray, JNIHaplotypeDataHolderClass[] haplotypeDataArray,
double[] likelihoodArray, int maxNumThreadsToUse);
/**
* {@inheritDoc}
*/
@Override
public PerReadAlleleLikelihoodMap computeLikelihoods(final List<GATKSAMRecord> reads, final Map<Allele, Haplotype> alleleHaplotypeMap, final Map<GATKSAMRecord, byte[]> GCPArrayMap)
{
// (re)initialize the pairHMM only if necessary
final int readMaxLength = verify ? findMaxReadLength(reads) : 0;
final int haplotypeMaxLength = verify ? findMaxHaplotypeLength(alleleHaplotypeMap) : 0;
if(verify)
{
if (!initialized || readMaxLength > maxReadLength || haplotypeMaxLength > maxHaplotypeLength)
{ initialize(readMaxLength, haplotypeMaxLength); }
if ( ! initialized )
throw new IllegalStateException("Must call initialize before calling jniComputeLikelihoods in debug/verify mode");
}
int readListSize = reads.size();
int numHaplotypes = alleleHaplotypeMap.size();
int numTestcases = readListSize*numHaplotypes;
if(debug0_1)
System.out.println("Java numReads "+readListSize+" numHaplotypes "+numHaplotypes);
JNIReadDataHolderClass[] readDataArray = new JNIReadDataHolderClass[readListSize];
int idx = 0;
for(GATKSAMRecord read : reads)
{
readDataArray[idx] = new JNIReadDataHolderClass();
readDataArray[idx].readBases = read.getReadBases();
readDataArray[idx].readQuals = read.getBaseQualities();
readDataArray[idx].insertionGOP = read.getBaseInsertionQualities();
readDataArray[idx].deletionGOP = read.getBaseDeletionQualities();
readDataArray[idx].overallGCP = GCPArrayMap.get(read);
if(debug0_1)
System.out.println("Java read length "+readDataArray[idx].readBases.length);
if(debug3)
{
for(int i=0;i<readDataArray[idx].readBases.length;++i)
{
debugDump("reads_java.txt",String.format("%d\n",(int)readDataArray[idx].readBases[i]),true);
debugDump("reads_java.txt",String.format("%d\n",(int)readDataArray[idx].readQuals[i]),true);
debugDump("reads_java.txt",String.format("%d\n",(int)readDataArray[idx].insertionGOP[i]),true);
debugDump("reads_java.txt",String.format("%d\n",(int)readDataArray[idx].deletionGOP[i]),true);
debugDump("reads_java.txt",String.format("%d\n",(int)readDataArray[idx].overallGCP[i]),true);
}
}
++idx;
}
JNIHaplotypeDataHolderClass[] haplotypeDataArray = new JNIHaplotypeDataHolderClass[numHaplotypes];
if(verify)
{
idx = 0;
for (Map.Entry<Allele, Haplotype> currEntry : alleleHaplotypeMap.entrySet()) //order is important - access in same order always
{
haplotypeDataArray[idx] = new JNIHaplotypeDataHolderClass();
haplotypeDataArray[idx].haplotypeBases = currEntry.getValue().getBases();
if(debug0_1)
System.out.println("Java haplotype length "+haplotypeDataArray[idx].haplotypeBases.length);
if(debug3)
{
for(int i=0;i<haplotypeDataArray[idx].haplotypeBases.length;++i)
debugDump("haplotype_bases_java.txt",String.format("%d\n",(int)haplotypeDataArray[idx].haplotypeBases[i]),true);
}
++idx;
}
}
double[] likelihoodArray = new double[readListSize*numHaplotypes]; //to store results
//for(reads)
// for(haplotypes)
// compute_full_prob()
jniComputeLikelihoods(readListSize, numHaplotypes, readDataArray, haplotypeDataArray, likelihoodArray, 12);
//final PerReadAlleleLikelihoodMap likelihoodMap = new PerReadAlleleLikelihoodMap();
PerReadAlleleLikelihoodMap likelihoodMap = new PerReadAlleleLikelihoodMap();
idx = 0;
int idxInsideHaplotypeList = 0;
int readIdx = 0;
for(GATKSAMRecord read : reads)
{
for (Map.Entry<Allele, Haplotype> currEntry : alleleHaplotypeMap.entrySet())//order is important - access in same order always
{
//Since the order of haplotypes in the List<Haplotype> and alleleHaplotypeMap is different,
//get idx of current haplotype in the list and use this idx to get the right likelihoodValue
idxInsideHaplotypeList = haplotypeToHaplotypeListIdxMap.get(currEntry.getValue());
likelihoodMap.add(read, currEntry.getKey(), likelihoodArray[readIdx + idxInsideHaplotypeList]);
++idx;
}
readIdx += numHaplotypes;
}
if(verify)
{
//re-order values in likelihoodArray
double[] tmpArray = new double[numHaplotypes];
idx = 0;
idxInsideHaplotypeList = 0;
readIdx = 0;
for(GATKSAMRecord read : reads)
{
for(int j=0;j<numHaplotypes;++j)
tmpArray[j] = likelihoodArray[readIdx+j];
for (Map.Entry<Allele, Haplotype> currEntry : alleleHaplotypeMap.entrySet())//order is important - access in same order always
{
idxInsideHaplotypeList = haplotypeToHaplotypeListIdxMap.get(currEntry.getValue());
likelihoodArray[idx] = tmpArray[idxInsideHaplotypeList];
++idx;
}
readIdx += numHaplotypes;
}
//to compare values
likelihoodMap = super.computeLikelihoods(reads, alleleHaplotypeMap, GCPArrayMap);
//for floating point values, no exact equality
//check whether numbers are close in terms of abs_error or relative_error
//For very large values, relative_error is relevant
//For very small values, abs_error is relevant
boolean toDump = false;
for(int i=0;i<likelihoodArray.length;++i)
{
double abs_error = Math.abs(likelihoodArray[i] - mLikelihoodArray[i]);
double relative_error = 0;
if(mLikelihoodArray[i] == 0)
relative_error = 0;
else
relative_error = Math.abs(abs_error/mLikelihoodArray[i]);
if(abs_error > 1e-5 && relative_error > 1e-5)
{
toDump = true;
break;
}
}
//if numbers are not close, then dump out the data that produced the inconsistency
if(toDump)
{
idx = 0;
System.out.println("Dump : Java numReads "+readListSize+" numHaplotypes "+numHaplotypes);
for(int i=0;i<readDataArray.length;++i)
{
for(int j=0;j<haplotypeDataArray.length;++j)
{
debugDump("debug_dump.txt",new String(haplotypeDataArray[j].haplotypeBases)+" ",true);
debugDump("debug_dump.txt",new String(readDataArray[i].readBases)+" ",true);
for(int k=0;k<readDataArray[i].readBases.length;++k)
debugDump("debug_dump.txt",String.format("%d",(int)(readDataArray[i].readQuals[k]))+" ",true);
for(int k=0;k<readDataArray[i].readBases.length;++k)
debugDump("debug_dump.txt",String.format("%d",(int)(readDataArray[i].insertionGOP[k]))+" ",true);
for(int k=0;k<readDataArray[i].readBases.length;++k)
debugDump("debug_dump.txt",String.format("%d",(int)(readDataArray[i].deletionGOP[k]))+" ",true);
for(int k=0;k<readDataArray[i].readBases.length;++k)
debugDump("debug_dump.txt",String.format("%d",(int)(readDataArray[i].overallGCP[k]))+" ",true);
debugDump("debug_dump.txt","\n",true);
debugDump("debug_results.txt",String.format("%e %e\n",mLikelihoodArray[idx],likelihoodArray[idx]),true);
++idx;
}
}
}
debugClose();
}
++numComputeLikelihoodCalls;
//if(numComputeLikelihoodCalls == 5)
//jniClose();
//System.exit(0);
return likelihoodMap;
}
/**
* {@inheritDoc}
*/
@Override
public double subComputeReadLikelihoodGivenHaplotypeLog10( final byte[] haplotypeBases,
final byte[] readBases,
final byte[] readQuals,
final byte[] insertionGOP,
final byte[] deletionGOP,
final byte[] overallGCP,
final int hapStartIndex,
final boolean recacheReadValues,
final int nextHapStartIndex) {
//System.out.println("#### START STACK TRACE ####");
//for (StackTraceElement ste : Thread.currentThread().getStackTrace()) {
//System.out.println(ste);
//}
//System.out.println("#### END STACK TRACE ####");
//
if(debug1)
jniSubComputeReadLikelihoodGivenHaplotypeLog10(readBases.length, haplotypeBases.length,
readBases, haplotypeBases, readQuals,
insertionGOP, deletionGOP, overallGCP,
hapStartIndex);
boolean doInitialization = (previousHaplotypeBases == null || previousHaplotypeBases.length != haplotypeBases.length);
if (doInitialization) {
final double initialValue = INITIAL_CONDITION / haplotypeBases.length;
// set the initial value (free deletions in the beginning) for the first row in the deletion matrix
for( int j = 0; j < paddedHaplotypeLength; j++ ) {
deletionMatrix[0][j] = initialValue;
}
}
if ( ! constantsAreInitialized || recacheReadValues ) {
initializeProbabilities(transition, insertionGOP, deletionGOP, overallGCP);
if(debug3)
{
System.out.println("Java: initializeProbabilities lengths : "+insertionGOP.length+" padded "+paddedReadLength+" "+paddedHaplotypeLength);
for(int i=0;i<insertionGOP.length;++i)
for(int j=0;j<6;++j)
debugDump("transitions_java.txt",String.format("%e\n",transition[i+1][j]),true);
}
if(debug2)
jniInitializeProbabilities(transition, insertionGOP, deletionGOP, overallGCP);
// note that we initialized the constants
constantsAreInitialized = true;
}
if(debug3)
System.out.println("Java: initializePriors : lengths "+readBases.length+" "+haplotypeBases.length+" padded "+paddedReadLength+" "+paddedHaplotypeLength + " doNotUseTristateCorrection "+doNotUseTristateCorrection);
initializePriors(haplotypeBases, readBases, readQuals, hapStartIndex);
for (int i = 1; i < paddedReadLength; i++) {
// +1 here is because hapStartIndex is 0-based, but our matrices are 1 based
for (int j = hapStartIndex+1; j < paddedHaplotypeLength; j++) {
updateCell(i, j, prior[i][j], transition[i]);
}
}
// final probability is the log10 sum of the last element in the Match and Insertion state arrays
// this way we ignore all paths that ended in deletions! (huge)
// but we have to sum all the paths ending in the M and I matrices, because they're no longer extended.
final int endI = paddedReadLength - 1;
double finalSumProbabilities = 0.0;
for (int j = 1; j < paddedHaplotypeLength; j++) {
finalSumProbabilities += matchMatrix[endI][j] + insertionMatrix[endI][j];
}
if(debug2)
jniInitializePriorsAndUpdateCells(doInitialization, paddedReadLength, paddedHaplotypeLength,
readBases, haplotypeBases, readQuals,
hapStartIndex);
if(debug)
debugDump("return_values_java.txt",String.format("%e\n",Math.log10(finalSumProbabilities) - INITIAL_CONDITION_LOG10),true);
++numCalls;
//if(numCalls == 100)
//{
//debugClose();
//System.exit(0);
//}
return Math.log10(finalSumProbabilities) - INITIAL_CONDITION_LOG10;
}
/**
* Initializes the matrix that holds all the constants related to the editing
* distance between the read and the haplotype.
*
* @param haplotypeBases the bases of the haplotype
* @param readBases the bases of the read
* @param readQuals the base quality scores of the read
* @param startIndex where to start updating the distanceMatrix (in case this read is similar to the previous read)
*/
protected void initializePriors(final byte[] haplotypeBases, final byte[] readBases, final byte[] readQuals, final int startIndex) {
// initialize the pBaseReadLog10 matrix for all combinations of read x haplotype bases
// Abusing the fact that java initializes arrays with 0.0, so no need to fill in rows and columns below 2.
if(debug3)
System.out.println("hapStartIndex "+startIndex);
for (int i = 0; i < readBases.length; i++) {
final byte x = readBases[i];
final byte qual = readQuals[i];
for (int j = startIndex; j < haplotypeBases.length; j++) {
final byte y = haplotypeBases[j];
prior[i+1][j+1] = ( x == y || x == (byte) 'N' || y == (byte) 'N' ?
QualityUtils.qualToProb(qual) : (QualityUtils.qualToErrorProb(qual) / (doNotUseTristateCorrection ? 1.0 : TRISTATE_CORRECTION)) );
if(debug3)
debugDump("priors_java.txt",String.format("%e\n",prior[i+1][j+1]),true);
}
}
}
/**
* Initializes the matrix that holds all the constants related to quality scores.
*
* @param insertionGOP insertion quality scores of the read
* @param deletionGOP deletion quality scores of the read
* @param overallGCP overall gap continuation penalty
*/
@Requires({
"insertionGOP != null",
"deletionGOP != null",
"overallGCP != null"
})
@Ensures("constantsAreInitialized")
protected static void initializeProbabilities(final double[][] transition, final byte[] insertionGOP, final byte[] deletionGOP, final byte[] overallGCP) {
for (int i = 0; i < insertionGOP.length; i++) {
final int qualIndexGOP = Math.min(insertionGOP[i] + deletionGOP[i], Byte.MAX_VALUE);
transition[i+1][matchToMatch] = QualityUtils.qualToProb((byte) qualIndexGOP);
transition[i+1][indelToMatch] = QualityUtils.qualToProb(overallGCP[i]);
transition[i+1][matchToInsertion] = QualityUtils.qualToErrorProb(insertionGOP[i]);
transition[i+1][insertionToInsertion] = QualityUtils.qualToErrorProb(overallGCP[i]);
transition[i+1][matchToDeletion] = QualityUtils.qualToErrorProb(deletionGOP[i]);
transition[i+1][deletionToDeletion] = QualityUtils.qualToErrorProb(overallGCP[i]);
//TODO it seems that it is not always the case that matchToMatch + matchToDeletion + matchToInsertion == 1.
//TODO We have detected cases of 1.00002 which can cause problems downstream. This are typically masked
//TODO by the fact that we always add a indelToMatch penalty to all PairHMM likelihoods (~ -0.1)
//TODO This is in fact not well justified and although it does not have any effect (since is equally added to all
//TODO haplotypes likelihoods) perhaps we should just remove it eventually and fix this != 1.0 issue here.
}
}
/**
* Updates a cell in the HMM matrix
*
* The read and haplotype indices are offset by one because the state arrays have an extra column to hold the
* initial conditions
* @param indI row index in the matrices to update
* @param indJ column index in the matrices to update
* @param prior the likelihood editing distance matrix for the read x haplotype
* @param transition an array with the six transition relevant to this location
*/
protected void updateCell( final int indI, final int indJ, final double prior, final double[] transition) {
matchMatrix[indI][indJ] = prior * ( matchMatrix[indI - 1][indJ - 1] * transition[matchToMatch] +
insertionMatrix[indI - 1][indJ - 1] * transition[indelToMatch] +
deletionMatrix[indI - 1][indJ - 1] * transition[indelToMatch] );
insertionMatrix[indI][indJ] = matchMatrix[indI - 1][indJ] * transition[matchToInsertion] + insertionMatrix[indI - 1][indJ] * transition[insertionToInsertion];
deletionMatrix[indI][indJ] = matchMatrix[indI][indJ - 1] * transition[matchToDeletion] + deletionMatrix[indI][indJ - 1] * transition[deletionToDeletion];
if(debug3)
{
debugDump("matrices_java.txt",String.format("%e\n",matchMatrix[indI][indJ]),true);
debugDump("matrices_java.txt",String.format("%e\n",insertionMatrix[indI][indJ]),true);
debugDump("matrices_java.txt",String.format("%e\n",deletionMatrix[indI][indJ]),true);
}
}
public abstract class JNILoglessPairHMM extends LoglessPairHMM {
public abstract HashMap<Haplotype, Integer> getHaplotypeToHaplotypeListIdxMap();
}

50
protected/java/src/org/broadinstitute/sting/utils/pairhmm/LoglessPairHMM.java
View File

@ -50,13 +50,6 @@ import com.google.java.contract.Ensures;
import com.google.java.contract.Requires;
import org.broadinstitute.sting.utils.QualityUtils;
import java.io.File;
import java.io.FileWriter;
import java.io.BufferedWriter;
import java.util.Map;
import java.util.HashMap;
import java.io.IOException;
/**
* Created with IntelliJ IDEA.
* User: rpoplin, carneiro
@ -76,48 +69,7 @@ public class LoglessPairHMM extends N2MemoryPairHMM {
protected static final int matchToDeletion = 4;
protected static final int deletionToDeletion = 5;
protected static int numCalls = 0;
protected HashMap<String, BufferedWriter> filenameToWriter = new HashMap<String, BufferedWriter>();
protected void debugDump(String filename, String s, boolean toAppend)
{
try
{
File file = new File(filename);
if (!file.exists())
file.createNewFile();
BufferedWriter currWriter = filenameToWriter.get(filename);
if(currWriter == null)
{
FileWriter fw = new FileWriter(file, toAppend);
currWriter = new BufferedWriter(fw);
filenameToWriter.put(filename, currWriter);
}
currWriter.write(s);
}
catch(IOException e)
{
e.printStackTrace();
}
}
protected void debugClose()
{
for(Map.Entry<String, BufferedWriter> currEntry : filenameToWriter.entrySet())
{
BufferedWriter currWriter = currEntry.getValue();
try
{
currWriter.flush();
currWriter.close();
}
catch(IOException e)
{
e.printStackTrace();
}
}
filenameToWriter.clear();
}
/**
* {@inheritDoc}
*/

198
protected/java/src/org/broadinstitute/sting/utils/pairhmm/VectorLoglessPairHMM.java 100644
View File

@ -0,0 +1,198 @@
/*
* By downloading the PROGRAM you agree to the following terms of use:
*
* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
*
* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
*
* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
*
* 1. DEFINITIONS
* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
*
* 2. LICENSE
* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
*
* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
* Copyright 2012 Broad Institute, Inc.
* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
*
* 4. INDEMNIFICATION
* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
*
* 5. NO REPRESENTATIONS OR WARRANTIES
* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
*
* 6. ASSIGNMENT
* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
*
* 7. MISCELLANEOUS
* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
*/
package org.broadinstitute.sting.utils.pairhmm;
import com.google.java.contract.Ensures;
import com.google.java.contract.Requires;
import org.broadinstitute.sting.utils.QualityUtils;
import org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap;
import org.broadinstitute.sting.utils.haplotype.Haplotype;
import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
import org.broadinstitute.variant.variantcontext.Allele;
import java.util.List;
import java.util.Map;
import java.util.HashMap;
/**
* Created with IntelliJ IDEA.
* User: rpoplin, carneiro
* Date: 10/16/12
*/
public class VectorLoglessPairHMM extends JNILoglessPairHMM {
//For machine capabilities
public static final long sse42Mask = 1;
public static final long avxMask = 2;
//Used to copy references to byteArrays to JNI from reads
protected class JNIReadDataHolderClass {
public byte[] readBases = null;
public byte[] readQuals = null;
public byte[] insertionGOP = null;
public byte[] deletionGOP = null;
public byte[] overallGCP = null;
}
//Used to copy references to byteArrays to JNI from haplotypes
protected class JNIHaplotypeDataHolderClass {
public byte[] haplotypeBases = null;
}
public native long jniGetMachineType();
private native void jniClose();
private native void jniGlobalInit(Class<?> readDataHolderClass, Class<?> haplotypeDataHolderClass, long mask);
private static boolean isVectorLoglessPairHMMLibraryLoaded = false;
public VectorLoglessPairHMM() {
super();
if(!isVectorLoglessPairHMMLibraryLoaded) {
System.loadLibrary("VectorLoglessPairHMM");
isVectorLoglessPairHMMLibraryLoaded = true;
jniGlobalInit(JNIReadDataHolderClass.class, JNIHaplotypeDataHolderClass.class, 0xFFFFFFFFFFFFFFFFl); //need to do this only once
}
}
@Override
public void close()
{
jniClose();
}
private native void jniInitializeHaplotypes(final int numHaplotypes, JNIHaplotypeDataHolderClass[] haplotypeDataArray);
//Hold the mapping between haplotype and index in the list of Haplotypes passed to initialize
//Use this mapping in computeLikelihoods to find the likelihood value corresponding to a given Haplotype
private HashMap<Haplotype,Integer> haplotypeToHaplotypeListIdxMap = new HashMap<Haplotype,Integer>();
@Override
public HashMap<Haplotype, Integer> getHaplotypeToHaplotypeListIdxMap() { return haplotypeToHaplotypeListIdxMap; }
//Used to transfer data to JNI
//Since the haplotypes are the same for all calls to computeLikelihoods within a region, transfer the haplotypes only once to the JNI per region
/**
* {@inheritDoc}
*/
@Override
public void initialize( final List<Haplotype> haplotypes, final Map<String, List<GATKSAMRecord>> perSampleReadList,
final int readMaxLength, final int haplotypeMaxLength ) {
int numHaplotypes = haplotypes.size();
JNIHaplotypeDataHolderClass[] haplotypeDataArray = new JNIHaplotypeDataHolderClass[numHaplotypes];
int idx = 0;
haplotypeToHaplotypeListIdxMap.clear();
for(final Haplotype currHaplotype : haplotypes)
{
haplotypeDataArray[idx] = new JNIHaplotypeDataHolderClass();
haplotypeDataArray[idx].haplotypeBases = currHaplotype.getBases();
haplotypeToHaplotypeListIdxMap.put(currHaplotype, idx);
++idx;
}
jniInitializeHaplotypes(numHaplotypes, haplotypeDataArray);
}
private native void jniFinalizeRegion();
//Tell JNI to release arrays - really important if native code is directly accessing Java memory, if not
//accessing Java memory directly, still important to release memory from C++
/**
* {@inheritDoc}
*/
@Override
public void finalizeRegion()
{
jniFinalizeRegion();
}
//Real compute kernel
private native void jniComputeLikelihoods(int numReads, int numHaplotypes, JNIReadDataHolderClass[] readDataArray,
JNIHaplotypeDataHolderClass[] haplotypeDataArray, double[] likelihoodArray, int maxNumThreadsToUse);
/**
* {@inheritDoc}
*/
@Override
public PerReadAlleleLikelihoodMap computeLikelihoods( final List<GATKSAMRecord> reads, final Map<Allele, Haplotype> alleleHaplotypeMap, final Map<GATKSAMRecord, byte[]> GCPArrayMap ) {
int readListSize = reads.size();
int numHaplotypes = alleleHaplotypeMap.size();
int numTestcases = readListSize*numHaplotypes;
JNIReadDataHolderClass[] readDataArray = new JNIReadDataHolderClass[readListSize];
int idx = 0;
for(GATKSAMRecord read : reads)
{
readDataArray[idx] = new JNIReadDataHolderClass();
readDataArray[idx].readBases = read.getReadBases();
readDataArray[idx].readQuals = read.getBaseQualities();
readDataArray[idx].insertionGOP = read.getBaseInsertionQualities();
readDataArray[idx].deletionGOP = read.getBaseDeletionQualities();
readDataArray[idx].overallGCP = GCPArrayMap.get(read);
++idx;
}
JNIHaplotypeDataHolderClass[] haplotypeDataArray = new JNIHaplotypeDataHolderClass[numHaplotypes];
mLikelihoodArray = new double[readListSize*numHaplotypes]; //to store results
//for(reads)
// for(haplotypes)
// compute_full_prob()
jniComputeLikelihoods(readListSize, numHaplotypes, readDataArray, haplotypeDataArray, mLikelihoodArray, 12);
final PerReadAlleleLikelihoodMap likelihoodMap = new PerReadAlleleLikelihoodMap();
idx = 0;
int idxInsideHaplotypeList = 0;
int readIdx = 0;
for(GATKSAMRecord read : reads)
{
for (Map.Entry<Allele, Haplotype> currEntry : alleleHaplotypeMap.entrySet())//order is important - access in same order always
{
//Since the order of haplotypes in the List<Haplotype> and alleleHaplotypeMap is different,
//get idx of current haplotype in the list and use this idx to get the right likelihoodValue
idxInsideHaplotypeList = haplotypeToHaplotypeListIdxMap.get(currEntry.getValue());
likelihoodMap.add(read, currEntry.getKey(), mLikelihoodArray[readIdx + idxInsideHaplotypeList]);
++idx;
}
readIdx += numHaplotypes;
}
return likelihoodMap;
}
}

5
public/java/src/org/broadinstitute/sting/utils/pairhmm/PairHMM.java
View File

@ -58,7 +58,9 @@ public abstract class PairHMM {
/* Optimized version of the PairHMM which caches per-read computations and operations in real space to avoid costly sums of log10'ed likelihoods */
LOGLESS_CACHING,
/* Optimized AVX implementation of LOGLESS_CACHING called through JNI */
JNI_LOGLESS_CACHING,
VECTOR_LOGLESS_CACHING,
/* Debugging for any JNI implementation of LOGLESS_CACHING */
DEBUG_JNI_LOGLESS_CACHING,
/* Logless caching PairHMM that stores computations in 1D arrays instead of matrices, and which proceeds diagonally over the (read x haplotype) intersection matrix */
ARRAY_LOGLESS
}
@ -301,6 +303,7 @@ public abstract class PairHMM {
return Math.min(haplotype1.length, haplotype2.length);
}
public double[] getLikelihoodArray() { return mLikelihoodArray; }
//Called at the end of all HC calls
public void close() { ; }
}