From e3b08f245fb707a41d9627f590a6f361deeeeb63 Mon Sep 17 00:00:00 2001 From: ebanks Date: Mon, 3 Aug 2009 17:00:20 +0000 Subject: [PATCH] Pull out RMS calculation into MathUtils for all to use git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1364 348d0f76-0448-11de-a6fe-93d51630548a --- .../gatk/walkers/variants/VECAlleleBalance.java | 2 +- .../gatk/walkers/variants/VECMappingQuality.java | 14 ++++++-------- .../walkers/variants/VECOnOffGenotypeRatio.java | 2 +- .../broadinstitute/sting/utils/MathUtils.java | 16 ++++++++++++++++ .../utils/genotype/calls/SSGGenotypeCall.java | 13 ++++++------- 5 files changed, 30 insertions(+), 17 deletions(-) diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/VECAlleleBalance.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/VECAlleleBalance.java index 68fafee61..21b94590a 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/VECAlleleBalance.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/VECAlleleBalance.java @@ -7,7 +7,7 @@ import org.broadinstitute.sting.utils.Utils; import org.broadinstitute.sting.utils.Pair; public class VECAlleleBalance implements VariantExclusionCriterion { //extends RatioFilter { - final private static GenotypeFeatureData.Tail tail = GenotypeFeatureData.Tail.TwoTailed; + //final private static GenotypeFeatureData.Tail tail = GenotypeFeatureData.Tail.TwoTailed; private boolean exclude; private double lowThreshold, highThreshold, ratio; diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/VECMappingQuality.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/VECMappingQuality.java index a016a779b..021adaf7f 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/VECMappingQuality.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/VECMappingQuality.java @@ -2,6 +2,7 @@ package org.broadinstitute.sting.playground.gatk.walkers.variants; import org.broadinstitute.sting.gatk.LocusContext; import org.broadinstitute.sting.gatk.refdata.rodVariants; +import org.broadinstitute.sting.utils.MathUtils; import net.sf.samtools.SAMRecord; import java.util.List; @@ -20,14 +21,11 @@ public class VECMappingQuality implements VariantExclusionCriterion { public void compute(char ref, LocusContext context, rodVariants variant) { List reads = context.getReads(); - - rms = 0.0; - for (int readIndex = 0; readIndex < reads.size(); readIndex++) { - int qual = reads.get(readIndex).getMappingQuality(); - rms += qual * qual; - } - rms /= reads.size(); - rms = Math.sqrt(rms); + int[] qualities = new int[reads.size()]; + for (int i=0; i < reads.size(); i++) + qualities[i] = reads.get(i).getMappingQuality(); + rms = MathUtils.rms(qualities); + exclude = rms < minQuality; } public boolean isExcludable() { diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/VECOnOffGenotypeRatio.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/VECOnOffGenotypeRatio.java index 21d42ca8a..8af1987f3 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/VECOnOffGenotypeRatio.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/VECOnOffGenotypeRatio.java @@ -5,7 +5,7 @@ import org.broadinstitute.sting.gatk.LocusContext; import org.broadinstitute.sting.utils.*; public class VECOnOffGenotypeRatio implements VariantExclusionCriterion { // extends RatioFilter { - final private static GenotypeFeatureData.Tail tail = GenotypeFeatureData.Tail.LeftTailed; + //final private static GenotypeFeatureData.Tail tail = GenotypeFeatureData.Tail.LeftTailed; private boolean exclude; private double lowThreshold, highThreshold, ratio; diff --git a/java/src/org/broadinstitute/sting/utils/MathUtils.java b/java/src/org/broadinstitute/sting/utils/MathUtils.java index 3a402c882..b93e8f0ae 100755 --- a/java/src/org/broadinstitute/sting/utils/MathUtils.java +++ b/java/src/org/broadinstitute/sting/utils/MathUtils.java @@ -142,4 +142,20 @@ public class MathUtils { return multinomialCoefficient*probs; } + + /** + * calculate the Root Mean Square of an array of integers + * @param x an int[] of numbers + * @return the RMS of the specified numbers. + */ + public static double rms(int[] x) { + if ( x.length == 0 ) + return 0.0; + + double rms = 0.0; + for (int i : x) + rms += i * i; + rms /= x.length; + return Math.sqrt(rms); + } } diff --git a/java/src/org/broadinstitute/sting/utils/genotype/calls/SSGGenotypeCall.java b/java/src/org/broadinstitute/sting/utils/genotype/calls/SSGGenotypeCall.java index 0760868cc..e5afd9d91 100644 --- a/java/src/org/broadinstitute/sting/utils/genotype/calls/SSGGenotypeCall.java +++ b/java/src/org/broadinstitute/sting/utils/genotype/calls/SSGGenotypeCall.java @@ -4,6 +4,7 @@ import net.sf.samtools.SAMRecord; import org.broadinstitute.sting.utils.GenomeLoc; import org.broadinstitute.sting.utils.ReadBackedPileup; import org.broadinstitute.sting.utils.Utils; +import org.broadinstitute.sting.utils.MathUtils; import org.broadinstitute.sting.utils.genotype.BasicGenotype; import org.broadinstitute.sting.utils.genotype.Genotype; import org.broadinstitute.sting.utils.genotype.GenotypeOutput; @@ -103,13 +104,11 @@ public class SSGGenotypeCall implements GenotypeCall, GenotypeOutput { * @return */ private double calculateRMS(ReadBackedPileup pileup) { - double rms = 0.0; - for (SAMRecord r : pileup.getReads()) { - rms += r.getMappingQuality() * r.getMappingQuality(); - } - rms /= pileup.getReads().size(); - rms = Math.sqrt(rms); - return rms; + List reads = pileup.getReads(); + int[] qualities = new int[reads.size()]; + for (int i=0; i < reads.size(); i++) + qualities[i] = reads.get(i).getMappingQuality(); + return MathUtils.rms(qualities); } /**