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

Bring back from the dead the old likelihoods model for indels, which has worse performance but is about 4x faster. Enabled with argument -GSA_PRODUCTION_ONLY in UG 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@5748 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
delangel 2011-05-03 22:38:33 +00:00
parent f35d955490
commit a19389528d
2 changed files with 36 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
java/src/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java
View File

@ -56,12 +56,33 @@ public class IndelGenotypeLikelihoodsCalculationModel extends GenotypeLikelihood
private boolean DEBUG = false; private boolean DEBUG = false;
private PairHMMIndelErrorModel pairModel; private PairHMMIndelErrorModel pairModel;
// gdebug removeme
// todo -cleanup
private HaplotypeIndelErrorModel model;
private boolean useOldWrongHorribleHackedUpLikelihoodModel = false;
//
private GenomeLoc lastSiteVisited; private GenomeLoc lastSiteVisited;
private List<Allele> alleleList; private List<Allele> alleleList;
protected IndelGenotypeLikelihoodsCalculationModel(UnifiedArgumentCollection UAC, Logger logger) { protected IndelGenotypeLikelihoodsCalculationModel(UnifiedArgumentCollection UAC, Logger logger) {
super(UAC, logger); super(UAC, logger);
if (UAC.GSA_PRODUCTION_ONLY == false) {
pairModel = new PairHMMIndelErrorModel(UAC.INDEL_GAP_OPEN_PENALTY,UAC.INDEL_GAP_CONTINUATION_PENALTY,
UAC.OUTPUT_DEBUG_INDEL_INFO, UAC.DO_CONTEXT_DEPENDENT_PENALTIES, UAC.dovit, UAC.GET_GAP_PENALTIES_FROM_DATA, UAC.INDEL_RECAL_FILE);
useOldWrongHorribleHackedUpLikelihoodModel = false;
}
else {
useOldWrongHorribleHackedUpLikelihoodModel = true;
double INSERTION_START_PROBABILITY = 1e-3;
double INSERTION_END_PROBABILITY = 0.5;
double ALPHA_DELETION_PROBABILITY = 1e-3;
model = new HaplotypeIndelErrorModel(3, INSERTION_START_PROBABILITY,
INSERTION_END_PROBABILITY,ALPHA_DELETION_PROBABILITY,UAC.INDEL_HAPLOTYPE_SIZE, false, UAC.OUTPUT_DEBUG_INDEL_INFO);
}
pairModel = new PairHMMIndelErrorModel(UAC.INDEL_GAP_OPEN_PENALTY,UAC.INDEL_GAP_CONTINUATION_PENALTY, pairModel = new PairHMMIndelErrorModel(UAC.INDEL_GAP_OPEN_PENALTY,UAC.INDEL_GAP_CONTINUATION_PENALTY,
UAC.OUTPUT_DEBUG_INDEL_INFO, UAC.DO_CONTEXT_DEPENDENT_PENALTIES, UAC.dovit, UAC.GET_GAP_PENALTIES_FROM_DATA, UAC.INDEL_RECAL_FILE); UAC.OUTPUT_DEBUG_INDEL_INFO, UAC.DO_CONTEXT_DEPENDENT_PENALTIES, UAC.dovit, UAC.GET_GAP_PENALTIES_FROM_DATA, UAC.INDEL_RECAL_FILE);
alleleList = new ArrayList<Allele>(); alleleList = new ArrayList<Allele>();
@ -322,7 +343,12 @@ public class IndelGenotypeLikelihoodsCalculationModel extends GenotypeLikelihood
List<Haplotype> haplotypesInVC; List<Haplotype> haplotypesInVC;
int hsize = (int)ref.getWindow().size()-Math.abs(eventLength)-1; int hsize = (int)ref.getWindow().size()-Math.abs(eventLength)-1;
int numPrefBases = ref.getLocus().getStart()-ref.getWindow().getStart()+1; int numPrefBases= ref.getLocus().getStart()-ref.getWindow().getStart()+1;
if (useOldWrongHorribleHackedUpLikelihoodModel) {
numPrefBases = 20;
hsize=80;
}
if (DEBUG) if (DEBUG)
System.out.format("hsize: %d eventLength: %d refSize: %d, locStart: %d numpr: %d\n",hsize,eventLength, System.out.format("hsize: %d eventLength: %d refSize: %d, locStart: %d numpr: %d\n",hsize,eventLength,
(int)ref.getWindow().size(), loc.getStart(), numPrefBases); (int)ref.getWindow().size(), loc.getStart(), numPrefBases);
@ -348,10 +374,13 @@ public class IndelGenotypeLikelihoodsCalculationModel extends GenotypeLikelihood
pileup = context.getBasePileup(); pileup = context.getBasePileup();
if (pileup != null ) { if (pileup != null ) {
if (useOldWrongHorribleHackedUpLikelihoodModel)
haplotypeLikehoodMatrix = model.computeReadHaplotypeLikelihoods( pileup, haplotypesInVC);
else
haplotypeLikehoodMatrix = pairModel.computeReadHaplotypeLikelihoods( pileup, haplotypesInVC, ref, HAPLOTYPE_SIZE, eventLength); haplotypeLikehoodMatrix = pairModel.computeReadHaplotypeLikelihoods( pileup, haplotypesInVC, ref, HAPLOTYPE_SIZE, eventLength);
double[] genotypeLikelihoods = HaplotypeIndelErrorModel.getHaplotypeLikelihoods( haplotypeLikehoodMatrix); double[] genotypeLikelihoods = HaplotypeIndelErrorModel.getHaplotypeLikelihoods( haplotypeLikehoodMatrix);
GLs.put(sample.getKey(), new BiallelicGenotypeLikelihoods(sample.getKey(), GLs.put(sample.getKey(), new BiallelicGenotypeLikelihoods(sample.getKey(),

4
java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedArgumentCollection.java
View File

@ -117,6 +117,9 @@ public class UnifiedArgumentCollection {
@Argument(fullName = "dovit", shortName = "dovit", doc = "Output indel debug info", required = false) @Argument(fullName = "dovit", shortName = "dovit", doc = "Output indel debug info", required = false)
public boolean dovit = false; public boolean dovit = false;
@Hidden @Hidden
@Argument(fullName = "GSA_PRODUCTION_ONLY", shortName = "GSA_PRODUCTION_ONLY", doc = "don't ever use me", required = false)
public boolean GSA_PRODUCTION_ONLY = false;
@Hidden
@Argument(fullName = "exactCalculation", shortName = "exactCalculation", doc = "expt", required = false) @Argument(fullName = "exactCalculation", shortName = "exactCalculation", doc = "expt", required = false)
public ExactAFCalculationModel.ExactCalculation EXACT_CALCULATION_TYPE = ExactAFCalculationModel.ExactCalculation.LINEAR_EXPERIMENTAL; public ExactAFCalculationModel.ExactCalculation EXACT_CALCULATION_TYPE = ExactAFCalculationModel.ExactCalculation.LINEAR_EXPERIMENTAL;
@ -161,6 +164,7 @@ public class UnifiedArgumentCollection {
// todo- arguments to remove // todo- arguments to remove
uac.COVERAGE_AT_WHICH_TO_ABORT = COVERAGE_AT_WHICH_TO_ABORT; uac.COVERAGE_AT_WHICH_TO_ABORT = COVERAGE_AT_WHICH_TO_ABORT;
uac.dovit = dovit; uac.dovit = dovit;
uac.GSA_PRODUCTION_ONLY = GSA_PRODUCTION_ONLY;
return uac; return uac;
} }