Skip to content

Commit

Permalink
Adding odds ratio filter
Browse files Browse the repository at this point in the history
  • Loading branch information
@lmose
lmose committed Aug 13, 2021
1 parent fe7ea29 commit b8e954f
Show file tree
Hide file tree
Showing 2 changed files with 27 additions and 10 deletions.
10 changes: 8 additions & 2 deletions src/main/java/abra/cadabra/CadabraOptions.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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() {
Expand All @@ -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;
Expand All @@ -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
Expand Down Expand Up @@ -162,4 +166,6 @@ public float getMinMapq() {
public float getMinVaf() {
return minVaf;
}

public float getOddsRatio() { return oddsr; }
}
27 changes: 19 additions & 8 deletions src/main/java/abra/cadabra/CadabraProcessor.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -329,6 +329,7 @@ private SampleCall processLocus(ReadsAtLocus reads, boolean isSomatic) {
Pair<Integer, String> 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());
Expand All @@ -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
Expand Down Expand Up @@ -400,6 +402,7 @@ public static class SampleCall {
int totalReads;
int usableDepth;
double qual;
double oddsRatio;
int repeatPeriod;
String repeatUnit;
int mapq0;
Expand All @@ -413,7 +416,7 @@ public static class SampleCall {
CadabraOptions options;

SampleCall(String chromosome, int position, Allele ref, Allele alt, Map<Allele, AlleleCounts> 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;
Expand All @@ -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;
Expand Down Expand Up @@ -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
Expand All @@ -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";
Expand Down Expand Up @@ -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;
Expand All @@ -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);

Expand All @@ -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);
}
Expand Down

0 comments on commit b8e954f

Please sign in to comment.