diff --git a/src/main/java/genepi/mut/objects/VariantLine.java b/src/main/java/genepi/mut/objects/VariantLine.java index f93a684..9f4828e 100644 --- a/src/main/java/genepi/mut/objects/VariantLine.java +++ b/src/main/java/genepi/mut/objects/VariantLine.java @@ -49,6 +49,8 @@ public class VariantLine implements Comparable { private double bayesProbability; private double bayesPercentageFWD; private double bayesPercentageREV; + + private double meanBaseQuality; private String insPosition; @@ -186,7 +188,7 @@ public void parseLine(BasePosition base, double level, HashMap f // set this.setCovFWD(totalFWD); this.setCovREV(totalREV); - + if (totalFWD > 0) { aFWDPercents = aFWD / (double) totalFWD; cFWDPercents = cFWD / (double) totalFWD; @@ -555,27 +557,41 @@ public void calcBayes(BasePosition base, HashMap freq) { } } + double meanQualityA = 0; + double meanQualityC = 0; + double meanQualityG = 0; + double meanQualityT = 0; + + double qualityA = 0; + double qualityC = 0; + double qualityG = 0; + double qualityT = 0; + for (int i = 0; i < base.getaFor(); i++) { byte err = base.getaForQ().get(i); double qualScore = Math.pow(10, (-err / 10)); + qualityA += err; probAFor += Math.log10(1 - qualScore); probCFor += Math.log10(qualScore / 3); probGFor += Math.log10(qualScore / 3); probTFor += Math.log10(qualScore / 3); } - + + for (int i = 0; i < base.getcFor(); i++) { byte err = base.getcForQ().get(i); double qualScore = Math.pow(10, (-err / 10)); + qualityC += err; probCFor += Math.log10(1 - qualScore); probAFor += Math.log10(qualScore / 3); probGFor += Math.log10(qualScore / 3); probTFor += Math.log10(qualScore / 3); } - + for (int i = 0; i < base.getgFor(); i++) { byte err = base.getgForQ().get(i); double qualScore = Math.pow(10, (-err / 10)); + qualityG += err; probGFor += Math.log10(1 - qualScore); probAFor += Math.log10(qualScore / 3); probCFor += Math.log10(qualScore / 3); @@ -586,6 +602,7 @@ public void calcBayes(BasePosition base, HashMap freq) { for (int i = 0; i < base.gettFor(); i++) { byte err = base.gettForQ().get(i); double qualScore = Math.pow(10, (-err / 10)); + qualityT += err; probTFor += Math.log10(1 - qualScore); probAFor += Math.log10(qualScore / 3); probCFor += Math.log10(qualScore / 3); @@ -595,6 +612,7 @@ public void calcBayes(BasePosition base, HashMap freq) { for (int i = 0; i < base.getaRev(); i++) { byte err = base.getaRevQ().get(i); double qualScore = Math.pow(10, (-err / 10)); + qualityA += err; probARev += Math.log10(1 - qualScore); probCRev += Math.log10(qualScore / 3); probGRev += Math.log10(qualScore / 3); @@ -604,6 +622,7 @@ public void calcBayes(BasePosition base, HashMap freq) { for (int i = 0; i < base.getcRev(); i++) { byte err = base.getcRevQ().get(i); double qualScore = Math.pow(10, (-err / 10)); + qualityC += err; probCRev += Math.log10((1 - qualScore)); probARev += Math.log10(qualScore / 3); probGRev += Math.log10(qualScore / 3); @@ -613,6 +632,7 @@ public void calcBayes(BasePosition base, HashMap freq) { for (int i = 0; i < base.getgRev(); i++) { byte err = base.getgRevQ().get(i); double qualScore = Math.pow(10, (-err / 10)); + qualityG += err; probGRev += Math.log10((1 - qualScore)); probARev += Math.log10(qualScore / 3); probCRev += Math.log10(qualScore / 3); @@ -622,11 +642,17 @@ public void calcBayes(BasePosition base, HashMap freq) { for (int i = 0; i < base.gettRev(); i++) { byte err = base.gettRevQ().get(i); double qualScore = Math.pow(10, (-err / 10)); + qualityT += err; probTRev += Math.log10((1 - qualScore)); probARev += Math.log10(qualScore / 3); probCRev += Math.log10(qualScore / 3); probGRev += Math.log10(qualScore / 3); } + + meanQualityA = qualityA / (base.getaForQ().size() + base.getaRevQ().size()); + meanQualityC = qualityC / (base.getcForQ().size() + base.getcRevQ().size()); + meanQualityG = qualityG / (base.getgForQ().size() + base.getgRevQ().size()); + meanQualityT = qualityT / (base.gettForQ().size() + base.gettRevQ().size()); // add prior double probA = (probAFor + probARev) + Math.log10(freqA); @@ -657,18 +683,22 @@ public void calcBayes(BasePosition base, HashMap freq) { if (finalBase == 'A') { this.bayesPercentageFWD = aPercentageFWD; this.bayesPercentageREV = aPercentageREV; + this.meanBaseQuality = meanQualityA; } if (finalBase == 'C') { this.bayesPercentageFWD = cPercentageFWD; this.bayesPercentageREV = cPercentageREV; + this.meanBaseQuality = meanQualityC; } if (finalBase == 'G') { this.bayesPercentageFWD = gPercentageFWD; this.bayesPercentageREV = gPercentageREV; + this.meanBaseQuality = meanQualityG; } if (finalBase == 'T') { this.bayesPercentageFWD = tPercentageFWD; this.bayesPercentageREV = tPercentageREV; + this.meanBaseQuality = meanQualityT; } this.setBayesProbability(Math.pow(10, bayesProb - d)); @@ -1400,4 +1430,12 @@ public void setFilter(Filter filter) { this.filter = filter; } + public double getMeanBaseQuality() { + return meanBaseQuality; + } + + public void setMeanBaseQuality(double meanBaseQuality) { + this.meanBaseQuality = meanBaseQuality; + } + } diff --git a/src/main/java/genepi/mut/objects/VariantResult.java b/src/main/java/genepi/mut/objects/VariantResult.java index 2e73e91..ef7a7ad 100644 --- a/src/main/java/genepi/mut/objects/VariantResult.java +++ b/src/main/java/genepi/mut/objects/VariantResult.java @@ -17,6 +17,7 @@ public class VariantResult { private int covREV; private int type; private Filter filter; + private String meanBaseQuality; public int getType() { return type; @@ -121,5 +122,14 @@ public Filter getFilter() { public void setFilter(Filter filter) { this.filter = filter; } + + public void setMeanBaseQuality(String meanBaseQuality) { + this.meanBaseQuality = meanBaseQuality; + + } + + public String getMeanBaseQuality() { + return meanBaseQuality; + } } diff --git a/src/main/java/genepi/mut/pileup/BamAnalyser.java b/src/main/java/genepi/mut/pileup/BamAnalyser.java index 52d907e..163fff9 100644 --- a/src/main/java/genepi/mut/pileup/BamAnalyser.java +++ b/src/main/java/genepi/mut/pileup/BamAnalyser.java @@ -19,7 +19,7 @@ public class BamAnalyser { public final static String headerRaw = "SAMPLE\tPOS\tREF\tTOP-FWD\tMINOR-FWD\tTOP-REV\tMINOR-REV\tCOV-FWD\tCOV-REV\tCOV-TOTAL\tTYPE\tLEVEL\t%A\t%C\t%G\t%T\t%D\t%N\t%a\t%c\t%g\t%t\t%d\t%n\tTOP-FWD-PERCENT\tTOP-REV-PERCENT\tMINOR-FWD-PERCENT\tMINOR-REV-PERCENT\tLLRFWD\tLLRREV\tLLRAFWD\tLLRCFWD\tLLRGFWD\tLLRTFWD\tLLRAREV\tLLRCREV\tLLRGREV\tLLRTREV\tLLRDFWD\tLLRDREV\tMINORS"; - public final static String headerVariants = "ID\tFilter\tPos\tRef\tVariant\tVariantLevel\tMajorBase\tMajorLevel\tMinorBase\tMinorLevel\tCoverage\tCoverageFWD\tCoverageREV\tType"; + public final static String headerVariants = "ID\tFilter\tPos\tRef\tVariant\tVariantLevel\tMajorBase\tMajorLevel\tMinorBase\tMinorLevel\tCoverage\tCoverageFWD\tCoverageREV\tType\tMeanBaseQuality"; HashMap counts; diff --git a/src/main/java/genepi/mut/pileup/VariantCaller.java b/src/main/java/genepi/mut/pileup/VariantCaller.java index 38075cf..dd8dd81 100644 --- a/src/main/java/genepi/mut/pileup/VariantCaller.java +++ b/src/main/java/genepi/mut/pileup/VariantCaller.java @@ -28,6 +28,8 @@ public enum Filter { public static int DELETION = 4; public static int INSERTION = 5; + + static DecimalFormat df = new DecimalFormat("#.##"); public static boolean isFinalVariant(VariantLine line) { @@ -152,6 +154,7 @@ private static VariantResult addVariantResult(VariantLine line, int type, Filter output.setCovFWD(line.getCovFWD()); output.setCovREV(line.getCovREV()); output.setType(type); + output.setMeanBaseQuality(df.format(line.getMeanBaseQuality())); if (filter == null) { filter = Filter.PASS; @@ -189,6 +192,7 @@ private static VariantResult addHomoplasmyResult(VariantLine line, int type) { output.setCovFWD(line.getCovFWD()); output.setCovREV(line.getCovREV()); output.setType(type); + output.setMeanBaseQuality(df.format(line.getMeanBaseQuality())); Filter filter = Filter.PASS; if (blacklist.contains(line.getPosition())) { @@ -450,7 +454,9 @@ public static String writeVariant(VariantResult result) throws IOException { build.append((result.getCovREV()) + "\t"); - build.append(result.getType()); + build.append(result.getType() + "\t"); + + build.append(result.getMeanBaseQuality()); return build.toString();