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

When merging MNPs, the phased flag and the phase quality (PQ) are determined simultaneously 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@4613 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
fromer 2010-11-02 14:44:26 +00:00
parent d496f2afde
commit a885ecf046
2 changed files with 24 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
java/src/org/broadinstitute/sting/gatk/contexts/variantcontext/VariantContextUtils.java
View File

@ -789,6 +789,9 @@ public class VariantContextUtils {
List<Allele> mergedAllelesForSample = new LinkedList<Allele>(); List<Allele> mergedAllelesForSample = new LinkedList<Allele>();
/* NOTE: Since merged alleles are added to mergedAllelesForSample in the SAME order as in the input VC records,
we preserve phase information (if any) relative to whatever precedes vc1:
*/
Iterator<Allele> all2It = site2Alleles.iterator(); Iterator<Allele> all2It = site2Alleles.iterator();
for (Allele all1 : site1Alleles) { for (Allele all1 : site1Alleles) {
Allele all2 = all2It.next(); // this is OK, since allSamplesAreMergeable() Allele all2 = all2It.next(); // this is OK, since allSamplesAreMergeable()
@ -801,11 +804,11 @@ public class VariantContextUtils {
Set<String> mergedGtFilters = new HashSet<String>(); // Since gt1 and gt2 were unfiltered, the Genotype remains unfiltered Set<String> mergedGtFilters = new HashSet<String>(); // Since gt1 and gt2 were unfiltered, the Genotype remains unfiltered
Map<String, Object> mergedGtAttribs = new HashMap<String, Object>(); Map<String, Object> mergedGtAttribs = new HashMap<String, Object>();
Object PQ = gt1.getAttribute(ReadBackedPhasingWalker.PQ_KEY); PhaseAndQuality phaseQual = calcPhaseForMergedGenotypes(gt1, gt2);
if (PQ != null) if (phaseQual.PQ != null)
mergedGtAttribs.put(ReadBackedPhasingWalker.PQ_KEY, PQ); mergedGtAttribs.put(ReadBackedPhasingWalker.PQ_KEY, phaseQual.PQ);
Genotype mergedGt = new Genotype(sample, mergedAllelesForSample, mergedGQ, mergedGtFilters, mergedGtAttribs, calcPhaseForMergedGenotypes(gt1, gt2)); Genotype mergedGt = new Genotype(sample, mergedAllelesForSample, mergedGQ, mergedGtFilters, mergedGtAttribs, phaseQual.isPhased);
mergedGenotypes.put(sample, mergedGt); mergedGenotypes.put(sample, mergedGt);
} }
@ -1029,10 +1032,21 @@ public class VariantContextUtils {
return (gt2.genotypesArePhased() || gt2.isHom() || gt1.isHom()); return (gt2.genotypesArePhased() || gt2.isHom() || gt1.isHom());
} }
private static class PhaseAndQuality {
public boolean isPhased;
public Double PQ = null;
public PhaseAndQuality(Genotype gt) {
this.isPhased = gt.genotypesArePhased();
if (this.isPhased)
this.PQ = VariantContextUtils.getDoubleAttribute(gt.getAttributes(), ReadBackedPhasingWalker.PQ_KEY);
}
}
// Assumes that alleleSegregationIsKnown(gt1, gt2): // Assumes that alleleSegregationIsKnown(gt1, gt2):
private static boolean calcPhaseForMergedGenotypes(Genotype gt1, Genotype gt2) { private static PhaseAndQuality calcPhaseForMergedGenotypes(Genotype gt1, Genotype gt2) {
if (gt2.genotypesArePhased() || gt2.isHom()) if (gt2.genotypesArePhased() || gt2.isHom())
return gt1.genotypesArePhased(); // maintain the phase of gt1 return new PhaseAndQuality(gt1); // maintain the phase of gt1
if (!gt1.isHom()) if (!gt1.isHom())
throw new ReviewedStingException("alleleSegregationIsKnown(gt1, gt2) implies: gt2.genotypesArePhased() || gt2.isHom() || gt1.isHom()"); throw new ReviewedStingException("alleleSegregationIsKnown(gt1, gt2) implies: gt2.genotypesArePhased() || gt2.isHom() || gt1.isHom()");
@ -1048,7 +1062,7 @@ public class VariantContextUtils {
0/1 0/1
1/2 1/2
*/ */
return false; // maintain the phase of gt2 [since !gt2.genotypesArePhased()] return new PhaseAndQuality(gt2); // maintain the phase of gt2 [since !gt2.genotypesArePhased()]
} }
/* Checks if any sample has a MNP of ALT alleles (segregating together): /* Checks if any sample has a MNP of ALT alleles (segregating together):

4
java/src/org/broadinstitute/sting/gatk/walkers/phasing/MergeSegregatingAlternateAllelesWalker.java
View File

@ -28,6 +28,7 @@ import org.broad.tribble.util.variantcontext.Allele;
import org.broad.tribble.util.variantcontext.VariantContext; import org.broad.tribble.util.variantcontext.VariantContext;
import org.broad.tribble.vcf.*; import org.broad.tribble.vcf.*;
import org.broadinstitute.sting.commandline.Argument; import org.broadinstitute.sting.commandline.Argument;
import org.broadinstitute.sting.commandline.Hidden;
import org.broadinstitute.sting.commandline.Output; import org.broadinstitute.sting.commandline.Output;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext; import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext; import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
@ -59,7 +60,8 @@ public class MergeSegregatingAlternateAllelesWalker extends RodWalker<Integer, I
@Argument(fullName = "useSingleSample", shortName = "useSample", doc = "Only output genotypes for the single sample given; [default:use all samples]", required = false) @Argument(fullName = "useSingleSample", shortName = "useSample", doc = "Only output genotypes for the single sample given; [default:use all samples]", required = false)
protected String useSingleSample = null; protected String useSingleSample = null;
@Argument(fullName = "emitOnlyMergedRecords", shortName = "emitOnlyMerged", doc = "Only output records that resulted from merging; [default:false]", required = false) @Hidden
@Argument(fullName = "emitOnlyMergedRecords", shortName = "emitOnlyMerged", doc = "Only output records that resulted from merging [For DEBUGGING purposes only - DO NOT USE, since it disregards the semantics of '|' as 'phased relative to previous non-filtered VC']; [default:false]", required = false)
protected boolean emitOnlyMergedRecords = false; protected boolean emitOnlyMergedRecords = false;
@Argument(fullName = "disablePrintAltAlleleStats", shortName = "noAlleleStats", doc = "Should the print-out of alternate allele statistics be disabled?; [default:false]", required = false) @Argument(fullName = "disablePrintAltAlleleStats", shortName = "noAlleleStats", doc = "Should the print-out of alternate allele statistics be disabled?; [default:false]", required = false)