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

First part of a two-stage commit. Removing old VariantEval to make room for VariantEval 3.0 in core. 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@5137 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
kiran 2011-01-31 17:03:41 +00:00
parent d406d9b3fc
commit b0432ee1e2
17 changed files with 0 additions and 4311 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

107
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/CompOverlap.java
View File

@ -1,107 +0,0 @@
package org.broadinstitute.sting.gatk.walkers.varianteval;
import org.broad.tribble.util.variantcontext.Allele;
import org.broad.tribble.util.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.refdata.*;
import org.broadinstitute.sting.utils.report.tags.Analysis;
import org.broadinstitute.sting.utils.report.tags.DataPoint;
/**
* The Broad Institute
* SOFTWARE COPYRIGHT NOTICE AGREEMENT
* This software and its documentation are copyright 2009 by the
* Broad Institute/Massachusetts Institute of Technology. All rights are reserved.
* <p/>
* This software is supplied without any warranty or guaranteed support whatsoever. Neither
* the Broad Institute nor MIT can be responsible for its use, misuse, or functionality.
*/
@Analysis(name = "Comp Overlap", description = "the overlap between eval and comp sites")
public class CompOverlap extends VariantEvaluator implements StandardEval {
@DataPoint(name = "eval sites", description = "number of eval SNP sites")
long nEvalSNPs = 0;
@DataPoint(name = "comp sites", description = "number of comp SNP sites")
long nCompSNPs = 0;
@DataPoint(name = "evals not at comp", description = "number of eval sites outside of comp sites")
long novelSites = 0;
@DataPoint(name = "evals at comp", description = "number of eval sites at comp sites")
long nSNPsAtComp = 0;
@DataPoint(name = "% evals at comp", description = "percentage of eval sites at comp sites")
double compRate = 0.0;
@DataPoint(name = "concordant", description = "number of concordant sites")
long nConcordant = 0;
@DataPoint(name = "% concordant", description = "the concordance rate")
double concordantRate = 0.0;
private boolean expectingIndels = false;
public CompOverlap(VariantEvalWalker parent) {
super(parent);
expectingIndels = parent.dels;
}
public String getName() {
return "compOverlap";
}
public int getComparisonOrder() {
return 2; // we need to see each eval track and each comp track
}
public long nNovelSites() { return nEvalSNPs - nSNPsAtComp; }
public double compRate() { return rate(nSNPsAtComp, nEvalSNPs); }
public double concordanceRate() { return rate(nConcordant, nSNPsAtComp); }
public void finalizeEvaluation() {
compRate = 100 * compRate();
concordantRate = 100 * concordanceRate();
novelSites = nNovelSites();
}
public boolean enabled() {
return true;
}
/**
* Returns true if every allele in eval is also in comp
*
* @param eval eval context
* @param comp db context
* @return true if eval and db are discordant
*/
public boolean discordantP(VariantContext eval, VariantContext comp) {
for (Allele a : eval.getAlleles()) {
if (!comp.hasAllele(a, true))
return true;
}
return false;
}
public String update2(VariantContext eval, VariantContext comp, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
boolean compIsGood = expectingIndels ? comp != null && comp.isNotFiltered() && comp.isIndel() : comp != null && comp.isNotFiltered() && comp.isSNP() ;
boolean evalIsGood = expectingIndels ? eval != null && eval.isIndel() : eval != null && eval.isSNP() ;
if ( compIsGood ) nCompSNPs++; // count the number of comp events
if (evalIsGood) nEvalSNPs++; // count the number of eval events
if (compIsGood && evalIsGood) {
nSNPsAtComp++;
if (!discordantP(eval, comp)) // count whether we're concordant or not with the comp value
nConcordant++;
}
return null; // we don't capture any interesting sites
}
}

77
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/CountFunctionalClasses.java
View File

@ -1,77 +0,0 @@
package org.broadinstitute.sting.gatk.walkers.varianteval;
import org.broad.tribble.util.variantcontext.VariantContext;
import org.broadinstitute.sting.utils.report.tags.Analysis;
import org.broadinstitute.sting.utils.report.tags.DataPoint;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
@Analysis(name = "Count Functional Classes", description = "Counts instances of different functional variant classes (provided the variants are annotated with that information)")
public class CountFunctionalClasses extends VariantEvaluator {
// the following fields are in output order:
@DataPoint(description = "miRNA")
long nMiRNA= 0;
@DataPoint(description = "3'-UTR")
long nUTR3 = 0;
@DataPoint(description = "Intron")
long nIntron = 0;
@DataPoint(description = "Splice-site")
long nSpliceSite= 0;
@DataPoint(description = "Read-through")
long nReadThrough = 0;
@DataPoint(description = "Nonsense")
long nNonsense = 0;
@DataPoint(description = "Missense")
long nMissense = 0;
@DataPoint(description = "Synonymous")
long nSynonymous = 0;
@DataPoint(description = "5'-UTR")
long nUTR5= 0;
@DataPoint(description = "Promoter")
long nPromoter = 0;
public CountFunctionalClasses(VariantEvalWalker parent) {
super(parent);
}
public String getName() {
return "functionalclasses";
}
public boolean enabled() {
return false;
}
public int getComparisonOrder() {
return 1;
}
public String update1(VariantContext vc1, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
String type = vc1.getAttributeAsString("type");
if (type != null) {
if (type.equalsIgnoreCase("miRNA")) { nMiRNA++; }
else if (type.equalsIgnoreCase("3'-UTR")) { nUTR3++; }
else if (type.equalsIgnoreCase("Intron")) { nIntron++; }
else if (type.equalsIgnoreCase("Splice_site")) { nSpliceSite++; }
else if (type.equalsIgnoreCase("Read-through")) { nReadThrough++; }
else if (type.equalsIgnoreCase("Nonsense")) { nNonsense++; }
else if (type.equalsIgnoreCase("Missense")) { nMissense++; }
else if (type.equalsIgnoreCase("Synonymous")) { nSynonymous++; }
else if (type.equalsIgnoreCase("5'-UTR")) { nUTR5++; }
else if (type.equalsIgnoreCase("Promoter")) { nPromoter++; }
}
return null;
}
}

150
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/CountVariants.java
View File

@ -1,150 +0,0 @@
package org.broadinstitute.sting.gatk.walkers.varianteval;
import org.broad.tribble.util.variantcontext.Genotype;
import org.broad.tribble.util.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.utils.report.tags.Analysis;
import org.broadinstitute.sting.utils.report.tags.DataPoint;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
@Analysis(name = "Count Variants", description = "Counts different classes of variants in the sample")
public class CountVariants extends VariantEvaluator implements StandardEval {
// the following fields are in output order:
// basic counts on various rates found
@DataPoint(description = "Number of processed loci")
long nProcessedLoci = 0;
@DataPoint(description = "Number of called loci")
long nCalledLoci = 0;
@DataPoint(description = "Number of reference loci")
long nRefLoci = 0;
@DataPoint(description = "Number of variant loci")
long nVariantLoci = 0;
// the following two calculations get set in the finalizeEvaluation
@DataPoint(description = "Variants per loci rate")
double variantRate = 0;
@DataPoint(description = "Number of variants per base")
double variantRatePerBp = 0;
@DataPoint(description = "Number of snp loci")
long nSNPs = 0;
@DataPoint(description = "Number of insertions")
long nInsertions = 0;
@DataPoint(description = "Number of deletions")
long nDeletions = 0;
@DataPoint(description = "Number of complex loci")
long nComplex = 0;
@DataPoint(description = "Number of no calls loci")
long nNoCalls = 0;
@DataPoint(description = "Number of het loci")
long nHets = 0;
@DataPoint(description = "Number of hom ref loci")
long nHomRef = 0;
@DataPoint(description = "Number of hom var loci")
long nHomVar = 0;
// calculations that get set in the finalizeEvaluation method
@DataPoint(description = "heterozygosity per locus rate")
double heterozygosity = 0;
@DataPoint(description = "heterozygosity per base pair")
double heterozygosityPerBp = 0;
@DataPoint(description = "heterozygosity to homozygosity ratio")
double hetHomRatio = 0;
@DataPoint(description = "indel rate (insertion count + deletion count)")
double indelRate = 0;
@DataPoint(description = "indel rate per base pair")
double indelRatePerBp = 0;
@DataPoint(description = "deletion to insertion ratio")
double deletionInsertionRatio = 0;
public CountVariants(VariantEvalWalker parent) {
super(parent);
}
private double perLocusRate(long n) {
return rate(n, nProcessedLoci);
}
private long perLocusRInverseRate(long n) {
return inverseRate(n, nProcessedLoci);
}
public String getName() {
return "counter";
}
public boolean enabled() {
return true;
}
public int getComparisonOrder() {
return 1; // we only need to see each eval track
}
public void update0(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
nProcessedLoci += context.getSkippedBases() + (ref == null ? 0 : 1);
}
public String update1(VariantContext vc1, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
//nProcessedLoci++;
nCalledLoci++;
if (vc1.isVariant()) nVariantLoci++;
switch (vc1.getType()) {
case NO_VARIATION:
nRefLoci++;
break;
case SNP:
nSNPs++;
break;
case INDEL:
if (vc1.isInsertion()) nInsertions++;
else nDeletions++;
break;
case MIXED:
nComplex++;
break;
default:
throw new ReviewedStingException("Unexpected VariantContext type " + vc1.getType());
}
for (Genotype g : vc1.getGenotypes().values()) {
switch (g.getType()) {
case NO_CALL:
nNoCalls++;
break;
case HOM_REF:
nHomRef++;
break;
case HET:
nHets++;
break;
case HOM_VAR:
nHomVar++;
break;
default:
throw new ReviewedStingException("BUG: Unexpected genotype type: " + g);
}
}
return null; // we don't capture any interesting sites
}
public void finalizeEvaluation() {
variantRate = perLocusRate(nVariantLoci);
variantRatePerBp = perLocusRInverseRate(nVariantLoci);
heterozygosity = perLocusRate(nHets);
heterozygosityPerBp = perLocusRInverseRate(nHets);
hetHomRatio = ratio(nHets, nHomVar);
indelRate = perLocusRate(nDeletions + nInsertions);
indelRatePerBp = perLocusRInverseRate(nDeletions + nInsertions);
deletionInsertionRatio = ratio(nDeletions, nInsertions);
}
}

789
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/GenotypeConcordance.java
View File

@ -1,789 +0,0 @@
package org.broadinstitute.sting.gatk.walkers.varianteval;
import org.apache.commons.lang.ArrayUtils;
import org.broad.tribble.util.variantcontext.Allele;
import org.broad.tribble.util.variantcontext.Genotype;
import org.broad.tribble.util.variantcontext.VariantContext;
import org.broad.tribble.vcf.VCFConstants;
import org.broadinstitute.sting.gatk.contexts.*;
import org.broadinstitute.sting.gatk.refdata.*;
import org.broadinstitute.sting.utils.exceptions.StingException;
import org.broadinstitute.sting.utils.report.tags.Analysis;
import org.broadinstitute.sting.utils.report.tags.DataPoint;
import org.broadinstitute.sting.utils.report.utils.TableType;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.utils.exceptions.UserException;
import java.util.*;
/*
* Copyright (c) 2010 The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
* THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
@Analysis(name = "Genotype Concordance", description = "Determine the genotype concordance between the genotypes in difference tracks")
public class GenotypeConcordance extends VariantEvaluator implements StandardEval {
private static final boolean PRINT_INTERESTING_SITES = true;
protected final static Logger logger = Logger.getLogger(GenotypeConcordance.class);
// a mapping from allele count to stats
@DataPoint(description = "the frequency statistics for each allele")
FrequencyStats alleleFreqStats = new FrequencyStats();
// a mapping from sample to stats
@DataPoint(name="samples", description = "the concordance statistics for each sample")
SampleStats sampleStats = null;
// a mapping from sample to stats summary
@DataPoint(name="summary", description = "the concordance statistics summary for each sample")
SampleSummaryStats sampleSummaryStats = null;
// two histograms of variant quality scores, for true positive and false positive calls
@DataPoint(description = "the variant quality score histograms for true positive and false positive calls")
QualityScoreHistograms qualityScoreHistograms = null;
@DataPoint(description = "the concordance statistics summary by allele count")
ACSummaryStats alleleCountSummary = null;
@DataPoint(description = "the concordance statistics by allele count")
ACStats alleleCountStats = null;
private static final int MAX_MISSED_VALIDATION_DATA = 100;
private boolean discordantInteresting = false;
private VariantEvalWalker.EvaluationContext group = null;
static class FrequencyStats implements TableType {
class Stats {
public Stats(int found, int missed) { nFound = found; nMissed = missed; }
public long nFound = 0;
public long nMissed = 0;
}
public HashMap<Integer, Stats> foundMissedByAC = new HashMap<Integer, Stats>();
public Object[] getRowKeys() {
String rows[] = new String[foundMissedByAC.size()];
int index = 0;
for (int i : foundMissedByAC.keySet()) rows[index++] = "Allele Count " + i;
return rows;
}
public Object[] getColumnKeys() {
return new String[]{"number_found", "number_missing"};
}
public String getName() {
return "FrequencyStats";
}
public String getCell(int x, int y) {
if (x >= foundMissedByAC.size()) throw new IllegalStateException(x + " is greater than the max index of " + (foundMissedByAC.size()-1));
if (y == 0) return String.valueOf(foundMissedByAC.get(foundMissedByAC.keySet().toArray(new Integer[foundMissedByAC.size()])[x]).nFound);
else return String.valueOf(foundMissedByAC.get(foundMissedByAC.keySet().toArray(new Integer[foundMissedByAC.size()])[x]).nMissed);
}
public void incrementFoundCount(int alleleFreq) {
if (!foundMissedByAC.containsKey(alleleFreq))
foundMissedByAC.put(alleleFreq,new Stats(1,0));
else
foundMissedByAC.get(alleleFreq).nFound++;
}
public void incrementMissedCount(int alleleFreq) {
if (!foundMissedByAC.containsKey(alleleFreq))
foundMissedByAC.put(alleleFreq,new Stats(0,1));
else
foundMissedByAC.get(alleleFreq).nMissed++;
}
}
static class QualityScoreHistograms implements TableType {
final static int NUM_BINS = 20;
final HashMap<Integer,Integer> truePositiveQualityScoreMap = new HashMap<Integer,Integer>(); // A HashMap holds all the quality scores until we are able to bin them appropriately
final HashMap<Integer,Integer> falsePositiveQualityScoreMap = new HashMap<Integer,Integer>();
final int truePositiveHist[] = new int[NUM_BINS]; // the final histograms that get reported out
final int falsePositiveHist[] = new int[NUM_BINS];
final String[] rowKeys = new String[]{"true_positive_hist", "false_positive_hist"};
public Object[] getRowKeys() {
return rowKeys;
}
public Object[] getColumnKeys() {
final String columnKeys[] = new String[NUM_BINS];
for( int iii = 0; iii < NUM_BINS; iii++ ) {
columnKeys[iii] = "histBin" + iii;
}
return columnKeys;
}
public String getName() {
return "QualityScoreHistogram";
}
public String getCell(int x, int y) {
if( x == 0 ) {
return String.valueOf(truePositiveHist[y]);
} else if ( x == 1 ) {
return String.valueOf(falsePositiveHist[y]);
} else {
throw new ReviewedStingException( "Unknown row in " + getName() + ", row = " + x );
}
}
public String toString() {
String returnString = "";
// output both histogram arrays
returnString += "TP: ";
for( int iii = 0; iii < NUM_BINS; iii++ ) {
returnString += truePositiveHist[iii] + " ";
}
returnString += "\nFP: ";
for( int iii = 0; iii < NUM_BINS; iii++ ) {
returnString += falsePositiveHist[iii] + " ";
}
return returnString;
}
public void incrValue( final double qual, final boolean isTruePositiveCall ) {
HashMap<Integer,Integer> qualScoreMap;
if( isTruePositiveCall ) {
qualScoreMap = truePositiveQualityScoreMap;
} else {
qualScoreMap = falsePositiveQualityScoreMap;
}
final Integer qualKey = Math.round((float) qual);
if( qualScoreMap.containsKey(qualKey) ) {
qualScoreMap.put(qualKey, qualScoreMap.get(qualKey) + 1);
} else {
qualScoreMap.put(qualKey, 1);
}
}
public void organizeHistogramTables() {
for( int iii = 0; iii < NUM_BINS; iii++ ) {
truePositiveHist[iii] = 0;
falsePositiveHist[iii] = 0;
}
int maxQual = 0;
// Calculate the maximum quality score for both TP and FP calls in order to normalize and histogram
for( final Integer qual : truePositiveQualityScoreMap.keySet()) {
if( qual > maxQual ) {
maxQual = qual;
}
}
for( final Integer qual : falsePositiveQualityScoreMap.keySet()) {
if( qual > maxQual ) {
maxQual = qual;
}
}
final double binSize = ((double)maxQual) / ((double) (NUM_BINS-1)); //BUGBUG: should be normalized max to min, not max to 0
for( final Integer qual : truePositiveQualityScoreMap.keySet()) {
final int index = (int)Math.floor( ((double)qual) / binSize );
if(index >= 0) { //BUGBUG: problem when maxQual is zero?
truePositiveHist[ index ] += truePositiveQualityScoreMap.get(qual);
}
}
for( final Integer qual : falsePositiveQualityScoreMap.keySet()) {
final int index = (int)Math.floor( ((double)qual) / binSize );
if(index >= 0) {
falsePositiveHist[ index ] += falsePositiveQualityScoreMap.get(qual);
}
}
}
}
// keep a list of the validation data we saw before the first eval data
private HashSet<VariantContext> missedValidationData = new HashSet<VariantContext>();
public GenotypeConcordance(VariantEvalWalker parent) {
super(parent);
discordantInteresting = parent.DISCORDANT_INTERESTING;
}
public String getName() {
return "genotypeConcordance";
}
public int getComparisonOrder() {
return 2; // we need to see each eval track and each comp track
}
public boolean enabled() {
return true;
}
public String toString() {
return getName() + ": <table>";
}
private boolean warnedAboutValidationData = false;
public String update2(VariantContext eval, VariantContext validation, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context, VariantEvalWalker.EvaluationContext group) {
this.group = group;
String interesting = null;
// sanity check that we at least have either eval or validation data
if (eval == null && !isValidVC(validation)) {
return interesting;
}
if( qualityScoreHistograms == null ) {
qualityScoreHistograms = new QualityScoreHistograms();
}
if ( alleleCountStats == null && eval != null && validation != null ) {
alleleCountStats = new ACStats(eval,validation,Genotype.Type.values().length);
alleleCountSummary = new ACSummaryStats(eval, validation);
}
if (sampleStats == null) {
if (eval != null) {
// initialize the concordance table
sampleStats = new SampleStats(eval,Genotype.Type.values().length);
sampleSummaryStats = new SampleSummaryStats(eval);
for (final VariantContext vc : missedValidationData) {
determineStats(null, vc);
}
missedValidationData = null;
} else {
// todo -- Eric, this results in a memory problem when eval is WEx data but you are using CG calls genome-wide
// todo -- perhaps you need should extend the evaluators with an initialize
// todo -- method that gets the header (or samples) for the first eval sites?
if (missedValidationData.size() > MAX_MISSED_VALIDATION_DATA) {
if (!warnedAboutValidationData) {
//logger.warn("Too many genotype sites missed before eval site appeared; ignoring");
warnedAboutValidationData = true;
}
} else {
missedValidationData.add(validation);
}
return interesting;
}
}
interesting = determineStats(eval, validation);
return interesting; // we don't capture any interesting sites
}
private String determineStats(final VariantContext eval, final VariantContext validation) {
String interesting = null;
final boolean validationIsValidVC = isValidVC(validation);
final String evalAC = ( vcHasGoodAC(eval) ) ? String.format("evalAC%d",getAC(eval)) : null ;
final String validationAC = ( vcHasGoodAC(validation) ) ? String.format("compAC%d",getAC(validation)) : null;
// determine concordance for eval data
if (eval != null) {
for (final String sample : eval.getGenotypes().keySet()) {
final Genotype.Type called = eval.getGenotype(sample).getType();
final Genotype.Type truth;
if (!validationIsValidVC || !validation.hasGenotype(sample)) {
truth = Genotype.Type.NO_CALL;
} else {
truth = validation.getGenotype(sample).getType();
// interesting = "ConcordanceStatus=FP";
if (discordantInteresting && truth.ordinal() != called.ordinal())
{
interesting = "ConcordanceStatus=" + truth + "/" + called;
}
}
sampleStats.incrValue(sample, truth, called);
if ( evalAC != null && validationAC != null) {
alleleCountStats.incrValue(evalAC,truth,called);
alleleCountStats.incrValue(validationAC,truth,called);
}
}
}
// otherwise, mark no-calls for all samples
else {
final Genotype.Type called = Genotype.Type.NO_CALL;
for (final String sample : validation.getGenotypes().keySet()) {
final Genotype.Type truth = validation.getGenotype(sample).getType();
sampleStats.incrValue(sample, truth, called);
if ( evalAC != null ) {
alleleCountStats.incrValue(evalAC,truth,called);
}
// print out interesting sites
if ( PRINT_INTERESTING_SITES && super.getVEWalker().gcLog != null ) {
if ( (truth == Genotype.Type.HOM_VAR || truth == Genotype.Type.HET) && called == Genotype.Type.NO_CALL ) {
super.getVEWalker().gcLog.printf("%s FN %s%n", group, validation);
}
if ( (called == Genotype.Type.HOM_VAR || called == Genotype.Type.HET) && truth == Genotype.Type.HOM_REF ) {
super.getVEWalker().gcLog.printf("%s FP %s%n", group, validation);
}
}
}
}
// determine allele count concordance () // this is really a FN rate estimate -- CH
if (validationIsValidVC && validation.isPolymorphic()) {
int trueAlleleCount = 0;
for (final Allele a : validation.getAlternateAlleles()) {
trueAlleleCount += validation.getChromosomeCount(a);
}
if (eval != null) {
alleleFreqStats.incrementFoundCount(trueAlleleCount);
} else {
alleleFreqStats.incrementMissedCount(trueAlleleCount);
}
}
// TP & FP quality score histograms
if( eval != null && eval.isPolymorphic() && validationIsValidVC ) {
if( eval.getGenotypes().keySet().size() == 1 ) { // single sample calls
for( final String sample : eval.getGenotypes().keySet() ) { // only one sample
if( validation.hasGenotype(sample) ) {
final Genotype truth = validation.getGenotype(sample);
qualityScoreHistograms.incrValue( eval.getPhredScaledQual(), !truth.isHomRef() );
}
}
} else { // multi sample calls
qualityScoreHistograms.incrValue( eval.getPhredScaledQual(), validation.isPolymorphic() );
}
}
return interesting;
}
private static boolean isValidVC(final VariantContext vc) {
return (vc != null && !vc.isFiltered());
}
public void finalizeEvaluation() {
if( qualityScoreHistograms != null ) {
qualityScoreHistograms.organizeHistogramTables();
}
if( sampleSummaryStats != null && sampleStats != null ) {
sampleSummaryStats.generateSampleSummaryStats( sampleStats );
}
if ( alleleCountSummary != null && alleleCountStats != null ) {
alleleCountSummary.generateSampleSummaryStats( alleleCountStats );
}
}
private boolean vcHasGoodAC(VariantContext vc) {
return ( vc != null && vc.getAlternateAlleles().size() == 1 && vc.hasAttribute(VCFConstants.ALLELE_COUNT_KEY) );
}
private int getAC(VariantContext vc) {
if ( List.class.isAssignableFrom(vc.getAttribute(VCFConstants.ALLELE_COUNT_KEY).getClass()) ) {
return ((List<Integer>) vc.getAttribute(VCFConstants.ALLELE_COUNT_KEY)).get(0);
} else if ( Integer.class.isAssignableFrom(vc.getAttribute(VCFConstants.ALLELE_COUNT_KEY).getClass())) {
return (Integer) vc.getAttribute(VCFConstants.ALLELE_COUNT_KEY);
} else if ( String.class.isAssignableFrom(vc.getAttribute(VCFConstants.ALLELE_COUNT_KEY).getClass()) ) {
// two ways of parsing
String ac = (String) vc.getAttribute(VCFConstants.ALLELE_COUNT_KEY);
if ( ac.startsWith("[") ) {
return Integer.parseInt(ac.replaceAll("\\[","").replaceAll("\\]",""));
} else {
try {
return Integer.parseInt(ac);
} catch ( NumberFormatException e ) {
throw new UserException(String.format("The format of the AC field is improperly formatted: AC=%s",ac));
}
}
} else {
throw new UserException(String.format("The format of the AC field does not appear to be of integer-list or String format, class was %s",vc.getAttribute(VCFConstants.ALLELE_COUNT_KEY).getClass()));
}
}
}
/**
* a table of sample names to genotype concordance figures
*/
class SampleStats implements TableType {
private final int nGenotypeTypes;
// sample to concordance stats object
public final HashMap<String, long[][]> concordanceStats = new HashMap<String, long[][]>();
/**
*
* @return one row per sample
*/
public Object[] getRowKeys() {
return concordanceStats.keySet().toArray(new String[concordanceStats.size()]);
}
/**
* increment the specified value
* @param sample the sample name
* @param truth the truth type
* @param called the called type
*/
public void incrValue(String sample, Genotype.Type truth, Genotype.Type called) {
if ( concordanceStats.containsKey(sample) )
concordanceStats.get(sample)[truth.ordinal()][called.ordinal()]++;
else if ( called != Genotype.Type.NO_CALL )
throw new UserException.CommandLineException("Sample " + sample + " has not been seen in a previous eval; this analysis module assumes that all samples are present in each variant context");
}
/**
* get the column keys
* @return a list of objects, in this case strings, that are the column names
*/
public Object[] getColumnKeys() {
return new String[]{"total_true_ref","%_ref/ref","n_ref/no-call",
"n_ref/ref","n_ref/het","n_ref/hom",
"total_true_het","%_het/het","n_het/no-call",
"n_het/ref","n_het/het","n_het/hom",
"total_true_hom","%_hom/hom","n_hom/no-call",
"n_hom/ref","n_hom/het","n_hom/hom"};
}
public SampleStats(VariantContext vc, int nGenotypeTypes) {
this.nGenotypeTypes = nGenotypeTypes;
for (String sample : vc.getGenotypes().keySet())
concordanceStats.put(sample, new long[nGenotypeTypes][nGenotypeTypes]);
}
public SampleStats(int genotypeTypes) {
nGenotypeTypes = genotypeTypes;
}
public Object getCell(int x, int y) {
// we have three rows of 6 right now for output (rows: ref, het, hom)
Genotype.Type type = Genotype.Type.values()[(y/6)+1]; // get the row type
// save some repeat work, get the total every time
long total = 0;
Object[] rowKeys = getRowKeys();
for (int called = 0; called < nGenotypeTypes; called++) {
total += concordanceStats.get(rowKeys[x])[type.ordinal()][called];
}
// now get the cell they're interested in
switch (y % 6) {
case (0): // get the total_true for this type
return total;
case (1):
return total == 0 ? 0.0 : (100.0 * (double) concordanceStats.get(rowKeys[x])[type.ordinal()][type.ordinal()] / (double) total);
default:
return concordanceStats.get(rowKeys[x])[type.ordinal()][(y % 6) - 2];
}
}
public String getName() {
return "Sample Statistics";
}
}
/**
* Sample stats, but for AC
*/
class ACStats extends SampleStats {
private String[] rowKeys;
public ACStats(VariantContext evalvc, VariantContext compvc, int nGenotypeTypes) {
super(nGenotypeTypes);
rowKeys = new String[1+2*evalvc.getGenotypes().size()+1+2*compvc.getGenotypes().size()];
for ( int i = 0; i <= 2*evalvc.getGenotypes().size(); i++ ) { // todo -- assuming ploidy 2 here...
concordanceStats.put(String.format("evalAC%d",i),new long[nGenotypeTypes][nGenotypeTypes]);
rowKeys[i] = String.format("evalAC%d",i);
}
for ( int i = 0; i <= 2*compvc.getGenotypes().size(); i++ ) {
concordanceStats.put(String.format("compAC%d",i), new long[nGenotypeTypes][nGenotypeTypes]);
rowKeys[1+2*evalvc.getGenotypes().size()+i] = String.format("compAC%d",i);
}
}
public String getName() {
return "Allele Count Statistics";
}
public Object[] getRowKeys() {
if ( rowKeys == null ) {
throw new StingException("RowKeys is null!");
}
return rowKeys;
}
}
/**
* a table of sample names to genotype concordance summary statistics
*/
class SampleSummaryStats implements TableType {
protected final static String ALL_SAMPLES_KEY = "allSamples";
protected final static String[] COLUMN_KEYS = new String[]{
"percent_comp_ref_called_var",
"percent_comp_het_called_het",
"percent_comp_het_called_var",
"percent_comp_hom_called_hom",
"percent_comp_hom_called_var",
"percent_non-reference_sensitivity",
"percent_overall_genotype_concordance",
"percent_non-reference_discrepancy_rate"};
// sample to concordance stats object
protected final HashMap<String, double[]> concordanceSummary = new HashMap<String, double[]>();
/**
*
* @return one row per sample
*/
public Object[] getRowKeys() {
return concordanceSummary.keySet().toArray(new String[concordanceSummary.size()]);
}
/**
* get the column keys
* @return a list of objects, in this case strings, that are the column names
*/
public Object[] getColumnKeys() {
return COLUMN_KEYS;
}
public SampleSummaryStats(final VariantContext vc) {
concordanceSummary.put(ALL_SAMPLES_KEY, new double[COLUMN_KEYS.length]);
for( final String sample : vc.getGenotypes().keySet() ) {
concordanceSummary.put(sample, new double[COLUMN_KEYS.length]);
}
}
public SampleSummaryStats() {
}
public Object getCell(int x, int y) {
final Object[] rowKeys = getRowKeys();
return String.format("%.2f",concordanceSummary.get(rowKeys[x])[y]);
}
/**
* Helper routine that sums up all columns / rows found in stats specified by all pairs in d1 x d2
*
* @param stats
* @param d1
* @param d2
* @return
*/
private long sumStatsAllPairs( final long[][] stats, EnumSet<Genotype.Type> d1, EnumSet<Genotype.Type> d2 ) {
long sum = 0L;
for(final Genotype.Type e1 : d1 ) {
for(final Genotype.Type e2 : d2 ) {
sum += stats[e1.ordinal()][e2.ordinal()];
}
}
return sum;
}
private long sumStatsDiag( final long[][] stats, EnumSet<Genotype.Type> d1) {
long sum = 0L;
for(final Genotype.Type e1 : d1 ) {
sum += stats[e1.ordinal()][e1.ordinal()];
}
return sum;
}
private double ratio(long numer, long denom) {
return denom != 0L ? 100.0 * ( ((double)numer) / ((double)denom) ) : 0.0;
}
final long[] allSamplesNumerators = new long[COLUMN_KEYS.length];
final long[] allSamplesDenominators = new long[COLUMN_KEYS.length];
private void updateSummaries(int i, double[] summary, long numer, long denom ) {
allSamplesNumerators[i] += numer;
allSamplesDenominators[i] += denom;
summary[i] = ratio(numer, denom);
}
/**
* Calculate the five summary stats per sample
* @param sampleStats The Map which holds concordance values per sample
*/
public void generateSampleSummaryStats( final SampleStats sampleStats ) {
EnumSet<Genotype.Type> allVariantGenotypes = EnumSet.of(Genotype.Type.HOM_VAR, Genotype.Type.HET);
EnumSet<Genotype.Type> allCalledGenotypes = EnumSet.of(Genotype.Type.HOM_VAR, Genotype.Type.HET, Genotype.Type.HOM_REF);
EnumSet<Genotype.Type> allGenotypes = EnumSet.allOf(Genotype.Type.class);
for( final String sample : concordanceSummary.keySet() ) {
if ( sample.equals(ALL_SAMPLES_KEY) ) continue;
final long[][] stats = sampleStats.concordanceStats.get(sample);
final double[] summary = concordanceSummary.get(sample);
if( stats == null ) { throw new ReviewedStingException( "SampleStats and SampleSummaryStats contain different samples! sample = " + sample ); }
long numer, denom;
// Summary 0: % ref called as var
numer = sumStatsAllPairs(stats, EnumSet.of(Genotype.Type.HOM_REF), allVariantGenotypes);
denom = sumStatsAllPairs(stats, EnumSet.of(Genotype.Type.HOM_REF), allGenotypes);
updateSummaries(0, summary, numer, denom);
// Summary 1: % het called as het
numer = stats[Genotype.Type.HET.ordinal()][Genotype.Type.HET.ordinal()];
denom = sumStatsAllPairs(stats, EnumSet.of(Genotype.Type.HET), allGenotypes);
updateSummaries(1, summary, numer, denom);
// Summary 2: % het called as var
numer = sumStatsAllPairs(stats, EnumSet.of(Genotype.Type.HET), allVariantGenotypes);
denom = sumStatsAllPairs(stats, EnumSet.of(Genotype.Type.HET), allGenotypes);
updateSummaries(2, summary, numer, denom);
// Summary 3: % homVar called as homVar
numer = stats[Genotype.Type.HOM_VAR.ordinal()][Genotype.Type.HOM_VAR.ordinal()];
denom = sumStatsAllPairs(stats, EnumSet.of(Genotype.Type.HOM_VAR), allGenotypes);
updateSummaries(3, summary, numer, denom);
// Summary 4: % homVars called as var
numer = sumStatsAllPairs(stats, EnumSet.of(Genotype.Type.HOM_VAR), allVariantGenotypes);
denom = sumStatsAllPairs(stats, EnumSet.of(Genotype.Type.HOM_VAR), allGenotypes);
updateSummaries(4, summary, numer, denom);
// Summary 5: % non-ref called as non-ref
// MAD: this is known as the non-reference sensitivity (# non-ref according to comp found in eval / # non-ref in comp)
numer = sumStatsAllPairs(stats, allVariantGenotypes, allVariantGenotypes);
denom = sumStatsAllPairs(stats, allVariantGenotypes, allGenotypes);
updateSummaries(5, summary, numer, denom);
// Summary 6: overall genotype concordance of sites called in eval track
// MAD: this is the tradition genotype concordance
numer = sumStatsDiag(stats, allCalledGenotypes);
denom = sumStatsAllPairs(stats, allCalledGenotypes, allCalledGenotypes);
updateSummaries(6, summary, numer, denom);
// Summary 7: overall genotype concordance of sites called non-ref in eval track
long homrefConcords = stats[Genotype.Type.HOM_REF.ordinal()][Genotype.Type.HOM_REF.ordinal()];
long diag = sumStatsDiag(stats, allVariantGenotypes);
long allNoHomRef = sumStatsAllPairs(stats, allCalledGenotypes, allCalledGenotypes) - homrefConcords;
numer = allNoHomRef - diag;
denom = allNoHomRef;
updateSummaries(7, summary, numer, denom);
}
// update the final summary stats
final double[] allSamplesSummary = concordanceSummary.get(ALL_SAMPLES_KEY);
for ( int i = 0; i < allSamplesSummary.length; i++) {
allSamplesSummary[i] = ratio(allSamplesNumerators[i], allSamplesDenominators[i]);
}
}
public String getName() {
return "Sample Summary Statistics";
}
}
/**
* SampleSummaryStats .. but for allele counts
*/
class ACSummaryStats extends SampleSummaryStats {
private String[] rowKeys;
public ACSummaryStats (final VariantContext evalvc, final VariantContext compvc) {
concordanceSummary.put(ALL_SAMPLES_KEY, new double[COLUMN_KEYS.length]);
rowKeys = new String[3+2*evalvc.getGenotypes().size() + 2*compvc.getGenotypes().size()];
rowKeys[0] = ALL_SAMPLES_KEY;
for( int i = 0; i <= 2*evalvc.getGenotypes().size() ; i ++ ) {
concordanceSummary.put(String.format("evalAC%d",i), new double[COLUMN_KEYS.length]);
rowKeys[i+1] = String.format("evalAC%d",i);
}
for( int i = 0; i <= 2*compvc.getGenotypes().size() ; i ++ ) {
concordanceSummary.put(String.format("compAC%d",i), new double[COLUMN_KEYS.length]);
rowKeys[2+2*evalvc.getGenotypes().size()+i] = String.format("compAC%d",i);
}
}
public String getName() {
return "Allele Count Summary Statistics";
}
public Object[] getRowKeys() {
if ( rowKeys == null) {
throw new StingException("rowKeys is null!!");
}
return rowKeys;
}
}
class CompACNames implements Comparator{
final Logger myLogger;
private boolean info = true;
public CompACNames(Logger l) {
myLogger = l;
}
public boolean equals(Object o) {
return ( o.getClass() == CompACNames.class );
}
public int compare(Object o1, Object o2) {
if ( info ) {
myLogger.info("Sorting AC names");
info = false;
}
//System.out.printf("Objects %s %s get ranks %d %d%n",o1.toString(),o2.toString(),getRank(o1),getRank(o2));
return getRank(o1) - getRank(o2);
}
public int getRank(Object o) {
if ( o.getClass() != String.class ) {
return Integer.MIN_VALUE/4;
} else {
String s = (String) o;
if ( s.startsWith("eval") ) {
return Integer.MIN_VALUE/4 + 1 + parseAC(s);
} else if ( s.startsWith("comp") ) {
return 1+ parseAC(s);
} else {
return Integer.MIN_VALUE/4;
}
}
}
public int parseAC(String s) {
String[] g = s.split("AC");
return Integer.parseInt(g[1]);
}
}

419
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/GenotypePhasingEvaluator.java
View File

@ -1,419 +0,0 @@
package org.broadinstitute.sting.gatk.walkers.varianteval;
import org.broad.tribble.util.variantcontext.Genotype;
import org.broad.tribble.util.variantcontext.VariantContext;
import org.broad.tribble.vcf.VCFConstants;
import org.broadinstitute.sting.gatk.contexts.*;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContextUtils;
import org.broadinstitute.sting.gatk.refdata.*;
import org.broadinstitute.sting.gatk.walkers.phasing.*;
import org.broadinstitute.sting.utils.report.tags.Analysis;
import org.broadinstitute.sting.utils.report.tags.DataPoint;
import org.broadinstitute.sting.utils.report.utils.TableType;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.utils.MathUtils;
import java.util.*;
/*
* Copyright (c) 2010 The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
* THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
@Analysis(name = "Genotype Phasing Evaluation", description = "Evaluates the phasing of genotypes in different tracks")
public class GenotypePhasingEvaluator extends VariantEvaluator {
protected final static Logger logger = Logger.getLogger(GenotypePhasingEvaluator.class);
// a mapping from sample to stats
@DataPoint(name = "samples", description = "the phasing statistics for each sample")
SamplePhasingStatistics samplePhasingStatistics = null;
SamplePreviousGenotypes samplePrevGenotypes = null;
public GenotypePhasingEvaluator(VariantEvalWalker parent) {
super(parent);
this.samplePhasingStatistics = new SamplePhasingStatistics(getVEWalker().minPhaseQuality);
this.samplePrevGenotypes = new SamplePreviousGenotypes();
}
public String getName() {
return "GenotypePhasingEvaluator";
}
public int getComparisonOrder() {
return 2; // we only need to see pairs of (comp, eval)
}
public boolean enabled() {
return true;
}
public String toString() {
return getName() + ": <table>";
}
public String update2(VariantContext eval, VariantContext comp, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context, VariantEvalWalker.EvaluationContext group) {
Reasons interesting = new Reasons();
if (ref == null)
return interesting.toString();
GenomeLoc curLocus = ref.getLocus();
logger.debug("update2() locus: " + curLocus);
logger.debug("comp = " + comp + " eval = " + eval);
Set<String> allSamples = new HashSet<String>();
Map<String, Genotype> compSampGenotypes = null;
if (isRelevantToPhasing(comp)) {
allSamples.addAll(comp.getSampleNames());
compSampGenotypes = comp.getGenotypes();
}
Map<String, Genotype> evalSampGenotypes = null;
if (isRelevantToPhasing(eval)) {
allSamples.addAll(eval.getSampleNames());
evalSampGenotypes = eval.getGenotypes();
}
for (String samp : allSamples) {
logger.debug("sample = " + samp);
Genotype compSampGt = null;
if (compSampGenotypes != null)
compSampGt = compSampGenotypes.get(samp);
Genotype evalSampGt = null;
if (evalSampGenotypes != null)
evalSampGt = evalSampGenotypes.get(samp);
if (compSampGt == null || evalSampGt == null) { // Since either comp or eval (or both) are missing the site, the best we can do is hope to preserve phase [if the non-missing one preserves phase]
// Having an unphased site breaks the phasing for the sample [does NOT permit "transitive phasing"] - hence, must reset phasing knowledge for both comp and eval [put a null CompEvalGenotypes]:
if (isNonNullButUnphased(compSampGt) || isNonNullButUnphased(evalSampGt))
samplePrevGenotypes.put(samp, null);
}
else { // Both comp and eval have a non-null Genotype at this site:
AllelePair compAllelePair = new AllelePair(compSampGt);
AllelePair evalAllelePair = new AllelePair(evalSampGt);
boolean breakPhasing = false;
if (compSampGt.isHet() != evalSampGt.isHet() || compSampGt.isHom() != evalSampGt.isHom())
breakPhasing = true; // since they are not both het or both hom
else { // both are het, or both are hom:
boolean topMatchesTopAndBottomMatchesBottom = (topMatchesTop(compAllelePair, evalAllelePair) && bottomMatchesBottom(compAllelePair, evalAllelePair));
boolean topMatchesBottomAndBottomMatchesTop = (topMatchesBottom(compAllelePair, evalAllelePair) && bottomMatchesTop(compAllelePair, evalAllelePair));
if (!topMatchesTopAndBottomMatchesBottom && !topMatchesBottomAndBottomMatchesTop)
breakPhasing = true; // since the 2 VCFs have different diploid genotypes for this sample
}
if (breakPhasing) {
samplePrevGenotypes.put(samp, null); // nothing to do for this site, AND must remove any history for the future
}
else if (compSampGt.isHet() && evalSampGt.isHet()) {
/* comp and eval have the HET same Genotype at this site:
[Note that if both are hom, then nothing is done here, but the het history IS preserved].
*/
CompEvalGenotypes prevCompAndEval = samplePrevGenotypes.get(samp);
if (prevCompAndEval != null && !prevCompAndEval.getLocus().onSameContig(curLocus)) // exclude curLocus if it is "phased" relative to a different chromosome
prevCompAndEval = null;
// Replace the previous hets with the current hets:
samplePrevGenotypes.put(samp, curLocus, compSampGt, evalSampGt);
if (prevCompAndEval != null) {
GenomeLoc prevLocus = prevCompAndEval.getLocus();
logger.debug("Potentially phaseable het locus: " + curLocus + " [relative to previous het locus: " + prevLocus + "]");
PhaseStats ps = samplePhasingStatistics.ensureSampleStats(samp);
boolean compSampIsPhased = genotypesArePhasedAboveThreshold(compSampGt);
boolean evalSampIsPhased = genotypesArePhasedAboveThreshold(evalSampGt);
if (compSampIsPhased || evalSampIsPhased) {
if (!evalSampIsPhased) {
ps.onlyCompPhased++;
interesting.addReason("ONLY_COMP", samp, group, prevLocus, "");
}
else if (!compSampIsPhased) {
ps.onlyEvalPhased++;
interesting.addReason("ONLY_EVAL", samp, group, prevLocus, "");
}
else { // both comp and eval are phased:
AllelePair prevCompAllelePair = new AllelePair(prevCompAndEval.getCompGenotpye());
AllelePair prevEvalAllelePair = new AllelePair(prevCompAndEval.getEvalGenotype());
// Sufficient to check only the top of comp, since we ensured that comp and eval have the same diploid genotypes for this sample:
boolean topsMatch = (topMatchesTop(prevCompAllelePair, prevEvalAllelePair) && topMatchesTop(compAllelePair, evalAllelePair));
boolean topMatchesBottom = (topMatchesBottom(prevCompAllelePair, prevEvalAllelePair) && topMatchesBottom(compAllelePair, evalAllelePair));
if (topsMatch || topMatchesBottom) {
ps.phasesAgree++;
Double compPQ = getPQ(compSampGt);
Double evalPQ = getPQ(evalSampGt);
if (compPQ != null && evalPQ != null && MathUtils.compareDoubles(compPQ, evalPQ) != 0)
interesting.addReason("PQ_CHANGE", samp, group, prevLocus, compPQ + " -> " + evalPQ);
}
else {
ps.phasesDisagree++;
logger.debug("SWITCHED locus: " + curLocus);
interesting.addReason("SWITCH", samp, group, prevLocus, toString(prevCompAllelePair, compAllelePair) + " -> " + toString(prevEvalAllelePair, evalAllelePair));
}
}
}
else {
ps.neitherPhased++;
}
}
}
}
}
logger.debug("\n" + samplePhasingStatistics + "\n");
return interesting.toString();
}
public static boolean isRelevantToPhasing(VariantContext vc) {
return (vc != null && !vc.isFiltered());
}
public boolean isNonNullButUnphased(Genotype gt) {
return (gt != null && !genotypesArePhasedAboveThreshold(gt));
}
public boolean genotypesArePhasedAboveThreshold(Genotype gt) {
if (gt.isHom()) // Can always consider a hom site to be phased to its predecessor, since its successor will only be phased to it if it's hom or "truly" phased
return true;
if (!gt.isPhased())
return false;
Double pq = getPQ(gt);
return (pq == null || pq >= getVEWalker().minPhaseQuality);
}
public static Double getPQ(Genotype gt) {
return gt.getAttributeAsDoubleNoException(ReadBackedPhasingWalker.PQ_KEY);
}
public boolean topMatchesTop(AllelePair b1, AllelePair b2) {
return b1.getTopAllele().equals(b2.getTopAllele());
}
public boolean topMatchesBottom(AllelePair b1, AllelePair b2) {
return b1.getTopAllele().equals(b2.getBottomAllele());
}
public boolean bottomMatchesTop(AllelePair b1, AllelePair b2) {
return topMatchesBottom(b2, b1);
}
public boolean bottomMatchesBottom(AllelePair b1, AllelePair b2) {
return b1.getBottomAllele().equals(b2.getBottomAllele());
}
public String toString(AllelePair prev, AllelePair cur) {
return prev.getTopAllele().getBaseString() + "+" + cur.getTopAllele().getBaseString() + "|" + prev.getBottomAllele().getBaseString() + "+" + cur.getBottomAllele().getBaseString();
}
public void finalizeEvaluation() {
}
private static class Reasons {
private StringBuilder sb;
public Reasons() {
sb = new StringBuilder();
}
public void addReason(String category, String sample, VariantEvalWalker.EvaluationContext evalGroup, GenomeLoc prevLoc, String reason) {
sb.append(category + "(" + sample + ", previous: " + prevLoc + " [" + evalGroup.compTrackName + ", " + evalGroup.evalTrackName + "]): " + reason + ";");
}
public String toString() {
if (sb.length() == 0)
return null;
return "reasons=" + sb.toString();
}
}
}
class CompEvalGenotypes {
private GenomeLoc loc;
private Genotype compGt;
private Genotype evalGt;
public CompEvalGenotypes(GenomeLoc loc, Genotype compGt, Genotype evalGt) {
this.loc = loc;
this.compGt = compGt;
this.evalGt = evalGt;
}
public GenomeLoc getLocus() {
return loc;
}
public Genotype getCompGenotpye() {
return compGt;
}
public Genotype getEvalGenotype() {
return evalGt;
}
public void setCompGenotype(Genotype compGt) {
this.compGt = compGt;
}
public void setEvalGenotype(Genotype evalGt) {
this.evalGt = evalGt;
}
}
class SamplePreviousGenotypes {
private HashMap<String, CompEvalGenotypes> sampleGenotypes = null;
public SamplePreviousGenotypes() {
this.sampleGenotypes = new HashMap<String, CompEvalGenotypes>();
}
public CompEvalGenotypes get(String sample) {
return sampleGenotypes.get(sample);
}
public void put(String sample, CompEvalGenotypes compEvalGts) {
sampleGenotypes.put(sample, compEvalGts);
}
public void put(String sample, GenomeLoc locus, Genotype compGt, Genotype evalGt) {
sampleGenotypes.put(sample, new CompEvalGenotypes(locus, compGt, evalGt));
}
}
class PhaseStats {
public int neitherPhased;
public int onlyCompPhased;
public int onlyEvalPhased;
public int phasesAgree;
public int phasesDisagree;
public PhaseStats() {
this.neitherPhased = 0;
this.onlyCompPhased = 0;
this.onlyEvalPhased = 0;
this.phasesAgree = 0;
this.phasesDisagree = 0;
}
public String toString() {
StringBuilder sb = new StringBuilder();
sb.append("Neither phased: " + neitherPhased + "\tOnly Comp: " + onlyCompPhased + "\tOnly Eval: " + onlyEvalPhased + "\tSame phase: " + phasesAgree + "\tOpposite phase: " + phasesDisagree);
return sb.toString();
}
public static String[] getFieldNamesArray() {
return new String[]{"total", "neither", "only_comp", "only_eval", "both", "match", "switch", "switch_rate"};
}
public Object getField(int index) {
switch (index) {
case (0):
return (neitherPhased + onlyCompPhased + onlyEvalPhased + phasesAgree + phasesDisagree);
case (1):
return neitherPhased;
case (2):
return onlyCompPhased;
case (3):
return onlyEvalPhased;
case (4):
return (phasesAgree + phasesDisagree);
case (5):
return phasesAgree;
case (6):
return phasesDisagree;
case (7):
return ((phasesDisagree == 0) ? 0 : ((double) phasesDisagree) / (phasesAgree + phasesDisagree));
default:
return -1;
}
}
}
/**
* a table of sample names to genotype phasing statistics
*/
class SamplePhasingStatistics implements TableType {
private HashMap<String, PhaseStats> sampleStats = null;
private double minPhaseQuality;
public SamplePhasingStatistics(double minPhaseQuality) {
this.sampleStats = new HashMap<String, PhaseStats>();
this.minPhaseQuality = minPhaseQuality;
}
public PhaseStats ensureSampleStats(String samp) {
PhaseStats ps = sampleStats.get(samp);
if (ps == null) {
ps = new PhaseStats();
sampleStats.put(samp, ps);
}
return ps;
}
/**
* @return one row per sample
*/
public String[] getRowKeys() {
return sampleStats.keySet().toArray(new String[sampleStats.size()]);
}
/**
* get the column keys
*
* @return a list of objects, in this case strings, that are the column names
*/
public String[] getColumnKeys() {
return PhaseStats.getFieldNamesArray();
}
public Object getCell(int x, int y) {
String[] rowKeys = getRowKeys();
PhaseStats ps = sampleStats.get(rowKeys[x]);
return ps.getField(y);
}
public String getName() {
return "Sample Phasing Statistics (for PQ >= " + minPhaseQuality + ")";
}
public String toString() {
StringBuilder sb = new StringBuilder();
for (Map.Entry<String, PhaseStats> sampPhaseStatsEnt : sampleStats.entrySet()) {
String sample = sampPhaseStatsEnt.getKey();
PhaseStats ps = sampPhaseStatsEnt.getValue();
sb.append(sample + "\t" + ps);
}
return sb.toString();
}
}

114
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/IndelLengthHistogram.java
View File

@ -1,114 +0,0 @@
package org.broadinstitute.sting.gatk.walkers.varianteval;
import org.broad.tribble.util.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.report.tags.Analysis;
import org.broadinstitute.sting.utils.report.tags.DataPoint;
import org.broadinstitute.sting.utils.report.utils.TableType;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
/**
* IF THERE IS NO JAVADOC RIGHT HERE, YELL AT chartl
*
* @Author chartl
* @Date May 26, 2010
*/
@Analysis(name = "Indel length histograms", description = "Shows the distrbution of insertion/deletion event lengths (negative for deletion, positive for insertion)")
public class IndelLengthHistogram extends VariantEvaluator {
private static final int SIZE_LIMIT = 50;
@DataPoint(name="indelLengthHistogram",description="Histogram of indel lengths")
IndelHistogram indelHistogram = new IndelHistogram(SIZE_LIMIT);
/*
* Indel length histogram table object
*/
static class IndelHistogram implements TableType {
private Integer[] colKeys;
private int limit;
private String[] rowKeys = {"EventLength"};
private Integer[] indelHistogram;
public IndelHistogram(int limit) {
colKeys = initColKeys(limit);
indelHistogram = initHistogram(limit);
this.limit = limit;
}
public Object[] getColumnKeys() {
return colKeys;
}
public Object[] getRowKeys() {
return rowKeys;
}
public Object getCell(int row, int col) {
return indelHistogram[col];
}
private Integer[] initColKeys(int size) {
Integer[] cK = new Integer[size*2+1];
int index = 0;
for ( int i = -size; i <= size; i ++ ) {
cK[index] = i;
index++;
}
return cK;
}
private Integer[] initHistogram(int size) {
Integer[] hist = new Integer[size*2+1];
for ( int i = 0; i < 2*size+1; i ++ ) {
hist[i] = 0;
}
return hist;
}
public String getName() { return "indelHistTable"; }
public void update(int eLength) {
indelHistogram[len2index(eLength)]++;
}
private int len2index(int len) {
if ( len > limit || len < -limit ) {
throw new ReviewedStingException("Indel length exceeds limit of "+limit+" please increase indel limit size");
}
return len + limit;
}
}
public IndelLengthHistogram(VariantEvalWalker parent) { super(parent); }
public boolean enabled() { return false; }
public String getName() { return "IndelLengthHistogram"; }
public int getComparisonOrder() { return 1; } // need only the evals
public String update1(VariantContext vc1, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
//System.out.println("Update1 called");
if ( ! vc1.isBiallelic() && vc1.isIndel() ) {
veWalker.getLogger().warn("[IndelLengthHistogram] Non-biallelic indel at "+ref.getLocus()+" ignored.");
return vc1.toString(); // biallelic sites are output
}
if ( vc1.isIndel() ) {
//System.out.println("Is indel");
if ( vc1.isInsertion() ) {
indelHistogram.update(vc1.getAlternateAllele(0).length());
} else if ( vc1.isDeletion() ) {
indelHistogram.update(-vc1.getReference().length());
} else {
throw new ReviewedStingException("Indel type that is not insertion or deletion.");
}
}
return null;
}
}

214
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/IndelMetricsByAC.java
View File

@ -1,214 +0,0 @@
package org.broadinstitute.sting.gatk.walkers.varianteval;
import org.broad.tribble.util.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.report.tags.Analysis;
import org.broadinstitute.sting.utils.report.tags.DataPoint;
import org.broadinstitute.sting.utils.report.utils.TableType;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import java.util.ArrayList;
/*
* Copyright (c) 2010 The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/
/**
* @author delangel
* @since Apr 11, 2010
*/
@Analysis(name = "Indel Metrics by allele count", description = "Shows various stats binned by allele count")
public class IndelMetricsByAC extends VariantEvaluator {
// a mapping from quality score histogram bin to Ti/Tv ratio
@DataPoint(name="Indel Metrics by AC", description = "Indel Metrics by allele count")
IndelMetricsByAc metrics = null;
//@DataPoint(name="Quality by Allele Count", description = "average variant quality for each allele count")
//AlleleCountStats alleleCountStats = null;
private static final int INDEL_SIZE_LIMIT = 100;
private static final int NUM_SCALAR_COLUMNS = 6;
static int len2Index(int ind) {
return ind+INDEL_SIZE_LIMIT+NUM_SCALAR_COLUMNS;
}
static int index2len(int ind) {
return ind-INDEL_SIZE_LIMIT-NUM_SCALAR_COLUMNS;
}
protected final static String[] METRIC_COLUMNS;
static {
METRIC_COLUMNS= new String[NUM_SCALAR_COLUMNS+2*INDEL_SIZE_LIMIT+1];
METRIC_COLUMNS[0] = "AC";
METRIC_COLUMNS[1] = "nIns";
METRIC_COLUMNS[2] = "nDels";
METRIC_COLUMNS[3] = "n";
METRIC_COLUMNS[4] = "nComplex";
METRIC_COLUMNS[5] = "nLong";
for (int k=NUM_SCALAR_COLUMNS; k < NUM_SCALAR_COLUMNS+ 2*INDEL_SIZE_LIMIT+1; k++)
METRIC_COLUMNS[k] = "indel_size_len"+Integer.valueOf(index2len(k));
}
class IndelMetricsAtAC {
public int ac = -1, nIns =0, nDel = 0, nComplex = 0, nLong;
public int sizeCount[] = new int[2*INDEL_SIZE_LIMIT+1];
public IndelMetricsAtAC(int ac) { this.ac = ac; }
public void update(VariantContext eval) {
int eventLength = 0;
if ( eval.isInsertion() ) {
eventLength = eval.getAlternateAllele(0).length();
nIns++;
} else if ( eval.isDeletion() ) {
eventLength = -eval.getReference().length();
nDel++;
}
else {
nComplex++;
}
if (Math.abs(eventLength) < INDEL_SIZE_LIMIT)
sizeCount[len2Index(eventLength)]++;
else
nLong++;
}
// corresponding to METRIC_COLUMNS
public String getColumn(int i) {
if (i >= NUM_SCALAR_COLUMNS && i <=NUM_SCALAR_COLUMNS+ 2*INDEL_SIZE_LIMIT)
return String.valueOf(sizeCount[i-NUM_SCALAR_COLUMNS]);
switch (i) {
case 0: return String.valueOf(ac);
case 1: return String.valueOf(nIns);
case 2: return String.valueOf(nDel);
case 3: return String.valueOf(nIns + nDel);
case 4: return String.valueOf(nComplex);
case 5: return String.valueOf(nLong);
default:
throw new ReviewedStingException("Unexpected column " + i);
}
}
}
class IndelMetricsByAc implements TableType {
ArrayList<IndelMetricsAtAC> metrics = new ArrayList<IndelMetricsAtAC>();
Object[] rows = null;
public IndelMetricsByAc( int nchromosomes ) {
rows = new Object[nchromosomes+1];
metrics = new ArrayList<IndelMetricsAtAC>(nchromosomes+1);
for ( int i = 0; i < nchromosomes + 1; i++ ) {
metrics.add(new IndelMetricsAtAC(i));
rows[i] = "ac" + i;
}
}
public Object[] getRowKeys() {
return rows;
}
public Object[] getColumnKeys() {
return METRIC_COLUMNS;
}
public String getName() {
return "IndelMetricsByAc";
}
//
public String getCell(int ac, int y) {
return metrics.get(ac).getColumn(y);
}
public String toString() {
String returnString = "";
return returnString;
}
public void incrValue( VariantContext eval ) {
int ac = -1;
if ( eval.hasGenotypes() )
ac = eval.getChromosomeCount(eval.getAlternateAllele(0));
else if ( eval.hasAttribute("AC") ) {
ac = Integer.valueOf(eval.getAttributeAsString("AC"));
}
if ( ac != -1 )
metrics.get(ac).update(eval);
}
}
public IndelMetricsByAC(VariantEvalWalker parent) {
super(parent);
// don't do anything
}
public String getName() {
return "IndelMetricsByAC";
}
public int getComparisonOrder() {
return 1; // we only need to see each eval track
}
public boolean enabled() {
return true;
}
public String toString() {
return getName();
}
public String update1(VariantContext eval, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
final String interesting = null;
if (eval != null ) {
if ( metrics == null ) {
int nSamples = this.getVEWalker().getNSamplesForEval(eval);
if ( nSamples != -1 )
metrics = new IndelMetricsByAc(2 * nSamples);
}
if ( eval.isIndel() && eval.isBiallelic() &&
metrics != null ) {
metrics.incrValue(eval);
}
}
return interesting; // This module doesn't capture any interesting sites, so return null
}
//public void finalizeEvaluation() {
//
//}
}

508
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/IndelStatistics.java
View File

@ -1,508 +0,0 @@
package org.broadinstitute.sting.gatk.walkers.varianteval;
import org.broad.tribble.util.variantcontext.Genotype;
import org.broad.tribble.util.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.report.tags.Analysis;
import org.broadinstitute.sting.utils.report.tags.DataPoint;
import org.broadinstitute.sting.utils.report.utils.TableType;
import java.util.Arrays;
import java.util.HashMap;
/*
* Copyright (c) 2010 The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/
@Analysis(name = "IndelStatistics", description = "Shows various indel metrics and statistics")
public class IndelStatistics extends VariantEvaluator {
@DataPoint(name="IndelStatistics", description = "Indel Statistics")
IndelStats indelStats = null;
@DataPoint(name="IndelClasses", description = "Indel Classification")
IndelClasses indelClasses = null;
private static final int INDEL_SIZE_LIMIT = 100;
private static final int NUM_SCALAR_COLUMNS = 10;
static int len2Index(int ind) {
return ind+INDEL_SIZE_LIMIT+NUM_SCALAR_COLUMNS;
}
static int index2len(int ind) {
return ind-INDEL_SIZE_LIMIT-NUM_SCALAR_COLUMNS;
}
static class IndelStats implements TableType {
protected final static String ALL_SAMPLES_KEY = "allSamples";
protected final static String[] COLUMN_KEYS;
static {
COLUMN_KEYS= new String[NUM_SCALAR_COLUMNS+2*INDEL_SIZE_LIMIT+1];
COLUMN_KEYS[0] = "heterozygosity";
COLUMN_KEYS[1] = "number_of_insertions";
COLUMN_KEYS[2] = "number_of_deletions";
COLUMN_KEYS[3] = "number_het_insertions";
COLUMN_KEYS[4] = "number_homozygous_insertions";
COLUMN_KEYS[5] = "number_het_deletions";
COLUMN_KEYS[6] = "number_homozygous_deletions";
COLUMN_KEYS[7] = "number of homozygous reference sites";
COLUMN_KEYS[8] = "number of complex events";
COLUMN_KEYS[9] = "number of long indels";
for (int k=NUM_SCALAR_COLUMNS; k < NUM_SCALAR_COLUMNS+ 2*INDEL_SIZE_LIMIT+1; k++)
COLUMN_KEYS[k] = "indel_size_len"+Integer.valueOf(index2len(k));
}
// map of sample to statistics
protected final HashMap<String, double[]> indelSummary = new HashMap<String, double[]>();
public IndelStats(final VariantContext vc) {
indelSummary.put(ALL_SAMPLES_KEY, new double[COLUMN_KEYS.length]);
for( final String sample : vc.getGenotypes().keySet() ) {
indelSummary.put(sample, new double[COLUMN_KEYS.length]);
}
}
/**
*
* @return one row per sample
*/
public Object[] getRowKeys() {
return indelSummary.keySet().toArray(new String[indelSummary.size()]);
}
public Object getCell(int x, int y) {
final Object[] rowKeys = getRowKeys();
return String.format("%4.2f",indelSummary.get(rowKeys[x])[y]);
}
/**
* get the column keys
* @return a list of objects, in this case strings, that are the column names
*/
public Object[] getColumnKeys() {
return COLUMN_KEYS;
}
public String getName() {
return "IndelStats";
}
public int getComparisonOrder() {
return 1; // we only need to see each eval track
}
public String toString() {
return getName();
}
/*
* increment the specified value
*/
public void incrValue(VariantContext vc) {
int eventLength = 0;
boolean isInsertion = false, isDeletion = false;
if ( vc.isInsertion() ) {
eventLength = vc.getAlternateAllele(0).length();
indelSummary.get(ALL_SAMPLES_KEY)[1]++;
isInsertion = true;
} else if ( vc.isDeletion() ) {
indelSummary.get(ALL_SAMPLES_KEY)[2]++;
eventLength = -vc.getReference().length();
isDeletion = true;
}
else {
indelSummary.get(ALL_SAMPLES_KEY)[8]++;
}
// make sure event doesn't overstep array boundaries
if (Math.abs(eventLength) < INDEL_SIZE_LIMIT)
indelSummary.get(ALL_SAMPLES_KEY)[len2Index(eventLength)]++;
else
indelSummary.get(ALL_SAMPLES_KEY)[9]++;
for( final String sample : vc.getGenotypes().keySet() ) {
if ( indelSummary.containsKey(sample) ) {
Genotype g = vc.getGenotype(sample);
boolean isVariant = (g.isCalled() && !g.isHomRef());
if (isVariant) {
// update ins/del count
if (isInsertion) {
indelSummary.get(sample)[1]++;
}
else if (isDeletion)
indelSummary.get(sample)[2]++;
else
indelSummary.get(sample)[8]++;
// update histogram
if (Math.abs(eventLength) < INDEL_SIZE_LIMIT)
indelSummary.get(sample)[len2Index(eventLength)]++;
else
indelSummary.get(sample)[9]++;
if (g.isHet())
if (isInsertion)
indelSummary.get(sample)[3]++;
else
indelSummary.get(sample)[5]++;
else
if (isInsertion)
indelSummary.get(sample)[4]++;
else
indelSummary.get(sample)[6]++;
}
else
indelSummary.get(sample)[7]++;
}
}
}
}
static class IndelClasses implements TableType {
protected final static String ALL_SAMPLES_KEY = "allSamples";
protected final static String[] COLUMN_KEYS;
static {
COLUMN_KEYS= new String[41];
COLUMN_KEYS[0] = "Novel_A";
COLUMN_KEYS[1] = "Novel_C";
COLUMN_KEYS[2] = "Novel_G";
COLUMN_KEYS[3] = "Novel_T";
COLUMN_KEYS[4] = "NOVEL_1";
COLUMN_KEYS[5] = "NOVEL_2";
COLUMN_KEYS[6] = "NOVEL_3";
COLUMN_KEYS[7] = "NOVEL_4";
COLUMN_KEYS[8] = "NOVEL_5";
COLUMN_KEYS[9] = "NOVEL_6";
COLUMN_KEYS[10] = "NOVEL_7";
COLUMN_KEYS[11] = "NOVEL_8";
COLUMN_KEYS[12] = "NOVEL_9";
COLUMN_KEYS[13] = "NOVEL_10orMore";
COLUMN_KEYS[14] = "RepeatExpansion_A";
COLUMN_KEYS[15] = "RepeatExpansion_C";
COLUMN_KEYS[16] = "RepeatExpansion_G";
COLUMN_KEYS[17] = "RepeatExpansion_T";
COLUMN_KEYS[18] = "RepeatExpansion_AC";
COLUMN_KEYS[19] = "RepeatExpansion_AG";
COLUMN_KEYS[20] = "RepeatExpansion_AT";
COLUMN_KEYS[21] = "RepeatExpansion_CA";
COLUMN_KEYS[22] = "RepeatExpansion_CG";
COLUMN_KEYS[23] = "RepeatExpansion_CT";
COLUMN_KEYS[24] = "RepeatExpansion_GA";
COLUMN_KEYS[25] = "RepeatExpansion_GC";
COLUMN_KEYS[26] = "RepeatExpansion_GT";
COLUMN_KEYS[27] = "RepeatExpansion_TA";
COLUMN_KEYS[28] = "RepeatExpansion_TC";
COLUMN_KEYS[29] = "RepeatExpansion_TG";
COLUMN_KEYS[30] = "RepeatExpansion_1";
COLUMN_KEYS[31] = "RepeatExpansion_2";
COLUMN_KEYS[32] = "RepeatExpansion_3";
COLUMN_KEYS[33] = "RepeatExpansion_4";
COLUMN_KEYS[34] = "RepeatExpansion_5";
COLUMN_KEYS[35] = "RepeatExpansion_6";
COLUMN_KEYS[36] = "RepeatExpansion_7";
COLUMN_KEYS[37] = "RepeatExpansion_8";
COLUMN_KEYS[38] = "RepeatExpansion_9";
COLUMN_KEYS[39] = "RepeatExpansion_10orMore";
COLUMN_KEYS[40] = "Other";
}
private static final int START_IND_NOVEL = 4;
private static final int STOP_IND_NOVEL = 13;
private static final int START_IND_FOR_REPEAT_EXPANSION_1 = 14;
private static final int STOP_IND_FOR_REPEAT_EXPANSION_2 = 29;
private static final int START_IND_FOR_REPEAT_EXPANSION_COUNTS = 30;
private static final int STOP_IND_FOR_REPEAT_EXPANSION_COUNTS = 39;
private static final int IND_FOR_OTHER_EVENT = 40;
private static final int START_IND_NOVEL_PER_BASE = 0;
private static final int STOP_IND_NOVEL_PER_BASE = 3;
// map of sample to statistics
protected final HashMap<String, int[]> indelClassSummary = new HashMap<String, int[]>();
public IndelClasses(final VariantContext vc) {
indelClassSummary.put(ALL_SAMPLES_KEY, new int[COLUMN_KEYS.length]);
for( final String sample : vc.getGenotypes().keySet() ) {
indelClassSummary.put(sample, new int[COLUMN_KEYS.length]);
}
}
/**
*
* @return one row per sample
*/
public Object[] getRowKeys() {
return indelClassSummary.keySet().toArray(new String[indelClassSummary.size()]);
}
public Object getCell(int x, int y) {
final Object[] rowKeys = getRowKeys();
return String.format("%d",indelClassSummary.get(rowKeys[x])[y]);
}
/**
* get the column keys
* @return a list of objects, in this case strings, that are the column names
*/
public Object[] getColumnKeys() {
return COLUMN_KEYS;
}
public String getName() {
return "IndelClasses";
}
public int getComparisonOrder() {
return 1; // we only need to see each eval track
}
public String toString() {
return getName();
}
private void incrementSampleStat(VariantContext vc, int index) {
indelClassSummary.get(ALL_SAMPLES_KEY)[index]++;
for( final String sample : vc.getGenotypes().keySet() ) {
if ( indelClassSummary.containsKey(sample) ) {
Genotype g = vc.getGenotype(sample);
boolean isVariant = (g.isCalled() && !g.isHomRef());
if (isVariant)
// update count
indelClassSummary.get(sample)[index]++;
}
}
}
/*
* increment the specified value
*/
private String findMinimalEvent(String eventString) {
// for each length up to given string length, see if event string is a repetition of units of size N
boolean foundSubstring = false;
String minEvent = eventString;
for (int k=1; k < eventString.length(); k++) {
if (eventString.length() % k > 0)
continue;
String str = eventString.substring(0,k);
// now see if event string is a repetition of str
int numReps = eventString.length() / k;
String r = "";
for (int j=0; j < numReps; j++)
r = r.concat(str);
if (r.matches(eventString)) {
foundSubstring = true;
minEvent = str;
break;
}
}
return minEvent;
}
public void incrValue(VariantContext vc, ReferenceContext ref) {
int eventLength = 0;
boolean isInsertion = false, isDeletion = false;
String indelAlleleString;
if ( vc.isInsertion() ) {
isInsertion = true;
indelAlleleString = vc.getAlternateAllele(0).getDisplayString();
} else if ( vc.isDeletion() ) {
isDeletion = true;
indelAlleleString = vc.getReference().getDisplayString();
}
else {
incrementSampleStat(vc, IND_FOR_OTHER_EVENT);
return;
}
byte[] refBases = ref.getBases();
indelAlleleString = findMinimalEvent(indelAlleleString);
eventLength = indelAlleleString.length();
// See first if indel is a repetition of bases before current
int indStart = refBases.length/2-eventLength+1;
boolean done = false;
int numRepetitions = 0;
while (!done) {
if (indStart < 0)
done = true;
else {
String refPiece = new String(Arrays.copyOfRange(refBases,indStart,indStart+eventLength));
if (refPiece.matches(indelAlleleString))
{
numRepetitions++;
indStart = indStart - eventLength;
}
else
done = true;
}
}
// now do it forward
done = false;
indStart = refBases.length/2+1;
while (!done) {
if (indStart + eventLength >= refBases.length)
break;
else {
String refPiece = new String(Arrays.copyOfRange(refBases,indStart,indStart+eventLength));
if (refPiece.matches(indelAlleleString))
{
numRepetitions++;
indStart = indStart + eventLength;
}
else
done = true;
}
}
if (numRepetitions == 0) {
//unrepeated sequence from surroundings
int ind = START_IND_NOVEL + (eventLength-1);
if (ind > STOP_IND_NOVEL)
ind = STOP_IND_NOVEL;
incrementSampleStat(vc, ind);
if (eventLength == 1) {
// log single base indels additionally by base
String keyStr = "Novel_" + indelAlleleString;
int k;
for (k=START_IND_NOVEL_PER_BASE; k <= STOP_IND_NOVEL_PER_BASE; k++) {
if (keyStr.matches(COLUMN_KEYS[k]))
break;
}
// log now event
incrementSampleStat(vc, k);
}
}
else {
int ind = START_IND_FOR_REPEAT_EXPANSION_COUNTS + (numRepetitions-1);
if (ind > STOP_IND_FOR_REPEAT_EXPANSION_COUNTS)
ind = STOP_IND_FOR_REPEAT_EXPANSION_COUNTS;
incrementSampleStat(vc, ind);
if (eventLength<=2) {
// for single or dinucleotide indels, we further log the base in which they occurred
String keyStr = "RepeatExpansion_" + indelAlleleString;
int k;
for (k=START_IND_FOR_REPEAT_EXPANSION_1; k <= STOP_IND_FOR_REPEAT_EXPANSION_2; k++) {
if (keyStr.matches(COLUMN_KEYS[k]))
break;
}
// log now event
incrementSampleStat(vc, k);
}
}
//g+
/*
System.out.format("RefBefore: %s\n",new String(refBefore));
System.out.format("RefAfter: %s\n",new String(refAfter));
System.out.format("Indel Allele: %s\n", indelAlleleString);
System.out.format("Num Repetitions: %d\n", numRepetitions);
*/
}
}
public IndelStatistics(VariantEvalWalker parent) {
super(parent);
// don't do anything
}
public String getName() {
return "IndelStatistics";
}
public int getComparisonOrder() {
return 1; // we only need to see each eval track
}
public boolean enabled() {
return true;
}
public String toString() {
return getName();
}
public String update2(VariantContext eval, VariantContext validation, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context, VariantEvalWalker.EvaluationContext group) {
return null;
}
public String update1(VariantContext eval, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
if (eval != null ) {
if ( indelStats == null ) {
int nSamples = this.getVEWalker().getNSamplesForEval(eval);
if ( nSamples != -1 )
indelStats = new IndelStats(eval);
}
if ( indelClasses == null ) {
indelClasses = new IndelClasses(eval);
}
if ( eval.isIndel() && eval.isBiallelic() ) {
if (indelStats != null )
indelStats.incrValue(eval);
if (indelClasses != null)
indelClasses.incrValue(eval, ref);
}
}
return null; // This module doesn't capture any interesting sites, so return null
}
public String update0(VariantContext eval, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
return null;
}
public void finalizeEvaluation() {
//
int k=0;
}
}

182
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/MendelianViolationEvaluator.java
View File

@ -1,182 +0,0 @@
package org.broadinstitute.sting.gatk.walkers.varianteval;
import org.broad.tribble.util.variantcontext.Allele;
import org.broad.tribble.util.variantcontext.Genotype;
import org.broad.tribble.util.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.utils.report.tags.Analysis;
import org.broadinstitute.sting.utils.report.tags.DataPoint;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import java.util.List;
import java.util.Arrays;
import java.util.regex.Pattern;
import java.util.regex.Matcher;
/**
* Mendelian violation detection and counting
* <p/>
* a violation looks like:
* Suppose dad = A/B and mom = C/D
* The child can be [A or B] / [C or D].
* If the child doesn't match this, the site is a violation
* <p/>
* Some examples:
* <p/>
* mom = A/A, dad = C/C
* child can be A/C only
* <p/>
* mom = A/C, dad = C/C
* child can be A/C or C/C
* <p/>
* mom = A/C, dad = A/C
* child can be A/A, A/C, C/C
* <p/>
* The easiest way to do this calculation is to:
* <p/>
* Get alleles for mom => A/B
* Get alleles for dad => C/D
* Make allowed genotypes for child: A/C, A/D, B/C, B/D
* Check that the child is one of these.
*/
@Analysis(name = "Mendelian Violation Evaluator", description = "Mendelian Violation Evaluator")
public class MendelianViolationEvaluator extends VariantEvaluator {
@DataPoint(name = "variants", description = "Number of mendelian variants found")
long nVariants;
@DataPoint(name = "violations", description = "Number of mendelian violations found")
long nViolations;
@DataPoint(name="KHR->PHV",description = "number of child hom ref calls where the parent was hom variant")
long KidHomRef_ParentHomVar;
@DataPoint(name="KHET->PHR",description = "number of child het calls where the parent was hom ref")
long KidHet_ParentsHomRef;
@DataPoint(name="KHET->PHV",description = "number of child het calls where the parent was hom variant")
long KidHet_ParentsHomVar;
@DataPoint(name="KHV->PHR",description = "number of child hom variant calls where the parent was hom ref")
long KidHomVar_ParentHomRef;
TrioStructure trio;
private static Pattern FAMILY_PATTERN = Pattern.compile("(.*)\\+(.*)=(.*)");
public static class TrioStructure {
public String mom, dad, child;
}
public static TrioStructure parseTrioDescription(String family) {
Matcher m = FAMILY_PATTERN.matcher(family);
if (m.matches()) {
TrioStructure trio = new TrioStructure();
//System.out.printf("Found a family pattern: %s%n", parent.FAMILY_STRUCTURE);
trio.mom = m.group(1);
trio.dad = m.group(2);
trio.child = m.group(3);
return trio;
} else {
throw new IllegalArgumentException("Malformatted family structure string: " + family + " required format is mom+dad=child");
}
}
public MendelianViolationEvaluator(VariantEvalWalker parent) {
super(parent);
if (enabled()) {
trio = parseTrioDescription(parent.FAMILY_STRUCTURE);
parent.getLogger().debug(String.format("Found a family pattern: %s mom=%s dad=%s child=%s",
parent.FAMILY_STRUCTURE, trio.mom, trio.dad, trio.child));
}
}
public boolean enabled() {
return getVEWalker().FAMILY_STRUCTURE != null;
}
private double getQThreshold() {
return getVEWalker().MENDELIAN_VIOLATION_QUAL_THRESHOLD / 10; // we aren't 10x scaled in the GATK a la phred
}
public String getName() {
return "mendelian_violations";
}
public int getComparisonOrder() {
return 1; // we only need to see each eval track
}
public String update1(VariantContext vc, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
if (vc.isBiallelic() && vc.hasGenotypes()) { // todo -- currently limited to biallelic loci
Genotype momG = vc.getGenotype(trio.mom);
Genotype dadG = vc.getGenotype(trio.dad);
Genotype childG = vc.getGenotype(trio.child);
if (includeGenotype(momG) && includeGenotype(dadG) && includeGenotype(childG)) {
nVariants++;
if (momG == null || dadG == null || childG == null)
throw new IllegalArgumentException(String.format("VariantContext didn't contain genotypes for expected trio members: mom=%s dad=%s child=%s", trio.mom, trio.dad, trio.child));
// all genotypes are good, so let's see if child is a violation
if (isViolation(vc, momG, dadG, childG)) {
nViolations++;
String label;
if (childG.isHomRef() && (momG.isHomVar() || dadG.isHomVar())) {
label = "KidHomRef_ParentHomVar";
KidHomRef_ParentHomVar++;
} else if (childG.isHet() && (momG.isHomRef() && dadG.isHomRef())) {
label = "KidHet_ParentsHomRef";
KidHet_ParentsHomRef++;
} else if (childG.isHet() && (momG.isHomVar() && dadG.isHomVar())) {
label = "KidHet_ParentsHomVar";
KidHet_ParentsHomVar++;
} else if (childG.isHomVar() && (momG.isHomRef() || dadG.isHomRef())) {
label = "KidHomVar_ParentHomRef";
KidHomVar_ParentHomRef++;
} else {
throw new ReviewedStingException("BUG: unexpected child genotype class " + childG);
}
return "MendelViolation=" + label;
}
}
}
return null; // we don't capture any intersting sites
}
private boolean includeGenotype(Genotype g) {
return g.getNegLog10PError() > getQThreshold() && g.isCalled();
}
public static boolean isViolation(VariantContext vc, Genotype momG, Genotype dadG, Genotype childG) {
return isViolation(vc, momG.getAlleles(), dadG.getAlleles(), childG.getAlleles());
}
public static boolean isViolation(VariantContext vc, TrioStructure trio ) {
return isViolation(vc, vc.getGenotype(trio.mom), vc.getGenotype(trio.dad), vc.getGenotype(trio.child) );
}
public static boolean isViolation(VariantContext vc, List<Allele> momA, List<Allele> dadA, List<Allele> childA) {
//VariantContext momVC = vc.subContextFromGenotypes(momG);
//VariantContext dadVC = vc.subContextFromGenotypes(dadG);
int i = 0;
Genotype childG = new Genotype("kidG", childA);
for (Allele momAllele : momA) {
for (Allele dadAllele : dadA) {
if (momAllele.isCalled() && dadAllele.isCalled()) {
Genotype possibleChild = new Genotype("possibleGenotype" + i, Arrays.asList(momAllele, dadAllele));
if (childG.sameGenotype(possibleChild, false)) {
return false;
}
}
}
}
return true;
}
}

67
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/PrintMissingComp.java
View File

@ -1,67 +0,0 @@
/*
* Copyright (c) 2010, The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.gatk.walkers.varianteval;
import org.broad.tribble.util.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.refdata.*;
import org.broadinstitute.sting.utils.report.tags.Analysis;
import org.broadinstitute.sting.utils.report.tags.DataPoint;
@Analysis(name = "PrintMissingComp", description = "the overlap between eval and comp sites")
public class PrintMissingComp extends VariantEvaluator {
@DataPoint(name = "evals not at comp", description = "number of eval sites outside of comp sites")
long nMissing = 0;
public PrintMissingComp(VariantEvalWalker parent) {
super(parent);
}
public String getName() {
return "PrintMissingComp";
}
public int getComparisonOrder() {
return 2; // we need to see each eval track and each comp track
}
public boolean enabled() {
return true;
}
public String update2(VariantContext eval, VariantContext comp, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
boolean compIsGood = comp != null && comp.isNotFiltered() && comp.isSNP();
boolean evalIsGood = eval != null && eval.isSNP();
if ( compIsGood & ! evalIsGood ) {
nMissing++;
return "MissingFrom" + comp.getSource();
} else {
return null;
}
}
}

175
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/SimpleMetricsByAC.java
View File

@ -1,175 +0,0 @@
package org.broadinstitute.sting.gatk.walkers.varianteval;
import org.broad.tribble.util.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContextUtils;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.report.tags.Analysis;
import org.broadinstitute.sting.utils.report.tags.DataPoint;
import org.broadinstitute.sting.utils.report.utils.TableType;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import java.util.ArrayList;
/*
* Copyright (c) 2010 The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/
/**
* @author depristo
* @since Apr 11, 2010
*/
@Analysis(name = "Quality Metrics by allele count", description = "Shows various stats binned by allele count")
public class SimpleMetricsByAC extends VariantEvaluator implements StandardEval {
// a mapping from quality score histogram bin to Ti/Tv ratio
@DataPoint(name="TiTv by AC", description = "TiTv by allele count")
MetricsByAc metrics = null;
//@DataPoint(name="Quality by Allele Count", description = "average variant quality for each allele count")
//AlleleCountStats alleleCountStats = null;
private final static Object[] METRIC_COLUMNS = {"AC", "nTi", "nTv", "n", "Ti/Tv"};
class MetricsAtAC {
public int ac = -1, nTi = 0, nTv = 0;
public MetricsAtAC(int ac) { this.ac = ac; }
public void update(VariantContext eval) {
if ( VariantContextUtils.isTransition(eval) )
nTi++;
else
nTv++;
}
// corresponding to METRIC_COLUMNS
public String getColumn(int i) {
switch (i) {
case 0: return String.valueOf(ac);
case 1: return String.valueOf(nTi);
case 2: return String.valueOf(nTv);
case 3: return String.valueOf(nTi + nTv);
case 4: return String.valueOf(ratio(nTi, nTv));
default:
throw new ReviewedStingException("Unexpected column " + i);
}
}
}
class MetricsByAc implements TableType {
ArrayList<MetricsAtAC> metrics = new ArrayList<MetricsAtAC>();
Object[] rows = null;
public MetricsByAc( int nchromosomes ) {
rows = new Object[nchromosomes+1];
metrics = new ArrayList<MetricsAtAC>(nchromosomes+1);
for ( int i = 0; i < nchromosomes + 1; i++ ) {
metrics.add(new MetricsAtAC(i));
rows[i] = "ac" + i;
}
}
public Object[] getRowKeys() {
return rows;
}
public Object[] getColumnKeys() {
return METRIC_COLUMNS;
}
public String getName() {
return "MetricsByAc";
}
//
public String getCell(int ac, int y) {
return metrics.get(ac).getColumn(y);
}
public String toString() {
String returnString = "";
return returnString;
}
public void incrValue( VariantContext eval ) {
int ac = -1;
if ( eval.hasGenotypes() )
ac = eval.getChromosomeCount(eval.getAlternateAllele(0));
else if ( eval.hasAttribute("AC") ) {
ac = Integer.valueOf(eval.getAttributeAsString("AC"));
}
if ( ac != -1 )
metrics.get(ac).update(eval);
}
}
public SimpleMetricsByAC(VariantEvalWalker parent) {
super(parent);
// don't do anything
}
public String getName() {
return "SimpleMetricsByAC";
}
public int getComparisonOrder() {
return 1; // we only need to see each eval track
}
public boolean enabled() {
return true;
}
public String toString() {
return getName();
}
public String update1(VariantContext eval, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
final String interesting = null;
if (eval != null ) {
if ( metrics == null ) {
int nSamples = this.getVEWalker().getNSamplesForEval(eval);
if ( nSamples != -1 )
metrics = new MetricsByAc(2 * nSamples);
}
if ( eval.isSNP() &&
eval.isBiallelic() &&
metrics != null ) {
metrics.incrValue(eval);
}
}
return interesting; // This module doesn't capture any interesting sites, so return null
}
//public void finalizeEvaluation() {
//
//}
}

28
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/StandardEval.java
View File

@ -1,28 +0,0 @@
/*
* Copyright (c) 2010.
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
* THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.gatk.walkers.varianteval;
public interface StandardEval {}

136
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/ThetaVariantEvaluator.java
View File

@ -1,136 +0,0 @@
package org.broadinstitute.sting.gatk.walkers.varianteval;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broad.tribble.util.variantcontext.VariantContext;
import org.broad.tribble.util.variantcontext.Genotype;
import org.broad.tribble.util.variantcontext.Allele;
import org.broadinstitute.sting.utils.report.tags.Analysis;
import org.broadinstitute.sting.utils.report.tags.DataPoint;
import java.util.concurrent.ConcurrentMap;
import java.util.concurrent.ConcurrentHashMap;
@Analysis(name = "Theta Variant Evaluator", description = "Computes different estimates of theta based on variant sites and genotypes")
public class ThetaVariantEvaluator extends VariantEvaluator {
@DataPoint(name = "avg_heterozygosity", description = "Average heterozygosity at variant sites; note that missing genotypes are ignored when computing this value")
double avgHet = 0.0;
@DataPoint(name = "avg_pairwise_diffs", description = "Average pairwise differences at aligned sequences; averaged over both number of sequeneces and number of variant sites; note that missing genotypes are ignored when computing this value")
double avgAvgDiffs = 0.0;
@DataPoint(name = "sum_heterozygosity", description = "Sum of heterozygosity over all variant sites; divide this by total target to get estimate of per base theta")
double totalHet = 0.0;
@DataPoint(name = "sum_pairwise_diffs", description = "Sum of pairwise diffs over all variant sites; divide this by total target to get estimate of per base theta")
double totalAvgDiffs = 0.0;
@DataPoint(name = "theta_region_num_sites", description = "Theta for entire region estimated based on number of segregating sites; divide ths by total target to get estimate of per base theta")
double thetaRegionNumSites = 0.0;
//helper variables
double numSites = 0;
public ThetaVariantEvaluator(VariantEvalWalker parent) {
super(parent);
}
public boolean enabled() {
return true;
}
public String getName() {
return "theta";
}
public int getComparisonOrder() {
return 1;
}
public String update1(VariantContext vc, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
if (vc == null || !vc.isSNP() || !vc.hasGenotypes()) {
return null; //no interesting sites
}
if (vc.hasGenotypes()) {
//this maps allele to a count
ConcurrentMap<String, Integer> alleleCounts = new ConcurrentHashMap<String, Integer>();
int numHetsHere = 0;
float numGenosHere = 0;
int numIndsHere = 0;
for (Genotype genotype : vc.getGenotypes().values()) {
numIndsHere++;
if (!genotype.isNoCall()) {
//increment stats for heterozygosity
if (genotype.isHet()) {
numHetsHere++;
}
numGenosHere++;
//increment stats for pairwise mismatches
for (Allele allele : genotype.getAlleles()) {
if (allele.isNonNull() && allele.isCalled()) {
String alleleString = allele.toString();
alleleCounts.putIfAbsent(alleleString, 0);
alleleCounts.put(alleleString, alleleCounts.get(alleleString) + 1);
}
}
}
}
if (numGenosHere > 0) {
//only if have one called genotype at least
this.numSites++;
this.totalHet += numHetsHere / numGenosHere;
//compute based on num sites
float harmonicFactor = 0;
for (int i = 1; i <= numIndsHere; i++) {
harmonicFactor += 1.0 / i;
}
this.thetaRegionNumSites += 1.0 / harmonicFactor;
//now compute pairwise mismatches
float numPairwise = 0;
float numDiffs = 0;
for (String allele1 : alleleCounts.keySet()) {
int allele1Count = alleleCounts.get(allele1);
for (String allele2 : alleleCounts.keySet()) {
if (allele1.compareTo(allele2) < 0) {
continue;
}
if (allele1 .compareTo(allele2) == 0) {
numPairwise += allele1Count * (allele1Count - 1) * .5;
}
else {
int allele2Count = alleleCounts.get(allele2);
numPairwise += allele1Count * allele2Count;
numDiffs += allele1Count * allele2Count;
}
}
}
if (numPairwise > 0) {
this.totalAvgDiffs += numDiffs / numPairwise;
}
}
}
return null;
}
@Override
public void finalizeEvaluation() {
if (this.numSites > 0) {
this.avgHet = this.totalHet / this.numSites;
this.avgAvgDiffs = this.totalAvgDiffs / this.numSites;
}
}
}

68
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/TiTvVariantEvaluator.java
View File

@ -1,68 +0,0 @@
package org.broadinstitute.sting.gatk.walkers.varianteval;
import org.broad.tribble.util.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContextUtils;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.utils.report.tags.Analysis;
import org.broadinstitute.sting.utils.report.tags.DataPoint;
@Analysis(name = "Ti/Tv Variant Evaluator", description = "Ti/Tv Variant Evaluator")
public class TiTvVariantEvaluator extends VariantEvaluator implements StandardEval {
@DataPoint(name = "ti_count", description = "number of transition loci")
long nTi = 0;
@DataPoint(name = "tv_count", description = "number of transversion loci")
long nTv = 0;
@DataPoint(name = "ti/tv_ratio", description = "the transition to transversion ratio")
double tiTvRatio = 0.0;
@DataPoint(name = "ti_count_comp", description = "number of comp transition sites")
long nTiInComp = 0;
@DataPoint(name = "tv_count_comp", description = "number of comp transversion sites")
long nTvInComp = 0;
@DataPoint(name = "ti/tv_ratio_standard", description = "the transition to transversion ratio for comp sites")
double TiTvRatioStandard = 0.0;
public TiTvVariantEvaluator(VariantEvalWalker parent) {
super(parent);
}
public boolean enabled() {
return true;
}
public String getName() {
return "titv";
}
public int getComparisonOrder() {
return 2; // we only need to see each eval track
}
public void updateTiTv(VariantContext vc, boolean updateStandard) {
if (vc != null && vc.isSNP() && vc.isBiallelic()) {
if (VariantContextUtils.isTransition(vc)) {
if (updateStandard) nTiInComp++;
else nTi++;
} else {
if (updateStandard) nTvInComp++;
else nTv++;
}
}
}
public String update2(VariantContext vc1, VariantContext vc2, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
if (vc1 != null) updateTiTv(vc1, false);
if (vc2 != null) updateTiTv(vc2, true);
return null; // we don't capture any intersting sites
}
@Override
public void finalizeEvaluation() {
// the ti/tv ratio needs to be set (it's not calculated per-variant).
this.tiTvRatio = rate(nTi,nTv);
this.TiTvRatioStandard = rate(nTiInComp, nTvInComp);
}
}

918
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/VariantEvalWalker.java
View File

@ -1,918 +0,0 @@
/*
* Copyright (c) 2010 The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
* THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.gatk.walkers.varianteval;
import org.apache.commons.jexl2.*;
import org.apache.commons.jexl2.introspection.*;
import org.apache.commons.logging.LogFactory;
import org.apache.log4j.Logger;
import org.broad.tribble.util.variantcontext.MutableVariantContext;
import org.broad.tribble.util.variantcontext.VariantContext;
import org.broad.tribble.vcf.VCFConstants;
import org.broad.tribble.vcf.VCFWriter;
import org.broad.tribble.vcf.VCFHeader;
import org.broad.tribble.vcf.VCFHeaderLine;
import org.broadinstitute.sting.commandline.Hidden;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContextUtils;
import org.broadinstitute.sting.gatk.datasources.simpleDataSources.ReferenceOrderedDataSource;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.refdata.utils.helpers.DbSNPHelper;
import org.broadinstitute.sting.gatk.walkers.Reference;
import org.broadinstitute.sting.gatk.walkers.RodWalker;
import org.broadinstitute.sting.gatk.walkers.Window;
import org.broadinstitute.sting.gatk.walkers.TreeReducible;
import org.broadinstitute.sting.gatk.walkers.variantrecalibration.Tranche;
import org.broadinstitute.sting.gatk.walkers.variantrecalibration.VariantRecalibrator;
import org.broadinstitute.sting.utils.SampleUtils;
import org.broadinstitute.sting.utils.report.ReportMarshaller;
import org.broadinstitute.sting.utils.report.VE2ReportFactory;
import org.broadinstitute.sting.utils.report.templates.ReportFormat;
import org.broadinstitute.sting.utils.report.utils.Node;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import org.broadinstitute.sting.utils.classloader.PluginManager;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.commandline.Argument;
import org.broadinstitute.sting.commandline.Output;
import org.broadinstitute.sting.utils.exceptions.DynamicClassResolutionException;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.text.XReadLines;
import java.io.File;
import java.io.FileNotFoundException;
import java.io.OutputStreamWriter;
import java.io.PrintStream;
import java.lang.reflect.Constructor;
import java.lang.reflect.Field;
import java.util.*;
// todo -- evalations should support comment lines
// todo -- add Mendelian variable explanations (nDeNovo and nMissingTransmissions)
// todo -- site frequency spectrum eval (freq. of variants in eval as a function of their AC and AN numbers)
// todo -- clustered SNP counter
// todo -- HWEs
// todo -- indel metrics [count of sizes in/del should be in CountVariants]
// todo -- port over SNP density walker:
// todo -- see walker for WG calc but will need to make it work with intervals correctly
// Todo -- should really include argument parsing @annotations from subclass in this walker. Very
// todo -- useful general capability. Right now you need to add arguments to VariantEval2 to handle new
// todo -- evaluation arguments (which is better than passing a string!)
// todo -- these really should be implemented as default select expression
// todo Extend VariantEval, our general-purpose tool for SNP evaluation, to differentiate Ti/Tv at CpG islands and also
// todo classify (and count) variants into coding, non-coding, synonomous/non-symonomous, 2/4 fold degenerate sites, etc.
// todo Assume that the incoming VCF has the annotations (you don't need to do this) but VE2 should split up results by
// todo these catogies automatically (using the default selects)
// todo -- this is really more a documentation issue. Really would be nice to have a pre-defined argument packet that
// todo -- can be provided to the system
// todo -- We agreed to report two standard values for variant evaluation from here out. One, we will continue to report
// todo -- the dbSNP 129 rate. Additionally, we will start to report the % of variants found that have already been seen in
// todo -- 1000 Genomes. This should be implemented as another standard comp_1kg binding, pointing to only variants
// todo -- discovered and released by 1KG. Might need to make this data set ourselves and keep it in GATK/data like
// todo -- dbsnp rod
//
// todo -- implement as select statment, but it's hard for multi-sample calls.
// todo -- Provide separate dbsnp rates for het only calls and any call where there is at least one hom-var genotype,
// todo -- since hets are much more likely to be errors
//
// todo -- Add Heng's hom run metrics -- single sample haplotype block lengths
/**
* General-purpose tool for variant evaluation (% in dbSNP, genotype concordance, Ts/Tv ratios, and a lot more)
*/
@Reference(window=@Window(start=-50,stop=50))
public class VariantEvalWalker extends RodWalker<Integer, Integer> implements TreeReducible<Integer> {
@Output
protected PrintStream out;
// --------------------------------------------------------------------------------------------------------------
//
// walker arguments
//
// --------------------------------------------------------------------------------------------------------------
@Argument(shortName="select", doc="One or more stratifications to use when evaluating the data", required=false)
protected ArrayList<String> SELECT_EXPS = new ArrayList<String>();
@Argument(shortName="selectName", doc="Names to use for the list of stratifications (must be a 1-to-1 mapping)", required=false)
protected ArrayList<String> SELECT_NAMES = new ArrayList<String>();
@Hidden
@Argument(shortName="summary", doc="One or more JEXL staments to log after evaluating the data", required=false)
protected ArrayList<String> SUMMARY_EXPS = new ArrayList<String>();
@Hidden
@Argument(shortName="validate", doc="One or more JEXL validations to use after evaluating the data", required=false)
protected ArrayList<String> VALIDATE_EXPS = new ArrayList<String>();
@Argument(shortName="known", doc="Name of ROD bindings containing variant sites that should be treated as known when splitting eval rods into known and novel subsets", required=false)
protected String[] KNOWN_NAMES = {DbSNPHelper.STANDARD_DBSNP_TRACK_NAME};
@Argument(shortName="sample", doc="Derive eval and comp contexts using only these sample genotypes, when genotypes are available in the original context", required=false)
protected String[] SAMPLES = {};
private List<String> SAMPLES_LIST = null;
//
// Arguments for choosing which modules to run
//
@Argument(fullName="evalModule", shortName="E", doc="One or more specific eval modules to apply to the eval track(s) (in addition to the standard modules, unless -noStandard is specified)", required=false)
protected String[] modulesToUse = {};
@Argument(fullName="doNotUseAllStandardModules", shortName="noStandard", doc="Do not use the standard modules by default (instead, only those that are specified with the -E option)")
protected Boolean NO_STANDARD = false;
@Argument(fullName="list", shortName="ls", doc="List the available eval modules and exit")
protected Boolean LIST = false;
//
// Arguments for Mendelian Violation calculations
//
@Argument(shortName="family", doc="If provided, genotypes in will be examined for mendelian violations: this argument is a string formatted as dad+mom=child where these parameters determine which sample names are examined", required=false)
protected String FAMILY_STRUCTURE;
@Argument(shortName="MVQ", fullName="MendelianViolationQualThreshold", doc="Minimum genotype QUAL score for each trio member required to accept a site as a violation", required=false)
protected double MENDELIAN_VIOLATION_QUAL_THRESHOLD = 50;
@Output(shortName="outputVCF", fullName="InterestingSitesVCF", doc="If provided, interesting sites emitted to this vcf and the INFO field annotated as to why they are interesting", required=false)
protected VCFWriter writer = null;
@Argument(shortName="gcLog", fullName="GenotypeCocordanceLog", doc="If provided, sites with genotype concordance problems (e.g., FP and FNs) will be emitted ot this file", required=false)
protected PrintStream gcLog = null;
private static double NO_MIN_QUAL_SCORE = -1.0;
@Argument(shortName = "Q", fullName="minPhredConfidenceScore", doc="Minimum confidence score to consider an evaluation SNP a variant", required=false)
public double minQualScore = NO_MIN_QUAL_SCORE;
@Argument(shortName = "Qcomp", fullName="minPhredConfidenceScoreForComp", doc="Minimum confidence score to consider a comp SNP a variant", required=false)
public double minCompQualScore = NO_MIN_QUAL_SCORE;
@Argument(shortName = "dels", fullName="indelCalls", doc="evaluate indels rather than SNPs", required = false)
public boolean dels = false;
// Right now we will only be looking at SNPS
// todo -- enable INDEL variant contexts, there's no reason not to but the integration tests need to be updated
EnumSet<VariantContext.Type> ALLOW_VARIANT_CONTEXT_TYPES = EnumSet.of(VariantContext.Type.SNP, VariantContext.Type.NO_VARIATION);
@Argument(shortName="rsID", fullName="rsID", doc="If provided, list of rsID and build number for capping known snps by their build date", required=false)
protected String rsIDFile = null;
@Argument(shortName="maxRsIDBuild", fullName="maxRsIDBuild", doc="If provided, only variants with rsIDs <= maxRsIDBuild will be included in the set of known snps", required=false)
protected int maxRsIDBuild = Integer.MAX_VALUE;
@Argument(shortName="reportType", fullName="reportType", doc="If provided, set the template type", required=false)
protected VE2ReportFactory.VE2TemplateType reportType = VE2ReportFactory.defaultReportFormat;
@Output(shortName="reportLocation", fullName="reportLocation", doc="If provided, set the base file for reports (Required for output formats with more than one file per analysis)", required=false)
protected File outputLocation = null;
@Argument(shortName="nSamples", fullName="nSamples", doc="If provided, analyses that need the number of samples in an eval track that has no genotype information will receive this number as the number of samples", required=false)
protected int nSamples = -1;
Set<String> rsIDsToExclude = null;
@Argument(shortName="aatk", fullName="aminoAcidTransitionKey", doc="required for the amino acid transition table; this is the key in the info field for the VCF for the transition", required = false)
public String aminoAcidTransitionKey = null;
@Argument(shortName="aats", fullName="aminoAcidTransitionSplit", doc="required for the amino acid transition table, this is the key on which to split the info field value to get the reference and alternate amino acids", required=false)
public String aminoAcidTransitionSplit = null;
@Argument(shortName="aatUseCodons", fullName="aminoAcidsRepresentedByCodons", doc="for the amino acid table, specifiy that the transitions are represented as codon changes, and not directly amino acid names", required = false)
public boolean aatUseCodons = false;
@Argument(shortName="disI", fullName="discordantInteresting", doc="If passed, write discordant sites as interesting", required=false)
protected boolean DISCORDANT_INTERESTING = false;
@Argument(fullName="tranchesFile", shortName="tf", doc="The input tranches file describing where to cut the data", required=false)
private String TRANCHE_FILENAME = null;
// For GenotypePhasingEvaluator:
@Argument(fullName = "minPhaseQuality", shortName = "minPQ", doc = "The minimum phasing quality (PQ) score required to consider phasing; [default:0]", required = false)
protected Double minPhaseQuality = 0.0; // accept any positive value of PQ
@Argument(shortName="min", fullName="minimalComparisons", doc="If passed, filters and raw site values won't be computed", required=false)
protected boolean MINIMAL = false;
// --------------------------------------------------------------------------------------------------------------
//
// private walker data
//
// --------------------------------------------------------------------------------------------------------------
/** private class holding all of the information about a single evaluation group (e.g., for eval ROD) */
public class EvaluationContext implements Comparable<EvaluationContext> {
// useful for typing
public String evalTrackName, compTrackName, novelty, filtered;
public boolean enableInterestingSiteCaptures = false;
VariantContextUtils.JexlVCMatchExp selectExp;
Set<VariantEvaluator> evaluations;
public boolean isIgnoringFilters() { return filtered.equals(RAW_SET_NAME); }
public boolean requiresFiltered() { return filtered.equals(FILTERED_SET_NAME); }
public boolean requiresNotFiltered() { return filtered.equals(RETAINED_SET_NAME); }
public boolean isIgnoringNovelty() { return novelty.equals(ALL_SET_NAME); }
public boolean requiresNovel() { return novelty.equals(NOVEL_SET_NAME); }
public boolean requiresKnown() { return novelty.equals(KNOWN_SET_NAME); }
public boolean isSelected() { return selectExp == null; }
public String getDisplayName() {
return getName(CONTEXT_SEPARATOR);
}
public String getJexlName() {
return getName(".");
}
private String getName(String separator) {
return Utils.join(separator, Arrays.asList(evalTrackName, compTrackName, selectExp == null ? "all" : selectExp.name, filtered, novelty));
}
public String toString() { return getDisplayName(); }
public int compareTo(EvaluationContext other) {
return this.getDisplayName().compareTo(other.getDisplayName());
}
public EvaluationContext( String evalName, String compName, String novelty, String filtered, VariantContextUtils.JexlVCMatchExp selectExp ) {
this.evalTrackName = evalName;
this.compTrackName = compName;
this.novelty = novelty;
this.filtered = filtered;
this.selectExp = selectExp;
this.enableInterestingSiteCaptures = selectExp == null;
this.evaluations = instantiateEvalationsSet();
}
}
private List<EvaluationContext> contexts = null;
// lists of all comp and eval ROD track names
private Set<String> compNames = new HashSet<String>();
private Set<String> evalNames = new HashSet<String>();
private List<String> variantEvaluationNames = new ArrayList<String>();
private static String RAW_SET_NAME = "raw";
private static String RETAINED_SET_NAME = "called";
private static String FILTERED_SET_NAME = "filtered";
private static String ALL_SET_NAME = "all";
private static String KNOWN_SET_NAME = "known";
private static String NOVEL_SET_NAME = "novel";
private static String NO_COMP_NAME = "N/A";
private final static String CONTEXT_SEPARATOR = "XXX";
//private final static String CONTEXT_SEPARATOR = "\\.";
private final static String CONTEXT_HEADER = Utils.join(CONTEXT_SEPARATOR, Arrays.asList("eval", "comp", "select", "filter", "novelty"));
private final static int N_CONTEXT_NAME_PARTS = CONTEXT_HEADER.split(CONTEXT_SEPARATOR).length;
private static int[] nameSizes = new int[N_CONTEXT_NAME_PARTS];
static {
int i = 0;
for ( String elt : CONTEXT_HEADER.split(CONTEXT_SEPARATOR) )
nameSizes[i++] = elt.length();
}
// Dynamically determined variantEvaluation classes
private Set<Class<? extends VariantEvaluator>> evaluationClasses = null;
// --------------------------------------------------------------------------------------------------------------
//
// initialize
//
// --------------------------------------------------------------------------------------------------------------
public boolean printInterestingSites() { return writer != null; }
public void initialize() {
if ( dels ) {
ALLOW_VARIANT_CONTEXT_TYPES = EnumSet.of(VariantContext.Type.INDEL, VariantContext.Type.NO_VARIATION);
}
ReportFormat.AcceptableOutputType type = (outputLocation == null) ? ReportFormat.AcceptableOutputType.STREAM : ReportFormat.AcceptableOutputType.FILE;
if (!VE2ReportFactory.isCompatibleWithOutputType(type,reportType))
throw new UserException.CommandLineException("The report format requested is not compatible with your output location. You specified a " + type + " output type which isn't an option for " + reportType);
if ( LIST )
listModulesAndExit();
SAMPLES_LIST = SampleUtils.getSamplesFromCommandLineInput(Arrays.asList(SAMPLES));
determineEvalations();
if ( TRANCHE_FILENAME != null ) {
// we are going to build a few select names automatically from the tranches file
for ( Tranche t : Tranche.readTraches(new File(TRANCHE_FILENAME)) ) {
logger.info("Adding select for all variant above the pCut of : " + t);
SELECT_EXPS.add(String.format(VariantRecalibrator.VQS_LOD_KEY + " >= %.2f", t.minVQSLod));
SELECT_NAMES.add(String.format("FDR-%.2f", t.fdr));
}
}
if ( SELECT_NAMES.size() > 0 ) {
logger.info("Selects: " + SELECT_NAMES);
logger.info("Selects: " + SELECT_EXPS);
}
List<VariantContextUtils.JexlVCMatchExp> selectExps = VariantContextUtils.initializeMatchExps(SELECT_NAMES, SELECT_EXPS);
for ( ReferenceOrderedDataSource d : this.getToolkit().getRodDataSources() ) {
if ( d.getName().startsWith("eval") ) {
evalNames.add(d.getName());
} else if ( d.getName().startsWith("comp") ) {
compNames.add(d.getName());
} else if ( d.getName().startsWith(DbSNPHelper.STANDARD_DBSNP_TRACK_NAME) || d.getName().startsWith("hapmap") ) {
compNames.add(d.getName());
} else {
logger.info("Not evaluating ROD binding " + d.getName());
}
}
// if no comp rod was provided, we still want to be able to do evaluations, so use a default comp name
if ( compNames.size() == 0 )
compNames.add(NO_COMP_NAME);
contexts = initializeEvaluationContexts(evalNames, compNames, selectExps);
determineContextNamePartSizes();
if ( rsIDFile != null ) {
if ( maxRsIDBuild == Integer.MAX_VALUE )
throw new IllegalArgumentException("rsIDFile " + rsIDFile + " was given but associated max RSID build parameter wasn't available");
rsIDsToExclude = getrsIDsToExclude(new File(rsIDFile), maxRsIDBuild);
}
if ( writer != null ) {
Set<String> samples = SampleUtils.getUniqueSamplesFromRods(getToolkit(), evalNames);
final VCFHeader vcfHeader = new VCFHeader(new HashSet<VCFHeaderLine>(), samples);
writer.writeHeader(vcfHeader);
}
}
private void listModulesAndExit() {
List<Class<? extends VariantEvaluator>> veClasses = new PluginManager<VariantEvaluator>( VariantEvaluator.class ).getPlugins();
out.println("\nAvailable eval modules:");
out.println("(Standard modules are starred)");
for (Class<? extends VariantEvaluator> veClass : veClasses)
out.println("\t" + veClass.getSimpleName() + (StandardEval.class.isAssignableFrom(veClass) ? "*" : ""));
out.println();
System.exit(0);
}
private static Set<String> getrsIDsToExclude(File rsIDFile, int maxRsIDBuild) {
List<String> toExclude = new LinkedList<String>();
int n = 1;
try {
for ( String line : new XReadLines(rsIDFile) ) {
String parts[] = line.split(" ");
if ( parts.length != 2 )
throw new UserException.MalformedFile(rsIDFile, "Invalid rsID / build pair at " + n + " line = " + line );
//System.out.printf("line %s %s %s%n", line, parts[0], parts[1]);
if ( Integer.valueOf(parts[1]) > maxRsIDBuild ) {
//System.out.printf("Excluding %s%n", line);
toExclude.add("rs"+parts[0]);
}
n++;
if ( n % 1000000 == 0 )
logger.info(String.format("Read %d rsIDs from rsID -> build file", n));
}
} catch (FileNotFoundException e) {
throw new UserException.CouldNotReadInputFile(rsIDFile, e);
}
logger.info(String.format("Excluding %d of %d (%.2f%%) rsIDs found from builds > %d",
toExclude.size(), n, ((100.0 * toExclude.size())/n), maxRsIDBuild));
return new HashSet<String>(toExclude);
}
private boolean excludeComp(VariantContext vc) {
String id = vc != null && vc.hasID() ? vc.getID() : null;
boolean ex = rsIDsToExclude != null && id != null && rsIDsToExclude.contains(id);
//System.out.printf("Testing id %s ex=%b against %s%n", id, ex, vc);
return ex;
}
private void determineEvalations() {
// create a map for all eval modules for easy lookup
HashMap<String, Class<? extends VariantEvaluator>> classMap = new HashMap<String, Class<? extends VariantEvaluator>>();
for ( Class<? extends VariantEvaluator> c : new PluginManager<VariantEvaluator>( VariantEvaluator.class ).getPlugins() )
classMap.put(c.getSimpleName(), c);
evaluationClasses = new HashSet<Class<? extends VariantEvaluator>>();
// by default, use standard eval modules
if ( !NO_STANDARD ) {
for ( Class<? extends StandardEval> myClass : new PluginManager<StandardEval>( StandardEval.class ).getPlugins() ) {
if ( classMap.containsKey(myClass.getSimpleName()) )
evaluationClasses.add(classMap.get(myClass.getSimpleName()));
}
}
// get the specific classes provided
for ( String module : modulesToUse ) {
if ( !classMap.containsKey(module) )
throw new UserException.CommandLineException("Module " + module + " could not be found; please check that you have specified the class name correctly");
evaluationClasses.add(classMap.get(module));
}
for ( VariantEvaluator e : instantiateEvalationsSet() ) {
// for collecting purposes
variantEvaluationNames.add(e.getName());
logger.debug("Including VariantEvaluator " + e.getName() + " of class " + e.getClass());
}
Collections.sort(variantEvaluationNames);
}
private <T> List<T> append(List<T> selectExps, T elt) {
List<T> l = new ArrayList<T>(selectExps);
l.add(elt);
return l;
}
private List<EvaluationContext> initializeEvaluationContexts(Set<String> evalNames, Set<String> compNames, List<VariantContextUtils.JexlVCMatchExp> selectExps) {
List<EvaluationContext> contexts = new ArrayList<EvaluationContext>();
// todo -- add another for loop for each sample (be smart about the selection here -
// honor specifications of just one or a few samples), and put an "all" in here so
// that we don't lose multi-sample evaluations
List<String> filterTypes = MINIMAL ? Arrays.asList(RETAINED_SET_NAME) : Arrays.asList(RAW_SET_NAME, RETAINED_SET_NAME, FILTERED_SET_NAME);
selectExps = append(selectExps, null);
for ( String evalName : evalNames ) {
for ( String compName : compNames ) {
for ( VariantContextUtils.JexlVCMatchExp e : selectExps ) {
for ( String filteredName : filterTypes ) {
for ( String novelty : Arrays.asList(ALL_SET_NAME, KNOWN_SET_NAME, NOVEL_SET_NAME) ) {
EvaluationContext context = new EvaluationContext(evalName, compName, novelty, filteredName, e);
contexts.add(context);
}
}
}
}
}
Collections.sort(contexts);
return contexts;
}
private Set<VariantEvaluator> instantiateEvalationsSet() {
Set<VariantEvaluator> evals = new HashSet<VariantEvaluator>();
Object[] args = new Object[]{this};
Class<?>[] argTypes = new Class<?>[]{VariantEvalWalker.class};
for ( Class<? extends VariantEvaluator> c : evaluationClasses ) {
try {
Constructor<? extends VariantEvaluator> constructor = c.getConstructor(argTypes);
VariantEvaluator eval = constructor.newInstance(args);
evals.add(eval);
} catch (Exception e) {
throw new DynamicClassResolutionException(c, e);
}
}
return evals;
}
private boolean captureInterestingSitesOfEvalSet(EvaluationContext group) {
//System.out.printf("checking %s%n", name);
return group.requiresNotFiltered() && group.isIgnoringNovelty();
}
// --------------------------------------------------------------------------------------------------------------
//
// map
//
// --------------------------------------------------------------------------------------------------------------
// todo -- call a single function to build a map from track name -> variant context / null for all
// -- eval + comp names. Use this data structure to get data throughout rest of the loops here
public Integer map(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
//System.out.printf("map at %s with %d skipped%n", context.getLocation(), context.getSkippedBases());
Map<String, VariantContext> vcs = getVariantContexts(ref, tracker, context);
//System.out.println("vcs has size "+vcs.size());
//Collection<VariantContext> comps = getCompVariantContexts(tracker, context);
// to enable walking over pairs where eval or comps have no elements
for ( EvaluationContext group : contexts ) {
VariantContext vc = vcs.get(group.evalTrackName);
//logger.debug(String.format("Updating %s with variant", vc));
Set<VariantEvaluator> evaluations = group.evaluations;
boolean evalWantsVC = applyVCtoEvaluation(vc, vcs, group);
VariantContext interestingVC = vc;
List<String> interestingReasons = new ArrayList<String>();
for ( VariantEvaluator evaluation : evaluations ) {
synchronized ( evaluation ) {
if ( evaluation.enabled() ) {
// we always call update0 in case the evaluation tracks things like number of bases covered
evaluation.update0(tracker, ref, context);
// the other updateN methods don't see a null context
if ( tracker == null )
continue;
// now call the single or paired update function
switch ( evaluation.getComparisonOrder() ) {
case 1:
if ( evalWantsVC && vc != null ) {
String interesting = evaluation.update1(vc, tracker, ref, context, group);
if ( interesting != null ) interestingReasons.add(interesting);
}
break;
case 2:
VariantContext comp = vcs.get(group.compTrackName);
if ( comp != null &&
minCompQualScore != NO_MIN_QUAL_SCORE &&
comp.hasNegLog10PError() &&
comp.getNegLog10PError() < (minCompQualScore / 10.0) )
comp = null;
String interesting = evaluation.update2( evalWantsVC ? vc : null, comp, tracker, ref, context, group );
/** TODO
-- for Eric: Fix me (current implementation causes GenotypeConcordance
to treat sites that don't match JEXL as no-calls)
String interesting = null;
if (evalWantsVC)
{
interesting = evaluation.update2( evalWantsVC ? vc : null, comp, tracker, ref, context, group );
}
**/
if ( interesting != null ) {
interestingVC = interestingVC == null ? ( vc == null ? comp : vc ) : interestingVC;
interestingReasons.add(interesting);
}
break;
default:
throw new ReviewedStingException("BUG: Unexpected evaluation order " + evaluation);
}
}
}
}
if ( tracker != null && group.enableInterestingSiteCaptures && captureInterestingSitesOfEvalSet(group) )
writeInterestingSite(interestingReasons, interestingVC, ref.getBase());
}
return 0;
}
private void writeInterestingSite(List<String> interestingReasons, VariantContext vc, byte ref) {
if ( vc != null && writer != null && interestingReasons.size() > 0 ) {
// todo -- the vc == null check is because you can be interesting because you are a FN, and so VC == null
MutableVariantContext mvc = new MutableVariantContext(vc);
for ( String why : interestingReasons ) {
String key, value;
String[] parts = why.split("=");
switch ( parts.length ) {
case 1:
key = parts[0];
value = "1";
break;
case 2:
key = parts[0];
value = parts[1];
break;
default:
throw new IllegalStateException("BUG: saw a interesting site reason sting with multiple = signs " + why);
}
mvc.putAttribute(key, value);
}
writer.add(mvc, ref);
//interestingReasons.clear();
}
}
private static Set<String> seenJEXLExceptions = new HashSet<String>();
private boolean applyVCtoEvaluation(VariantContext vc, Map<String, VariantContext> vcs, EvaluationContext group) {
if ( vc == null )
return true;
if ( minQualScore != NO_MIN_QUAL_SCORE &&
vc.hasNegLog10PError() &&
vc.getNegLog10PError() < (minQualScore / 10.0) ) {
//System.out.printf("exclude %s%n", vc);
return false;
}
if ( group.requiresFiltered() && vc.isNotFiltered() )
return false;
if ( group.requiresNotFiltered() && vc.isFiltered() )
return false;
boolean vcKnown = vcIsKnown(vc, vcs, KNOWN_NAMES);
if ( group.requiresKnown() && ! vcKnown )
return false;
else if ( group.requiresNovel() && vcKnown )
return false;
if ( group.selectExp != null ) {
try {
if ( ! VariantContextUtils.match(vc, group.selectExp) )
return false;
} catch ( RuntimeException e ) {
if ( ! seenJEXLExceptions.contains(group.selectExp.name) ) {
seenJEXLExceptions.add(group.selectExp.name);
logger.warn("JEXL evaluation error for SELECT " + group.selectExp.name + ": " + e.getMessage() +
"; this may be an error or may simply result from some variants not having INFO fields keys " +
"referenced in the JEXL expressions. Variants generating exceptions will *NOT* be matched " +
"by the expression. Occurred with variant " + vc);
}
return false;
}
}
// nothing invalidated our membership in this set
return true;
}
private boolean vcIsKnown(VariantContext vc, Map<String, VariantContext> vcs, String[] knownNames ) {
for ( String knownName : knownNames ) {
VariantContext known = vcs.get(knownName);
if ( known != null && known.isNotFiltered() && known.getType() == vc.getType() ) {
return true;
}
}
return false;
}
// can't handle this situation
// todo -- warning, this leads to some missing SNPs at complex loci, such as:
// todo -- 591 1 841619 841620 rs4970464 0 - A A -/C/T genomic mixed unknown 0 0 near-gene-3 exact 1
// todo -- 591 1 841619 841620 rs62677860 0 + A A C/T genomic single unknown 0 0 near-gene-3 exact 1
//
//logger.info(String.format("Ignore second+ events at locus %s in rod %s => rec is %s", context.getLocation(), rodList.getName(), rec));
private Map<String, VariantContext> getVariantContexts(ReferenceContext ref, RefMetaDataTracker tracker, AlignmentContext context) {
// todo -- we need to deal with dbSNP where there can be multiple records at the same start site. A potential solution is to
// todo -- allow the variant evaluation to specify the type of variants it wants to see and only take the first such record at a site
Map<String, VariantContext> bindings = new HashMap<String, VariantContext>();
if ( tracker != null ) {
//System.out.println("Tracker is not null");
bindVariantContexts(ref, bindings, evalNames, tracker, context, false);
bindVariantContexts(ref, bindings, compNames, tracker, context, true);
}
return bindings;
}
private void bindVariantContexts(ReferenceContext ref, Map<String, VariantContext> map, Collection<String> names,
RefMetaDataTracker tracker, AlignmentContext context, boolean allowExcludes ) {
for ( String name : names ) {
Collection<VariantContext> contexts = tracker.getVariantContexts(ref, name, ALLOW_VARIANT_CONTEXT_TYPES, context.getLocation(), true, true);
if ( contexts.size() > 1 )
throw new UserException.CommandLineException("Found multiple variant contexts found in " + name + " at " + context.getLocation() + "; VariantEval assumes one variant per position");
VariantContext vc = contexts.size() == 1 ? contexts.iterator().next() : null;
if ( SAMPLES_LIST.size() > 0 && vc != null ) {
boolean hasGenotypes = vc.hasGenotypes(SAMPLES_LIST);
if ( hasGenotypes ) {
//if ( ! name.equals("eval") ) logger.info(String.format("subsetting VC %s", vc));
vc = vc.subContextFromGenotypes(vc.getGenotypes(SAMPLES_LIST).values());
HashMap<String,Object> newAts = new HashMap<String,Object>(vc.getAttributes());
VariantContextUtils.calculateChromosomeCounts(vc,newAts,true);
vc = VariantContext.modifyAttributes(vc,newAts);
logger.debug(String.format("VC %s subset to %s AC%n",vc.getSource(),vc.getAttributeAsString(VCFConstants.ALLELE_COUNT_KEY)));
//if ( ! name.equals("eval") ) logger.info(String.format(" => VC %s", vc));
} else if ( !hasGenotypes && !name.equals("dbsnp") ) {
throw new UserException(String.format("Genotypes for the variant context %s do not contain all the provided samples %s",vc.getSource(), getMissingSamples(SAMPLES_LIST,vc)));
}
}
map.put(name, allowExcludes && excludeComp(vc) ? null : vc);
}
}
private static String getMissingSamples(Collection<String> soughtSamples, VariantContext vc) {
StringBuffer buf = new StringBuffer();
buf.append("Missing samples are:");
for ( String s : soughtSamples ) {
if ( ! vc.getGenotypes().keySet().contains(s) ) {
buf.append(String.format("%n%s",s));
}
}
return buf.toString();
}
// --------------------------------------------------------------------------------------------------------------
//
// reduce
//
// --------------------------------------------------------------------------------------------------------------
public Integer reduceInit() {
return 0;
}
public Integer reduce(Integer point, Integer sum) {
return point + sum;
}
public Integer treeReduce(Integer point, Integer sum) {
return point + sum;
}
public VariantEvaluator getEvalByName(String name, Set<VariantEvaluator> s) {
for ( VariantEvaluator e : s )
if ( e.getName().equals(name) )
return e;
return null;
}
private void determineContextNamePartSizes() {
for ( EvaluationContext group : contexts ) {
String[] parts = group.getDisplayName().split(CONTEXT_SEPARATOR);
if ( parts.length != N_CONTEXT_NAME_PARTS ) {
throw new ReviewedStingException("Unexpected number of eval name parts " + group.getDisplayName() + " length = " + parts.length + ", expected " + N_CONTEXT_NAME_PARTS);
} else {
for ( int i = 0; i < parts.length; i++ )
nameSizes[i] = Math.max(nameSizes[i], parts[i].length());
}
}
}
public void onTraversalDone(Integer result) {
writeReport();
validateContext();
}
/**
* Writes the report out to disk.
*/
private void writeReport() {
// our report mashaller
ReportMarshaller marshaller;
// create the report marshaller early, so that we can fail-fast if something is wrong with the output sources
if (outputLocation == null)
marshaller = VE2ReportFactory.createMarhsaller(new OutputStreamWriter(out), reportType, createExtraOutputTags());
else
marshaller = VE2ReportFactory.createMarhsaller(outputLocation, reportType, createExtraOutputTags());
for ( String evalName : variantEvaluationNames ) {
for ( EvaluationContext group : contexts ) {
VariantEvaluator eval = getEvalByName(evalName, group.evaluations);
// finalize the evaluation
eval.finalizeEvaluation();
if ( eval.enabled() )
marshaller.write(createPrependNodeList(group),eval);
}
}
marshaller.close();
}
/**
* Validates the JEXL expressions and throws an exception if they do not all return true.
*/
private void validateContext() {
if (SUMMARY_EXPS.size() + VALIDATE_EXPS.size() == 0)
return;
JexlContext jc = new MapContext();
for (EvaluationContext context : contexts)
for (VariantEvaluator eval: context.evaluations)
jc.set(context.getJexlName() + "." + eval.getName(), eval);
Uberspect uberspect = new UberspectImpl(LogFactory.getLog(JexlEngine.class)) {
/** Gets the field, even if the field was non-public. */
@Override
public Field getField(Object obj, String name, JexlInfo info) {
Field result = super.getField(obj, name, info);
if (result == null && obj != null) {
Class<?> clazz = obj instanceof Class<?> ? (Class<?>)obj : obj.getClass();
try {
// TODO: Default UberspectImpl uses an internal field cache by class type
result = clazz.getDeclaredField(name);
result.setAccessible(true);
} catch (NoSuchFieldException nsfe) {
/* ignore */
}
}
return result;
}
};
JexlEngine jexl = new JexlEngine(uberspect, null, null, null);
for (String expression: SUMMARY_EXPS) {
Object jexlResult = jexl.createExpression(expression).evaluate(jc);
logger.info("Summary: " + expression + " = " + jexlResult);
}
List<String> failedExpressions = new ArrayList<String>();
for (String expression: VALIDATE_EXPS) {
// ex: evalYRI.compYRI.all.called.novel.titv.tiTvRatio > 1.0
Object jexlResult = jexl.createExpression(expression).evaluate(jc);
boolean pass = Boolean.TRUE.equals(jexlResult);
if (!pass) {
logger.error("FAIL: " + expression);
failedExpressions.add(expression);
} else if (logger.isDebugEnabled()) {
logger.debug("PASS: " + expression);
}
}
int failed = failedExpressions.size();
int total = VALIDATE_EXPS.size();
logger.info(String.format("Validations: Total %s, Passed %s, Failed %s", total, (total-failed), failed));
if (failed > 0) {
StringBuilder message = new StringBuilder("The validation expressions below did not return true. Please check the report output for more info.");
for (String expression: failedExpressions)
message.append(String.format("%n ")).append(expression);
throw new UserException(message.toString());
}
}
/**
* create some additional output lines about the analysis
* @return a list of nodes to attach to the report as tags
*/
private List<Node> createExtraOutputTags() {
List<Node> list = new ArrayList<Node>();
list.add(new Node("reference file",getToolkit().getArguments().referenceFile.getName(),"The reference sequence file"));
for (String binding : getToolkit().getArguments().RODBindings)
list.add(new Node("ROD binding",binding,"The reference sequence file"));
return list;
}
/**
* given the evaluation name, and the context, create the list of pre-pended nodes for the output system.
* Currently it expects the the following list: jexl_expression, evaluation_name, comparison_name, filter_name,
* novelty_name
* @param group the evaluation context
* @return a list of Nodes to prepend the analysis module output with
*/
private List<Node> createPrependNodeList(EvaluationContext group) {
// add the branching nodes: jexl expression, comparison track, eval track etc
Node jexlNode = new Node("jexl_expression",(group.selectExp != null) ? group.selectExp.name : "none","The jexl filtering expression");
Node compNode = new Node("comparison_name",group.compTrackName,"The comparison track name");
Node evalNode = new Node("evaluation_name",group.evalTrackName,"The evaluation name");
Node filterNode = new Node("filter_name",group.filtered,"The filter name");
Node noveltyNode = new Node("novelty_name",group.novelty,"The novelty name");
// the ordering is important below, this is the order the columns will appear in any output format
return Arrays.asList(evalNode,compNode,jexlNode,filterNode,noveltyNode);
}
//
// utility functions
//
/**
* Takes an eval generated VariantContext and attempts to return the number of samples in the
* VC. If there are genotypes, it returns that value first. Otherwise it returns nSamples, if
* provided. Returns -1 if no sample counts can be obtained.
*
* @param eval
* @return
*/
public int getNSamplesForEval( VariantContext eval ) {
return eval.hasGenotypes() ? eval.getNSamples() : nSamples;
}
public Logger getLogger() { return logger; }
}

101
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/VariantEvaluator.java
View File

@ -1,101 +0,0 @@
package org.broadinstitute.sting.gatk.walkers.varianteval;
import org.apache.log4j.Logger;
import org.broad.tribble.util.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
/**
* The Broad Institute
* SOFTWARE COPYRIGHT NOTICE AGREEMENT
* This software and its documentation are copyright 2009 by the
* Broad Institute/Massachusetts Institute of Technology. All rights are reserved.
* <p/>
* This software is supplied without any warranty or guaranteed support whatsoever. Neither
* the Broad Institute nor MIT can be responsible for its use, misuse, or functionality.
*/
public abstract class VariantEvaluator {
// protected boolean accumulateInterestingSites = false, printInterestingSites = false;
// protected String interestingSitePrefix = null;
// protected boolean processedASite = false;
// protected List<VariantContext> interestingSites = new ArrayList<VariantContext>();
// do we want to keep track of things that are interesting
// public void accumulateInterestingSites(boolean enable) { accumulateInterestingSites = enable; }
// public void printInterestingSites(String prefix) { printInterestingSites = true; interestingSitePrefix = prefix; }
// public boolean isAccumulatingInterestingSites() { return accumulateInterestingSites; }
// public List<VariantContext> getInterestingSites() { return interestingSites; }
// protected void addInterestingSite(String why, VariantContext vc) {
// if ( accumulateInterestingSites )
// interestingSites.add(vc);
// if ( printInterestingSites )
// System.out.printf("%40s %s%n", interestingSitePrefix, why);
// }
public abstract String getName();
protected VariantEvalWalker veWalker = null;
public VariantEvaluator(VariantEvalWalker parent) {
veWalker = parent;
// don't do anything
}
protected VariantEvalWalker getVEWalker() {
return veWalker;
}
protected Logger getLogger() {
return veWalker.getLogger();
}
public abstract boolean enabled();
//public boolean processedAnySites() { return processedASite; }
//protected void markSiteAsProcessed() { processedASite = true; }
// Should return the number of VariantContexts expected as inputs to update. Can be 1 or 2
public abstract int getComparisonOrder();
// called at all sites, regardless of eval context itself; useful for counting processed bases
public void update0(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) { }
public String update1(VariantContext vc1, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
return null;
}
public String update1(VariantContext vc1, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context, VariantEvalWalker.EvaluationContext group) {
return update1(vc1, tracker, ref, context);
}
public String update2(VariantContext vc1, VariantContext vc2, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
return null;
}
public String update2(VariantContext vc1, VariantContext vc2, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context, VariantEvalWalker.EvaluationContext group) {
return update2(vc1, vc2, tracker, ref, context);
}
/**
* override this method for any finalization of calculations after the analysis is completed
*/
public void finalizeEvaluation() {}
//
// useful common utility routines
//
protected double rate(long n, long d) {
return n / (1.0 * Math.max(d, 1));
}
protected long inverseRate(long n, long d) {
return n == 0 ? 0 : d / Math.max(n, 1);
}
protected double ratio(long num, long denom) {
return ((double)num) / (Math.max(denom, 1));
}
}

258
java/src/org/broadinstitute/sting/gatk/walkers/varianteval/VariantQualityScore.java
View File

@ -1,258 +0,0 @@
/*
* Copyright (c) 2010 The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
* THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.gatk.walkers.varianteval;
import org.broad.tribble.util.variantcontext.Allele;
import org.broad.tribble.util.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.variantcontext.VariantContextUtils;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.report.tags.Analysis;
import org.broadinstitute.sting.utils.report.tags.DataPoint;
import org.broadinstitute.sting.utils.report.utils.TableType;
import org.broadinstitute.sting.utils.collections.Pair;
import java.util.ArrayList;
import java.util.HashMap;
/**
* @author rpoplin
* @since Apr 6, 2010
*/
@Analysis(name = "Variant Quality Score", description = "Shows various stats of sets of variants binned by variant quality score")
public class VariantQualityScore extends VariantEvaluator {
// a mapping from quality score histogram bin to Ti/Tv ratio
@DataPoint(name="TiTv by Quality", description = "the Ti/Tv ratio broken out by variant quality")
TiTvStats titvStats = null;
@DataPoint(name="Quality by Allele Count", description = "average variant quality for each allele count")
AlleleCountStats alleleCountStats = null;
static class TiTvStats implements TableType {
final static int NUM_BINS = 20;
final HashMap<Integer, Pair<Long,Long>> qualByIsTransition = new HashMap<Integer, Pair<Long,Long>>(); // A hashMap holds all the qualities until we are able to bin them appropriately
final long transitionByQuality[] = new long[NUM_BINS];
final long transversionByQuality[] = new long[NUM_BINS];
final double titvByQuality[] = new double[NUM_BINS]; // the final ti/tv sets that get reported out
public Object[] getRowKeys() {
return new String[]{"sample"};
}
public Object[] getColumnKeys() {
final String columnKeys[] = new String[NUM_BINS];
for( int iii = 0; iii < NUM_BINS; iii++ ) {
columnKeys[iii] = "titvBin" + iii;
}
return columnKeys;
}
public String getName() {
return "TiTvStats";
}
public String getCell(int x, int y) {
return String.valueOf(titvByQuality[y]);
}
public String toString() {
StringBuffer returnString = new StringBuffer();
// output the ti/tv array
returnString.append("titvByQuality: ");
for( int iii = 0; iii < NUM_BINS; iii++ ) {
returnString.append(titvByQuality[iii]);
returnString.append(" ");
}
return returnString.toString();
}
public void incrValue( final double qual, final boolean isTransition ) {
final Integer qualKey = Math.round((float) qual);
final long numTransition = (isTransition ? 1L : 0L);
final long numTransversion = (isTransition ? 0L : 1L);
if( qualByIsTransition.containsKey(qualKey) ) {
Pair<Long,Long> transitionPair = qualByIsTransition.get(qualKey);
transitionPair.set(transitionPair.getFirst() + numTransition, transitionPair.getSecond() + numTransversion);
qualByIsTransition.put(qualKey, transitionPair);
} else {
qualByIsTransition.put(qualKey, new Pair<Long,Long>(numTransition,numTransversion));
}
}
public void organizeTiTvTables() {
for( int iii = 0; iii < NUM_BINS; iii++ ) {
transitionByQuality[iii] = 0L;
transversionByQuality[iii] = 0L;
titvByQuality[iii] = 0.0;
}
int maxQual = 0;
// Calculate the maximum quality score in order to normalize and histogram
for( final Integer qual : qualByIsTransition.keySet() ) {
if( qual > maxQual ) {
maxQual = qual;
}
}
final double binSize = ((double)maxQual) / ((double) (NUM_BINS-1));
for( final Integer qual : qualByIsTransition.keySet() ) {
final int index = (int)Math.floor( ((double) qual) / binSize );
if( index >= 0 ) { // BUGBUG: why is there overflow here?
Pair<Long,Long> transitionPair = qualByIsTransition.get(qual);
transitionByQuality[index] += transitionPair.getFirst();
transversionByQuality[index] += transitionPair.getSecond();
}
}
for( int iii = 0; iii < NUM_BINS; iii++ ) {
if( transitionByQuality[iii] + transversionByQuality[iii] > 800L ) { // need to have a sufficient number of variants to get a useful Ti/Tv ratio
titvByQuality[iii] = ((double) transitionByQuality[iii]) / ((double) transversionByQuality[iii]);
} else {
titvByQuality[iii] = 0.0;
}
}
}
}
class AlleleCountStats implements TableType {
final HashMap<Integer, ArrayList<Double>> qualityListMap = new HashMap<Integer, ArrayList<Double>>();
final HashMap<Integer, Double> qualityMap = new HashMap<Integer, Double>();
public Object[] getRowKeys() {
final int NUM_BINS = qualityListMap.keySet().size();
final String rowKeys[] = new String[NUM_BINS];
int iii = 0;
for( final Integer key : qualityListMap.keySet() ) {
rowKeys[iii] = "AC" + key;
iii++;
}
return rowKeys;
}
public Object[] getColumnKeys() {
return new String[]{"alleleCount","avgQual"};
}
public String getName() {
return "AlleleCountStats";
}
public String getCell(int x, int y) {
int iii = 0;
for( final Integer key : qualityListMap.keySet() ) {
if(iii == x) {
if(y == 0) { return String.valueOf(key); }
else { return String.valueOf(qualityMap.get(key)); }
}
iii++;
}
return null;
}
public String toString() {
String returnString = "";
// output the quality map
returnString += "AlleleCountStats: ";
//for( int iii = 0; iii < NUM_BINS; iii++ ) {
// returnString += titvByQuality[iii] + " ";
//}
return returnString;
}
public void incrValue( final double qual, final int alleleCount ) {
ArrayList<Double> list = qualityListMap.get(alleleCount);
if(list==null) { list = new ArrayList<Double>(); }
list.add(qual);
qualityListMap.put(alleleCount, list);
}
public void organizeAlleleCountTables() {
for( final Integer key : qualityListMap.keySet() ) {
final ArrayList<Double> list = qualityListMap.get(key);
double meanQual = 0.0;
final double numQuals = (double)list.size();
for( Double qual : list ) {
meanQual += qual / numQuals;
}
qualityMap.put(key, meanQual);
}
}
}
public VariantQualityScore(VariantEvalWalker parent) {
super(parent);
}
public String getName() {
return "VariantQualityScore";
}
public int getComparisonOrder() {
return 1; // we only need to see each eval track
}
public boolean enabled() {
return true;
}
public String toString() {
return getName();
}
public String update1(VariantContext eval, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
final String interesting = null;
if( eval != null && eval.isSNP() && eval.isBiallelic() ) { //BUGBUG: only counting biallelic sites (revisit what to do with triallelic sites)
if( titvStats == null ) { titvStats = new TiTvStats(); }
titvStats.incrValue(eval.getPhredScaledQual(), VariantContextUtils.isTransition(eval));
if( alleleCountStats == null ) { alleleCountStats = new AlleleCountStats(); }
int alternateAlleleCount = 0;
for (final Allele a : eval.getAlternateAlleles()) {
alternateAlleleCount += eval.getChromosomeCount(a);
}
alleleCountStats.incrValue(eval.getPhredScaledQual(), alternateAlleleCount);
}
return interesting; // This module doesn't capture any interesting sites, so return null
}
public void finalizeEvaluation() {
if( titvStats != null ) {
titvStats.organizeTiTvTables();
}
if( alleleCountStats != null ) {
alleleCountStats.organizeAlleleCountTables();
}
}
}