From 969c995298f2b6cd17872f6887464e8709614f45 Mon Sep 17 00:00:00 2001 From: Ami Levy-Moonshine Date: Sat, 1 Dec 2012 00:08:19 -0500 Subject: [PATCH] work under development - catVariants. Changes to AssessRRQuals based on Eric todo comments. bug fix in CombineVariants --- .../gatk/walkers/qc/AssessReducedQuals.java | 44 ++++++++++++------- .../variantcontext/VariantContextUtils.java | 8 ++-- 2 files changed, 31 insertions(+), 21 deletions(-) diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/qc/AssessReducedQuals.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/qc/AssessReducedQuals.java index 78bcf1228..2c70d44e2 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/qc/AssessReducedQuals.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/qc/AssessReducedQuals.java @@ -45,7 +45,6 @@ import java.util.List; public class AssessReducedQuals extends LocusWalker implements TreeReducible { private static final String reduced = "reduced"; - private static final String original = "original"; private static final int originalQualsIndex = 0; private static final int reducedQualsIndex = 1; @@ -55,14 +54,14 @@ public class AssessReducedQuals extends LocusWalker implem @Argument(fullName = "qual_epsilon", shortName = "epsilon", doc = "when |Quals_reduced_bam - Quals_original_bam| > epsilon*Quals_original_bam we output this interval", required = false) public int qual_epsilon = 0; - @Argument(fullName = "debugLevel", shortName = "debug", doc = "debug mode on") // TODO -- best to make this optional - public int debugLevel = 0; // TODO -- best to make this an enum or boolean + @Argument(fullName = "debugLevel", shortName = "debug", doc = "debug level: NO_DEBUG, PRINT_LOCI,PRINT_PILEUPS", required = false) + public DebugLevel debugLevel = DebugLevel.NO_DEBUG; @Output protected PrintStream out; public void initialize() { - if (debugLevel != 0) + if (debugLevel != DebugLevel.NO_DEBUG) out.println(" Debug mode" + "Debug:\tsufficientQualSum: "+sufficientQualSum+ "\n " + "Debug:\tqual_epsilon: "+qual_epsilon); @@ -78,20 +77,20 @@ public class AssessReducedQuals extends LocusWalker implem boolean reportLocus; final int[] quals = getPileupQuals(context.getBasePileup()); - if (debugLevel != 0) + if (debugLevel != DebugLevel.NO_DEBUG) out.println("Debug:\tLocus Quals\t"+ref.getLocus()+"\toriginal\t"+quals[originalQualsIndex]+"\treduced\t"+quals[reducedQualsIndex]); final int epsilon = MathUtils.fastRound(quals[originalQualsIndex]*qual_epsilon); final int calcOriginalQuals = Math.min(quals[originalQualsIndex],sufficientQualSum); final int calcReducedQuals = Math.min(quals[reducedQualsIndex],sufficientQualSum); final int OriginalReducedQualDiff = calcOriginalQuals - calcReducedQuals; reportLocus = OriginalReducedQualDiff > epsilon || OriginalReducedQualDiff < -1*epsilon; - if(debugLevel != 0 && reportLocus) + if(debugLevel != DebugLevel.NO_DEBUG && reportLocus) out.println("Debug:\tEmited Locus\t"+ref.getLocus()+"\toriginal\t"+quals[originalQualsIndex]+"\treduced\t"+quals[reducedQualsIndex]+"\tepsilon\t"+epsilon+"\tdiff\t"+OriginalReducedQualDiff); return reportLocus ? ref.getLocus() : null; } - private final int[] getPileupQuals(final ReadBackedPileup readPileup) { + private int[] getPileupQuals(final ReadBackedPileup readPileup) { final int[] quals = new int[2]; String[] printPileup = {"Debug 2:\toriginal pileup:\t"+readPileup.getLocation()+"\nDebug 2:----------------------------------\n", @@ -99,30 +98,29 @@ public class AssessReducedQuals extends LocusWalker implem for( PileupElement p : readPileup ){ final List tags = getToolkit().getReaderIDForRead(p.getRead()).getTags().getPositionalTags(); - if ( isGoodRead(p,tags) ){ + if ( isGoodRead(p) ){ final int tempQual = (int)(p.getQual()) * p.getRepresentativeCount(); final int tagIndex = getTagIndex(tags); quals[tagIndex] += tempQual; - if(debugLevel == 2) + if(debugLevel == DebugLevel.PRINT_PILEUPS) printPileup[tagIndex] += "\tDebug 2: ("+p+")\tMQ="+p.getMappingQual()+":QU="+p.getQual()+":RC="+p.getRepresentativeCount()+":OS="+p.getOffset()+"\n"; } } - if(debugLevel == 2){ + if(debugLevel == DebugLevel.PRINT_PILEUPS){ out.println(printPileup[originalQualsIndex]); out.println(printPileup[reducedQualsIndex]); } return quals; } - // TODO -- arguments/variables should be final, not method declaration - private final boolean isGoodRead(PileupElement p, List tags){ - // TODO -- this isn't quite right. You don't need the tags here. Instead, you want to check whether the read itself (which - // TODO -- you can get from the PileupElement) is a reduced read (not all reads from the reduced bam are reduced) and only - // TODO -- for them do you want to ignore that min mapping quality cutoff (but you *do* still want the min base cutoff). - return !p.isDeletion() && (tags.contains(reduced) || (tags.contains(original) && (int)p.getQual() >= 20 && p.getMappingQual() >= 20)); + + private boolean isGoodRead(PileupElement p){ + // TODO -- You want to check whether the read itself is a reduced read and only + // TODO -- for them you want to ignore that min mapping quality cutoff (but you *do* still want the min base cutoff). + return !p.isDeletion() && ((p.getRead().isReducedRead()) || (!p.getRead().isReducedRead() && (int)p.getQual() >= 20 && p.getMappingQual() >= 20)); } - private final int getTagIndex(List tags){ + private int getTagIndex(List tags){ return tags.contains(reduced) ? 1 : 0; } @@ -170,4 +168,16 @@ public class AssessReducedQuals extends LocusWalker implem out.println(sum); return value; } + + public enum DebugLevel { + NO_DEBUG, + /** + * Print locus level information (such as locus quals) and loci with unmatch quals + */ + PRINT_LOCI, + /** + * Print the pileup infomarion of the reduced bam files and the original bam files + */ + PRINT_PILEUPS + } } diff --git a/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextUtils.java b/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextUtils.java index 614b234e9..3b6bd0182 100755 --- a/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextUtils.java +++ b/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextUtils.java @@ -451,7 +451,7 @@ public class VariantContextUtils { if ( unsortedVCs == null || unsortedVCs.size() == 0 ) return null; - if ( annotateOrigin && priorityListOfVCs == null ) + if ( annotateOrigin && priorityListOfVCs == null && genotypeMergeOptions == GenotypeMergeType.PRIORITIZE) throw new IllegalArgumentException("Cannot merge calls and annotate their origins without a complete priority list of VariantContexts"); if ( genotypeMergeOptions == GenotypeMergeType.REQUIRE_UNIQUE ) @@ -597,7 +597,7 @@ public class VariantContextUtils { if ( annotateOrigin ) { // we care about where the call came from String setValue; - if ( nFiltered == 0 && variantSources.size() == priorityListOfVCs.size() ) // nothing was unfiltered + if ( nFiltered == 0 && variantSources.size() == preFilteredVCs.size() ) // nothing was unfiltered setValue = MERGE_INTERSECTION; else if ( nFiltered == VCs.size() ) // everything was filtered out setValue = MERGE_FILTER_IN_ALL; @@ -840,9 +840,9 @@ public class VariantContextUtils { if ( mergeOption == GenotypeMergeType.PRIORITIZE && priorityListOfVCs == null ) throw new IllegalArgumentException("Cannot merge calls by priority with a null priority list"); - if ( mergeOption == GenotypeMergeType.UNSORTED ){ + if ( mergeOption != GenotypeMergeType.PRIORITIZE ){ if (priorityListOfVCs != null ) - logger.info("Priority string was provided but is not used since GenotypeMergeType is UNSORTED"); + logger.info("Priority string was provided but is not used since GenotypeMergeType is not PRIORITIZE"); return new ArrayList(unsortedVCs); } else {