From e199562c2521032abd003bd315cbfad294208e93 Mon Sep 17 00:00:00 2001 From: Eric Banks Date: Tue, 27 Nov 2012 10:26:17 -0500 Subject: [PATCH 01/52] I have pulled out all of the documentation URLs and put them into the HelpUtils class as static variables; this way, Appistry can change links as needed to point commercial users to their own internal forum without having to muck things up all over our source. Added some TODOs for Geraldine to update links in the GATK docs that still point to the old wiki. Sorry that I am pushing into stable, but that's what Appistry is pulling from for their release next week (and unstable has been failing forever). --- .../sting/commandline/CommandLineProgram.java | 3 ++- .../sting/gatk/CommandLineGATK.java | 24 ++++++++++++------- .../sting/gatk/filters/FilterManager.java | 6 ++--- .../gatk/walkers/diffengine/DiffEngine.java | 3 ++- .../gatk/walkers/diffengine/DiffObjects.java | 3 ++- .../stratifications/JexlExpression.java | 2 +- .../VariantDataManager.java | 3 ++- .../VariantValidationAssessor.java | 3 +-- .../sting/utils/exceptions/UserException.java | 5 ++-- .../sting/utils/help/ForumAPIUtils.java | 8 +++---- .../sting/utils/help/GATKDocUtils.java | 2 +- .../sting/utils/help/HelpUtils.java | 9 +++++++ .../SortingVariantContextWriterBase.java | 1 - .../utils/codecs/vcf/VCFIntegrationTest.java | 1 - 14 files changed, 43 insertions(+), 30 deletions(-) diff --git a/public/java/src/org/broadinstitute/sting/commandline/CommandLineProgram.java b/public/java/src/org/broadinstitute/sting/commandline/CommandLineProgram.java index d77ae67cf..fb15a3722 100644 --- a/public/java/src/org/broadinstitute/sting/commandline/CommandLineProgram.java +++ b/public/java/src/org/broadinstitute/sting/commandline/CommandLineProgram.java @@ -33,6 +33,7 @@ import org.broadinstitute.sting.gatk.CommandLineGATK; import org.broadinstitute.sting.utils.exceptions.ReviewedStingException; import org.broadinstitute.sting.utils.help.ApplicationDetails; import org.broadinstitute.sting.utils.help.HelpFormatter; +import org.broadinstitute.sting.utils.help.HelpUtils; import java.io.IOException; import java.util.*; @@ -288,7 +289,7 @@ public abstract class CommandLineProgram { */ private static void printDocumentationReference() { errorPrintf("Visit our website and forum for extensive documentation and answers to %n"); - errorPrintf("commonly asked questions http://www.broadinstitute.org/gatk%n"); + errorPrintf("commonly asked questions " + HelpUtils.BASE_GATK_URL + "%n"); } diff --git a/public/java/src/org/broadinstitute/sting/gatk/CommandLineGATK.java b/public/java/src/org/broadinstitute/sting/gatk/CommandLineGATK.java index 0daad2c2b..d1711ba4c 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/CommandLineGATK.java +++ b/public/java/src/org/broadinstitute/sting/gatk/CommandLineGATK.java @@ -39,6 +39,7 @@ import org.broadinstitute.sting.utils.exceptions.UserException; import org.broadinstitute.sting.utils.help.ApplicationDetails; import org.broadinstitute.sting.utils.help.DocumentedGATKFeature; import org.broadinstitute.sting.utils.help.GATKDocUtils; +import org.broadinstitute.sting.utils.help.HelpUtils; import org.broadinstitute.sting.utils.text.TextFormattingUtils; import java.util.*; @@ -118,17 +119,24 @@ public class CommandLineGATK extends CommandLineExecutable { public static final String DISK_QUOTA_EXCEEDED_ERROR = "Disk quota exceeded"; private static void checkForMaskedUserErrors(final Throwable t) { + // masked out of memory error + if ( t instanceof OutOfMemoryError ) + exitSystemWithUserError(new UserException.NotEnoughMemory()); + // masked user error + if ( t instanceof UserException || t instanceof TribbleException ) + exitSystemWithUserError(new UserException(t.getMessage())); + + // no message means no masked error final String message = t.getMessage(); if ( message == null ) return; - // we know what to do about the common "Too many open files" error + // too many open files error if ( message.contains("Too many open files") ) exitSystemWithUserError(new UserException.TooManyOpenFiles()); // malformed BAM looks like a SAM file - if ( message.contains(PICARD_TEXT_SAM_FILE_ERROR_1) || - message.contains(PICARD_TEXT_SAM_FILE_ERROR_2) ) + if ( message.contains(PICARD_TEXT_SAM_FILE_ERROR_1) || message.contains(PICARD_TEXT_SAM_FILE_ERROR_2) ) exitSystemWithSamError(t); // can't close tribble index when writing @@ -138,12 +146,10 @@ public class CommandLineGATK extends CommandLineExecutable { // disk is full if ( message.contains(NO_SPACE_LEFT_ON_DEVICE_ERROR) || message.contains(DISK_QUOTA_EXCEEDED_ERROR) ) exitSystemWithUserError(new UserException.NoSpaceOnDevice()); - if ( t.getCause() != null && (t.getCause().getMessage().contains(NO_SPACE_LEFT_ON_DEVICE_ERROR) || t.getCause().getMessage().contains(DISK_QUOTA_EXCEEDED_ERROR)) ) - exitSystemWithUserError(new UserException.NoSpaceOnDevice()); - // masked out of memory error - if ( t.getCause() != null && t.getCause() instanceof OutOfMemoryError ) - exitSystemWithUserError(new UserException.NotEnoughMemory()); + // masked error wrapped in another one + if ( t.getCause() != null ) + checkForMaskedUserErrors(t.getCause()); } /** @@ -155,7 +161,7 @@ public class CommandLineGATK extends CommandLineExecutable { List header = new ArrayList(); header.add(String.format("The Genome Analysis Toolkit (GATK) v%s, Compiled %s",getVersionNumber(), getBuildTime())); header.add("Copyright (c) 2010 The Broad Institute"); - header.add("For support and documentation go to http://www.broadinstitute.org/gatk"); + header.add("For support and documentation go to " + HelpUtils.BASE_GATK_URL); return header; } diff --git a/public/java/src/org/broadinstitute/sting/gatk/filters/FilterManager.java b/public/java/src/org/broadinstitute/sting/gatk/filters/FilterManager.java index 5ca8a1779..89099c587 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/filters/FilterManager.java +++ b/public/java/src/org/broadinstitute/sting/gatk/filters/FilterManager.java @@ -25,11 +25,9 @@ package org.broadinstitute.sting.gatk.filters; -import com.google.common.base.Function; -import com.google.common.collect.Collections2; -import org.broadinstitute.sting.utils.Utils; import org.broadinstitute.sting.utils.classloader.PluginManager; import org.broadinstitute.sting.utils.help.GATKDocUtils; +import org.broadinstitute.sting.utils.help.HelpUtils; import java.util.Collection; import java.util.List; @@ -73,7 +71,7 @@ public class FilterManager extends PluginManager { return String.format("Read filter %s not found. Available read filters:%n%n%s%n%n%s",pluginName, userFriendlyListofReadFilters(availableFilters), - "Please consult the GATK Documentation (http://www.broadinstitute.org/gatk/gatkdocs/) for more information."); + "Please consult the GATK Documentation (" + HelpUtils.GATK_DOCS_URL + ") for more information."); } private String userFriendlyListofReadFilters(List> filters) { diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/diffengine/DiffEngine.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/diffengine/DiffEngine.java index 40ed26608..579b84d96 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/diffengine/DiffEngine.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/diffengine/DiffEngine.java @@ -282,7 +282,8 @@ public class DiffEngine { // now that we have a specific list of values we want to show, display them GATKReport report = new GATKReport(); final String tableName = "differences"; - report.addTable(tableName, "Summarized differences between the master and test files. See http://www.broadinstitute.org/gsa/wiki/index.php/DiffEngine for more information", 3); + // TODO for Geraldine -- link needs to be updated below + report.addTable(tableName, "Summarized differences between the master and test files. See [ask Geraldine to fix link to DiffEngine wiki] for more information", 3); final GATKReportTable table = report.getTable(tableName); table.addColumn("Difference"); table.addColumn("NumberOfOccurrences"); diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/diffengine/DiffObjects.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/diffengine/DiffObjects.java index 92e2e2dc4..5951ee7d0 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/diffengine/DiffObjects.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/diffengine/DiffObjects.java @@ -138,7 +138,8 @@ public class DiffObjects extends RodWalker { /** * Writes out a file of the DiffEngine format: * - * http://www.broadinstitute.org/gsa/wiki/index.php/DiffEngine + * TODO for Geraldine -- link needs to be updated below (and also in SelectVariants and RefSeqCodec GATK docs) + * http://www.broadinstitute.org/gsa/wiki/index.php/DiffEngine */ @Output(doc="File to which results should be written",required=true) protected PrintStream out; diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/varianteval/stratifications/JexlExpression.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/varianteval/stratifications/JexlExpression.java index dc5438358..c89c4be66 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/varianteval/stratifications/JexlExpression.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/varianteval/stratifications/JexlExpression.java @@ -13,7 +13,7 @@ import java.util.Set; /** * Stratifies the eval RODs by user-supplied JEXL expressions * - * See http://www.broadinstitute.org/gsa/wiki/index.php/Using_JEXL_expressions for more details + * See http://gatkforums.broadinstitute.org/discussion/1255/what-are-jexl-expressions-and-how-can-i-use-them-with-the-gatk for more details */ public class JexlExpression extends VariantStratifier implements StandardStratification { // needs to know the jexl expressions diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantDataManager.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantDataManager.java index aacd987d5..3382a1d9b 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantDataManager.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantDataManager.java @@ -32,6 +32,7 @@ import org.broadinstitute.sting.utils.GenomeLoc; import org.broadinstitute.sting.utils.MathUtils; import org.broadinstitute.sting.utils.codecs.vcf.VCFConstants; import org.broadinstitute.sting.utils.exceptions.ReviewedStingException; +import org.broadinstitute.sting.utils.help.HelpUtils; import org.broadinstitute.sting.utils.variantcontext.writer.VariantContextWriter; import org.broadinstitute.sting.utils.collections.ExpandingArrayList; import org.broadinstitute.sting.utils.exceptions.UserException; @@ -80,7 +81,7 @@ public class VariantDataManager { final double theSTD = standardDeviation(theMean, iii); logger.info( annotationKeys.get(iii) + String.format(": \t mean = %.2f\t standard deviation = %.2f", theMean, theSTD) ); if( Double.isNaN(theMean) ) { - throw new UserException.BadInput("Values for " + annotationKeys.get(iii) + " annotation not detected for ANY training variant in the input callset. VariantAnnotator may be used to add these annotations. See http://www.broadinstitute.org/gsa/wiki/index.php/VariantAnnotator"); + throw new UserException.BadInput("Values for " + annotationKeys.get(iii) + " annotation not detected for ANY training variant in the input callset. VariantAnnotator may be used to add these annotations. See " + HelpUtils.GATK_FORUM_URL + "discussion/49/using-variant-annotator"); } foundZeroVarianceAnnotation = foundZeroVarianceAnnotation || (theSTD < 1E-6); diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/VariantValidationAssessor.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/VariantValidationAssessor.java index a301867fc..9236247f1 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/VariantValidationAssessor.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/VariantValidationAssessor.java @@ -51,8 +51,7 @@ import java.util.*; * The Variant Validation Assessor is a tool for vetting/assessing validation data (containing genotypes). * The tool produces a VCF that is annotated with information pertaining to plate quality control and by * default is soft-filtered by high no-call rate or low Hardy-Weinberg probability. - * If you have .ped files, please first convert them to VCF format - * (see http://www.broadinstitute.org/gsa/wiki/index.php/Converting_ped_to_vcf). + * If you have .ped files, please first convert them to VCF format. * *

Input

*

diff --git a/public/java/src/org/broadinstitute/sting/utils/exceptions/UserException.java b/public/java/src/org/broadinstitute/sting/utils/exceptions/UserException.java index a49a12292..a2ec35ae2 100755 --- a/public/java/src/org/broadinstitute/sting/utils/exceptions/UserException.java +++ b/public/java/src/org/broadinstitute/sting/utils/exceptions/UserException.java @@ -30,6 +30,7 @@ import net.sf.samtools.SAMSequenceDictionary; import org.broadinstitute.sting.gatk.phonehome.GATKRunReport; import org.broadinstitute.sting.utils.GenomeLoc; import org.broadinstitute.sting.utils.help.DocumentedGATKFeature; +import org.broadinstitute.sting.utils.help.HelpUtils; import org.broadinstitute.sting.utils.sam.ReadUtils; import org.broadinstitute.sting.utils.variantcontext.VariantContext; @@ -267,7 +268,7 @@ public class UserException extends ReviewedStingException { public static class ReadMissingReadGroup extends MalformedBAM { public ReadMissingReadGroup(SAMRecord read) { - super(read, String.format("Read %s is either missing the read group or its read group is not defined in the BAM header, both of which are required by the GATK. Please use http://www.broadinstitute.org/gsa/wiki/index.php/ReplaceReadGroups to fix this problem", read.getReadName())); + super(read, String.format("Read %s is either missing the read group or its read group is not defined in the BAM header, both of which are required by the GATK. Please use " + HelpUtils.GATK_FORUM_URL + "discussion/59/companion-utilities-replacereadgroups to fix this problem", read.getReadName())); } } @@ -343,7 +344,7 @@ public class UserException extends ReviewedStingException { super(String.format("Lexicographically sorted human genome sequence detected in %s." + "\nFor safety's sake the GATK requires human contigs in karyotypic order: 1, 2, ..., 10, 11, ..., 20, 21, 22, X, Y with M either leading or trailing these contigs." + "\nThis is because all distributed GATK resources are sorted in karyotypic order, and your processing will fail when you need to use these files." - + "\nYou can use the ReorderSam utility to fix this problem: http://www.broadinstitute.org/gsa/wiki/index.php/ReorderSam" + + "\nYou can use the ReorderSam utility to fix this problem: " + HelpUtils.GATK_FORUM_URL + "discussion/58/companion-utilities-reordersam" + "\n %s contigs = %s", name, name, ReadUtils.prettyPrintSequenceRecords(dict))); } diff --git a/public/java/src/org/broadinstitute/sting/utils/help/ForumAPIUtils.java b/public/java/src/org/broadinstitute/sting/utils/help/ForumAPIUtils.java index fe5f48a48..64238dc73 100644 --- a/public/java/src/org/broadinstitute/sting/utils/help/ForumAPIUtils.java +++ b/public/java/src/org/broadinstitute/sting/utils/help/ForumAPIUtils.java @@ -44,14 +44,13 @@ public class ForumAPIUtils { /** * How we post to the forum */ - private final static String API_URL = "https://gatkforums.broadinstitute.org/api/v1/"; final private static String ACCESS_TOKEN = "access_token="; public static List getPostedTools(String forumKey) { Gson gson = new Gson(); List output = new ArrayList(); - String text = httpGet(API_URL + "categories.json?CategoryIdentifier=tool-bulletin&page=1-100000&" + ACCESS_TOKEN + forumKey); + String text = httpGet(HelpUtils.GATK_FORUM_API_URL + "categories.json?CategoryIdentifier=tool-bulletin&page=1-100000&" + ACCESS_TOKEN + forumKey); APIQuery details = gson.fromJson(text, APIQuery.class); ForumDiscussion[] discussions = details.Discussions; @@ -159,7 +158,7 @@ public class ForumAPIUtils { Gson gson = new Gson(); String data = gson.toJson(post.getPostData()); - httpPost(data, API_URL + "post/discussion.json?" + ACCESS_TOKEN + forumKey); + httpPost(data, HelpUtils.GATK_FORUM_API_URL + "post/discussion.json?" + ACCESS_TOKEN + forumKey); } @@ -167,8 +166,7 @@ public class ForumAPIUtils { class APIQuery { ForumDiscussion[] Discussions; - public APIQuery() { - } + public APIQuery() {} } } diff --git a/public/java/src/org/broadinstitute/sting/utils/help/GATKDocUtils.java b/public/java/src/org/broadinstitute/sting/utils/help/GATKDocUtils.java index 4ec2ac6d7..21054a794 100644 --- a/public/java/src/org/broadinstitute/sting/utils/help/GATKDocUtils.java +++ b/public/java/src/org/broadinstitute/sting/utils/help/GATKDocUtils.java @@ -28,7 +28,7 @@ public class GATKDocUtils { /** * The URL root for RELEASED GATKDOC units */ - public final static String URL_ROOT_FOR_RELEASE_GATKDOCS = "http://www.broadinstitute.org/gatk/gatkdocs/"; + public final static String URL_ROOT_FOR_RELEASE_GATKDOCS = HelpUtils.GATK_DOCS_URL; /** * The URL root for STABLE GATKDOC units */ diff --git a/public/java/src/org/broadinstitute/sting/utils/help/HelpUtils.java b/public/java/src/org/broadinstitute/sting/utils/help/HelpUtils.java index 645ab34c1..1bc20d5a0 100644 --- a/public/java/src/org/broadinstitute/sting/utils/help/HelpUtils.java +++ b/public/java/src/org/broadinstitute/sting/utils/help/HelpUtils.java @@ -32,6 +32,15 @@ import org.broadinstitute.sting.utils.classloader.JVMUtils; import java.lang.reflect.Field; public class HelpUtils { + + public final static String BASE_GATK_URL = "http://www.broadinstitute.org/gatk"; + public final static String GATK_DOCS_URL = BASE_GATK_URL + "/gatkdocs/"; + public final static String GATK_FORUM_URL = "http://gatkforums.broadinstitute.org/"; + public final static String GATK_FORUM_API_URL = "https://gatkforums.broadinstitute.org/api/v1/"; + + + + protected static boolean assignableToClass(ProgramElementDoc classDoc, Class lhsClass, boolean requireConcrete) { try { Class type = getClassForDoc(classDoc); diff --git a/public/java/src/org/broadinstitute/sting/utils/variantcontext/writer/SortingVariantContextWriterBase.java b/public/java/src/org/broadinstitute/sting/utils/variantcontext/writer/SortingVariantContextWriterBase.java index 18d91ef3f..1f3cdd0fe 100755 --- a/public/java/src/org/broadinstitute/sting/utils/variantcontext/writer/SortingVariantContextWriterBase.java +++ b/public/java/src/org/broadinstitute/sting/utils/variantcontext/writer/SortingVariantContextWriterBase.java @@ -71,7 +71,6 @@ abstract class SortingVariantContextWriterBase implements VariantContextWriter { this.takeOwnershipOfInner = takeOwnershipOfInner; // has to be PriorityBlockingQueue to be thread-safe - // see http://getsatisfaction.com/gsa/topics/missing_loci_output_in_multi_thread_mode_when_implement_sortingvcfwriterbase?utm_content=topic_link&utm_medium=email&utm_source=new_topic this.queue = new PriorityBlockingQueue(50, new VariantContextComparator()); this.mostUpstreamWritableLoc = BEFORE_MOST_UPSTREAM_LOC; diff --git a/public/java/test/org/broadinstitute/sting/utils/codecs/vcf/VCFIntegrationTest.java b/public/java/test/org/broadinstitute/sting/utils/codecs/vcf/VCFIntegrationTest.java index b2a4ac2da..b9ce58992 100644 --- a/public/java/test/org/broadinstitute/sting/utils/codecs/vcf/VCFIntegrationTest.java +++ b/public/java/test/org/broadinstitute/sting/utils/codecs/vcf/VCFIntegrationTest.java @@ -28,7 +28,6 @@ public class VCFIntegrationTest extends WalkerTest { } @Test(enabled = true) - // See https://getsatisfaction.com/gsa/topics/support_vcf_4_1_structural_variation_breakend_alleles?utm_content=topic_link&utm_medium=email&utm_source=new_topic public void testReadingAndWritingBreakpointAlleles() { String testVCF = privateTestDir + "breakpoint-example.vcf"; From 4543ece0889cc07acc9d2a9cbefda7882ae340ef Mon Sep 17 00:00:00 2001 From: Eric Banks Date: Tue, 27 Nov 2012 11:00:33 -0500 Subject: [PATCH 02/52] Fixing parsing of genomelocs that contain colons in the contig names (which is allowed by the spec) as reported on the forum. Added unit test for this case. --- .../sting/utils/GenomeLocParser.java | 2 +- .../sting/utils/GenomeLocParserUnitTest.java | 21 ++++++++++++++++++- 2 files changed, 21 insertions(+), 2 deletions(-) diff --git a/public/java/src/org/broadinstitute/sting/utils/GenomeLocParser.java b/public/java/src/org/broadinstitute/sting/utils/GenomeLocParser.java index a3ffe708c..bf60b4a80 100644 --- a/public/java/src/org/broadinstitute/sting/utils/GenomeLocParser.java +++ b/public/java/src/org/broadinstitute/sting/utils/GenomeLocParser.java @@ -374,7 +374,7 @@ public final class GenomeLocParser { int start = 1; int stop = -1; - final int colonIndex = str.indexOf(":"); + final int colonIndex = str.lastIndexOf(":"); if(colonIndex == -1) { contig = str.substring(0, str.length()); // chr1 stop = Integer.MAX_VALUE; diff --git a/public/java/test/org/broadinstitute/sting/utils/GenomeLocParserUnitTest.java b/public/java/test/org/broadinstitute/sting/utils/GenomeLocParserUnitTest.java index e9f138a0e..e4313b30a 100644 --- a/public/java/test/org/broadinstitute/sting/utils/GenomeLocParserUnitTest.java +++ b/public/java/test/org/broadinstitute/sting/utils/GenomeLocParserUnitTest.java @@ -2,6 +2,8 @@ package org.broadinstitute.sting.utils; import net.sf.samtools.SAMFileHeader; +import net.sf.samtools.SAMSequenceDictionary; +import net.sf.samtools.SAMSequenceRecord; import org.broadinstitute.sting.BaseTest; import org.broadinstitute.sting.utils.exceptions.ReviewedStingException; import org.broadinstitute.sting.utils.exceptions.UserException; @@ -74,6 +76,23 @@ public class GenomeLocParserUnitTest extends BaseTest { genomeLocParser.parseGenomeLoc("Bad:0-1"); } + @Test + public void testContigHasColon() { + SAMFileHeader header = new SAMFileHeader(); + header.setSortOrder(net.sf.samtools.SAMFileHeader.SortOrder.coordinate); + SAMSequenceDictionary dict = new SAMSequenceDictionary(); + SAMSequenceRecord rec = new SAMSequenceRecord("c:h:r1", 10); + rec.setSequenceLength(10); + dict.addSequence(rec); + header.setSequenceDictionary(dict); + + final GenomeLocParser myGenomeLocParser = new GenomeLocParser(header.getSequenceDictionary()); + GenomeLoc loc = myGenomeLocParser.parseGenomeLoc("c:h:r1:4-5"); + assertEquals(0, loc.getContigIndex()); + assertEquals(loc.getStart(), 4); + assertEquals(loc.getStop(), 5); + } + @Test public void testParseGoodString() { GenomeLoc loc = genomeLocParser.parseGenomeLoc("chr1:1-10"); @@ -81,7 +100,7 @@ public class GenomeLocParserUnitTest extends BaseTest { assertEquals(loc.getStop(), 10); assertEquals(loc.getStart(), 1); } - + @Test public void testCreateGenomeLoc1() { GenomeLoc loc = genomeLocParser.createGenomeLoc("chr1", 1, 100); From cc550b4145bd8f439a46e407ba015eaeea36edea Mon Sep 17 00:00:00 2001 From: Joel Thibault Date: Mon, 26 Nov 2012 11:48:05 -0500 Subject: [PATCH 03/52] Add a read and interval on a different contig --- .../traversals/TraverseActiveRegionsTest.java | 27 ++++++++++++++----- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java b/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java index e4c7b2db0..018e92d84 100644 --- a/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java +++ b/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java @@ -96,8 +96,7 @@ public class TraverseActiveRegionsTest extends BaseTest { intervals.add(genomeLocParser.createGenomeLoc("1", 1000, 1999)); intervals.add(genomeLocParser.createGenomeLoc("1", 2000, 2999)); intervals.add(genomeLocParser.createGenomeLoc("1", 10000, 20000)); - // TODO: this fails! - //intervals.add(genomeLocParser.createGenomeLoc("20", 10000, 20000)); + intervals.add(genomeLocParser.createGenomeLoc("20", 10000, 10100)); intervals = IntervalUtils.sortAndMergeIntervals(genomeLocParser, intervals, IntervalMergingRule.OVERLAPPING_ONLY).toList(); reads = new ArrayList(); @@ -109,8 +108,7 @@ public class TraverseActiveRegionsTest extends BaseTest { reads.add(buildSAMRecord("extended_only", "1", 3000, 3100)); reads.add(buildSAMRecord("extended_and_np", "1", 990, 1990)); reads.add(buildSAMRecord("outside_intervals", "1", 5000, 6000)); - // TODO - //reads.add(buildSAMRecord("simple20", "20", 10100, 10150)); + reads.add(buildSAMRecord("simple20", "20", 10025, 10075)); } @Test @@ -204,9 +202,7 @@ public class TraverseActiveRegionsTest extends BaseTest { // extended_only: Extended in 1:2000-2999 // extended_and_np: Non-Primary in 1:1-999, Primary in 1:1000-1999, Extended in 1:2000-2999 // outside_intervals: none - - // TODO - // simple20: Primary in 20:10000-20000 + // simple20: Primary in 20:10000-10100 Map activeRegions = getActiveRegions(walker, intervals); ActiveRegion region; @@ -221,6 +217,7 @@ public class TraverseActiveRegionsTest extends BaseTest { verifyReadNotPlaced(region, "extended_only"); // TODO: fail verifyReadNonPrimary(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); + verifyReadNotPlaced(region, "simple20"); region = activeRegions.get(genomeLocParser.createGenomeLoc("1", 1000, 1999)); @@ -232,6 +229,7 @@ public class TraverseActiveRegionsTest extends BaseTest { verifyReadNotPlaced(region, "extended_only"); // TODO: fail verifyReadPrimary(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); + verifyReadNotPlaced(region, "simple20"); region = activeRegions.get(genomeLocParser.createGenomeLoc("1", 2000, 2999)); @@ -243,6 +241,19 @@ public class TraverseActiveRegionsTest extends BaseTest { // TODO: fail verifyReadExtended(region, "extended_only"); // TODO: fail verifyReadExtended(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); + verifyReadNotPlaced(region, "simple20"); + + region = activeRegions.get(genomeLocParser.createGenomeLoc("20", 10000, 10100)); + + verifyReadNotPlaced(region, "simple"); + verifyReadNotPlaced(region, "overlap_equal"); + verifyReadNotPlaced(region, "overlap_unequal"); + verifyReadNotPlaced(region, "boundary_equal"); + verifyReadNotPlaced(region, "boundary_unequal"); + verifyReadNotPlaced(region, "extended_only"); + verifyReadNotPlaced(region, "extended_and_np"); + verifyReadNotPlaced(region, "outside_intervals"); + verifyReadPrimary(region, "simple20"); // TODO: more tests and edge cases } @@ -282,6 +293,8 @@ public class TraverseActiveRegionsTest extends BaseTest { for (LocusShardDataProvider dataProvider : createDataProviders(intervals)) t.traverse(walker, dataProvider, 0); + t.endTraversal(walker, 0); + return walker.mappedActiveRegions; } From d83ad906eff14027908376eceb37dd502a6fdd78 Mon Sep 17 00:00:00 2001 From: Joel Thibault Date: Mon, 26 Nov 2012 13:44:13 -0500 Subject: [PATCH 04/52] Add profile range contract --- .../sting/gatk/walkers/ActiveRegionWalker.java | 2 ++ .../traversals/TraverseActiveRegionsTest.java | 17 +++++++++++++++++ 2 files changed, 19 insertions(+) diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/ActiveRegionWalker.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/ActiveRegionWalker.java index fed2c995e..c6e28df05 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/ActiveRegionWalker.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/ActiveRegionWalker.java @@ -1,5 +1,6 @@ package org.broadinstitute.sting.gatk.walkers; +import com.google.java.contract.Ensures; import net.sf.picard.reference.IndexedFastaSequenceFile; import org.broad.tribble.Feature; import org.broadinstitute.sting.commandline.Input; @@ -75,6 +76,7 @@ public abstract class ActiveRegionWalker extends Walker= 0.0", "result.isActiveProb <= 1.0"}) public abstract ActivityProfileResult isActive(final RefMetaDataTracker tracker, final ReferenceContext ref, final AlignmentContext context); // Map over the ActiveRegion diff --git a/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java b/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java index 018e92d84..8a4be48be 100644 --- a/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java +++ b/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java @@ -1,5 +1,6 @@ package org.broadinstitute.sting.gatk.traversals; +import com.google.java.contract.PreconditionError; import net.sf.samtools.*; import org.broadinstitute.sting.utils.interval.IntervalMergingRule; import org.broadinstitute.sting.utils.interval.IntervalUtils; @@ -52,6 +53,10 @@ public class TraverseActiveRegionsTest extends BaseTest { this.prob = 1.0; } + public DummyActiveRegionWalker(double constProb) { + this.prob = constProb; + } + @Override public ActivityProfileResult isActive(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) { isActiveCalls.add(ref.getLocus()); @@ -132,6 +137,18 @@ public class TraverseActiveRegionsTest extends BaseTest { return activeIntervals; } + @Test (expectedExceptions = PreconditionError.class) + public void testIsActiveRangeLow () { + DummyActiveRegionWalker walker = new DummyActiveRegionWalker(-0.1); + getActiveRegions(walker, intervals).values(); + } + + @Test (expectedExceptions = PreconditionError.class) + public void testIsActiveRangeHigh () { + DummyActiveRegionWalker walker = new DummyActiveRegionWalker(1.1); + getActiveRegions(walker, intervals).values(); + } + @Test public void testActiveRegionCoverage() { DummyActiveRegionWalker walker = new DummyActiveRegionWalker(); From 9bfe39411ee4465860d6cf1a1cb1f0fe32d0a1b3 Mon Sep 17 00:00:00 2001 From: Joel Thibault Date: Mon, 26 Nov 2012 14:29:22 -0500 Subject: [PATCH 05/52] Equal overlap should match right/later region --- .../gatk/traversals/TraverseActiveRegionsTest.java | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java b/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java index 8a4be48be..66504da11 100644 --- a/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java +++ b/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java @@ -109,7 +109,7 @@ public class TraverseActiveRegionsTest extends BaseTest { reads.add(buildSAMRecord("overlap_equal", "1", 10, 20)); reads.add(buildSAMRecord("overlap_unequal", "1", 10, 21)); reads.add(buildSAMRecord("boundary_equal", "1", 1990, 2009)); - reads.add(buildSAMRecord("boundary_unequal", "1", 1995, 2050)); + reads.add(buildSAMRecord("boundary_unequal", "1", 1990, 2008)); reads.add(buildSAMRecord("extended_only", "1", 3000, 3100)); reads.add(buildSAMRecord("extended_and_np", "1", 990, 1990)); reads.add(buildSAMRecord("outside_intervals", "1", 5000, 6000)); @@ -214,8 +214,8 @@ public class TraverseActiveRegionsTest extends BaseTest { // simple: Primary in 1:1-999 // overlap_equal: Primary in 1:1-999 // overlap_unequal: Primary in 1:1-999 - // boundary_equal: Primary in 1:1000-1999, Non-Primary in 1:2000-2999 - // boundary_unequal: Non-Primary in 1:1000-1999, Primary in 1:2000-2999 + // boundary_equal: Non-Primary in 1:1000-1999, Primary in 1:2000-2999 + // boundary_unequal: Primary in 1:1000-1999, Non-Primary in 1:2000-2999 // extended_only: Extended in 1:2000-2999 // extended_and_np: Non-Primary in 1:1-999, Primary in 1:1000-1999, Extended in 1:2000-2999 // outside_intervals: none @@ -241,8 +241,8 @@ public class TraverseActiveRegionsTest extends BaseTest { verifyReadNotPlaced(region, "simple"); verifyReadNotPlaced(region, "overlap_equal"); verifyReadNotPlaced(region, "overlap_unequal"); - // TODO: fail verifyReadPrimary(region, "boundary_equal"); - // TODO: fail verifyReadNonPrimary(region, "boundary_unequal"); + // TODO: fail verifyReadNonPrimary(region, "boundary_equal"); + verifyReadPrimary(region, "boundary_unequal"); verifyReadNotPlaced(region, "extended_only"); // TODO: fail verifyReadPrimary(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); @@ -253,8 +253,8 @@ public class TraverseActiveRegionsTest extends BaseTest { verifyReadNotPlaced(region, "simple"); verifyReadNotPlaced(region, "overlap_equal"); verifyReadNotPlaced(region, "overlap_unequal"); - // TODO: fail verifyReadNonPrimary(region, "boundary_equal"); - verifyReadPrimary(region, "boundary_unequal"); + verifyReadPrimary(region, "boundary_equal"); + // TODO: fail verifyReadNonPrimary(region, "boundary_unequal"); // TODO: fail verifyReadExtended(region, "extended_only"); // TODO: fail verifyReadExtended(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); From 7e4b9c9e6e38a1f20999fa0b6b48a5ce2e313c5f Mon Sep 17 00:00:00 2001 From: Mark DePristo Date: Tue, 27 Nov 2012 10:12:39 -0500 Subject: [PATCH 07/52] Fix failing unit tests for VariantContextUtilsUnitTest -- Previous version was adding multiple samples with the same name to the variant context --- .../sting/utils/variantcontext/VariantContextUtilsUnitTest.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/public/java/test/org/broadinstitute/sting/utils/variantcontext/VariantContextUtilsUnitTest.java b/public/java/test/org/broadinstitute/sting/utils/variantcontext/VariantContextUtilsUnitTest.java index f3daa9e4c..3ad438b26 100644 --- a/public/java/test/org/broadinstitute/sting/utils/variantcontext/VariantContextUtilsUnitTest.java +++ b/public/java/test/org/broadinstitute/sting/utils/variantcontext/VariantContextUtilsUnitTest.java @@ -793,7 +793,7 @@ public class VariantContextUtilsUnitTest extends BaseTest { int sampleI = 0; for ( final List alleles : Utils.makePermutations(vc.getAlleles(), 2, true) ) { - genotypes.add(GenotypeBuilder.create("sample" + sampleI, alleles)); + genotypes.add(GenotypeBuilder.create("sample" + sampleI++, alleles)); } genotypes.add(GenotypeBuilder.createMissing("missing", 2)); From 01abcc3e0f718a2ca17e8b182bb642a7926cf2ef Mon Sep 17 00:00:00 2001 From: Eric Banks Date: Tue, 27 Nov 2012 14:40:49 -0500 Subject: [PATCH 09/52] Tests didn't like my note to Geraldine in the output logs; apparently it's tested in integration tests --- .../sting/gatk/walkers/diffengine/DiffEngine.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/diffengine/DiffEngine.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/diffengine/DiffEngine.java index 579b84d96..0da23077d 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/diffengine/DiffEngine.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/diffengine/DiffEngine.java @@ -283,7 +283,7 @@ public class DiffEngine { GATKReport report = new GATKReport(); final String tableName = "differences"; // TODO for Geraldine -- link needs to be updated below - report.addTable(tableName, "Summarized differences between the master and test files. See [ask Geraldine to fix link to DiffEngine wiki] for more information", 3); + report.addTable(tableName, "Summarized differences between the master and test files. See http://www.broadinstitute.org/gsa/wiki/index.php/DiffEngine for more information", 3); final GATKReportTable table = report.getTable(tableName); table.addColumn("Difference"); table.addColumn("NumberOfOccurrences"); From 1cc0b48caab07426a3d54b34db3043ca96a28a4e Mon Sep 17 00:00:00 2001 From: Jacob Silterra Date: Tue, 27 Nov 2012 17:44:35 -0500 Subject: [PATCH 13/52] Abstract connection to MongoDB so we can specify it through JSON file. Include 2 JSON spec files in GenomeAnalysisTK.jar Create MongoDBManager, which keeps track of connections based on Locator class. Locators can be instantiated directly, or read from JSON files (NA12878DBArgumentCollection uses the GSon library) --- build.xml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/build.xml b/build.xml index a93918ec8..3a264b476 100644 --- a/build.xml +++ b/build.xml @@ -681,6 +681,9 @@ + + + From 79bc878e6a9a280fd7873eca7b2861690dbf2628 Mon Sep 17 00:00:00 2001 From: Menachem Fromer Date: Tue, 27 Nov 2012 22:37:41 -0500 Subject: [PATCH 14/52] Allow debugging to be set from the command line --- .../sting/gatk/walkers/phasing/ReadBackedPhasing.java | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/phasing/ReadBackedPhasing.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/phasing/ReadBackedPhasing.java index d8ae6b28b..eda43e6a5 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/phasing/ReadBackedPhasing.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/phasing/ReadBackedPhasing.java @@ -95,7 +95,8 @@ import static org.broadinstitute.sting.utils.codecs.vcf.VCFUtils.getVCFHeadersFr @DocumentedGATKFeature( groupName = "Variant Discovery Tools", extraDocs = {CommandLineGATK.class} ) public class ReadBackedPhasing extends RodWalker { - private static final boolean DEBUG = false; + @Argument(fullName="debug", shortName="debug", doc="If specified, print out very verbose debug information (if -l DEBUG is also specified)", required = false) + protected boolean DEBUG = false; /** * The VCF file we are phasing variants from. * @@ -949,7 +950,7 @@ public class ReadBackedPhasing extends RodWalker Date: Wed, 28 Nov 2012 11:35:41 -0500 Subject: [PATCH 16/52] Critical bugfix to AFCalcResult affecting UG/HC quality score emission thresholds As reported by Menachem Fromer: a critical bug in AFCalcResult: Specifically, the implementation: public boolean isPolymorphic(final Allele allele, final double log10minPNonRef) { return getLog10PosteriorOfAFGt0ForAllele(allele) >= log10minPNonRef; } seems incorrect and should probably be: getLog10PosteriorOfAFEq0ForAllele(allele) <= log10minPNonRef The issue here is that the 30 represents a Phred-scaled probability of *error* and it's currently being compared to a log probability of *non-error*. Instead, we need to require that our probability of error be less than the error threshold. This bug has only a minor impact on the calls -- hardly any sites change -- which is good. But the inverted logic effects multi-allelic sites significantly. Basically you only hit this logic with multiple alleles, and in that case it'\s including extra alt alleles incorrectly, and throwing out good ones. Change was to create a new function that properly handles thresholds that are PhredScaled quality scores: /** * Same as #isPolymorphic but takes a phred-scaled quality score as input */ public boolean isPolymorphicPhredScaledQual(final Allele allele, final double minPNonRefPhredScaledQual) { if ( minPNonRefPhredScaledQual < 0 ) throw new IllegalArgumentException("phredScaledQual " + minPNonRefPhredScaledQual + " < 0 "); final double log10Threshold = Math.log10(QualityUtils.qualToProb(minPNonRefPhredScaledQual)); return isPolymorphic(allele, log10Threshold); } --- .../UnifiedGenotyperIntegrationTest.java | 8 +-- .../afcalc/AFCalcResultUnitTest.java | 56 +++++++++++++++++-- .../genotyper/UnifiedGenotyperEngine.java | 8 +-- .../genotyper/afcalc/AFCalcResult.java | 11 ++++ 4 files changed, 70 insertions(+), 13 deletions(-) diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java index 9212d0e53..88f2ab3ea 100755 --- a/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java +++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java @@ -62,7 +62,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest { public void testMultipleSNPAlleles() { WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec( "-T UnifiedGenotyper -R " + b37KGReference + " --no_cmdline_in_header -glm BOTH --dbsnp " + b37dbSNP129 + " -I " + privateTestDir + "multiallelic.snps.bam -o %s -L " + privateTestDir + "multiallelic.snps.intervals", 1, - Arrays.asList("94dc17d76d841f1d3a36160767ffa034")); + Arrays.asList("a373979d01c3a3fb20159235d27eb92c")); executeTest("test Multiple SNP alleles", spec); } @@ -197,7 +197,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest { @Test public void testOutputParameterAllSites() { - testOutputParameters("--output_mode EMIT_ALL_SITES", "8a263fd0a94463ce1de9990f2b8ec841"); + testOutputParameters("--output_mode EMIT_ALL_SITES", "5c75cecb523cac988beecd59186289ff"); } private void testOutputParameters(final String args, final String md5) { @@ -345,13 +345,13 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest { public void testMultiSampleIndels1() { WalkerTest.WalkerTestSpec spec1 = new WalkerTest.WalkerTestSpec( baseCommandIndels + " -I " + validationDataLocation + "low_coverage_CEU.chr1.10k-11k.bam -o %s -L 1:10450700-10551000", 1, - Arrays.asList("f7d0d0aee603df25c1f0525bb8df189e")); + Arrays.asList("a5a81bf1b10be860a6a5272fb928e8eb")); List result = executeTest("test MultiSample Pilot1 CEU indels", spec1).getFirst(); WalkerTest.WalkerTestSpec spec2 = new WalkerTest.WalkerTestSpec( baseCommandIndels + " --genotyping_mode GENOTYPE_GIVEN_ALLELES -alleles " + result.get(0).getAbsolutePath() + " -I " + validationDataLocation + "low_coverage_CEU.chr1.10k-11k.bam -o %s -L 1:10450700-10551000", 1, - Arrays.asList("fc91d457a16b4ca994959c2b5f3f0352")); + Arrays.asList("ad52814cd6c45df424fc992699feead6")); executeTest("test MultiSample Pilot1 CEU indels using GENOTYPE_GIVEN_ALLELES", spec2); } diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcResultUnitTest.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcResultUnitTest.java index cbe2eb268..96e055e92 100644 --- a/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcResultUnitTest.java +++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcResultUnitTest.java @@ -2,16 +2,14 @@ package org.broadinstitute.sting.gatk.walkers.genotyper.afcalc; import org.broadinstitute.sting.BaseTest; import org.broadinstitute.sting.utils.MathUtils; +import org.broadinstitute.sting.utils.QualityUtils; import org.broadinstitute.sting.utils.Utils; import org.broadinstitute.sting.utils.variantcontext.Allele; import org.testng.Assert; import org.testng.annotations.DataProvider; import org.testng.annotations.Test; -import java.util.ArrayList; -import java.util.Arrays; -import java.util.Collections; -import java.util.List; +import java.util.*; public class AFCalcResultUnitTest extends BaseTest { private static class MyTest { @@ -79,4 +77,54 @@ public class AFCalcResultUnitTest extends BaseTest { final double[] actualPosteriors = new double[]{result.getLog10PosteriorOfAFEq0(), result.getLog10PosteriorOfAFGT0()}; Assert.assertEquals(MathUtils.sumLog10(actualPosteriors), 1.0, 1e-3, "Posteriors don't sum to 1 with 1e-3 precision"); } + + @DataProvider(name = "TestIsPolymorphic") + public Object[][] makeTestIsPolymorphic() { + List tests = new ArrayList(); + + final List pValues = new LinkedList(); + for ( final double p : Arrays.asList(0.01, 0.1, 0.9, 0.99, 0.999) ) + for ( final double espilon : Arrays.asList(-1e-5, 0.0, 1e-5) ) + pValues.add(p + espilon); + + for ( final double pNonRef : pValues ) { + for ( final double pThreshold : pValues ) { + final boolean shouldBePoly = pNonRef >= pThreshold; + if ( pNonRef != pThreshold) + // let's not deal with numerical instability + tests.add(new Object[]{ pNonRef, pThreshold, shouldBePoly }); + } + } + + return tests.toArray(new Object[][]{}); + } + + private AFCalcResult makePolymorphicTestData(final double pNonRef) { + return new AFCalcResult( + new int[]{0}, + 1, + alleles, + MathUtils.normalizeFromLog10(new double[]{1 - pNonRef, pNonRef}, true, false), + log10Even, + Collections.singletonMap(C, Math.log10(pNonRef))); + } + + @Test(enabled = true, dataProvider = "TestIsPolymorphic") + private void testIsPolymorphic(final double pNonRef, final double pThreshold, final boolean shouldBePoly) { + final AFCalcResult result = makePolymorphicTestData(pNonRef); + final boolean actualIsPoly = result.isPolymorphic(C, Math.log10(pThreshold)); + Assert.assertEquals(actualIsPoly, shouldBePoly, + "isPolymorphic with pNonRef " + pNonRef + " and threshold " + pThreshold + " returned " + + actualIsPoly + " but the expected result is " + shouldBePoly); + } + + @Test(enabled = true, dataProvider = "TestIsPolymorphic") + private void testIsPolymorphicQual(final double pNonRef, final double pThreshold, final boolean shouldBePoly) { + final AFCalcResult result = makePolymorphicTestData(pNonRef); + final double qual = QualityUtils.phredScaleCorrectRate(pThreshold); + final boolean actualIsPoly = result.isPolymorphicPhredScaledQual(C, qual); + Assert.assertEquals(actualIsPoly, shouldBePoly, + "isPolymorphic with pNonRef " + pNonRef + " and threshold " + pThreshold + " returned " + + actualIsPoly + " but the expected result is " + shouldBePoly); + } } diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java index 97254c478..f22187363 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java @@ -378,11 +378,9 @@ public class UnifiedGenotyperEngine { if ( alternateAllele.isReference() ) continue; - // we are non-ref if the probability of being non-ref > the emit confidence. - // the emit confidence is phred-scaled, say 30 => 10^-3. - // the posterior AF > 0 is log10: -5 => 10^-5 - // we are non-ref if 10^-5 < 10^-3 => -5 < -3 - final boolean isNonRef = AFresult.isPolymorphic(alternateAllele, UAC.STANDARD_CONFIDENCE_FOR_EMITTING / -10.0); + // Compute if the site is considered polymorphic with sufficient confidence relative to our + // phred-scaled emission QUAL + final boolean isNonRef = AFresult.isPolymorphicPhredScaledQual(alternateAllele, UAC.STANDARD_CONFIDENCE_FOR_EMITTING); // if the most likely AC is not 0, then this is a good alternate allele to use if ( isNonRef ) { diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcResult.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcResult.java index a65772444..dbb0e8cdd 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcResult.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcResult.java @@ -28,6 +28,7 @@ package org.broadinstitute.sting.gatk.walkers.genotyper.afcalc; import com.google.java.contract.Ensures; import com.google.java.contract.Requires; import org.broadinstitute.sting.utils.MathUtils; +import org.broadinstitute.sting.utils.QualityUtils; import org.broadinstitute.sting.utils.Utils; import org.broadinstitute.sting.utils.variantcontext.Allele; @@ -234,10 +235,20 @@ public class AFCalcResult { * * @return true if there's enough confidence (relative to log10minPNonRef) to reject AF == 0 */ + @Requires("MathUtils.goodLog10Probability(log10minPNonRef)") public boolean isPolymorphic(final Allele allele, final double log10minPNonRef) { return getLog10PosteriorOfAFGt0ForAllele(allele) >= log10minPNonRef; } + /** + * Same as #isPolymorphic but takes a phred-scaled quality score as input + */ + public boolean isPolymorphicPhredScaledQual(final Allele allele, final double minPNonRefPhredScaledQual) { + if ( minPNonRefPhredScaledQual < 0 ) throw new IllegalArgumentException("phredScaledQual " + minPNonRefPhredScaledQual + " < 0 "); + final double log10Threshold = Math.log10(QualityUtils.qualToProb(minPNonRefPhredScaledQual)); + return isPolymorphic(allele, log10Threshold); + } + /** * Are any of the alleles polymorphic w.r.t. #isPolymorphic? * From 6030605242386c3f04313b797b7d9f6719c67475 Mon Sep 17 00:00:00 2001 From: Eric Banks Date: Wed, 28 Nov 2012 13:26:31 -0500 Subject: [PATCH 17/52] Added quick check for creation of bad BAQ values associated with badly encoded base qualities; hopefully this can help us debug the non-reproducible issue seen by many users. --- .../broadinstitute/sting/utils/QualityUtils.java | 2 +- .../src/org/broadinstitute/sting/utils/baq/BAQ.java | 13 +++++++++---- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/public/java/src/org/broadinstitute/sting/utils/QualityUtils.java b/public/java/src/org/broadinstitute/sting/utils/QualityUtils.java index 1242e5b00..848beccb8 100755 --- a/public/java/src/org/broadinstitute/sting/utils/QualityUtils.java +++ b/public/java/src/org/broadinstitute/sting/utils/QualityUtils.java @@ -9,7 +9,7 @@ import net.sf.samtools.SAMUtils; * @author Kiran Garimella */ public class QualityUtils { - public final static byte MAX_RECALIBRATED_Q_SCORE = 93; + public final static byte MAX_RECALIBRATED_Q_SCORE = SAMUtils.MAX_PHRED_SCORE; public final static byte MAX_QUAL_SCORE = SAMUtils.MAX_PHRED_SCORE; public final static double ERROR_RATE_OF_MAX_QUAL_SCORE = qualToErrorProbRaw(MAX_QUAL_SCORE); diff --git a/public/java/src/org/broadinstitute/sting/utils/baq/BAQ.java b/public/java/src/org/broadinstitute/sting/utils/baq/BAQ.java index cf4d699ee..9ad1bf773 100644 --- a/public/java/src/org/broadinstitute/sting/utils/baq/BAQ.java +++ b/public/java/src/org/broadinstitute/sting/utils/baq/BAQ.java @@ -406,10 +406,15 @@ public class BAQ { // so BQi = Qi - BAQi + 64 byte[] bqTag = new byte[baq.length]; for ( int i = 0; i < bqTag.length; i++) { - int bq = (int)read.getBaseQualities()[i] + 64; - int baq_i = (int)baq[i]; - int tag = bq - baq_i; - if ( tag < 0 ) throw new ReviewedStingException("BAQ tag calculation error. BAQ value above base quality at " + read); + final int bq = (int)read.getBaseQualities()[i] + 64; + final int baq_i = (int)baq[i]; + final int tag = bq - baq_i; + // problem with the calculation of the correction factor; this is our problem + if ( tag < 0 ) + throw new ReviewedStingException("BAQ tag calculation error. BAQ value above base quality at " + read); + // the original quality is too high, almost certainly due to using the wrong encoding in the BAM file + if ( tag > Byte.MAX_VALUE ) + throw new UserException.MalformedBAM(read, "we encountered an extremely high quality score (" + (bq - 64) + ") with BAQ correction factor of " + baq_i + "; the BAM file appears to be using the wrong encoding for quality scores"); bqTag[i] = (byte)tag; } return new String(bqTag); From f0395b457ac4a2c7e40a573f6a46d8e0065b33d4 Mon Sep 17 00:00:00 2001 From: Ryan Poplin Date: Wed, 28 Nov 2012 13:56:32 -0500 Subject: [PATCH 18/52] Adding the work-in-progress, experimental RepeatLengthCovariate to the BQSR so Chris can continue the development. --- .../covariates/RepeatLengthCovariate.java | 64 +++++++++++++++++++ .../variantcontext/VariantContextUtils.java | 2 +- 2 files changed, 65 insertions(+), 1 deletion(-) create mode 100644 public/java/src/org/broadinstitute/sting/utils/recalibration/covariates/RepeatLengthCovariate.java diff --git a/public/java/src/org/broadinstitute/sting/utils/recalibration/covariates/RepeatLengthCovariate.java b/public/java/src/org/broadinstitute/sting/utils/recalibration/covariates/RepeatLengthCovariate.java new file mode 100644 index 000000000..d4e4ab65e --- /dev/null +++ b/public/java/src/org/broadinstitute/sting/utils/recalibration/covariates/RepeatLengthCovariate.java @@ -0,0 +1,64 @@ +package org.broadinstitute.sting.utils.recalibration.covariates; + +import org.broadinstitute.sting.gatk.walkers.bqsr.RecalibrationArgumentCollection; +import org.broadinstitute.sting.gatk.walkers.varianteval.stratifications.TandemRepeat; +import org.broadinstitute.sting.utils.QualityUtils; +import org.broadinstitute.sting.utils.recalibration.ReadCovariates; +import org.broadinstitute.sting.utils.sam.GATKSAMRecord; +import org.broadinstitute.sting.utils.variantcontext.VariantContextUtils; + +import java.util.Arrays; + +/** + * Created with IntelliJ IDEA. + * User: rpoplin + * Date: 11/3/12 + */ + +public class RepeatLengthCovariate implements ExperimentalCovariate { + final int MAX_REPEAT_LENGTH = 20; + + // Initialize any member variables using the command-line arguments passed to the walkers + @Override + public void initialize(final RecalibrationArgumentCollection RAC) {} + + @Override + public void recordValues(final GATKSAMRecord read, final ReadCovariates values) { + byte[] readBytes = read.getReadBases(); + for (int i = 0; i < readBytes.length; i++) { + int maxRL = 0; + for (int str = 1; str <= 8; str++) { + if (i + str <= readBytes.length) { + maxRL = Math.max(maxRL, VariantContextUtils.findNumberofRepetitions( + Arrays.copyOfRange(readBytes,i,i + str), + Arrays.copyOfRange(readBytes,i,readBytes.length) + )); + } + } + if(maxRL > MAX_REPEAT_LENGTH) { maxRL = MAX_REPEAT_LENGTH; } + values.addCovariate(maxRL, maxRL, maxRL, i); + } + } + + // Used to get the covariate's value from input csv file during on-the-fly recalibration + @Override + public final Object getValue(final String str) { + return Byte.parseByte(str); + } + + @Override + public String formatKey(final int key) { + return String.format("%d", key); + } + + @Override + public int keyFromValue(final Object value) { + return (value instanceof String) ? Integer.parseInt((String) value) : (Integer) value; + } + + @Override + public int maximumKeyValue() { + return MAX_REPEAT_LENGTH + 1; + } + +} diff --git a/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextUtils.java b/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextUtils.java index 1f1867f75..b3e3cf8df 100755 --- a/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextUtils.java +++ b/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextUtils.java @@ -1267,7 +1267,7 @@ public class VariantContextUtils { * @param testString String to test * @return Number of repetitions (0 if testString is not a concatenation of n repeatUnit's */ - protected static int findNumberofRepetitions(byte[] repeatUnit, byte[] testString) { + public static int findNumberofRepetitions(byte[] repeatUnit, byte[] testString) { int numRepeats = 0; for (int start = 0; start < testString.length; start += repeatUnit.length) { int end = start + repeatUnit.length; From 198923b597e3635ec5f71aaa192ac4bcaca36ddd Mon Sep 17 00:00:00 2001 From: Joel Thibault Date: Mon, 26 Nov 2012 15:01:13 -0500 Subject: [PATCH 19/52] Add ActiveRegionReadState handling --- .../haplotypecaller/HaplotypeCaller.java | 11 +- .../traversals/TraverseActiveRegions.java | 16 +- .../gatk/walkers/ActiveRegionWalker.java | 20 +- .../activeregion/ActiveRegionReadState.java | 16 ++ .../traversals/TraverseActiveRegionsTest.java | 207 +++++++++++++++--- 5 files changed, 230 insertions(+), 40 deletions(-) create mode 100644 public/java/src/org/broadinstitute/sting/utils/activeregion/ActiveRegionReadState.java diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java index 24b3309f1..d194e2620 100755 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java @@ -45,6 +45,7 @@ import org.broadinstitute.sting.gatk.walkers.genotyper.UnifiedArgumentCollection import org.broadinstitute.sting.gatk.walkers.genotyper.UnifiedGenotyperEngine; import org.broadinstitute.sting.gatk.walkers.genotyper.VariantCallContext; import org.broadinstitute.sting.utils.*; +import org.broadinstitute.sting.utils.activeregion.ActiveRegionReadState; import org.broadinstitute.sting.utils.activeregion.ActivityProfileResult; import org.broadinstitute.sting.utils.clipping.ReadClipper; import org.broadinstitute.sting.utils.codecs.vcf.*; @@ -295,9 +296,15 @@ public class HaplotypeCaller extends ActiveRegionWalker implem @Override public boolean includeReadsWithDeletionAtLoci() { return true; } - // enable non primary reads in the active region + // enable non primary and extended reads in the active region @Override - public boolean wantsNonPrimaryReads() { return true; } + public EnumSet desiredReadStates() { + return EnumSet.of( + ActiveRegionReadState.PRIMARY, + ActiveRegionReadState.NONPRIMARY, + ActiveRegionReadState.EXTENDED + ); + } @Override @Ensures({"result.isActiveProb >= 0.0", "result.isActiveProb <= 1.0"}) diff --git a/public/java/src/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegions.java b/public/java/src/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegions.java index 3f20db0af..06fc01232 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegions.java +++ b/public/java/src/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegions.java @@ -258,13 +258,23 @@ public class TraverseActiveRegions extends TraversalEngine extends Walker desiredReadStates() { + return EnumSet.of(ActiveRegionReadState.PRIMARY); + } + + public final boolean wantsNonPrimaryReads() { + return desiredReadStates().contains(ActiveRegionReadState.NONPRIMARY); + } + + public boolean wantsExtendedReads() { + return desiredReadStates().contains(ActiveRegionReadState.EXTENDED); + } + + public boolean wantsUnmappedReads() { + return desiredReadStates().contains(ActiveRegionReadState.UNMAPPED); } // Determine probability of active status over the AlignmentContext diff --git a/public/java/src/org/broadinstitute/sting/utils/activeregion/ActiveRegionReadState.java b/public/java/src/org/broadinstitute/sting/utils/activeregion/ActiveRegionReadState.java new file mode 100644 index 000000000..00e491eb0 --- /dev/null +++ b/public/java/src/org/broadinstitute/sting/utils/activeregion/ActiveRegionReadState.java @@ -0,0 +1,16 @@ +package org.broadinstitute.sting.utils.activeregion; + +/** + * Created with IntelliJ IDEA. + * User: thibault + * Date: 11/26/12 + * Time: 2:35 PM + * + * Describes how a read relates to an assigned ActiveRegion + */ +public enum ActiveRegionReadState { + PRIMARY, // This is the read's primary region + NONPRIMARY, // This region overlaps the read, but it is not primary + EXTENDED, // This region would overlap the read if it were extended + UNMAPPED // This read is not mapped +} diff --git a/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java b/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java index 66504da11..b70085eff 100644 --- a/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java +++ b/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java @@ -2,6 +2,7 @@ package org.broadinstitute.sting.gatk.traversals; import com.google.java.contract.PreconditionError; import net.sf.samtools.*; +import org.broadinstitute.sting.utils.activeregion.ActiveRegionReadState; import org.broadinstitute.sting.utils.interval.IntervalMergingRule; import org.broadinstitute.sting.utils.interval.IntervalUtils; import org.broadinstitute.sting.utils.sam.GATKSAMRecord; @@ -46,6 +47,8 @@ public class TraverseActiveRegionsTest extends BaseTest { private class DummyActiveRegionWalker extends ActiveRegionWalker { private final double prob; + private EnumSet states = super.desiredReadStates(); + protected List isActiveCalls = new ArrayList(); protected Map mappedActiveRegions = new HashMap(); @@ -57,6 +60,16 @@ public class TraverseActiveRegionsTest extends BaseTest { this.prob = constProb; } + public DummyActiveRegionWalker(EnumSet wantStates) { + this.prob = 1.0; + this.states = wantStates; + } + + @Override + public EnumSet desiredReadStates() { + return states; + } + @Override public ActivityProfileResult isActive(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) { isActiveCalls.add(ref.getLocus()); @@ -202,12 +215,158 @@ public class TraverseActiveRegionsTest extends BaseTest { } @Test - public void testReadMapping() { + public void testPrimaryReadMapping() { DummyActiveRegionWalker walker = new DummyActiveRegionWalker(); // Contract: Each read has the Primary state in a single region (or none) // This is the region of maximum overlap for the read (earlier if tied) + // simple: Primary in 1:1-999 + // overlap_equal: Primary in 1:1-999 + // overlap_unequal: Primary in 1:1-999 + // boundary_equal: Non-Primary in 1:1000-1999, Primary in 1:2000-2999 + // boundary_unequal: Primary in 1:1000-1999, Non-Primary in 1:2000-2999 + // extended_only: Extended in 1:2000-2999 + // extended_and_np: Non-Primary in 1:1-999, Primary in 1:1000-1999, Extended in 1:2000-2999 + // outside_intervals: none + // simple20: Primary in 20:10000-10100 + + Map activeRegions = getActiveRegions(walker, intervals); + ActiveRegion region; + + region = activeRegions.get(genomeLocParser.createGenomeLoc("1", 1, 999)); + + getRead(region, "simple"); + getRead(region, "overlap_equal"); + getRead(region, "overlap_unequal"); + verifyReadNotPlaced(region, "boundary_equal"); + verifyReadNotPlaced(region, "boundary_unequal"); + verifyReadNotPlaced(region, "extended_only"); + verifyReadNotPlaced(region, "extended_and_np"); + verifyReadNotPlaced(region, "outside_intervals"); + verifyReadNotPlaced(region, "simple20"); + + region = activeRegions.get(genomeLocParser.createGenomeLoc("1", 1000, 1999)); + + verifyReadNotPlaced(region, "simple"); + verifyReadNotPlaced(region, "overlap_equal"); + verifyReadNotPlaced(region, "overlap_unequal"); + verifyReadNotPlaced(region, "boundary_equal"); + getRead(region, "boundary_unequal"); + verifyReadNotPlaced(region, "extended_only"); + // TODO: fail getRead(region, "extended_and_np"); + verifyReadNotPlaced(region, "outside_intervals"); + verifyReadNotPlaced(region, "simple20"); + + region = activeRegions.get(genomeLocParser.createGenomeLoc("1", 2000, 2999)); + + verifyReadNotPlaced(region, "simple"); + verifyReadNotPlaced(region, "overlap_equal"); + verifyReadNotPlaced(region, "overlap_unequal"); + getRead(region, "boundary_equal"); + verifyReadNotPlaced(region, "boundary_unequal"); + verifyReadNotPlaced(region, "extended_only"); + verifyReadNotPlaced(region, "extended_and_np"); + verifyReadNotPlaced(region, "outside_intervals"); + verifyReadNotPlaced(region, "simple20"); + + region = activeRegions.get(genomeLocParser.createGenomeLoc("20", 10000, 10100)); + + verifyReadNotPlaced(region, "simple"); + verifyReadNotPlaced(region, "overlap_equal"); + verifyReadNotPlaced(region, "overlap_unequal"); + verifyReadNotPlaced(region, "boundary_equal"); + verifyReadNotPlaced(region, "boundary_unequal"); + verifyReadNotPlaced(region, "extended_only"); + verifyReadNotPlaced(region, "extended_and_np"); + verifyReadNotPlaced(region, "outside_intervals"); + getRead(region, "simple20"); + + // TODO: more tests and edge cases + } + + @Test + public void testNonPrimaryReadMapping() { + DummyActiveRegionWalker walker = new DummyActiveRegionWalker( + EnumSet.of(ActiveRegionReadState.PRIMARY, ActiveRegionReadState.NONPRIMARY)); + + // Contract: Each read has the Primary state in a single region (or none) + // This is the region of maximum overlap for the read (earlier if tied) + + // Contract: Each read has the Non-Primary state in all other regions it overlaps + + // simple: Primary in 1:1-999 + // overlap_equal: Primary in 1:1-999 + // overlap_unequal: Primary in 1:1-999 + // boundary_equal: Non-Primary in 1:1000-1999, Primary in 1:2000-2999 + // boundary_unequal: Primary in 1:1000-1999, Non-Primary in 1:2000-2999 + // extended_only: Extended in 1:2000-2999 + // extended_and_np: Non-Primary in 1:1-999, Primary in 1:1000-1999, Extended in 1:2000-2999 + // outside_intervals: none + // simple20: Primary in 20:10000-10100 + + Map activeRegions = getActiveRegions(walker, intervals); + ActiveRegion region; + + region = activeRegions.get(genomeLocParser.createGenomeLoc("1", 1, 999)); + + getRead(region, "simple"); + getRead(region, "overlap_equal"); + getRead(region, "overlap_unequal"); + verifyReadNotPlaced(region, "boundary_equal"); + verifyReadNotPlaced(region, "boundary_unequal"); + verifyReadNotPlaced(region, "extended_only"); + // TODO: fail getRead(region, "extended_and_np"); + verifyReadNotPlaced(region, "outside_intervals"); + verifyReadNotPlaced(region, "simple20"); + + region = activeRegions.get(genomeLocParser.createGenomeLoc("1", 1000, 1999)); + + verifyReadNotPlaced(region, "simple"); + verifyReadNotPlaced(region, "overlap_equal"); + verifyReadNotPlaced(region, "overlap_unequal"); + getRead(region, "boundary_equal"); + getRead(region, "boundary_unequal"); + verifyReadNotPlaced(region, "extended_only"); + // TODO: fail getRead(region, "extended_and_np"); + verifyReadNotPlaced(region, "outside_intervals"); + verifyReadNotPlaced(region, "simple20"); + + region = activeRegions.get(genomeLocParser.createGenomeLoc("1", 2000, 2999)); + + verifyReadNotPlaced(region, "simple"); + verifyReadNotPlaced(region, "overlap_equal"); + verifyReadNotPlaced(region, "overlap_unequal"); + getRead(region, "boundary_equal"); + getRead(region, "boundary_unequal"); + verifyReadNotPlaced(region, "extended_only"); + verifyReadNotPlaced(region, "extended_and_np"); + verifyReadNotPlaced(region, "outside_intervals"); + verifyReadNotPlaced(region, "simple20"); + + region = activeRegions.get(genomeLocParser.createGenomeLoc("20", 10000, 10100)); + + verifyReadNotPlaced(region, "simple"); + verifyReadNotPlaced(region, "overlap_equal"); + verifyReadNotPlaced(region, "overlap_unequal"); + verifyReadNotPlaced(region, "boundary_equal"); + verifyReadNotPlaced(region, "boundary_unequal"); + verifyReadNotPlaced(region, "extended_only"); + verifyReadNotPlaced(region, "extended_and_np"); + verifyReadNotPlaced(region, "outside_intervals"); + getRead(region, "simple20"); + + // TODO: more tests and edge cases + } + + @Test + public void testExtendedReadMapping() { + DummyActiveRegionWalker walker = new DummyActiveRegionWalker( + EnumSet.of(ActiveRegionReadState.PRIMARY, ActiveRegionReadState.NONPRIMARY, ActiveRegionReadState.EXTENDED)); + + // Contract: Each read has the Primary state in a single region (or none) + // This is the region of maximum overlap for the read (earlier if tied) + // Contract: Each read has the Non-Primary state in all other regions it overlaps // Contract: Each read has the Extended state in regions where it only overlaps if the region is extended @@ -226,13 +385,13 @@ public class TraverseActiveRegionsTest extends BaseTest { region = activeRegions.get(genomeLocParser.createGenomeLoc("1", 1, 999)); - verifyReadPrimary(region, "simple"); - verifyReadPrimary(region, "overlap_equal"); - verifyReadPrimary(region, "overlap_unequal"); + getRead(region, "simple"); + getRead(region, "overlap_equal"); + getRead(region, "overlap_unequal"); verifyReadNotPlaced(region, "boundary_equal"); verifyReadNotPlaced(region, "boundary_unequal"); verifyReadNotPlaced(region, "extended_only"); - // TODO: fail verifyReadNonPrimary(region, "extended_and_np"); + // TODO: fail getRead(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); verifyReadNotPlaced(region, "simple20"); @@ -241,10 +400,10 @@ public class TraverseActiveRegionsTest extends BaseTest { verifyReadNotPlaced(region, "simple"); verifyReadNotPlaced(region, "overlap_equal"); verifyReadNotPlaced(region, "overlap_unequal"); - // TODO: fail verifyReadNonPrimary(region, "boundary_equal"); - verifyReadPrimary(region, "boundary_unequal"); + getRead(region, "boundary_equal"); + getRead(region, "boundary_unequal"); verifyReadNotPlaced(region, "extended_only"); - // TODO: fail verifyReadPrimary(region, "extended_and_np"); + // TODO: fail getRead(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); verifyReadNotPlaced(region, "simple20"); @@ -253,10 +412,10 @@ public class TraverseActiveRegionsTest extends BaseTest { verifyReadNotPlaced(region, "simple"); verifyReadNotPlaced(region, "overlap_equal"); verifyReadNotPlaced(region, "overlap_unequal"); - verifyReadPrimary(region, "boundary_equal"); - // TODO: fail verifyReadNonPrimary(region, "boundary_unequal"); - // TODO: fail verifyReadExtended(region, "extended_only"); - // TODO: fail verifyReadExtended(region, "extended_and_np"); + getRead(region, "boundary_equal"); + getRead(region, "boundary_unequal"); + verifyReadNotPlaced(region, "extended_only"); + verifyReadNotPlaced(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); verifyReadNotPlaced(region, "simple20"); @@ -267,33 +426,19 @@ public class TraverseActiveRegionsTest extends BaseTest { verifyReadNotPlaced(region, "overlap_unequal"); verifyReadNotPlaced(region, "boundary_equal"); verifyReadNotPlaced(region, "boundary_unequal"); - verifyReadNotPlaced(region, "extended_only"); - verifyReadNotPlaced(region, "extended_and_np"); + // TODO: fail getRead(region, "extended_only"); + // TODO: fail getRead(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); - verifyReadPrimary(region, "simple20"); + getRead(region, "simple20"); // TODO: more tests and edge cases } - private void verifyReadPrimary(ActiveRegion region, String readName) { - SAMRecord read = getRead(region, readName); - Assert.assertFalse(read.getNotPrimaryAlignmentFlag(), "Read " + read + " not primary in active region " + region); - } - - private void verifyReadNonPrimary(ActiveRegion region, String readName) { - SAMRecord read = getRead(region, readName); - Assert.assertTrue(read.getNotPrimaryAlignmentFlag(), "Read " + read + " primary in active region " + region); - } - - private void verifyReadExtended(ActiveRegion region, String readName) { - Assert.fail("The Extended read state has not been implemented"); - } - private void verifyReadNotPlaced(ActiveRegion region, String readName) { for (SAMRecord read : region.getReads()) { if (read.getReadName().equals(readName)) Assert.fail("Read " + readName + " found in active region " + region); - } + } } private SAMRecord getRead(ActiveRegion region, String readName) { @@ -302,7 +447,7 @@ public class TraverseActiveRegionsTest extends BaseTest { return read; } - Assert.fail("Read " + readName + " not found in active region " + region); + Assert.fail("Read " + readName + " not assigned to active region " + region); return null; } From c76c808268596f12d81ca247c826653f4ffa2f56 Mon Sep 17 00:00:00 2001 From: Joel Thibault Date: Wed, 28 Nov 2012 11:09:12 -0500 Subject: [PATCH 20/52] Reads are required to be sorted - Remove the extended_only case because it's outside intervals --- .../traversals/TraverseActiveRegionsTest.java | 42 +++++++------------ 1 file changed, 15 insertions(+), 27 deletions(-) diff --git a/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java b/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java index b70085eff..a65b0cb45 100644 --- a/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java +++ b/public/java/test/org/broadinstitute/sting/gatk/traversals/TraverseActiveRegionsTest.java @@ -25,6 +25,7 @@ import org.broadinstitute.sting.utils.activeregion.ActiveRegion; import org.broadinstitute.sting.utils.activeregion.ActivityProfileResult; import org.broadinstitute.sting.utils.fasta.CachingIndexedFastaSequenceFile; import org.broadinstitute.sting.utils.sam.ArtificialSAMUtils; +import org.broadinstitute.sting.utils.sam.ReadUtils; import org.testng.Assert; import org.testng.annotations.BeforeClass; import org.testng.annotations.Test; @@ -100,7 +101,7 @@ public class TraverseActiveRegionsTest extends BaseTest { private GenomeLocParser genomeLocParser; private List intervals; - private List reads; + private List reads; @BeforeClass private void init() throws FileNotFoundException { @@ -117,16 +118,18 @@ public class TraverseActiveRegionsTest extends BaseTest { intervals.add(genomeLocParser.createGenomeLoc("20", 10000, 10100)); intervals = IntervalUtils.sortAndMergeIntervals(genomeLocParser, intervals, IntervalMergingRule.OVERLAPPING_ONLY).toList(); - reads = new ArrayList(); + reads = new ArrayList(); reads.add(buildSAMRecord("simple", "1", 100, 200)); reads.add(buildSAMRecord("overlap_equal", "1", 10, 20)); reads.add(buildSAMRecord("overlap_unequal", "1", 10, 21)); reads.add(buildSAMRecord("boundary_equal", "1", 1990, 2009)); reads.add(buildSAMRecord("boundary_unequal", "1", 1990, 2008)); - reads.add(buildSAMRecord("extended_only", "1", 3000, 3100)); reads.add(buildSAMRecord("extended_and_np", "1", 990, 1990)); reads.add(buildSAMRecord("outside_intervals", "1", 5000, 6000)); reads.add(buildSAMRecord("simple20", "20", 10025, 10075)); + + // required by LocusIteratorByState, and I prefer to list them in test case order above + ReadUtils.sortReadsByCoordinate(reads); } @Test @@ -226,7 +229,6 @@ public class TraverseActiveRegionsTest extends BaseTest { // overlap_unequal: Primary in 1:1-999 // boundary_equal: Non-Primary in 1:1000-1999, Primary in 1:2000-2999 // boundary_unequal: Primary in 1:1000-1999, Non-Primary in 1:2000-2999 - // extended_only: Extended in 1:2000-2999 // extended_and_np: Non-Primary in 1:1-999, Primary in 1:1000-1999, Extended in 1:2000-2999 // outside_intervals: none // simple20: Primary in 20:10000-10100 @@ -241,7 +243,6 @@ public class TraverseActiveRegionsTest extends BaseTest { getRead(region, "overlap_unequal"); verifyReadNotPlaced(region, "boundary_equal"); verifyReadNotPlaced(region, "boundary_unequal"); - verifyReadNotPlaced(region, "extended_only"); verifyReadNotPlaced(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); verifyReadNotPlaced(region, "simple20"); @@ -253,8 +254,7 @@ public class TraverseActiveRegionsTest extends BaseTest { verifyReadNotPlaced(region, "overlap_unequal"); verifyReadNotPlaced(region, "boundary_equal"); getRead(region, "boundary_unequal"); - verifyReadNotPlaced(region, "extended_only"); - // TODO: fail getRead(region, "extended_and_np"); + getRead(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); verifyReadNotPlaced(region, "simple20"); @@ -265,7 +265,6 @@ public class TraverseActiveRegionsTest extends BaseTest { verifyReadNotPlaced(region, "overlap_unequal"); getRead(region, "boundary_equal"); verifyReadNotPlaced(region, "boundary_unequal"); - verifyReadNotPlaced(region, "extended_only"); verifyReadNotPlaced(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); verifyReadNotPlaced(region, "simple20"); @@ -277,7 +276,6 @@ public class TraverseActiveRegionsTest extends BaseTest { verifyReadNotPlaced(region, "overlap_unequal"); verifyReadNotPlaced(region, "boundary_equal"); verifyReadNotPlaced(region, "boundary_unequal"); - verifyReadNotPlaced(region, "extended_only"); verifyReadNotPlaced(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); getRead(region, "simple20"); @@ -300,7 +298,6 @@ public class TraverseActiveRegionsTest extends BaseTest { // overlap_unequal: Primary in 1:1-999 // boundary_equal: Non-Primary in 1:1000-1999, Primary in 1:2000-2999 // boundary_unequal: Primary in 1:1000-1999, Non-Primary in 1:2000-2999 - // extended_only: Extended in 1:2000-2999 // extended_and_np: Non-Primary in 1:1-999, Primary in 1:1000-1999, Extended in 1:2000-2999 // outside_intervals: none // simple20: Primary in 20:10000-10100 @@ -315,8 +312,7 @@ public class TraverseActiveRegionsTest extends BaseTest { getRead(region, "overlap_unequal"); verifyReadNotPlaced(region, "boundary_equal"); verifyReadNotPlaced(region, "boundary_unequal"); - verifyReadNotPlaced(region, "extended_only"); - // TODO: fail getRead(region, "extended_and_np"); + getRead(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); verifyReadNotPlaced(region, "simple20"); @@ -327,8 +323,7 @@ public class TraverseActiveRegionsTest extends BaseTest { verifyReadNotPlaced(region, "overlap_unequal"); getRead(region, "boundary_equal"); getRead(region, "boundary_unequal"); - verifyReadNotPlaced(region, "extended_only"); - // TODO: fail getRead(region, "extended_and_np"); + getRead(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); verifyReadNotPlaced(region, "simple20"); @@ -339,7 +334,6 @@ public class TraverseActiveRegionsTest extends BaseTest { verifyReadNotPlaced(region, "overlap_unequal"); getRead(region, "boundary_equal"); getRead(region, "boundary_unequal"); - verifyReadNotPlaced(region, "extended_only"); verifyReadNotPlaced(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); verifyReadNotPlaced(region, "simple20"); @@ -351,7 +345,6 @@ public class TraverseActiveRegionsTest extends BaseTest { verifyReadNotPlaced(region, "overlap_unequal"); verifyReadNotPlaced(region, "boundary_equal"); verifyReadNotPlaced(region, "boundary_unequal"); - verifyReadNotPlaced(region, "extended_only"); verifyReadNotPlaced(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); getRead(region, "simple20"); @@ -375,7 +368,6 @@ public class TraverseActiveRegionsTest extends BaseTest { // overlap_unequal: Primary in 1:1-999 // boundary_equal: Non-Primary in 1:1000-1999, Primary in 1:2000-2999 // boundary_unequal: Primary in 1:1000-1999, Non-Primary in 1:2000-2999 - // extended_only: Extended in 1:2000-2999 // extended_and_np: Non-Primary in 1:1-999, Primary in 1:1000-1999, Extended in 1:2000-2999 // outside_intervals: none // simple20: Primary in 20:10000-10100 @@ -390,8 +382,7 @@ public class TraverseActiveRegionsTest extends BaseTest { getRead(region, "overlap_unequal"); verifyReadNotPlaced(region, "boundary_equal"); verifyReadNotPlaced(region, "boundary_unequal"); - verifyReadNotPlaced(region, "extended_only"); - // TODO: fail getRead(region, "extended_and_np"); + getRead(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); verifyReadNotPlaced(region, "simple20"); @@ -402,8 +393,7 @@ public class TraverseActiveRegionsTest extends BaseTest { verifyReadNotPlaced(region, "overlap_unequal"); getRead(region, "boundary_equal"); getRead(region, "boundary_unequal"); - verifyReadNotPlaced(region, "extended_only"); - // TODO: fail getRead(region, "extended_and_np"); + getRead(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); verifyReadNotPlaced(region, "simple20"); @@ -414,8 +404,7 @@ public class TraverseActiveRegionsTest extends BaseTest { verifyReadNotPlaced(region, "overlap_unequal"); getRead(region, "boundary_equal"); getRead(region, "boundary_unequal"); - verifyReadNotPlaced(region, "extended_only"); - verifyReadNotPlaced(region, "extended_and_np"); + getRead(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); verifyReadNotPlaced(region, "simple20"); @@ -426,8 +415,7 @@ public class TraverseActiveRegionsTest extends BaseTest { verifyReadNotPlaced(region, "overlap_unequal"); verifyReadNotPlaced(region, "boundary_equal"); verifyReadNotPlaced(region, "boundary_unequal"); - // TODO: fail getRead(region, "extended_only"); - // TODO: fail getRead(region, "extended_and_np"); + verifyReadNotPlaced(region, "extended_and_np"); verifyReadNotPlaced(region, "outside_intervals"); getRead(region, "simple20"); @@ -438,7 +426,7 @@ public class TraverseActiveRegionsTest extends BaseTest { for (SAMRecord read : region.getReads()) { if (read.getReadName().equals(readName)) Assert.fail("Read " + readName + " found in active region " + region); - } + } } private SAMRecord getRead(ActiveRegion region, String readName) { @@ -520,7 +508,7 @@ public class TraverseActiveRegionsTest extends BaseTest { engine.setGenomeLocParser(genomeLocParser); t.initialize(engine); - StingSAMIterator iterator = ArtificialSAMUtils.createReadIterator(reads); + StingSAMIterator iterator = ArtificialSAMUtils.createReadIterator(new ArrayList(reads)); Shard shard = new MockLocusShard(genomeLocParser, intervals); List providers = new ArrayList(); From 26d9c41615ccd502b75f23a42110af925995beba Mon Sep 17 00:00:00 2001 From: David Roazen Date: Wed, 28 Nov 2012 14:06:58 -0500 Subject: [PATCH 21/52] Allow arbitrary resources to be packaged in the GATK jar, selecting among public/private/protected appropriately -Resources must be in a "resources" or "templates" subdirectory within the Java package hierarchy -Remove direct inclusion of private resources from the main jar packaging target added in Jacob's patch: this would break builds where the private directory was absent, and did not respect build settings (include.private, etc.) --- build.xml | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/build.xml b/build.xml index 3a264b476..4db71a9ab 100644 --- a/build.xml +++ b/build.xml @@ -226,10 +226,17 @@ - - - - + + + + + + + + + + + @@ -681,9 +688,6 @@ - - - From b2e699169cf1c545e92db2536be2cab656c472fe Mon Sep 17 00:00:00 2001 From: David Roazen Date: Wed, 28 Nov 2012 15:26:05 -0500 Subject: [PATCH 22/52] Update GATK packaging settings to package arbitrary resources With the newly-added support for packaging arbitrary resources, the resources were getting packaged in a normal build but not when creating a standalone GATK jar. This corrects this oversight. --- public/packages/GATKEngine.xml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/public/packages/GATKEngine.xml b/public/packages/GATKEngine.xml index 2de0273f3..d0b4a52b5 100644 --- a/public/packages/GATKEngine.xml +++ b/public/packages/GATKEngine.xml @@ -36,6 +36,9 @@

+ + + From b06e71cedf057e3848a46a9002cf83c496d6b8ef Mon Sep 17 00:00:00 2001 From: David Roazen Date: Wed, 28 Nov 2012 20:44:09 -0500 Subject: [PATCH 30/52] Use build jars in test classpaths by default -Allows packaged resource files to be accessed within tests -Guards against packaging errors in dist/ jars by testing the jars that actually get run rather than unpackaged class files. Previously we were only protected against packaging errors in the monolithic jars posted to our website, not the dist/ jars used in everyday runs. -"ant fasttest" still uses the unpackaged class files for speed (don't want to have to rebuild the jars in fasttest). Relies on dubious methods to get at the resource files that would end up in the jars. -Eliminated the stupid separate "test" ivy config. Now we only invoke ivy ONCE during an ant build that includes tests. --- build.xml | 64 +++++++++++++++++++++++++++---------------------------- ivy.xml | 11 ++++------ 2 files changed, 36 insertions(+), 39 deletions(-) diff --git a/build.xml b/build.xml index 4db71a9ab..cc45467d8 100644 --- a/build.xml +++ b/build.xml @@ -185,10 +185,7 @@ - - - - + @@ -205,6 +202,16 @@ + + + + + + + + + + @@ -240,13 +247,6 @@ - - - - - - - @@ -1110,15 +1110,10 @@ - - + - - - - @@ -1126,9 +1121,6 @@ - - - @@ -1136,10 +1128,8 @@ - - + - @@ -1150,11 +1140,9 @@ - + - - @@ -1167,9 +1155,8 @@ - + - @@ -1376,14 +1363,13 @@ - + - @@ -1401,13 +1387,27 @@ - - + + + + + + + + + + + + + + + + diff --git a/ivy.xml b/ivy.xml index 1d2f95dc1..b7ca65406 100644 --- a/ivy.xml +++ b/ivy.xml @@ -24,11 +24,8 @@ - + - - - @@ -83,9 +80,9 @@ - - - + + + From f837e6ced7ac1fc187eb98f3f2c707a18e7f4a8c Mon Sep 17 00:00:00 2001 From: Mark DePristo Date: Thu, 29 Nov 2012 14:38:09 -0500 Subject: [PATCH 34/52] Refactored entire NA12878KB to allow us to easily build a na12878kb.jar for IGV integration -- Just separated infrastructure into core package, away from the walkers themselves. -- Added na12878kb.jar target that builds a jar that can run a test main function (see testNA12878kbJar.csh) --- build.xml | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/build.xml b/build.xml index cc45467d8..834aef3cd 100644 --- a/build.xml +++ b/build.xml @@ -679,6 +679,24 @@ + + + + + + + + + + + + + + + + + + From daf6269b6503055b5f1a932ec4debb0a43a41bd3 Mon Sep 17 00:00:00 2001 From: Johan Dahlberg Date: Mon, 1 Oct 2012 11:28:46 +0200 Subject: [PATCH 35/52] Setting the walltime Signed-off-by: Joel Thibault --- .../src/org/broadinstitute/sting/queue/QSettings.scala | 4 ++++ .../sting/queue/engine/drmaa/DrmaaJobRunner.scala | 3 +++ .../sting/queue/function/CommandLineFunction.scala | 10 ++++++++++ 3 files changed, 17 insertions(+) diff --git a/public/scala/src/org/broadinstitute/sting/queue/QSettings.scala b/public/scala/src/org/broadinstitute/sting/queue/QSettings.scala index 2c0f43bac..fb21700ac 100644 --- a/public/scala/src/org/broadinstitute/sting/queue/QSettings.scala +++ b/public/scala/src/org/broadinstitute/sting/queue/QSettings.scala @@ -31,6 +31,10 @@ import org.broadinstitute.sting.commandline.{ClassType, Argument} * Default settings settable on the command line and passed to CommandLineFunctions. */ class QSettings { + + @Argument(fullName="job_walltime", shortName="wallTime", doc="Setting the required walltime when using the drmaa job runner.", required=false) + var jobWalltime: Option[Long] = None + @Argument(fullName="run_name", shortName="runName", doc="A name for this run used for various status messages.", required=false) var runName: String = _ diff --git a/public/scala/src/org/broadinstitute/sting/queue/engine/drmaa/DrmaaJobRunner.scala b/public/scala/src/org/broadinstitute/sting/queue/engine/drmaa/DrmaaJobRunner.scala index 2aae2fc6b..31b314c79 100644 --- a/public/scala/src/org/broadinstitute/sting/queue/engine/drmaa/DrmaaJobRunner.scala +++ b/public/scala/src/org/broadinstitute/sting/queue/engine/drmaa/DrmaaJobRunner.scala @@ -65,6 +65,9 @@ class DrmaaJobRunner(val session: Session, val function: CommandLineFunction) ex drmaaJob.setJoinFiles(true) } + if(function.wallTime != null) + drmaaJob.setHardWallclockTimeLimit(function.wallTime.get) + drmaaJob.setNativeSpecification(functionNativeSpec) // Instead of running the function.commandLine, run "sh " diff --git a/public/scala/src/org/broadinstitute/sting/queue/function/CommandLineFunction.scala b/public/scala/src/org/broadinstitute/sting/queue/function/CommandLineFunction.scala index eb426d301..d5870a6c3 100644 --- a/public/scala/src/org/broadinstitute/sting/queue/function/CommandLineFunction.scala +++ b/public/scala/src/org/broadinstitute/sting/queue/function/CommandLineFunction.scala @@ -33,6 +33,9 @@ import org.broadinstitute.sting.commandline.Argument trait CommandLineFunction extends QFunction with Logging { def commandLine: String + /** Setting the wall time request for drmaa job*/ + var wallTime: Option[Long] = None + /** Upper memory limit */ @Argument(doc="Memory limit", required=false) var memoryLimit: Option[Double] = None @@ -67,6 +70,9 @@ trait CommandLineFunction extends QFunction with Logging { super.copySettingsTo(function) function match { case commandLineFunction: CommandLineFunction => + if(commandLineFunction.wallTime.isEmpty) + commandLineFunction.wallTime = this.wallTime + if (commandLineFunction.memoryLimit.isEmpty) commandLineFunction.memoryLimit = this.memoryLimit @@ -110,6 +116,10 @@ trait CommandLineFunction extends QFunction with Logging { * Sets all field values. */ override def freezeFieldValues() { + + if(wallTime.isEmpty) + wallTime = qSettings.jobWalltime + if (jobQueue == null) jobQueue = qSettings.jobQueue From 97d29f203e35fd98393f61a28e78134da2b84755 Mon Sep 17 00:00:00 2001 From: Joel Thibault Date: Fri, 12 Oct 2012 17:16:56 -0400 Subject: [PATCH 36/52] Add walltime changes to LSF - Check whether the specified attribute is available - Add pipeline test (disabled due to missing attribute) --- .../sting/jna/drmaa/v1_0/JnaSession.java | 18 ++++++++++++++---- .../broadinstitute/sting/queue/QSettings.scala | 8 ++++---- .../queue/engine/drmaa/DrmaaJobRunner.scala | 2 +- .../queue/engine/lsf/Lsf706JobRunner.scala | 7 +++++-- .../queue/function/CommandLineFunction.scala | 2 +- .../examples/HelloWorldPipelineTest.scala | 11 +++++++++++ 6 files changed, 36 insertions(+), 12 deletions(-) diff --git a/public/java/src/org/broadinstitute/sting/jna/drmaa/v1_0/JnaSession.java b/public/java/src/org/broadinstitute/sting/jna/drmaa/v1_0/JnaSession.java index 480113e1e..830c6590d 100644 --- a/public/java/src/org/broadinstitute/sting/jna/drmaa/v1_0/JnaSession.java +++ b/public/java/src/org/broadinstitute/sting/jna/drmaa/v1_0/JnaSession.java @@ -210,13 +210,23 @@ public class JnaSession implements Session { } public static void setAttribute(Pointer jt, String name, String value) throws DrmaaException { - checkError(LibDrmaa.drmaa_set_attribute(jt, name, value, getError(), LibDrmaa.DRMAA_ERROR_STRING_BUFFER_LEN)); + if (getAttrNames().contains(name)) { + checkError(LibDrmaa.drmaa_set_attribute(jt, name, value, getError(), LibDrmaa.DRMAA_ERROR_STRING_BUFFER_LEN)); + } + else { + throw new InvalidAttributeValueException("Attribute " + name + " is not supported by this implementation of DRMAA"); + } } public static String getAttribute(Pointer jt, String name) throws DrmaaException { - Memory attrBuffer = new Memory(LibDrmaa.DRMAA_ATTR_BUFFER); - checkError(LibDrmaa.drmaa_get_attribute(jt, name, attrBuffer, LibDrmaa.DRMAA_ATTR_BUFFER_LEN, getError(), LibDrmaa.DRMAA_ERROR_STRING_BUFFER_LEN)); - return attrBuffer.getString(0); + if (getAttrNames().contains(name)) { + Memory attrBuffer = new Memory(LibDrmaa.DRMAA_ATTR_BUFFER); + checkError(LibDrmaa.drmaa_get_attribute(jt, name, attrBuffer, LibDrmaa.DRMAA_ATTR_BUFFER_LEN, getError(), LibDrmaa.DRMAA_ERROR_STRING_BUFFER_LEN)); + return attrBuffer.getString(0); + } + else { + throw new InvalidAttributeValueException("Attribute " + name + " is not supported by this implementation of DRMAA"); + } } public static void setVectorAttribute(Pointer jt, String name, Collection values) throws DrmaaException { diff --git a/public/scala/src/org/broadinstitute/sting/queue/QSettings.scala b/public/scala/src/org/broadinstitute/sting/queue/QSettings.scala index fb21700ac..b1e98a0e2 100644 --- a/public/scala/src/org/broadinstitute/sting/queue/QSettings.scala +++ b/public/scala/src/org/broadinstitute/sting/queue/QSettings.scala @@ -31,10 +31,6 @@ import org.broadinstitute.sting.commandline.{ClassType, Argument} * Default settings settable on the command line and passed to CommandLineFunctions. */ class QSettings { - - @Argument(fullName="job_walltime", shortName="wallTime", doc="Setting the required walltime when using the drmaa job runner.", required=false) - var jobWalltime: Option[Long] = None - @Argument(fullName="run_name", shortName="runName", doc="A name for this run used for various status messages.", required=false) var runName: String = _ @@ -76,6 +72,10 @@ class QSettings { @Argument(fullName="resident_memory_request_parameter", shortName="resMemReqParam", doc="Parameter for resident memory requests. By default not requested.", required=false) var residentRequestParameter: String = _ + @Argument(fullName="job_walltime", shortName="wallTime", doc="Setting the required DRMAA walltime or LSF run limit.", required=false) + @ClassType(classOf[Long]) + var jobWalltime: Option[Long] = None + /** The name of the parallel environment (required for SGE, for example) */ @Argument(fullName="job_parallel_env", shortName="jobParaEnv", doc="An SGE style parallel environment to use for jobs requesting more than 1 core. Equivalent to submitting jobs with -pe ARG nt for jobs with nt > 1", required=false) var parallelEnvironmentName: String = "smp_pe" // Broad default diff --git a/public/scala/src/org/broadinstitute/sting/queue/engine/drmaa/DrmaaJobRunner.scala b/public/scala/src/org/broadinstitute/sting/queue/engine/drmaa/DrmaaJobRunner.scala index 31b314c79..1dca22981 100644 --- a/public/scala/src/org/broadinstitute/sting/queue/engine/drmaa/DrmaaJobRunner.scala +++ b/public/scala/src/org/broadinstitute/sting/queue/engine/drmaa/DrmaaJobRunner.scala @@ -65,7 +65,7 @@ class DrmaaJobRunner(val session: Session, val function: CommandLineFunction) ex drmaaJob.setJoinFiles(true) } - if(function.wallTime != null) + if(!function.wallTime.isEmpty) drmaaJob.setHardWallclockTimeLimit(function.wallTime.get) drmaaJob.setNativeSpecification(functionNativeSpec) diff --git a/public/scala/src/org/broadinstitute/sting/queue/engine/lsf/Lsf706JobRunner.scala b/public/scala/src/org/broadinstitute/sting/queue/engine/lsf/Lsf706JobRunner.scala index 2fbea1497..5dc126e49 100644 --- a/public/scala/src/org/broadinstitute/sting/queue/engine/lsf/Lsf706JobRunner.scala +++ b/public/scala/src/org/broadinstitute/sting/queue/engine/lsf/Lsf706JobRunner.scala @@ -151,8 +151,11 @@ class Lsf706JobRunner(val function: CommandLineFunction) extends CommandLineJobR throw new QException("setOption_() returned -1 while setting esub"); } - // LSF specific: get the max runtime for the jobQueue and pass it for this job - request.rLimits(LibLsf.LSF_RLIMIT_RUN) = Lsf706JobRunner.getRlimitRun(function.jobQueue) + if(!function.wallTime.isEmpty) + request.rLimits(LibLsf.LSF_RLIMIT_RUN) = function.wallTime.get.toInt + else + // LSF specific: get the max runtime for the jobQueue and pass it for this job + request.rLimits(LibLsf.LSF_RLIMIT_RUN) = Lsf706JobRunner.getRlimitRun(function.jobQueue) // Run the command as sh request.command = "sh " + jobScript diff --git a/public/scala/src/org/broadinstitute/sting/queue/function/CommandLineFunction.scala b/public/scala/src/org/broadinstitute/sting/queue/function/CommandLineFunction.scala index d5870a6c3..2453cc50a 100644 --- a/public/scala/src/org/broadinstitute/sting/queue/function/CommandLineFunction.scala +++ b/public/scala/src/org/broadinstitute/sting/queue/function/CommandLineFunction.scala @@ -33,7 +33,7 @@ import org.broadinstitute.sting.commandline.Argument trait CommandLineFunction extends QFunction with Logging { def commandLine: String - /** Setting the wall time request for drmaa job*/ + /** Setting the wall time request for DRMAA / run limit for LSF */ var wallTime: Option[Long] = None /** Upper memory limit */ diff --git a/public/scala/test/org/broadinstitute/sting/queue/pipeline/examples/HelloWorldPipelineTest.scala b/public/scala/test/org/broadinstitute/sting/queue/pipeline/examples/HelloWorldPipelineTest.scala index 50fc529dd..c8085784d 100644 --- a/public/scala/test/org/broadinstitute/sting/queue/pipeline/examples/HelloWorldPipelineTest.scala +++ b/public/scala/test/org/broadinstitute/sting/queue/pipeline/examples/HelloWorldPipelineTest.scala @@ -126,4 +126,15 @@ class HelloWorldPipelineTest { spec.jobRunners = Seq("GridEngine") PipelineTest.executeTest(spec) } + + // disabled because our DRMAA implementation doesn't support wallTime + @Test(enabled=false, timeOut=36000000) + def testHelloWorldWithWalltime() { + val spec = new PipelineTestSpec + spec.name = "HelloWorldWithWalltime" + spec.args = "-S public/scala/qscript/org/broadinstitute/sting/queue/qscripts/examples/HelloWorld.scala" + + " -wallTime 100" + spec.jobRunners = PipelineTest.allJobRunners + PipelineTest.executeTest(spec) + } } From fc7fab5f3b0798671d10b3371c87e62eb50c0ffe Mon Sep 17 00:00:00 2001 From: Mauricio Carneiro Date: Thu, 29 Nov 2012 22:10:25 -0500 Subject: [PATCH 38/52] Fixed ReadBackedPileup downsampling Downsampling in the PerSampleReadBackedPileup was broken, it didn't downsample anything, always returning a copy the original pileup. --- .../utils/pileup/AbstractReadBackedPileup.java | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/public/java/src/org/broadinstitute/sting/utils/pileup/AbstractReadBackedPileup.java b/public/java/src/org/broadinstitute/sting/utils/pileup/AbstractReadBackedPileup.java index 42938d2a6..25f0bfa6d 100644 --- a/public/java/src/org/broadinstitute/sting/utils/pileup/AbstractReadBackedPileup.java +++ b/public/java/src/org/broadinstitute/sting/utils/pileup/AbstractReadBackedPileup.java @@ -652,23 +652,19 @@ public abstract class AbstractReadBackedPileup tracker = (PerSamplePileupElementTracker) pileupElementTracker; PerSamplePileupElementTracker filteredTracker = new PerSamplePileupElementTracker(); - int current = 0; for (final String sample : tracker.getSamples()) { PileupElementTracker perSampleElements = tracker.getElements(sample); - List filteredPileup = new ArrayList(); - for (PileupElement p : perSampleElements) { + int current = 0; + UnifiedPileupElementTracker filteredPileup = new UnifiedPileupElementTracker(); + for (PE p : perSampleElements) { if (positions.contains(current)) filteredPileup.add(p); - } + current++; - if (!filteredPileup.isEmpty()) { - AbstractReadBackedPileup pileup = createNewPileup(loc, perSampleElements); - filteredTracker.addElements(sample, pileup.pileupElementTracker); } - - current++; + filteredTracker.addElements(sample, filteredPileup); } return (RBP) createNewPileup(loc, filteredTracker); From 8020ba14db0bb38422510729d50efcf044a71543 Mon Sep 17 00:00:00 2001 From: Mark DePristo Date: Fri, 30 Nov 2012 15:04:33 -0500 Subject: [PATCH 44/52] Minor cleanup of SAMDataSource as part of my system review -- Changed a few function from public to protected, as they are only used by the package contents, to simplify the SAMDataSource interface --- .../sting/gatk/datasources/reads/SAMDataSource.java | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/SAMDataSource.java b/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/SAMDataSource.java index bb788c89f..88de3ac9b 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/SAMDataSource.java +++ b/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/SAMDataSource.java @@ -30,12 +30,10 @@ import net.sf.samtools.*; import net.sf.samtools.util.CloseableIterator; import net.sf.samtools.util.RuntimeIOException; import org.apache.log4j.Logger; -import org.broadinstitute.sting.gatk.downsampling.*; -import org.broadinstitute.sting.gatk.downsampling.DownsampleType; -import org.broadinstitute.sting.gatk.downsampling.DownsamplingMethod; import org.broadinstitute.sting.gatk.ReadMetrics; import org.broadinstitute.sting.gatk.ReadProperties; import org.broadinstitute.sting.gatk.arguments.ValidationExclusion; +import org.broadinstitute.sting.gatk.downsampling.*; import org.broadinstitute.sting.gatk.filters.CountingFilteringIterator; import org.broadinstitute.sting.gatk.filters.ReadFilter; import org.broadinstitute.sting.gatk.iterators.*; @@ -567,7 +565,7 @@ public class SAMDataSource { * * @return the start positions of the first chunk of reads for all BAM files */ - public Map getInitialReaderPositions() { + protected Map getInitialReaderPositions() { Map initialPositions = new HashMap(); SAMReaders readers = resourcePool.getAvailableReaders(); @@ -585,7 +583,7 @@ public class SAMDataSource { * @param shard The shard specifying the data limits. * @return An iterator over the selected data. */ - public StingSAMIterator getIterator( Shard shard ) { + protected StingSAMIterator getIterator( Shard shard ) { return getIterator(resourcePool.getAvailableReaders(), shard, shard instanceof ReadShard); } From 2849889af5ccb850f297c644c23ca0e40e887f2f Mon Sep 17 00:00:00 2001 From: Mark DePristo Date: Sat, 1 Dec 2012 14:23:57 -0500 Subject: [PATCH 47/52] Updating md5 for UG --- .../gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java index 9d12b0ded..8ded61af8 100755 --- a/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java +++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java @@ -351,7 +351,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest { WalkerTest.WalkerTestSpec spec2 = new WalkerTest.WalkerTestSpec( baseCommandIndels + " --genotyping_mode GENOTYPE_GIVEN_ALLELES -alleles " + result.get(0).getAbsolutePath() + " -I " + validationDataLocation + "low_coverage_CEU.chr1.10k-11k.bam -o %s -L 1:10450700-10551000", 1, - Arrays.asList("c526c234947482d1cd2ffc5102083a08")); + Arrays.asList("1256a7eceff2c2374c231ff981df486d")); executeTest("test MultiSample Pilot1 CEU indels using GENOTYPE_GIVEN_ALLELES", spec2); }