From d79b5f0931db2a65138c636c23c67d602688e480 Mon Sep 17 00:00:00 2001 From: sathibault Date: Wed, 8 May 2013 11:01:20 -0500 Subject: [PATCH] Adding Convey HC-1 HMM acceleration --- .../LikelihoodCalculationEngine.java | 38 ++++++++---- .../sting/utils/pairhmm/CnyPairHMM.java | 59 +++++++++++++++++++ .../sting/utils/pairhmm/BatchPairHMM.java | 16 +++++ 3 files changed, 103 insertions(+), 10 deletions(-) create mode 100644 protected/java/src/org/broadinstitute/sting/utils/pairhmm/CnyPairHMM.java create mode 100644 public/java/src/org/broadinstitute/sting/utils/pairhmm/BatchPairHMM.java diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java index 8697833a6..62d4d17fd 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java @@ -90,7 +90,10 @@ public class LikelihoodCalculationEngine { pairHMM = new Log10PairHMM(false); break; case LOGLESS_CACHING: - pairHMM = new LoglessPairHMM(); + if (CnyPairHMM.isAvailable()) + pairHMM = new CnyPairHMM(); + else + pairHMM = new LoglessPairHMM(); break; default: throw new UserException.BadArgumentValue("pairHMM", "Specified pairHMM implementation is unrecognized or incompatible with the HaplotypeCaller. Acceptable options are ORIGINAL, EXACT, CACHING, and LOGLESS_CACHING."); @@ -151,6 +154,8 @@ public class LikelihoodCalculationEngine { private PerReadAlleleLikelihoodMap computeReadLikelihoods( final List haplotypes, final List reads) { // first, a little set up to get copies of the Haplotypes that are Alleles (more efficient than creating them each time) + final BatchPairHMM batchPairHMM = (pairHMM instanceof BatchPairHMM) ? (BatchPairHMM)pairHMM : null; + final Vector batchedReads = new Vector(reads.size()); final int numHaplotypes = haplotypes.size(); final Map alleleVersions = new HashMap(numHaplotypes); for ( final Haplotype haplotype : haplotypes ) { @@ -177,16 +182,29 @@ public class LikelihoodCalculationEngine { readQuals[kkk] = ( readQuals[kkk] < (byte) 18 ? QualityUtils.MIN_USABLE_Q_SCORE : readQuals[kkk] ); } - for( int jjj = 0; jjj < numHaplotypes; jjj++ ) { - final Haplotype haplotype = haplotypes.get(jjj); - final boolean isFirstHaplotype = jjj == 0; - final double log10l = pairHMM.computeReadLikelihoodGivenHaplotypeLog10(haplotype.getBases(), - read.getReadBases(), readQuals, readInsQuals, readDelQuals, overallGCP, isFirstHaplotype); - - perReadAlleleLikelihoodMap.add(read, alleleVersions.get(haplotype), log10l); - } + if ( batchPairHMM != null ) { + batchPairHMM.batchAdd(haplotypes, read.getReadBases(), readQuals, readInsQuals, readDelQuals, overallGCP); + batchedReads.add(read); + } else { + for( int jjj = 0; jjj < numHaplotypes; jjj++ ) { + final Haplotype haplotype = haplotypes.get(jjj); + final boolean isFirstHaplotype = jjj == 0; + final double log10l = pairHMM.computeReadLikelihoodGivenHaplotypeLog10(haplotype.getBases(), + read.getReadBases(), readQuals, readInsQuals, readDelQuals, overallGCP, isFirstHaplotype); + + perReadAlleleLikelihoodMap.add(read, alleleVersions.get(haplotype), log10l); + } + } } - + if ( batchPairHMM != null ) { + for( final GATKSAMRecord read : batchedReads ) { + final double[] likelihoods = batchPairHMM.batchResult(); + for( int jjj = 0; jjj < numHaplotypes; jjj++ ) { + final Haplotype haplotype = haplotypes.get(jjj); + perReadAlleleLikelihoodMap.add(read, alleleVersions.get(haplotype), likelihoods[jjj]); + } + } + } return perReadAlleleLikelihoodMap; } diff --git a/protected/java/src/org/broadinstitute/sting/utils/pairhmm/CnyPairHMM.java b/protected/java/src/org/broadinstitute/sting/utils/pairhmm/CnyPairHMM.java new file mode 100644 index 000000000..51611bb08 --- /dev/null +++ b/protected/java/src/org/broadinstitute/sting/utils/pairhmm/CnyPairHMM.java @@ -0,0 +1,59 @@ +package org.broadinstitute.sting.utils.pairhmm; + +import java.util.*; + +import org.broadinstitute.sting.utils.haplotype.Haplotype; + +public final class CnyPairHMM extends PairHMM implements BatchPairHMM { + private static class HmmInput { + public List haplotypes; + public byte[] readBases; + public byte[] readQuals; + public byte[] insertionGOP; + public byte[] deletionGOP; + public byte[] overallGCP; + }; + + private static boolean loaded = false; + private List pending = new LinkedList(); + + static public boolean isAvailable() { + return true; + } + + public void batchAdd(final List haplotypes, + final byte[] readBases, + final byte[] readQuals, + final byte[] insertionGOP, + final byte[] deletionGOP, + final byte[] overallGCP) { + HmmInput test=new HmmInput(); + test.haplotypes=haplotypes; + test.readBases=readBases; + test.readQuals=readQuals; + test.insertionGOP=insertionGOP; + test.deletionGOP=deletionGOP; + test.overallGCP=overallGCP; + pending.add(test); + } + + public double[] batchResult() { + HmmInput test=pending.remove(0); + double[] results=new double[test.haplotypes.size()]; + for (int i=0; i haplotypes, + final byte[] readBases, + final byte[] readQuals, + final byte[] insertionGOP, + final byte[] deletionGOP, + final byte[] overallGCP); + + public double[] batchResult(); +}