Reorg of UG args; checking in first before upcoming changes that will break integration tests.
git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@3274 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
parent
ba78d146ec
commit
1714c322c2
|
|
@ -74,7 +74,7 @@ public class DiploidGenotypeCalculationModel extends JointEstimateGenotypeCalcul
|
||||||
ReadBackedPileup pileup = context.getContext(contextType).getBasePileup();
|
ReadBackedPileup pileup = context.getContext(contextType).getBasePileup();
|
||||||
|
|
||||||
// create the GenotypeLikelihoods object
|
// create the GenotypeLikelihoods object
|
||||||
GenotypeLikelihoods GL = new GenotypeLikelihoods(baseModel, priors, defaultPlatform);
|
GenotypeLikelihoods GL = new GenotypeLikelihoods(UAC.baseModel, priors, UAC.defaultPlatform);
|
||||||
|
|
||||||
GL.add(pileup, true);
|
GL.add(pileup, true);
|
||||||
GLs.put(sample, GL);
|
GLs.put(sample, GL);
|
||||||
|
|
|
||||||
|
|
@ -25,7 +25,6 @@
|
||||||
|
|
||||||
package org.broadinstitute.sting.gatk.walkers.genotyper;
|
package org.broadinstitute.sting.gatk.walkers.genotyper;
|
||||||
|
|
||||||
import org.broadinstitute.sting.utils.Utils;
|
|
||||||
import org.broadinstitute.sting.utils.MathUtils;
|
import org.broadinstitute.sting.utils.MathUtils;
|
||||||
import org.broadinstitute.sting.utils.genotype.DiploidGenotype;
|
import org.broadinstitute.sting.utils.genotype.DiploidGenotype;
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -25,18 +25,10 @@ public abstract class GenotypeCalculationModel implements Cloneable {
|
||||||
INDELS
|
INDELS
|
||||||
}
|
}
|
||||||
|
|
||||||
protected BaseMismatchModel baseModel;
|
protected UnifiedArgumentCollection UAC;
|
||||||
protected Set<String> samples;
|
protected Set<String> samples;
|
||||||
protected Logger logger;
|
protected Logger logger;
|
||||||
protected double heterozygosity;
|
|
||||||
protected EmpiricalSubstitutionProbabilities.SequencerPlatform defaultPlatform;
|
|
||||||
protected GenotypeWriterFactory.GENOTYPE_FORMAT OUTPUT_FORMAT;
|
protected GenotypeWriterFactory.GENOTYPE_FORMAT OUTPUT_FORMAT;
|
||||||
protected boolean ALL_BASE_MODE;
|
|
||||||
protected boolean GENOTYPE_MODE;
|
|
||||||
protected int POOL_SIZE;
|
|
||||||
protected double CONFIDENCE_THRESHOLD;
|
|
||||||
protected double MINIMUM_ALLELE_FREQUENCY;
|
|
||||||
protected boolean REPORT_SLOD;
|
|
||||||
protected PrintStream verboseWriter;
|
protected PrintStream verboseWriter;
|
||||||
protected PrintStream beagleWriter;
|
protected PrintStream beagleWriter;
|
||||||
|
|
||||||
|
|
@ -62,17 +54,10 @@ public abstract class GenotypeCalculationModel implements Cloneable {
|
||||||
GenotypeWriterFactory.GENOTYPE_FORMAT outputFormat,
|
GenotypeWriterFactory.GENOTYPE_FORMAT outputFormat,
|
||||||
PrintStream verboseWriter,
|
PrintStream verboseWriter,
|
||||||
PrintStream beagleWriter) {
|
PrintStream beagleWriter) {
|
||||||
|
this.UAC = UAC.clone();
|
||||||
this.samples = new TreeSet<String>(samples);
|
this.samples = new TreeSet<String>(samples);
|
||||||
this.logger = logger;
|
|
||||||
baseModel = UAC.baseModel;
|
|
||||||
heterozygosity = UAC.heterozygosity;
|
|
||||||
defaultPlatform = UAC.defaultPlatform;
|
|
||||||
OUTPUT_FORMAT = outputFormat;
|
OUTPUT_FORMAT = outputFormat;
|
||||||
ALL_BASE_MODE = UAC.ALL_BASES;
|
this.logger = logger;
|
||||||
GENOTYPE_MODE = UAC.GENOTYPE;
|
|
||||||
POOL_SIZE = UAC.POOLSIZE;
|
|
||||||
CONFIDENCE_THRESHOLD = UAC.CONFIDENCE_THRESHOLD;
|
|
||||||
REPORT_SLOD = ! UAC.NO_SLOD;
|
|
||||||
this.verboseWriter = verboseWriter;
|
this.verboseWriter = verboseWriter;
|
||||||
this.beagleWriter = beagleWriter;
|
this.beagleWriter = beagleWriter;
|
||||||
}
|
}
|
||||||
|
|
|
||||||
|
|
@ -60,7 +60,7 @@ public abstract class JointEstimateGenotypeCalculationModel extends GenotypeCalc
|
||||||
// if there are no non-ref bases...
|
// if there are no non-ref bases...
|
||||||
if ( bestAlternateAllele == null ) {
|
if ( bestAlternateAllele == null ) {
|
||||||
// if we don't want all bases, then we don't need to calculate genotype likelihoods
|
// if we don't want all bases, then we don't need to calculate genotype likelihoods
|
||||||
if ( !triggerTrack && !ALL_BASE_MODE && !GENOTYPE_MODE ) {
|
if ( !triggerTrack && !UAC.ALL_BASES_MODE && !UAC.GENOTYPE_MODE ) {
|
||||||
VariantCallContext vcc = new VariantCallContext(false);
|
VariantCallContext vcc = new VariantCallContext(false);
|
||||||
estimateReferenceConfidence(vcc, contexts, DiploidGenotypePriors.HUMAN_HETEROZYGOSITY, false);
|
estimateReferenceConfidence(vcc, contexts, DiploidGenotypePriors.HUMAN_HETEROZYGOSITY, false);
|
||||||
return vcc;
|
return vcc;
|
||||||
|
|
@ -152,7 +152,7 @@ public abstract class JointEstimateGenotypeCalculationModel extends GenotypeCalc
|
||||||
sigma_1_over_I += 1.0 / (double)i;
|
sigma_1_over_I += 1.0 / (double)i;
|
||||||
|
|
||||||
// delta = theta / sum(1/i)
|
// delta = theta / sum(1/i)
|
||||||
double delta = heterozygosity / sigma_1_over_I;
|
double delta = UAC.heterozygosity / sigma_1_over_I;
|
||||||
|
|
||||||
// calculate the null allele frequencies for 1-N
|
// calculate the null allele frequencies for 1-N
|
||||||
AFs = new double[N];
|
AFs = new double[N];
|
||||||
|
|
@ -190,7 +190,7 @@ public abstract class JointEstimateGenotypeCalculationModel extends GenotypeCalc
|
||||||
P_of_ref *= 1.0 - (theta / 2.0) * MathUtils.binomialProbability(0, depth, 0.5);
|
P_of_ref *= 1.0 - (theta / 2.0) * MathUtils.binomialProbability(0, depth, 0.5);
|
||||||
}
|
}
|
||||||
|
|
||||||
vcc.confidentlyCalled = QualityUtils.phredScaleErrorRate(1.0 - P_of_ref) >= CONFIDENCE_THRESHOLD;
|
vcc.confidentlyCalled = QualityUtils.phredScaleErrorRate(1.0 - P_of_ref) >= UAC.CONFIDENCE_THRESHOLD;
|
||||||
}
|
}
|
||||||
|
|
||||||
protected void calculateAlleleFrequencyPosteriors(char ref, int frequencyEstimationPoints, Map<String, StratifiedAlignmentContext> contexts, StratifiedAlignmentContext.StratifiedContextType contextType) {
|
protected void calculateAlleleFrequencyPosteriors(char ref, int frequencyEstimationPoints, Map<String, StratifiedAlignmentContext> contexts, StratifiedAlignmentContext.StratifiedContextType contextType) {
|
||||||
|
|
@ -350,8 +350,8 @@ public abstract class JointEstimateGenotypeCalculationModel extends GenotypeCalc
|
||||||
}
|
}
|
||||||
|
|
||||||
// return a null call if we don't pass the confidence cutoff or the most likely allele frequency is zero
|
// return a null call if we don't pass the confidence cutoff or the most likely allele frequency is zero
|
||||||
if ( !triggerTrack && !ALL_BASE_MODE && ((!GENOTYPE_MODE && bestAFguess == 0) || phredScaledConfidence < CONFIDENCE_THRESHOLD) )
|
if ( !triggerTrack && !UAC.ALL_BASES_MODE && ((!UAC.GENOTYPE_MODE && bestAFguess == 0) || phredScaledConfidence < UAC.CONFIDENCE_THRESHOLD) )
|
||||||
return new VariantCallContext(phredScaledConfidence >= CONFIDENCE_THRESHOLD);
|
return new VariantCallContext(phredScaledConfidence >= UAC.CONFIDENCE_THRESHOLD);
|
||||||
|
|
||||||
// output to beagle file if requested
|
// output to beagle file if requested
|
||||||
if ( beagleWriter != null ) {
|
if ( beagleWriter != null ) {
|
||||||
|
|
@ -384,7 +384,7 @@ public abstract class JointEstimateGenotypeCalculationModel extends GenotypeCalc
|
||||||
if ( dbsnp != null )
|
if ( dbsnp != null )
|
||||||
attributes.put("ID", dbsnp.getRS_ID());
|
attributes.put("ID", dbsnp.getRS_ID());
|
||||||
|
|
||||||
if ( REPORT_SLOD ) {
|
if ( !UAC.NO_SLOD ) {
|
||||||
// the overall lod
|
// the overall lod
|
||||||
double overallLog10PofNull = log10AlleleFrequencyPriors[0] + log10PofDgivenAFi[indexOfMax][0];
|
double overallLog10PofNull = log10AlleleFrequencyPriors[0] + log10PofDgivenAFi[indexOfMax][0];
|
||||||
double overallLog10PofF = log10AlleleFrequencyPriors[bestAFguess] + log10PofDgivenAFi[indexOfMax][bestAFguess];
|
double overallLog10PofF = log10AlleleFrequencyPriors[bestAFguess] + log10PofDgivenAFi[indexOfMax][bestAFguess];
|
||||||
|
|
@ -422,6 +422,6 @@ public abstract class JointEstimateGenotypeCalculationModel extends GenotypeCalc
|
||||||
|
|
||||||
VariantContext vc = new VariantContext("UG_SNP_call", loc, alleles, genotypes, phredScaledConfidence/10.0, null, attributes);
|
VariantContext vc = new VariantContext("UG_SNP_call", loc, alleles, genotypes, phredScaledConfidence/10.0, null, attributes);
|
||||||
|
|
||||||
return new VariantCallContext(vc, phredScaledConfidence >= CONFIDENCE_THRESHOLD);
|
return new VariantCallContext(vc, phredScaledConfidence >= UAC.CONFIDENCE_THRESHOLD);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
@ -1,11 +1,9 @@
|
||||||
package org.broadinstitute.sting.gatk.walkers.genotyper;
|
package org.broadinstitute.sting.gatk.walkers.genotyper;
|
||||||
|
|
||||||
import org.broadinstitute.sting.utils.*;
|
import org.broadinstitute.sting.utils.*;
|
||||||
import org.broadinstitute.sting.utils.genotype.vcf.VCFRecord;
|
|
||||||
import org.broadinstitute.sting.utils.collections.Pair;
|
import org.broadinstitute.sting.utils.collections.Pair;
|
||||||
import org.broadinstitute.sting.utils.pileup.*;
|
import org.broadinstitute.sting.utils.pileup.*;
|
||||||
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
|
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
|
||||||
import org.broadinstitute.sting.gatk.refdata.rodDbSNP;
|
|
||||||
import org.broadinstitute.sting.gatk.contexts.*;
|
import org.broadinstitute.sting.gatk.contexts.*;
|
||||||
import org.broadinstitute.sting.gatk.contexts.variantcontext.*;
|
import org.broadinstitute.sting.gatk.contexts.variantcontext.*;
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -45,10 +45,10 @@ public class UnifiedArgumentCollection {
|
||||||
|
|
||||||
// control the output
|
// control the output
|
||||||
@Argument(fullName = "genotype", shortName = "genotype", doc = "Should we output confident genotypes (i.e. including ref calls) or just the variants?", required = false)
|
@Argument(fullName = "genotype", shortName = "genotype", doc = "Should we output confident genotypes (i.e. including ref calls) or just the variants?", required = false)
|
||||||
public boolean GENOTYPE = false;
|
public boolean GENOTYPE_MODE = false;
|
||||||
|
|
||||||
@Argument(fullName = "output_all_callable_bases", shortName = "all_bases", doc = "Should we output all callable bases?", required = false)
|
@Argument(fullName = "output_all_callable_bases", shortName = "all_bases", doc = "Should we output all callable bases?", required = false)
|
||||||
public boolean ALL_BASES = false;
|
public boolean ALL_BASES_MODE = false;
|
||||||
|
|
||||||
@Argument(fullName = "noSLOD", shortName = "nsl", doc = "If provided, we will not calculate the SLOD", required = false)
|
@Argument(fullName = "noSLOD", shortName = "nsl", doc = "If provided, we will not calculate the SLOD", required = false)
|
||||||
public boolean NO_SLOD = false;
|
public boolean NO_SLOD = false;
|
||||||
|
|
@ -80,4 +80,28 @@ public class UnifiedArgumentCollection {
|
||||||
|
|
||||||
@Argument(fullName = "max_deletion_fraction", shortName = "deletions", doc = "Maximum fraction of reads with deletions spanning this locus for it to be callable [to disable, set to < 0 or > 1; default:0.05]", required = false)
|
@Argument(fullName = "max_deletion_fraction", shortName = "deletions", doc = "Maximum fraction of reads with deletions spanning this locus for it to be callable [to disable, set to < 0 or > 1; default:0.05]", required = false)
|
||||||
public Double MAX_DELETION_FRACTION = 0.05;
|
public Double MAX_DELETION_FRACTION = 0.05;
|
||||||
|
|
||||||
|
|
||||||
|
public UnifiedArgumentCollection clone() {
|
||||||
|
UnifiedArgumentCollection uac = new UnifiedArgumentCollection();
|
||||||
|
|
||||||
|
uac.genotypeModel = genotypeModel;
|
||||||
|
uac.baseModel = baseModel;
|
||||||
|
uac.heterozygosity = heterozygosity;
|
||||||
|
uac.GENOTYPE_MODE = GENOTYPE_MODE;
|
||||||
|
uac.ALL_BASES_MODE = ALL_BASES_MODE;
|
||||||
|
uac.NO_SLOD = NO_SLOD;
|
||||||
|
uac.ASSUME_SINGLE_SAMPLE = ASSUME_SINGLE_SAMPLE;
|
||||||
|
uac.defaultPlatform = defaultPlatform;
|
||||||
|
uac.CONFIDENCE_THRESHOLD = CONFIDENCE_THRESHOLD;
|
||||||
|
uac.MIN_BASE_QUALTY_SCORE = MIN_BASE_QUALTY_SCORE;
|
||||||
|
uac.MIN_MAPPING_QUALTY_SCORE = MIN_MAPPING_QUALTY_SCORE;
|
||||||
|
uac.MAX_MISMATCHES = MAX_MISMATCHES;
|
||||||
|
uac.USE_BADLY_MATED_READS = USE_BADLY_MATED_READS;
|
||||||
|
uac.MAX_DELETION_FRACTION = MAX_DELETION_FRACTION;
|
||||||
|
|
||||||
|
return uac;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
|
||||||
|
|
@ -85,7 +85,7 @@ public class BaseTransitionTableCalculatorJavaWalker extends LocusWalker<Set<Bas
|
||||||
}
|
}
|
||||||
UnifiedArgumentCollection uac = new UnifiedArgumentCollection();
|
UnifiedArgumentCollection uac = new UnifiedArgumentCollection();
|
||||||
uac.baseModel = BaseMismatchModel.THREE_STATE;
|
uac.baseModel = BaseMismatchModel.THREE_STATE;
|
||||||
uac.ALL_BASES = true;
|
uac.ALL_BASES_MODE = true;
|
||||||
ug = new UnifiedGenotyperEngine(getToolkit(), uac);
|
ug = new UnifiedGenotyperEngine(getToolkit(), uac);
|
||||||
// refWindow = new ReferenceContextWindow(nPreviousBases);
|
// refWindow = new ReferenceContextWindow(nPreviousBases);
|
||||||
usePreviousBases = new ArrayList<Boolean>();
|
usePreviousBases = new ArrayList<Boolean>();
|
||||||
|
|
|
||||||
|
|
@ -77,7 +77,7 @@ public class LocusMismatchWalker extends LocusWalker<String,Integer> implements
|
||||||
public void initialize() {
|
public void initialize() {
|
||||||
UnifiedArgumentCollection uac = new UnifiedArgumentCollection();
|
UnifiedArgumentCollection uac = new UnifiedArgumentCollection();
|
||||||
uac.baseModel = BaseMismatchModel.THREE_STATE;
|
uac.baseModel = BaseMismatchModel.THREE_STATE;
|
||||||
uac.ALL_BASES = true;
|
uac.ALL_BASES_MODE = true;
|
||||||
ug = new UnifiedGenotyperEngine(getToolkit(), uac);
|
ug = new UnifiedGenotyperEngine(getToolkit(), uac);
|
||||||
|
|
||||||
// print the header
|
// print the header
|
||||||
|
|
|
||||||
|
|
@ -61,7 +61,7 @@ public class SnpCallRateByCoverageWalker extends LocusWalker<List<String>, Strin
|
||||||
public void initialize() {
|
public void initialize() {
|
||||||
UnifiedArgumentCollection uac = new UnifiedArgumentCollection();
|
UnifiedArgumentCollection uac = new UnifiedArgumentCollection();
|
||||||
uac.CONFIDENCE_THRESHOLD = confidence;
|
uac.CONFIDENCE_THRESHOLD = confidence;
|
||||||
uac.ALL_BASES = true;
|
uac.ALL_BASES_MODE = true;
|
||||||
UG = new UnifiedGenotyperEngine(getToolkit(), uac);
|
UG = new UnifiedGenotyperEngine(getToolkit(), uac);
|
||||||
|
|
||||||
out.println("#locus\tid\tdownsampled_coverage\tpct_coverage\titeration\tref\teval_call\tcomp_call\tvariant_concordance\tgenotype_concordance");
|
out.println("#locus\tid\tdownsampled_coverage\tpct_coverage\titeration\tref\teval_call\tcomp_call\tvariant_concordance\tgenotype_concordance");
|
||||||
|
|
|
||||||
|
|
@ -32,7 +32,7 @@ public class SecondaryBaseTransitionTableWalker extends LocusWalker<Integer, Int
|
||||||
public void initialize() {
|
public void initialize() {
|
||||||
UnifiedArgumentCollection uac = new UnifiedArgumentCollection();
|
UnifiedArgumentCollection uac = new UnifiedArgumentCollection();
|
||||||
uac.CONFIDENCE_THRESHOLD = 50;
|
uac.CONFIDENCE_THRESHOLD = 50;
|
||||||
uac.ALL_BASES = true;
|
uac.ALL_BASES_MODE = true;
|
||||||
ug = new UnifiedGenotyperEngine(getToolkit(), uac);
|
ug = new UnifiedGenotyperEngine(getToolkit(), uac);
|
||||||
|
|
||||||
altTable = new NamedTable();
|
altTable = new NamedTable();
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue