Skip to content

Commit

Permalink
Write mean base quality to variants file
Browse files Browse the repository at this point in the history
  • Loading branch information
@seppinho
seppinho committed Mar 2, 2024
1 parent 125f4b8 commit e08a03b
Show file tree
Hide file tree
Showing 4 changed files with 59 additions and 5 deletions.
44 changes: 41 additions & 3 deletions src/main/java/genepi/mut/objects/VariantLine.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,8 @@ public class VariantLine implements Comparable<VariantLine> {
private double bayesProbability;
private double bayesPercentageFWD;
private double bayesPercentageREV;

private double meanBaseQuality;

private String insPosition;

Expand Down Expand Up @@ -186,7 +188,7 @@ public void parseLine(BasePosition base, double level, HashMap<String, Double> f
// set
this.setCovFWD(totalFWD);
this.setCovREV(totalREV);

if (totalFWD > 0) {
aFWDPercents = aFWD / (double) totalFWD;
cFWDPercents = cFWD / (double) totalFWD;
Expand Down Expand Up @@ -555,27 +557,41 @@ public void calcBayes(BasePosition base, HashMap<String, Double> 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);
Expand All @@ -586,6 +602,7 @@ public void calcBayes(BasePosition base, HashMap<String, Double> 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);
Expand All @@ -595,6 +612,7 @@ public void calcBayes(BasePosition base, HashMap<String, Double> 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);
Expand All @@ -604,6 +622,7 @@ public void calcBayes(BasePosition base, HashMap<String, Double> 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);
Expand All @@ -613,6 +632,7 @@ public void calcBayes(BasePosition base, HashMap<String, Double> 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);
Expand All @@ -622,11 +642,17 @@ public void calcBayes(BasePosition base, HashMap<String, Double> 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);
Expand Down Expand Up @@ -657,18 +683,22 @@ public void calcBayes(BasePosition base, HashMap<String, Double> 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));
Expand Down Expand Up @@ -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;
}

}
10 changes: 10 additions & 0 deletions src/main/java/genepi/mut/objects/VariantResult.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ public class VariantResult {
private int covREV;
private int type;
private Filter filter;
private String meanBaseQuality;

public int getType() {
return type;
Expand Down Expand Up @@ -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;
}

}
2 changes: 1 addition & 1 deletion src/main/java/genepi/mut/pileup/BamAnalyser.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Integer, BasePosition> counts;

Expand Down
8 changes: 7 additions & 1 deletion src/main/java/genepi/mut/pileup/VariantCaller.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {

Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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())) {
Expand Down Expand Up @@ -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();

Expand Down

0 comments on commit e08a03b

Please sign in to comment.