Handle null VCs (which can arise when indels are present in the file)

This commit is contained in:
Eric Banks 2012-04-23 10:30:00 -04:00
parent adad76b36f
commit 2761da975e
1 changed files with 49 additions and 44 deletions

View File

@ -744,62 +744,67 @@ public class PhaseByTransmission extends RodWalker<HashMap<Byte,Integer>, HashMa
String mvfLine;
if (tracker != null) {
VariantContext vc = tracker.getFirstValue(variantCollection.variants, context.getLocation());
VariantContextBuilder builder = new VariantContextBuilder(vc);
if (tracker == null)
return metricsCounters;
GenotypesContext genotypesContext = GenotypesContext.copy(vc.getGenotypes());
for (Sample sample : trios) {
Genotype mother = vc.getGenotype(sample.getMaternalID());
Genotype father = vc.getGenotype(sample.getPaternalID());
Genotype child = vc.getGenotype(sample.getID());
final VariantContext vc = tracker.getFirstValue(variantCollection.variants, context.getLocation());
if (vc == null)
return metricsCounters;
//Keep only trios and parent/child pairs
if(mother == null && father == null || child == null)
continue;
final VariantContextBuilder builder = new VariantContextBuilder(vc);
ArrayList<Genotype> trioGenotypes = new ArrayList<Genotype>(3);
final int mvCount = phaseTrioGenotypes(vc.getReference(), vc.getAltAlleleWithHighestAlleleCount(), mother, father, child,trioGenotypes);
final GenotypesContext genotypesContext = GenotypesContext.copy(vc.getGenotypes());
for (Sample sample : trios) {
Genotype mother = vc.getGenotype(sample.getMaternalID());
Genotype father = vc.getGenotype(sample.getPaternalID());
Genotype child = vc.getGenotype(sample.getID());
Genotype phasedMother = trioGenotypes.get(0);
Genotype phasedFather = trioGenotypes.get(1);
Genotype phasedChild = trioGenotypes.get(2);
//Keep only trios and parent/child pairs
if(mother == null && father == null || child == null)
continue;
//Fill the genotype map with the new genotypes and increment metrics counters
genotypesContext.replace(phasedChild);
if(mother != null){
genotypesContext.replace(phasedMother);
if(father != null){
genotypesContext.replace(phasedFather);
updateTrioMetricsCounters(phasedMother,phasedFather,phasedChild,mvCount,metricsCounters);
mvfLine = String.format("%s\t%d\t%s\t%s\t%s\t%s\t%s:%s:%s:%s\t%s:%s:%s:%s\t%s:%s:%s:%s",vc.getChr(),vc.getStart(),vc.getFilters(),vc.getAttribute(VCFConstants.ALLELE_COUNT_KEY),sample.toString(),phasedMother.getAttribute(TRANSMISSION_PROBABILITY_TAG_NAME),phasedMother.getGenotypeString(),phasedMother.getAttribute(VCFConstants.DEPTH_KEY),phasedMother.getAttribute("AD"),phasedMother.getLikelihoods().toString(),phasedFather.getGenotypeString(),phasedFather.getAttribute(VCFConstants.DEPTH_KEY),phasedFather.getAttribute("AD"),phasedFather.getLikelihoods().toString(),phasedChild.getGenotypeString(),phasedChild.getAttribute(VCFConstants.DEPTH_KEY),phasedChild.getAttribute("AD"),phasedChild.getLikelihoods().toString());
if(!(phasedMother.getType()==mother.getType() && phasedFather.getType()==father.getType() && phasedChild.getType()==child.getType()))
metricsCounters.put(NUM_GENOTYPES_MODIFIED,metricsCounters.get(NUM_GENOTYPES_MODIFIED)+1);
}
else{
updatePairMetricsCounters(phasedMother,phasedChild,mvCount,metricsCounters);
if(!(phasedMother.getType()==mother.getType() && phasedChild.getType()==child.getType()))
metricsCounters.put(NUM_GENOTYPES_MODIFIED,metricsCounters.get(NUM_GENOTYPES_MODIFIED)+1);
mvfLine = String.format("%s\t%d\t%s\t%s\t%s\t%s\t%s:%s:%s:%s\t.:.:.:.\t%s:%s:%s:%s",vc.getChr(),vc.getStart(),vc.getFilters(),vc.getAttribute(VCFConstants.ALLELE_COUNT_KEY),sample.toString(),phasedMother.getAttribute(TRANSMISSION_PROBABILITY_TAG_NAME),phasedMother.getGenotypeString(),phasedMother.getAttribute(VCFConstants.DEPTH_KEY),phasedMother.getAttribute("AD"),phasedMother.getLikelihoods().toString(),phasedChild.getGenotypeString(),phasedChild.getAttribute(VCFConstants.DEPTH_KEY),phasedChild.getAttribute("AD"),phasedChild.getLikelihoods().toString());
}
ArrayList<Genotype> trioGenotypes = new ArrayList<Genotype>(3);
final int mvCount = phaseTrioGenotypes(vc.getReference(), vc.getAltAlleleWithHighestAlleleCount(), mother, father, child,trioGenotypes);
Genotype phasedMother = trioGenotypes.get(0);
Genotype phasedFather = trioGenotypes.get(1);
Genotype phasedChild = trioGenotypes.get(2);
//Fill the genotype map with the new genotypes and increment metrics counters
genotypesContext.replace(phasedChild);
if(mother != null){
genotypesContext.replace(phasedMother);
if(father != null){
genotypesContext.replace(phasedFather);
updateTrioMetricsCounters(phasedMother,phasedFather,phasedChild,mvCount,metricsCounters);
mvfLine = String.format("%s\t%d\t%s\t%s\t%s\t%s\t%s:%s:%s:%s\t%s:%s:%s:%s\t%s:%s:%s:%s",vc.getChr(),vc.getStart(),vc.getFilters(),vc.getAttribute(VCFConstants.ALLELE_COUNT_KEY),sample.toString(),phasedMother.getAttribute(TRANSMISSION_PROBABILITY_TAG_NAME),phasedMother.getGenotypeString(),phasedMother.getAttribute(VCFConstants.DEPTH_KEY),phasedMother.getAttribute("AD"),phasedMother.getLikelihoods().toString(),phasedFather.getGenotypeString(),phasedFather.getAttribute(VCFConstants.DEPTH_KEY),phasedFather.getAttribute("AD"),phasedFather.getLikelihoods().toString(),phasedChild.getGenotypeString(),phasedChild.getAttribute(VCFConstants.DEPTH_KEY),phasedChild.getAttribute("AD"),phasedChild.getLikelihoods().toString());
if(!(phasedMother.getType()==mother.getType() && phasedFather.getType()==father.getType() && phasedChild.getType()==child.getType()))
metricsCounters.put(NUM_GENOTYPES_MODIFIED,metricsCounters.get(NUM_GENOTYPES_MODIFIED)+1);
}
else{
genotypesContext.replace(phasedFather);
updatePairMetricsCounters(phasedFather,phasedChild,mvCount,metricsCounters);
if(!(phasedFather.getType()==father.getType() && phasedChild.getType()==child.getType()))
updatePairMetricsCounters(phasedMother,phasedChild,mvCount,metricsCounters);
if(!(phasedMother.getType()==mother.getType() && phasedChild.getType()==child.getType()))
metricsCounters.put(NUM_GENOTYPES_MODIFIED,metricsCounters.get(NUM_GENOTYPES_MODIFIED)+1);
mvfLine = String.format("%s\t%d\t%s\t%s\t%s\t%s\t.:.:.:.\t%s:%s:%s:%s\t%s:%s:%s:%s",vc.getChr(),vc.getStart(),vc.getFilters(),vc.getAttribute(VCFConstants.ALLELE_COUNT_KEY),sample.toString(),phasedFather.getAttribute(TRANSMISSION_PROBABILITY_TAG_NAME),phasedFather.getGenotypeString(),phasedFather.getAttribute(VCFConstants.DEPTH_KEY),phasedFather.getAttribute("AD"),phasedFather.getLikelihoods().toString(),phasedChild.getGenotypeString(),phasedChild.getAttribute(VCFConstants.DEPTH_KEY),phasedChild.getAttribute("AD"),phasedChild.getLikelihoods().toString());
mvfLine = String.format("%s\t%d\t%s\t%s\t%s\t%s\t%s:%s:%s:%s\t.:.:.:.\t%s:%s:%s:%s",vc.getChr(),vc.getStart(),vc.getFilters(),vc.getAttribute(VCFConstants.ALLELE_COUNT_KEY),sample.toString(),phasedMother.getAttribute(TRANSMISSION_PROBABILITY_TAG_NAME),phasedMother.getGenotypeString(),phasedMother.getAttribute(VCFConstants.DEPTH_KEY),phasedMother.getAttribute("AD"),phasedMother.getLikelihoods().toString(),phasedChild.getGenotypeString(),phasedChild.getAttribute(VCFConstants.DEPTH_KEY),phasedChild.getAttribute("AD"),phasedChild.getLikelihoods().toString());
}
//Report violation if set so
//TODO: ADAPT FOR PAIRS TOO!!
if(mvCount>0 && mvFile != null)
mvFile.println(mvfLine);
}
else{
genotypesContext.replace(phasedFather);
updatePairMetricsCounters(phasedFather,phasedChild,mvCount,metricsCounters);
if(!(phasedFather.getType()==father.getType() && phasedChild.getType()==child.getType()))
metricsCounters.put(NUM_GENOTYPES_MODIFIED,metricsCounters.get(NUM_GENOTYPES_MODIFIED)+1);
mvfLine = String.format("%s\t%d\t%s\t%s\t%s\t%s\t.:.:.:.\t%s:%s:%s:%s\t%s:%s:%s:%s",vc.getChr(),vc.getStart(),vc.getFilters(),vc.getAttribute(VCFConstants.ALLELE_COUNT_KEY),sample.toString(),phasedFather.getAttribute(TRANSMISSION_PROBABILITY_TAG_NAME),phasedFather.getGenotypeString(),phasedFather.getAttribute(VCFConstants.DEPTH_KEY),phasedFather.getAttribute("AD"),phasedFather.getLikelihoods().toString(),phasedChild.getGenotypeString(),phasedChild.getAttribute(VCFConstants.DEPTH_KEY),phasedChild.getAttribute("AD"),phasedChild.getLikelihoods().toString());
}
builder.genotypes(genotypesContext);
vcfWriter.add(builder.make());
//Report violation if set so
//TODO: ADAPT FOR PAIRS TOO!!
if(mvCount>0 && mvFile != null)
mvFile.println(mvfLine);
}
builder.genotypes(genotypesContext);
vcfWriter.add(builder.make());
return metricsCounters;
}