From b8e954f83f68cb971fba4a833c2be0d85c77f38c Mon Sep 17 00:00:00 2001 From: mozack Date: Fri, 13 Aug 2021 13:07:46 -0400 Subject: [PATCH] Adding odds ratio filter --- .../java/abra/cadabra/CadabraOptions.java | 10 +++++-- .../java/abra/cadabra/CadabraProcessor.java | 27 +++++++++++++------ 2 files changed, 27 insertions(+), 10 deletions(-) diff --git a/src/main/java/abra/cadabra/CadabraOptions.java b/src/main/java/abra/cadabra/CadabraOptions.java index 6ac4896..a089771 100644 --- a/src/main/java/abra/cadabra/CadabraOptions.java +++ b/src/main/java/abra/cadabra/CadabraOptions.java @@ -23,7 +23,8 @@ public class CadabraOptions extends Options { private static final String MIN_MAPQ = "mapq"; private static final String MIN_VAF = "mf"; private static final String PCR_PENALTY = "pen"; - + private static final String ODDS_RATIO = "oddsr"; + private OptionParser parser; private boolean isValid; @@ -42,6 +43,7 @@ public class CadabraOptions extends Options { private float minQual; private float minMapq; private float minVaf; + private float oddsr; @Override protected OptionParser getOptionParser() { @@ -61,7 +63,8 @@ protected OptionParser getOptionParser() { parser.accepts(LOW_MQ_FILTER, "Filter variants with fraction of low quality reads greater than specified value").withRequiredArg().ofType(Float.class).defaultsTo(.5f); parser.accepts(MIN_QUAL, "Variants with quality below specified threshold are not output").withRequiredArg().ofType(Float.class).defaultsTo(5.0f); parser.accepts(MIN_MAPQ, "Reads with mapping quality below specified value are excluded from processing (except unmapped reads)").withRequiredArg().ofType(Integer.class).defaultsTo(20); - parser.accepts(MIN_VAF, "Do not output variants with frequency below specified value").withRequiredArg().ofType(Float.class).defaultsTo(0.0f); + parser.accepts(MIN_VAF, "Do not output variants with frequency below specified value").withRequiredArg().ofType(Float.class).defaultsTo(0.0f); + parser.accepts(ODDS_RATIO, "Filter variants with odds ratio below specified value").withRequiredArg().ofType(Float.class).defaultsTo(5.0f); } return parser; @@ -85,6 +88,7 @@ public void init() { this.minMapq = (Integer) getOptions().valueOf(MIN_MAPQ); this.minVaf = (Float) getOptions().valueOf(MIN_VAF); this.pcrPenalty = (Integer) getOptions().valueOf(PCR_PENALTY); + this.oddsr = (Float) getOptions().valueOf(ODDS_RATIO); } @Override @@ -162,4 +166,6 @@ public float getMinMapq() { public float getMinVaf() { return minVaf; } + + public float getOddsRatio() { return oddsr; } } diff --git a/src/main/java/abra/cadabra/CadabraProcessor.java b/src/main/java/abra/cadabra/CadabraProcessor.java index e1808ed..4df3e61 100644 --- a/src/main/java/abra/cadabra/CadabraProcessor.java +++ b/src/main/java/abra/cadabra/CadabraProcessor.java @@ -329,6 +329,7 @@ private SampleCall processLocus(ReadsAtLocus reads, boolean isSomatic) { Pair repeat = getRepeatPeriod(chromosome, position, alt, altCounts); double qual = 0; + double oddsRatio = 0; int usableDepth = 0; if (!isSomatic) { int repeatLength = getRepeatLength(repeat.getFirst(), repeat.getSecond(), alt.getType()); @@ -348,16 +349,17 @@ private SampleCall processLocus(ReadsAtLocus reads, boolean isSomatic) { } call = new SampleCall(chromosome, position, refAllele, alt, alleleCounts, totalDepth, - usableDepth, qual, repeat.getFirst(), repeat.getSecond(), tumorMapq0, refField, altField, mismatchExceededReads, refSeq, options); + usableDepth, qual, oddsRatio, repeat.getFirst(), repeat.getSecond(), tumorMapq0, refField, altField, mismatchExceededReads, refSeq, options); } else { String refField = getInsRefField(chromosome, position); String altField = "."; double qual = 0; + double oddsRatio = 0; int rp = 0; String ru = ""; call = new SampleCall(chromosome, position, refAllele, Allele.UNK, alleleCounts, totalDepth, - 0, qual, rp, ru, tumorMapq0, refField, altField, mismatchExceededReads, refSeq, options); + 0, qual, oddsRatio, rp, ru, tumorMapq0, refField, altField, mismatchExceededReads, refSeq, options); if (!isSomatic) { // Adjust qual score for PCR slippage @@ -400,6 +402,7 @@ public static class SampleCall { int totalReads; int usableDepth; double qual; + double oddsRatio; int repeatPeriod; String repeatUnit; int mapq0; @@ -413,7 +416,7 @@ public static class SampleCall { CadabraOptions options; SampleCall(String chromosome, int position, Allele ref, Allele alt, Map alleleCounts, - int totalReads, int usableDepth, double qual, int repeatPeriod, String repeatUnit, int mapq0, String refField, String altField, + int totalReads, int usableDepth, double qual, double oddsRatio, int repeatPeriod, String repeatUnit, int mapq0, String refField, String altField, int mismatchExceededReads, String context, CadabraOptions options) { this.chromosome = chromosome; this.position = position; @@ -423,6 +426,7 @@ public static class SampleCall { this.totalReads = totalReads; this.usableDepth = usableDepth; this.qual = qual; + this.oddsRatio = oddsRatio; this.repeatPeriod = repeatPeriod; this.repeatUnit = repeatUnit; this.mapq0 = mapq0; @@ -500,13 +504,13 @@ public String toString() { String sampleInfo = getSampleInfo(ref, alt); - String filter = CadabraProcessor.applyFilters(this, null, options, hrunLen, qual); + String filter = CadabraProcessor.applyFilters(this, null, options, hrunLen, qual, oddsRatio); return String.join("\t", chromosome, pos, ".", refField, altField, qualStr, filter, info, SampleCall.FORMAT, sampleInfo); } } - static String applyFilters(SampleCall tumor, SampleCall normal, CadabraOptions options, int hrunLen, double qual) { + static String applyFilters(SampleCall tumor, SampleCall normal, CadabraOptions options, int hrunLen, double qual, double oddsRatio) { String filter = ""; // Filter variants that do not appear in sufficiently varying read positions @@ -528,6 +532,10 @@ static String applyFilters(SampleCall tumor, SampleCall normal, CadabraOptions o if (qual < options.getQualFilter()) { filter += "LOW_QUAL;"; } + + if (oddsRatio < options.getOddsRatio()) { + filter += "ODDS_R"; + } if (filter.equals("")) { filter = "PASS"; @@ -571,6 +579,7 @@ public static class SomaticCall { HomopolymerRun hrun; String context; CadabraOptions options; + double oddsRatio; public SomaticCall(SampleCall normal, SampleCall tumor, String context, CadabraOptions options) { this.normal = normal; @@ -583,6 +592,8 @@ public SomaticCall(SampleCall normal, SampleCall tumor, String context, CadabraO int tumorAlt = tumor.alleleCounts.get(tumor.alt).getCount(); this.qual = calcFisherExactPhredScaledQuality(normalRef, normalAlt, tumorRef, tumorAlt); + + this.oddsRatio = ((float) tumorAlt / ((float) normalAlt + .0000001)) / ((float) tumorRef / ((float) normalRef + .0000001)); this.hrun = HomopolymerRun.find(context); @@ -607,13 +618,13 @@ public String toString() { char hrunBase = hrun != null ? hrun.getBase() : 'N'; int hrunPos = hrun != null ? hrun.getPos() : 0; - String info = String.format("RP=%d;RU=%s;HRUN=%d,%d;CTX=%s", tumor.repeatPeriod, tumor.repeatUnit, - hrunLen, hrunPos, context); + String info = String.format("RP=%d;RU=%s;HRUN=%d,%d;CTX=%s,ODDSR=%.2f", tumor.repeatPeriod, tumor.repeatUnit, + hrunLen, hrunPos, context, oddsRatio); String normalInfo = normal.getSampleInfo(tumor.ref, tumor.alt); String tumorInfo = tumor.getSampleInfo(tumor.ref, tumor.alt); - String filter = CadabraProcessor.applyFilters(tumor, normal, options, hrunLen, qual); + String filter = CadabraProcessor.applyFilters(tumor, normal, options, hrunLen, qual, oddsRatio); return String.join("\t", tumor.chromosome, pos, ".", tumor.refField, tumor.altField, qualStr, filter, info, SampleCall.FORMAT, normalInfo, tumorInfo); }