PercentOfBasesCovered functionality moved to DepthOfCoverageWalker. Added integration tests.

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@2244 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
ebanks 2009-12-03 16:11:09 +00:00
parent 126d1eca35
commit 44b9f60735
3 changed files with 22 additions and 67 deletions

View File

@ -51,6 +51,9 @@ public class DepthOfCoverageWalker extends LocusWalker<DepthOfCoverageWalker.DoC
@Argument(fullName="minMAPQ", shortName = "minMAPQ", doc="If provided, we will also list read counts with MAPQ >= this value at a locus in coverage",required=false)
protected int excludeMAPQBelowThis = -1;
@Argument(fullName="minDepth", shortName = "minDepth", doc="If provided, we will also list the percentage of loci with depth >= this value per interval",required=false)
protected int minDepthForPercentage = -1;
@Argument(fullName="byReadGroup", shortName="byRG", doc="List read depths for each read group")
protected boolean byReadGroup = false;
@ -153,6 +156,9 @@ public class DepthOfCoverageWalker extends LocusWalker<DepthOfCoverageWalker.DoC
sum.totalCoverage += value.totalCoverage;
sum.numDeletions += value.numDeletions;
sum.numBadMQReads += value.numBadMQReads;
if ( value.totalCoverage >= minDepthForPercentage ) {
sum.minDepthCoveredLoci++;
}
if ( byReadGroup ) {
for ( String rg : readGroupNames ) {
int oldDepth = sum.depthByReadGroup.get(rg);
@ -179,6 +185,10 @@ public class DepthOfCoverageWalker extends LocusWalker<DepthOfCoverageWalker.DoC
header.append("\taverage_coverage_atleast_MQ");
header.append(excludeMAPQBelowThis);
}
if ( minDepthForPercentage >= 0 ) {
header.append("\tpercent_loci_covered_atleast_depth");
header.append(minDepthForPercentage);
}
if ( byReadGroup ) {
for ( String rg : readGroupNames ) {
header.append("\tcoverage_for_");
@ -199,28 +209,34 @@ public class DepthOfCoverageWalker extends LocusWalker<DepthOfCoverageWalker.DoC
private void printDoCInfo(GenomeLoc loc, DoCInfo info, boolean printAverageCoverage) {
double totalBases = (double)(loc.getStop() - loc.getStart() + 1);
StringBuilder sb = new StringBuilder();
sb.append(loc);
sb.append("\t");
sb.append(info.totalCoverage);
sb.append("\t");
if ( printAverageCoverage ) {
sb.append(String.format("%.2f", ((double)info.totalCoverage) / (double)(loc.getStop() - loc.getStart() + 1)));
sb.append(String.format("%.2f", ((double)info.totalCoverage) / totalBases));
sb.append("\t");
}
sb.append((info.totalCoverage - info.numDeletions));
if ( printAverageCoverage ) {
sb.append("\t");
sb.append(String.format("%.2f", ((double)(info.totalCoverage - info.numDeletions)) / (double)(loc.getStop() - loc.getStart() + 1)));
sb.append(String.format("%.2f", ((double)(info.totalCoverage - info.numDeletions)) / totalBases));
}
if ( excludeMAPQBelowThis > 0 ) {
sb.append("\t");
sb.append((info.totalCoverage - info.numBadMQReads));
if ( printAverageCoverage ) {
sb.append("\t");
sb.append(String.format("%.2f", ((double)(info.totalCoverage - info.numBadMQReads)) / (double)(loc.getStop() - loc.getStart() + 1)));
sb.append(String.format("%.2f", ((double)(info.totalCoverage - info.numBadMQReads)) / totalBases));
}
}
if ( minDepthForPercentage >= 0 ) {
sb.append("\t");
sb.append(String.format("%.2f", ((double)info.minDepthCoveredLoci) / totalBases));
}
if ( byReadGroup ) {
for ( String rg : readGroupNames ) {
sb.append("\t");
@ -242,6 +258,7 @@ public class DepthOfCoverageWalker extends LocusWalker<DepthOfCoverageWalker.DoC
public int totalCoverage = 0;
public int numDeletions = 0;
public int numBadMQReads = 0;
public int minDepthCoveredLoci = 0;
public HashMap<String, Integer> depthByReadGroup = null;
public HashMap<String, Integer> depthBySample = null;

View File

@ -1,64 +0,0 @@
package org.broadinstitute.sting.playground.gatk.walkers;
import org.broadinstitute.sting.gatk.walkers.LocusWalker;
import org.broadinstitute.sting.gatk.walkers.TreeReducible;
import org.broadinstitute.sting.gatk.walkers.DataSource;
import org.broadinstitute.sting.gatk.walkers.By;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.utils.Pair;
import org.broadinstitute.sting.utils.cmdLine.Argument;
/**
* Created by IntelliJ IDEA.
* User: Ghost
* Date: Nov 25, 2009
* Time: 7:54:42 PM
* To change this template use File | Settings | File Templates.
*/
@By(DataSource.REFERENCE)
public class PercentOfBasesCoveredWalker extends LocusWalker<Boolean, Pair<Integer,Integer>> implements TreeReducible<Pair<Integer,Integer>> {
@Argument(fullName="minimumDepth", shortName="minDepth", doc="Minimum depth beyond which to discount the base as uncovered", required=false)
int minDepth = 0;
@Argument(fullName="includeName",doc="include this name in the output of this walker (e.g. NAME: percent of bases covered . . .) ", required=false)
String name = null;
public void initialize() {
// do nothing
}
public Pair<Integer,Integer> reduceInit() {
return new Pair<Integer,Integer>(0,1); // robust initialization -- seen one base with zero coverage
}
public Boolean map(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
return context.size() > minDepth;
}
public Pair<Integer,Integer> reduce(Boolean map, Pair<Integer,Integer> prevReduce) {
if ( map.booleanValue() ) {
return new Pair<Integer,Integer>(prevReduce.first+1,prevReduce.second+1);
} else {
return new Pair<Integer,Integer>(prevReduce.first,prevReduce.second+1);
}
}
public Pair<Integer,Integer> treeReduce( Pair<Integer,Integer> left, Pair<Integer,Integer> right) {
right.set(right.first+left.first,right.second+left.second);
return right;
}
public void onTraversalDone(Pair<Integer,Integer> reduceFinal) {
String msg;
if (name == null)
msg = "Percent of bases covered at "+minDepth+"x=";
else
msg = name+": Percent of bases covered at "+minDepth+"x=";
out.printf("%s%s",msg,pair2percent(reduceFinal));
}
private String pair2percent(Pair<Integer,Integer> p ) {
return String.format("%d%s", (int) Math.floor( ( ( double) p.first*100.0) / p.second ), "%" );
}
}

View File

@ -15,6 +15,8 @@ public class DepthOfCoverageIntegrationTest extends WalkerTest {
static {
expectations.put("-minMAPQ 1", "59c6071105a598e19f460640c35768c6");
expectations.put("-minMAPQ 100", "e997fb5d61eaec21518722b0de90af20");
expectations.put("-minDepth 8", "3e50afef0e751119cd27c324bdfae544");
expectations.put("-minDepth 10", "d4c336d9e748347e1082bbc92d2489a3");
expectations.put("-bySample", "160ffa185dbfa8b0d2dc57f60f5b1e48");
expectations.put("-byRG", "dd3b4d040df7325dad4760ac6fa5252d");
expectations.put("-minMAPQ 1 -bySample -byRG", "bd2a07ef548b86e82ac6cce534225612");