diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/diagnostics/FindCoveredIntervals.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/diagnostics/FindCoveredIntervals.java index ad6023579..bd69cbdbd 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/diagnostics/FindCoveredIntervals.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/diagnostics/FindCoveredIntervals.java @@ -101,11 +101,24 @@ public class FindCoveredIntervals extends ActiveRegionWalker { @Argument(fullName = "coverage_threshold", shortName = "cov", doc = "The minimum allowable coverage to be considered covered", required = false) private int coverageThreshold = 20; + @Argument(fullName = "minBaseQuality", shortName = "minBQ", doc = "The minimum allowable base quality score to be counted for coverage",required = false) + private int minBaseQuality = 0; + + @Argument(fullName = "minMappingQuality", shortName = "minMQ", doc = "The minimum allowable mapping quality score to be counted for coverage",required = false) + private int minMappingQuality = 0; + + + + @Override // Look to see if the region has sufficient coverage public ActivityProfileState isActive(final RefMetaDataTracker tracker, final ReferenceContext ref, final AlignmentContext context) { - int depth = context.getBasePileup().getBaseFilteredPileup(coverageThreshold).depthOfCoverage(); + int depth; + if(minBaseQuality == 0 && minMappingQuality == 0) + depth = context.getBasePileup().getBaseFilteredPileup(coverageThreshold).depthOfCoverage(); + else + depth = context.getBasePileup().getBaseAndMappingFilteredPileup(minBaseQuality,minMappingQuality).depthOfCoverage(); // note the linear probability scale return new ActivityProfileState(ref.getLocus(), Math.min(depth / coverageThreshold, 1)); diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/diagnostics/ReadLengthDistribution.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/diagnostics/ReadLengthDistribution.java index a269a94bc..796c817ff 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/diagnostics/ReadLengthDistribution.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/diagnostics/ReadLengthDistribution.java @@ -38,7 +38,10 @@ import org.broadinstitute.sting.utils.help.HelpConstants; import org.broadinstitute.sting.utils.sam.GATKSAMRecord; import java.io.PrintStream; +import java.util.HashMap; import java.util.List; +import java.util.Map; +import java.util.TreeMap; /** * Outputs the read lengths of all the reads in a file. @@ -77,51 +80,101 @@ public class ReadLengthDistribution extends ReadWalker { @Output public PrintStream out; - private GATKReport report; + //A map from RG to its column number (its index in an int[] array) + private Map readGroupsLocation; + //Each line in the table is a read length and each column it the number of reads of a specific RG with that length. Thus a table is a map between read lengths to array of values (one for each RG). + private Map table; + private List readGroups; public void initialize() { - final List readGroups = getToolkit().getSAMFileHeader().getReadGroups(); + readGroups = getToolkit().getSAMFileHeader().getReadGroups(); + readGroupsLocation = new HashMap<>(); + table = new TreeMap<>(); + int readGroupsNum = 0; - report = new GATKReport(); - report.addTable("ReadLengthDistribution", "Table of read length distributions", 1 + (readGroups.isEmpty() ? 1 : readGroups.size())); - GATKReportTable table = report.getTable("ReadLengthDistribution"); - - table.addColumn("readLength"); - - if (readGroups.isEmpty()) - table.addColumn("SINGLE_SAMPLE"); - else - for (SAMReadGroupRecord rg : readGroups) - table.addColumn(rg.getSample()); - } - - public boolean filter(ReferenceContext ref, GATKSAMRecord read) { - return ( !read.getReadPairedFlag() || read.getReadPairedFlag() && read.getFirstOfPairFlag()); + if (!readGroups.isEmpty()){ + for (SAMReadGroupRecord rg : readGroups){ + readGroupsLocation.put(rg,readGroupsNum); + readGroupsNum++; + } + } } @Override - public Integer map(ReferenceContext referenceContext, GATKSAMRecord samRecord, RefMetaDataTracker RefMetaDataTracker) { - GATKReportTable table = report.getTable("ReadLengthDistribution"); + public Integer map(final ReferenceContext referenceContext,final GATKSAMRecord samRecord,final RefMetaDataTracker RefMetaDataTracker) { - int length = Math.abs(samRecord.getReadLength()); - String sample = samRecord.getReadGroup().getSample(); + final int length = Math.abs(samRecord.getReadLength()); + final SAMReadGroupRecord rg = samRecord.getReadGroup(); - table.increment(length, sample); + increment(table,length, rg); return null; } + final private void increment(final Map table,final int length,final SAMReadGroupRecord rg){ + if(readGroupsLocation.isEmpty()){ + if(table.containsKey(length)) + table.get(length)[0]++; + else{ + final int[] newLength = {1}; + table.put(length,newLength); + } + } + else{ + final int rgLocation = readGroupsLocation.get(rg); + if(table.containsKey(length)) + table.get(length)[rgLocation]++; + else{ + table.put(length,new int[readGroupsLocation.size()]); + table.get(length)[rgLocation]++; + } + } + } + @Override public Integer reduceInit() { return null; } @Override - public Integer reduce(Integer integer, Integer integer1) { + public Integer reduce(final Integer integer,final Integer integer1) { return null; } - public void onTraversalDone(Integer sum) { + public void onTraversalDone(final Integer sum) { + final GATKReport report = createGATKReport(); report.print(out); } + + final private GATKReport createGATKReport(){ + final GATKReport report = new GATKReport(); + report.addTable("ReadLengthDistribution", "Table of read length distributions", 1 + (readGroupsLocation.isEmpty() ? 1 : readGroupsLocation.size())); + final GATKReportTable tableReport = report.getTable("ReadLengthDistribution"); + + tableReport.addColumn("readLength"); + + if (readGroupsLocation.isEmpty()){ + tableReport.addColumn("SINGLE_SAMPLE"); + int rowIndex = 0; + for (Integer length : table.keySet()){ + tableReport.set(rowIndex,0,length); + tableReport.set(rowIndex,1,table.get(length)[0]); + rowIndex++; + } + } + else{ + for (SAMReadGroupRecord rg : readGroups) + tableReport.addColumn(rg.getSample()); + int rowIndex = 0; + for (Integer length : table.keySet()){ + tableReport.set(rowIndex,0,length); + for (int i=0; i < readGroupsLocation.size(); i++) + tableReport.set(rowIndex,i+1,table.get(length)[i]); + rowIndex++; + } + + } + + return report; + } }