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

work under development - catVariants. Changes to AssessRRQuals based on Eric todo comments. bug fix in CombineVariants

This commit is contained in:
Ami Levy-Moonshine 2012-12-01 00:08:19 -05:00
parent 4714ccc284
commit 969c995298
2 changed files with 31 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
public/java/src/org/broadinstitute/sting/gatk/walkers/qc/AssessReducedQuals.java
View File

@ -45,7 +45,6 @@ import java.util.List;
public class AssessReducedQuals extends LocusWalker<GenomeLoc, GenomeLoc> implements TreeReducible<GenomeLoc> {
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<GenomeLoc, GenomeLoc> 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<GenomeLoc, GenomeLoc> 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<GenomeLoc, GenomeLoc> implem
for( PileupElement p : readPileup ){
final List<String> 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<String> 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<String> tags){
private int getTagIndex(List<String> tags){
return tags.contains(reduced) ? 1 : 0;
}
@ -170,4 +168,16 @@ public class AssessReducedQuals extends LocusWalker<GenomeLoc, GenomeLoc> 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
}
}

8
public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextUtils.java
View File

@ -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<VariantContext>(unsortedVCs);
}
else {