From 3a341b2f06d1a39c0efef5a5991566f39c0649da Mon Sep 17 00:00:00 2001 From: depristo Date: Fri, 18 Sep 2009 21:01:43 +0000 Subject: [PATCH] Fixes for VariantEval for genotyping mode git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1659 348d0f76-0448-11de-a6fe-93d51630548a --- .../broadinstitute/sting/gatk/refdata/RodGeliText.java | 9 +++++++-- .../walkers/varianteval/CallableBasesAnalysis.java | 10 +++++----- .../gatk/walkers/varianteval/VariantCounter.java | 9 +++++++-- python/MergeBAMBatch.py | 6 ++++++ python/MergeBAMsUtils.py | 3 ++- 5 files changed, 27 insertions(+), 10 deletions(-) diff --git a/java/src/org/broadinstitute/sting/gatk/refdata/RodGeliText.java b/java/src/org/broadinstitute/sting/gatk/refdata/RodGeliText.java index 9fc892573..43e9dbe30 100755 --- a/java/src/org/broadinstitute/sting/gatk/refdata/RodGeliText.java +++ b/java/src/org/broadinstitute/sting/gatk/refdata/RodGeliText.java @@ -68,7 +68,12 @@ public class RodGeliText extends BasicReferenceOrderedDatum implements Variation refBase = parts[2].charAt(0); depth = Integer.valueOf(parts[3]); maxMappingQuality = Integer.valueOf(parts[4]); - bestGenotype = parts[5]; + + System.out.printf("%s%n", parts[5]); + char[] x = parts[5].toUpperCase().toCharArray(); + Arrays.sort(x); + bestGenotype = new String(x); + lodBtr = Double.valueOf(parts[6]); lodBtnb = Double.valueOf(parts[7]); @@ -349,7 +354,7 @@ public class RodGeliText extends BasicReferenceOrderedDatum implements Variation */ @Override public Genotype getGenotype(DiploidGenotype x) { - if (x.toString() != this.getAltBasesFWD()) throw new IllegalStateException("We don't contain that genotype!"); + if (x.toString() != this.getAltBasesFWD()) throw new IllegalStateException("We don't contain genotype " + x); return new BasicGenotype(getLocation(), x.toString(), refBase, lodBtnb); } diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/varianteval/CallableBasesAnalysis.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/varianteval/CallableBasesAnalysis.java index c11825254..2ccf4da4d 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/varianteval/CallableBasesAnalysis.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/varianteval/CallableBasesAnalysis.java @@ -69,14 +69,14 @@ public class CallableBasesAnalysis extends BasicVariantAnalysis implements Genot // For every threshold, updated discoverable and callable for (int i = 0; i < thresholds.length; i++) { double threshold = thresholds[i]; - DiploidGenotype g = DiploidGenotype.createHomGenotype(ref); - Genotype genotype = ((VariantBackedByGenotype) eval).getGenotype(g); + Genotype genotype = ((VariantBackedByGenotype)eval).getGenotypes().get(0); + // update discoverable - if (eval.isSNP() && eval.getNegLog10PError() >= threshold) + if ( eval.isSNP() && eval.getNegLog10PError() >= threshold) discoverable_bases[i]++; - if (!eval.isSNP() && genotype.getNegLog10PError() >= threshold) + if ( eval.isReference() && genotype.getNegLog10PError() >= threshold) discoverable_bases[i]++; - if (genotype.getNegLog10PError() >= threshold) + if ( genotype.getNegLog10PError() >= threshold) genotypable_bases[i]++; } diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/varianteval/VariantCounter.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/varianteval/VariantCounter.java index ca182da84..d4b50f884 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/varianteval/VariantCounter.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/varianteval/VariantCounter.java @@ -5,6 +5,7 @@ import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker; import org.broadinstitute.sting.utils.genotype.DiploidGenotype; import org.broadinstitute.sting.utils.genotype.VariantBackedByGenotype; import org.broadinstitute.sting.utils.genotype.Variation; +import org.broadinstitute.sting.utils.genotype.Genotype; import java.util.ArrayList; import java.util.List; @@ -31,8 +32,12 @@ public class VariantCounter extends BasicVariantAnalysis implements GenotypeAnal public String update(Variation eval, RefMetaDataTracker tracker, char ref, AlignmentContext context) { nSNPs += eval == null ? 0 : 1; - if ( this.getMaster().evalContainsGenotypes && eval != null && eval.isSNP() && ((VariantBackedByGenotype)eval).getGenotype( DiploidGenotype.valueOf(eval.getAlternateBase())).isHet() ) - nHets++; + if ( this.getMaster().evalContainsGenotypes && eval != null ) { + List genotypes = ((VariantBackedByGenotype)eval).getGenotypes(); + if ( eval.isSNP() && eval.isBiallelic() && genotypes.get(0).isHet() ) { + nHets++; + } + } return null; } diff --git a/python/MergeBAMBatch.py b/python/MergeBAMBatch.py index 5dfe55de5..e2936a84a 100755 --- a/python/MergeBAMBatch.py +++ b/python/MergeBAMBatch.py @@ -37,6 +37,9 @@ if __name__ == "__main__": parser.add_option("-n", "--naIDPops", dest="NAIDS2POP", type="string", default=None, help="Path to file contains NAID POP names. If provided, input files will be merged by population") + parser.add_option("-p", "--pop", dest="onlyPop", + type="string", default=None, + help="If provided, only this population will be processed") (OPTIONS, args) = parser.parse_args() if len(args) != 1: @@ -58,6 +61,9 @@ if __name__ == "__main__": allSources = reduce( operator.__add__, map( glob.glob, s[1:] ), [] ) print 'Merging info:' for spec in splitSourcesByPopulation(allSources, merged_filename_base, NAID2Pop): + if OPTIONS.onlyPop <> None and spec.group() <> OPTIONS.onlyPop: + continue + spec.setPath(directory) spec.pprint() diff --git a/python/MergeBAMsUtils.py b/python/MergeBAMsUtils.py index d9004568b..039f23d8d 100755 --- a/python/MergeBAMsUtils.py +++ b/python/MergeBAMsUtils.py @@ -68,7 +68,8 @@ class MergeFilesSpec: print ' Population: ', self.group() print ' Merged filename: ', self.getMergedBAM() print ' N sources: ', len(self.sources()) - print ' Sources: ', self.sources() + for source in sorted(self.sources()): + print ' Source: ', source print ' Sizes: ', self.sourceSizes(humanReadable=True) print ' Est. merged size: ', greek(reduce(operator.__add__, self.sourceSizes(), 0))