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

Fixed precision issues with PQ (phasing quality) 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@4068 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
fromer 2010-08-19 20:34:47 +00:00
parent 88ca1fb22c
commit 3af4e618cc
1 changed files with 10 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
java/src/org/broadinstitute/sting/playground/gatk/walkers/ReadBackedPhasingWalker.java
View File

@ -87,7 +87,7 @@ public class ReadBackedPhasingWalker extends LocusWalker<PhasingStatsAndOutput,
Set<VCFHeaderLine> hInfo = new HashSet<VCFHeaderLine>();
hInfo.addAll(VCFUtils.getHeaderFields(getToolkit()));
hInfo.add(new VCFHeaderLine("reference", getToolkit().getArguments().referenceFile.getName()));
hInfo.add(new VCFFormatHeaderLine("PQ", 1, VCFHeaderLineType.Integer, "Read-backed phasing quality"));
hInfo.add(new VCFFormatHeaderLine("PQ", 1, VCFHeaderLineType.Float, "Read-backed phasing quality"));
writer = new VCFWriterImpl(new File(phasedVCFFile));
writer.writeHeader(new VCFHeader(hInfo, new TreeSet<String>(vc.getSampleNames())));
@ -282,8 +282,15 @@ public class ReadBackedPhasingWalker extends LocusWalker<PhasingStatsAndOutput,
logger.debug("\nPhasing table [AFTER NORMALIZATION]:\n" + sampleHaps + "\n");
PhasingTable.PhasingTableEntry maxEntry = sampleHaps.maxEntry();
PreciseNonNegativeDouble sumErrorProbs = new PreciseNonNegativeDouble(ZERO);
for (PhasingTable.PhasingTableEntry pte : sampleHaps) {
if (pte != maxEntry)
sumErrorProbs.plusEqual(pte.getScore());
}
// log10(x) = log(x) / log(10):
double phaseQuality = -10.0 * (sumErrorProbs.getLogValue() / Math.log(10.0)); // convert to PHRED scale, do NOT cap the quality as in QualityUtils.probToQual(posteriorProb)
double posteriorProb = maxEntry.getScore().getValue();
int phaseQuality = new Integer(QualityUtils.probToQual(posteriorProb, 0.0)); // 0.0 <=> do NOT cap the quality!
logger.debug("MAX hap:\t" + maxEntry.getHaplotypeClass() + "\tposteriorProb:\t" + posteriorProb + "\tphaseQuality:\t" + phaseQuality);
if (statsWriter != null)
@ -1111,7 +1118,7 @@ class PhasingQualityStatsWriter {
this.variantStatsFilePrefix = variantStatsFilePrefix;
}
public void addStat(String sample, int distanceFromPrevious, int phasingQuality) {
public void addStat(String sample, int distanceFromPrevious, double phasingQuality) {
BufferedWriter sampWriter = sampleToStatsWriter.get(sample);
if (sampWriter == null) {
String fileName = variantStatsFilePrefix + "." + sample + ".distance_PQ.txt";