{
/**
* 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/fasta/FastaAlternateReferenceMaker.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/fasta/FastaAlternateReferenceMaker.java
index 2b9744b89..22c6097cf 100755
--- a/public/java/src/org/broadinstitute/sting/gatk/walkers/fasta/FastaAlternateReferenceMaker.java
+++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/fasta/FastaAlternateReferenceMaker.java
@@ -47,6 +47,12 @@ import java.util.List;
*
* Given variant tracks, it replaces the reference bases at variation sites with the bases supplied by the ROD(s).
* Additionally, allows for one or more "snpmask" VCFs to set overlapping bases to 'N'.
+ *
+ * The output format can be partially controlled using the provided command-line arguments.
+ * Specify intervals with the usual -L argument to output only the reference bases within your intervals.
+ * Overlapping intervals are automatically merged; reference bases for each disjoint interval will be output as a
+ * separate fasta sequence (named numerically in order).
+ *
* Several important notes:
* 1) if there are multiple variants that start at a site, it chooses one of them randomly.
* 2) when there are overlapping indels (but with different start positions) only the first will be chosen.
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 80bc04845..cc086b148 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
@@ -382,11 +382,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?
*
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 RodWalkerInput
*
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/java/src/org/broadinstitute/sting/utils/GenomeLoc.java b/public/java/src/org/broadinstitute/sting/utils/GenomeLoc.java
index 6df9c9f1d..4d2c26a79 100644
--- a/public/java/src/org/broadinstitute/sting/utils/GenomeLoc.java
+++ b/public/java/src/org/broadinstitute/sting/utils/GenomeLoc.java
@@ -315,6 +315,20 @@ public class GenomeLoc implements Comparable, Serializable, HasGenome
return ( comparison == -1 || ( comparison == 0 && this.getStop() < that.getStart() ));
}
+ /**
+ * Tests whether this genome loc starts at the same position as that.
+ *
+ * i.e., do this and that have the same contig and the same start position
+ *
+ * @param that genome loc to compare to
+ * @return true if this and that have the same contig and the same start position
+ */
+ @Requires("that != null")
+ public final boolean startsAt( GenomeLoc that ) {
+ int comparison = this.compareContigs(that);
+ return comparison == 0 && this.getStart() == that.getStart();
+ }
+
/**
* Tests whether any portion of this contig is before that contig.
* @param that Other contig to test.
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/src/org/broadinstitute/sting/utils/GenomeLocSortedSet.java b/public/java/src/org/broadinstitute/sting/utils/GenomeLocSortedSet.java
index d11adf9e3..394220106 100755
--- a/public/java/src/org/broadinstitute/sting/utils/GenomeLocSortedSet.java
+++ b/public/java/src/org/broadinstitute/sting/utils/GenomeLocSortedSet.java
@@ -43,6 +43,9 @@ public class GenomeLocSortedSet extends AbstractSet {
// our private storage for the GenomeLoc's
private List mArray = new ArrayList();
+ // cache this to make overlap checking much more efficient
+ private int previousOverlapSearchIndex = -1;
+
/** default constructor */
public GenomeLocSortedSet(GenomeLocParser parser) {
this.genomeLocParser = parser;
@@ -101,7 +104,7 @@ public class GenomeLocSortedSet extends AbstractSet {
* Return the number of bps before loc in the sorted set
*
* @param loc the location before which we are counting bases
- * @return
+ * @return the number of base pairs over all previous intervals
*/
public long sizeBeforeLoc(GenomeLoc loc) {
long s = 0;
@@ -110,7 +113,7 @@ public class GenomeLocSortedSet extends AbstractSet {
if ( e.isBefore(loc) )
s += e.size();
else if ( e.isPast(loc) )
- ; // don't do anything
+ break; // we are done
else // loc is inside of s
s += loc.getStart() - e.getStart();
}
@@ -131,15 +134,43 @@ public class GenomeLocSortedSet extends AbstractSet {
* Determine if the given loc overlaps any loc in the sorted set
*
* @param loc the location to test
- * @return
+ * @return trip if the location overlaps any loc
*/
public boolean overlaps(final GenomeLoc loc) {
- for(final GenomeLoc e : mArray) {
- if(e.overlapsP(loc)) {
- return true;
- }
+ // edge condition
+ if ( mArray.isEmpty() )
+ return false;
+
+ // use the cached version first
+ if ( previousOverlapSearchIndex != -1 && overlapsAtOrImmediatelyAfterCachedIndex(loc, true) )
+ return true;
+
+ // update the cached index
+ previousOverlapSearchIndex = Collections.binarySearch(mArray, loc);
+
+ // if it matches an interval exactly, we are done
+ if ( previousOverlapSearchIndex >= 0 )
+ return true;
+
+ // check whether it overlaps the interval before or after the insertion point
+ previousOverlapSearchIndex = Math.max(0, -1 * previousOverlapSearchIndex - 2);
+ return overlapsAtOrImmediatelyAfterCachedIndex(loc, false);
+ }
+
+ private boolean overlapsAtOrImmediatelyAfterCachedIndex(final GenomeLoc loc, final boolean updateCachedIndex) {
+ // check the cached entry
+ if ( mArray.get(previousOverlapSearchIndex).overlapsP(loc) )
+ return true;
+
+ // check the entry after the cached entry since we may have moved to it
+ boolean returnValue = false;
+ if ( previousOverlapSearchIndex < mArray.size() - 1 ) {
+ returnValue = mArray.get(previousOverlapSearchIndex + 1).overlapsP(loc);
+ if ( updateCachedIndex )
+ previousOverlapSearchIndex++;
}
- return false;
+
+ return returnValue;
}
/**
@@ -155,7 +186,7 @@ public class GenomeLocSortedSet extends AbstractSet {
mArray.add(e);
return true;
} else {
- int loc = Collections.binarySearch(mArray,e);
+ final int loc = Collections.binarySearch(mArray,e);
if (loc >= 0) {
throw new ReviewedStingException("Genome Loc Sorted Set already contains the GenomicLoc " + e.toString());
} else {
diff --git a/public/java/src/org/broadinstitute/sting/utils/Haplotype.java b/public/java/src/org/broadinstitute/sting/utils/Haplotype.java
index b30d47074..30fdce75d 100755
--- a/public/java/src/org/broadinstitute/sting/utils/Haplotype.java
+++ b/public/java/src/org/broadinstitute/sting/utils/Haplotype.java
@@ -49,7 +49,9 @@ public class Haplotype {
private int alignmentStartHapwrtRef;
public int leftBreakPoint = 0;
public int rightBreakPoint = 0;
-
+ private Allele artificialAllele = null;
+ private int artificialAllelePosition = -1;
+
/**
* Create a simple consensus sequence with provided bases and a uniform quality over all bases of qual
*
@@ -71,6 +73,12 @@ public class Haplotype {
this(bases, 0);
}
+ protected Haplotype( final byte[] bases, final Allele artificialAllele, final int artificialAllelePosition ) {
+ this(bases, 0);
+ this.artificialAllele = artificialAllele;
+ this.artificialAllelePosition = artificialAllelePosition;
+ }
+
public Haplotype( final byte[] bases, final GenomeLoc loc ) {
this(bases);
this.genomeLocation = loc;
@@ -171,8 +179,25 @@ public class Haplotype {
this.cigar = cigar;
}
+ public boolean isArtificialHaplotype() {
+ return artificialAllele != null;
+ }
+
+ public Allele getArtificialAllele() {
+ return artificialAllele;
+ }
+
+ public int getArtificialAllelePosition() {
+ return artificialAllelePosition;
+ }
+
+ public void setArtificialAllele(final Allele artificialAllele, final int artificialAllelePosition) {
+ this.artificialAllele = artificialAllele;
+ this.artificialAllelePosition = artificialAllelePosition;
+ }
+
@Requires({"refInsertLocation >= 0"})
- public Haplotype insertAllele( final Allele refAllele, final Allele altAllele, final int refInsertLocation ) {
+ public Haplotype insertAllele( final Allele refAllele, final Allele altAllele, final int refInsertLocation, final int genomicInsertLocation ) {
// refInsertLocation is in ref haplotype offset coordinates NOT genomic coordinates
final int haplotypeInsertLocation = ReadUtils.getReadCoordinateForReferenceCoordinate(alignmentStartHapwrtRef, cigar, refInsertLocation, ReadUtils.ClippingTail.RIGHT_TAIL, true);
if( haplotypeInsertLocation == -1 || haplotypeInsertLocation + refAllele.length() >= bases.length ) { // desired change falls inside deletion so don't bother creating a new haplotype
@@ -182,7 +207,7 @@ public class Haplotype {
newHaplotypeBases = ArrayUtils.addAll(newHaplotypeBases, ArrayUtils.subarray(bases, 0, haplotypeInsertLocation)); // bases before the variant
newHaplotypeBases = ArrayUtils.addAll(newHaplotypeBases, altAllele.getBases()); // the alt allele of the variant
newHaplotypeBases = ArrayUtils.addAll(newHaplotypeBases, ArrayUtils.subarray(bases, haplotypeInsertLocation + refAllele.length(), bases.length)); // bases after the variant
- return new Haplotype(newHaplotypeBases);
+ return new Haplotype(newHaplotypeBases, altAllele, genomicInsertLocation);
}
public static class HaplotypeBaseComparator implements Comparator, Serializable {
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/Utils.java b/public/java/src/org/broadinstitute/sting/utils/Utils.java
index b780d0966..e4d6f6233 100755
--- a/public/java/src/org/broadinstitute/sting/utils/Utils.java
+++ b/public/java/src/org/broadinstitute/sting/utils/Utils.java
@@ -293,6 +293,10 @@ public class Utils {
}
}
+ public static String join(final String separator, final T ... objects) {
+ return join(separator, Arrays.asList(objects));
+ }
+
public static String dupString(char c, int nCopies) {
char[] chars = new char[nCopies];
Arrays.fill(chars, c);
@@ -701,11 +705,13 @@ public class Utils {
List oldRecords = header.getProgramRecords();
List newRecords = new ArrayList(oldRecords.size()+1);
for ( SAMProgramRecord record : oldRecords )
- if ( !record.getId().startsWith(programRecord.getId()) || KEEP_ALL_PG_RECORDS )
+ if ( (programRecord != null && !record.getId().startsWith(programRecord.getId())) || KEEP_ALL_PG_RECORDS )
newRecords.add(record);
- newRecords.add(programRecord);
- header.setProgramRecords(newRecords);
+ if (programRecord != null) {
+ newRecords.add(programRecord);
+ header.setProgramRecords(newRecords);
+ }
return header;
}
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/src/org/broadinstitute/sting/utils/activeregion/ActivityProfile.java b/public/java/src/org/broadinstitute/sting/utils/activeregion/ActivityProfile.java
index 38cfbb38d..e96eb843d 100644
--- a/public/java/src/org/broadinstitute/sting/utils/activeregion/ActivityProfile.java
+++ b/public/java/src/org/broadinstitute/sting/utils/activeregion/ActivityProfile.java
@@ -103,11 +103,6 @@ public class ActivityProfile {
isActiveList.add(result);
}
- // for unit testing
- public List getActiveList() {
- return isActiveList;
- }
-
public int size() {
return isActiveList.size();
}
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);
diff --git a/public/java/src/org/broadinstitute/sting/utils/codecs/vcf/VCFUtils.java b/public/java/src/org/broadinstitute/sting/utils/codecs/vcf/VCFUtils.java
index be87e7306..a8aefb703 100755
--- a/public/java/src/org/broadinstitute/sting/utils/codecs/vcf/VCFUtils.java
+++ b/public/java/src/org/broadinstitute/sting/utils/codecs/vcf/VCFUtils.java
@@ -30,12 +30,17 @@ import net.sf.samtools.SAMSequenceRecord;
import org.apache.commons.io.FilenameUtils;
import org.apache.log4j.Logger;
import org.broad.tribble.Feature;
+import org.broad.tribble.FeatureCodecHeader;
+import org.broad.tribble.readers.PositionalBufferedStream;
import org.broadinstitute.sting.commandline.RodBinding;
import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
import org.broadinstitute.sting.gatk.datasources.rmd.ReferenceOrderedDataSource;
+import org.broadinstitute.sting.utils.collections.Pair;
import org.broadinstitute.sting.utils.variantcontext.VariantContext;
import java.io.File;
+import java.io.FileInputStream;
+import java.io.IOException;
import java.util.*;
/**
@@ -317,4 +322,33 @@ public class VCFUtils {
assembly = "hg19";
return assembly;
}
+
+ /**
+ * Read all of the VCF records from source into memory, returning the header and the VariantContexts
+ *
+ * @param source the file to read, must be in VCF4 format
+ * @return
+ * @throws IOException
+ */
+ public static Pair> readVCF(final File source) throws IOException {
+ // read in the features
+ final List vcs = new ArrayList();
+ final VCFCodec codec = new VCFCodec();
+ PositionalBufferedStream pbs = new PositionalBufferedStream(new FileInputStream(source));
+ FeatureCodecHeader header = codec.readHeader(pbs);
+ pbs.close();
+
+ pbs = new PositionalBufferedStream(new FileInputStream(source));
+ pbs.skip(header.getHeaderEnd());
+
+ final VCFHeader vcfHeader = (VCFHeader)header.getHeaderValue();
+
+ while ( ! pbs.isDone() ) {
+ final VariantContext vc = codec.decode(pbs);
+ if ( vc != null )
+ vcs.add(vc);
+ }
+
+ return new Pair>(vcfHeader, vcs);
+ }
}
\ No newline at end of file
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/pileup/AbstractReadBackedPileup.java b/public/java/src/org/broadinstitute/sting/utils/pileup/AbstractReadBackedPileup.java
index ed6fc46bb..25f0bfa6d 100644
--- a/public/java/src/org/broadinstitute/sting/utils/pileup/AbstractReadBackedPileup.java
+++ b/public/java/src/org/broadinstitute/sting/utils/pileup/AbstractReadBackedPileup.java
@@ -254,19 +254,32 @@ public abstract class AbstractReadBackedPileup tracker = (PerSamplePileupElementTracker) pileupElementTracker;
PerSamplePileupElementTracker filteredTracker = new PerSamplePileupElementTracker();
for (final String sample : tracker.getSamples()) {
PileupElementTracker perSampleElements = tracker.getElements(sample);
- AbstractReadBackedPileup pileup = createNewPileup(loc, perSampleElements).getOverlappingFragmentFilteredPileup();
+ AbstractReadBackedPileup pileup = createNewPileup(loc, perSampleElements).getOverlappingFragmentFilteredPileup(discardDiscordant, baseQualNotMapQual);
filteredTracker.addElements(sample, pileup.pileupElementTracker);
}
return (RBP) createNewPileup(loc, filteredTracker);
@@ -284,11 +297,16 @@ public abstract class AbstractReadBackedPileup tracker = (PerSamplePileupElementTracker) pileupElementTracker;
PerSamplePileupElementTracker filteredTracker = new PerSamplePileupElementTracker