(new AlignmentStartWithNoTiesComparator());
- if (!regions.isEmpty()) {
- for (SingleSampleCompressor sample : compressorsPerSample.values()) {
- reads.addAll(sample.closeVariantRegions(regions, knownSnpPositions));
- }
- }
- return reads;
- }
-}
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReads.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReads.java
deleted file mode 100644
index 383ba5ee9..000000000
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReads.java
+++ /dev/null
@@ -1,782 +0,0 @@
-/*
-* By downloading the PROGRAM you agree to the following terms of use:
-*
-* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
-*
-* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
-*
-* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
-* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
-* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
-*
-* 1. DEFINITIONS
-* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
-*
-* 2. LICENSE
-* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
-* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
-* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
-* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
-*
-* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
-* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
-* Copyright 2012 Broad Institute, Inc.
-* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
-* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
-*
-* 4. INDEMNIFICATION
-* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
-*
-* 5. NO REPRESENTATIONS OR WARRANTIES
-* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
-* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
-*
-* 6. ASSIGNMENT
-* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
-*
-* 7. MISCELLANEOUS
-* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
-* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
-* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
-* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
-* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
-* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
-* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
-*/
-
-package org.broadinstitute.sting.gatk.walkers.compression.reducereads;
-
-import it.unimi.dsi.fastutil.objects.Object2LongOpenHashMap;
-import it.unimi.dsi.fastutil.objects.ObjectAVLTreeSet;
-import it.unimi.dsi.fastutil.objects.ObjectArrayList;
-import it.unimi.dsi.fastutil.objects.ObjectSortedSet;
-import net.sf.samtools.SAMFileHeader;
-import net.sf.samtools.SAMFileWriter;
-import net.sf.samtools.SAMProgramRecord;
-import net.sf.samtools.util.SequenceUtil;
-import org.broadinstitute.sting.commandline.*;
-import org.broadinstitute.sting.gatk.CommandLineGATK;
-import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
-import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
-import org.broadinstitute.sting.gatk.downsampling.DownsampleType;
-import org.broadinstitute.sting.gatk.filters.*;
-import org.broadinstitute.sting.gatk.io.StingSAMFileWriter;
-import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
-import org.broadinstitute.sting.gatk.walkers.*;
-import org.broadinstitute.sting.utils.GenomeLoc;
-import org.broadinstitute.sting.utils.SampleUtils;
-import org.broadinstitute.sting.utils.Utils;
-import org.broadinstitute.sting.utils.clipping.ReadClipper;
-import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
-import org.broadinstitute.sting.utils.exceptions.UserException;
-import org.broadinstitute.sting.utils.help.DocumentedGATKFeature;
-import org.broadinstitute.sting.utils.help.HelpConstants;
-import org.broadinstitute.sting.utils.sam.BySampleSAMFileWriter;
-import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
-import org.broadinstitute.sting.utils.sam.ReadUtils;
-import org.broadinstitute.variant.variantcontext.VariantContext;
-
-import java.util.Collections;
-import java.util.List;
-
-
-/**
- * Reduces the BAM file using read based compression that keeps only essential information for variant calling
- *
- *
- * This tool will generate reduced versions of the BAM files that still follow the BAM specification
- * and contain all the information necessary to call variants according to the GATK Best Practices recommendations.
- * Some options allow you to tune how much compression you want to achieve. The default values have been
- * shown to reduce a typical whole exome BAM file by 100x. The higher the coverage, the bigger the
- * savings in file size and performance of the downstream tools.
- *
- *
Input
- *
- * The BAM file to be compressed
- *
- *
- * Output
- *
- * The compressed (reduced) BAM file.
- *
- *
- * Examples
- *
- * java -Xmx4g -jar GenomeAnalysisTK.jar \
- * -R ref.fasta \
- * -T ReduceReads \
- * -I myData.bam \
- * -o myData.reduced.bam
- *
- */
-
-@DocumentedGATKFeature( groupName = HelpConstants.DOCS_CAT_DATA, extraDocs = {CommandLineGATK.class} )
-@PartitionBy(PartitionType.CONTIG)
-@ReadFilters({UnmappedReadFilter.class, NotPrimaryAlignmentFilter.class, DuplicateReadFilter.class, FailsVendorQualityCheckFilter.class, BadCigarFilter.class})
-@Downsample(by=DownsampleType.BY_SAMPLE, toCoverage=40)
-public class ReduceReads extends ReadWalker, ReduceReadsStash> {
-
- @Output(required = false, defaultToStdout = false)
- private StingSAMFileWriter out = null;
- private SAMFileWriter writerToUse = null;
-
- /**
- *
- */
- @Argument(fullName = "context_size", shortName = "cs", doc = "The number of bases to keep around mismatches (potential variation)", required = false)
- public int contextSize = 10;
-
- /**
- * Reads that have
- * mapping quality below this threshold will not be counted towards consensus, but are still counted
- * towards variable regions.
- */
- @Argument(fullName = "minimum_mapping_quality", shortName = "minmap", doc = "The minimum mapping quality to be considered for the consensus synthetic read", required = false)
- public int minMappingQuality = 20;
-
- /**
- * Reads that have
- * base quality below this threshold will not be counted towards consensus, but are still counted
- * towards variable regions.
- */
- @Argument(fullName = "minimum_base_quality_to_consider", shortName = "minqual", doc = "The minimum base quality to be considered for the consensus synthetic read", required = false)
- public byte minBaseQual = 15;
-
- /**
- * Reads have notoriously low quality bases on the tails (left and right). Consecutive bases at the tails with
- * quality at or lower than this threshold will be hard clipped off before entering the reduce reads algorithm.
- */
- @Argument(fullName = "minimum_tail_qualities", shortName = "mintail", doc = "", required = false)
- public byte minTailQuality = 2;
-
- /**
- * Any number of VCF files representing known SNPs to be used for the polyploid-based reduction.
- * Could be e.g. dbSNP and/or official 1000 Genomes SNP calls. Non-SNP variants in these files will be ignored.
- * If provided, the polyploid ("het") compression will work only when a single SNP from the known set is present
- * in a consensus window (otherwise there will be no reduction); if not provided then polyploid compression will
- * be triggered anywhere there is a single SNP present in a consensus window.
- */
- @Input(fullName="known_sites_for_polyploid_reduction", shortName = "known", doc="Input VCF file(s) with known SNPs", required=false)
- public List> known = Collections.emptyList();
-
- /**
- * This strips away all extra information of the read -- anything other than bases, quals
- * and read group.
- */
- @Argument(fullName = "dont_simplify_reads", shortName = "nosimplify", doc = "Do not simplify read", required = false)
- public boolean DONT_SIMPLIFY_READS = false;
-
- /**
- * Note that it is not necessary to turn this on for reads that are not mate paired.
- * The program will behave correctly by default in those cases.
- */
- @Argument(fullName = "dont_hardclip_adaptor_sequences", shortName = "noclip_ad", doc = "Do not hard clip adaptor sequences", required = false)
- public boolean DONT_CLIP_ADAPTOR_SEQUENCES = false;
-
- /**
- * This option overrides the argument of minimum tail
- * quality.
- */
- @Argument(fullName = "dont_hardclip_low_qual_tails", shortName = "noclip_tail", doc = "Do not hard clip the low quality tails of the reads", required = false)
- public boolean DONT_CLIP_LOW_QUAL_TAILS = false;
-
- /**
- * By default, ReduceReads will hard clip away any low quality soft clipped
- * base left by the aligner and use the high quality soft clipped bases in it's traversal algorithm to identify variant
- * regions. The minimum quality for soft clipped bases is the same as the minimum base quality to consider (minqual)
- */
- @Argument(fullName = "dont_use_softclipped_bases", shortName = "no_soft", doc = "Do not use high quality soft-clipped bases", required = false)
- public boolean DONT_USE_SOFTCLIPPED_BASES = false;
-
- /**
- * By default, ReduceReads will compress read names to numbers and guarantee
- * uniqueness and reads with similar name will still have similar compressed names. Note: If you scatter/gather
- * there is no guarantee that read name uniqueness will be maintained -- in this case we recommend not compressing.
- */
- @Argument(fullName = "dont_compress_read_names", shortName = "nocmp_names", doc = "Do not compress read names", required = false)
- public boolean DONT_COMPRESS_READ_NAMES = false;
-
- /**
- * The hard clips will happen exactly at the interval border.
- */
- @Argument(fullName = "hard_clip_to_interval", shortName = "clip_int", doc = "Hard clip all incoming reads to the desired intervals", required = false)
- public boolean HARD_CLIP_TO_INTERVAL = false;
-
- /**
- * Anything below this will be
- * considered consensus and reduced (otherwise we will try to trigger polyploid compression). Note that
- * this value is used only regions with high coverage.
- */
- @Advanced
- @Argument(fullName = "minimum_alt_proportion_to_trigger_variant", shortName = "minvar", doc = "Minimum proportion of mismatches in a site to trigger a variant region", required = false)
- public double minAltProportionToTriggerVariant = 0.05;
-
- /**
- * Any site with a value falling below this will be considered consensus and reduced (otherwise we will try to
- * trigger polyploid compression). Note that this value is used only regions with low coverage.
- */
- @Advanced
- @Argument(fullName = "minimum_alt_pvalue_to_trigger_variant", shortName = "min_pvalue", doc = "Minimum p-value from binomial distribution of mismatches in a site to trigger a variant region", required = false)
- public double minAltPValueToTriggerVariant = 0.01;
-
- /**
- * Anything below this will be considered consensus.
- */
- @Argument(fullName = "minimum_del_proportion_to_trigger_variant", shortName = "mindel", doc = "Minimum proportion of indels in a site to trigger a variant region", required = false)
- public double minIndelProportionToTriggerVariant = 0.05;
-
- /**
- * This level of downsampling only happens after the region has been evaluated, therefore it can
- * be combined with the engine level downsampling.
- * A value of 0 turns downsampling off.
- */
- @Argument(fullName = "downsample_coverage", shortName = "ds", doc = "Downsample the number of reads emitted per sample in a variant region for better compression", required = false)
- public int downsampleCoverage = 250;
-
- /**
- * Generally, this tool is not meant to be run for more than 1 sample at a time. The one valid exception
- * brought to our attention by colleagues is the specific case of tumor/normal pairs in cancer analysis.
- * To prevent users from unintentionally running the tool in a less than ideal manner, we require them
- * to explicitly enable multi-sample analysis with this argument.
- */
- @Argument(fullName = "cancer_mode", shortName = "cancer_mode", doc = "Enable multi-sample reduction for cancer analysis", required = false)
- public boolean ALLOW_MULTIPLE_SAMPLES = false;
-
- @Hidden
- @Argument(fullName = "nwayout", shortName = "nw", doc = "Generate separate output files per input file", required = false)
- public boolean nwayout = false;
-
- @Hidden
- @Argument(fullName = "", shortName = "dl", doc = "Debug level", required = false)
- public int debugLevel = 0;
-
- @Hidden
- @Argument(fullName = "", shortName = "dr", doc = "Debug read", required = false)
- public String debugRead = "";
-
- @Hidden
- @Argument(fullName = "downsample_strategy", shortName = "dm", doc = "Downsampling strategy", required = false)
- public DownsampleStrategy downsampleStrategy = DownsampleStrategy.Normal;
-
- @Hidden
- @Argument(fullName = "no_pg_tag", shortName = "npt", doc ="Discard program tags", required = false)
- public boolean NO_PG_TAG = false;
-
- public enum DownsampleStrategy {
- Normal,
- Adaptive
- }
-
- int nCompressedReads = 0;
-
- private static int READ_NAME_HASH_DEFAULT_SIZE = 1000;
- Long nextReadNumber = 1L; // The next number to use for the compressed read name.
- Object2LongOpenHashMap readNameHash; // This hash will keep the name of the original read the new compressed name (a number).
-
- ObjectSortedSet intervalList;
-
- ObjectSortedSet knownSnpPositions;
-
- // IMPORTANT: DO NOT CHANGE THE VALUE OF THIS CONSTANT VARIABLE; IT IS NOW PERMANENTLY THE @PG NAME THAT EXTERNAL TOOLS LOOK FOR IN THE BAM HEADER
- public static final String PROGRAM_RECORD_NAME = "GATK ReduceReads"; // The name that will go in the @PG tag
- private static final String PROGRAM_FILENAME_EXTENSION = ".reduced.bam";
-
- /**
- * Basic generic initialization of the readNameHash and the intervalList. Output initialization
- * is done at the reduceInit method
- */
- @Override
- public void initialize() {
- super.initialize();
-
- if ( !nwayout && out == null )
- throw new UserException.MissingArgument("out", "the output must be provided and is optional only for certain debugging modes");
-
- if ( nwayout && out != null )
- throw new UserException.CommandLineException("--out and --nwayout cannot be used simultaneously; please use one or the other");
-
- if ( minAltPValueToTriggerVariant < 0.0 || minAltPValueToTriggerVariant > 1.0 )
- throw new UserException.BadArgumentValue("--minimum_alt_pvalue_to_trigger_variant", "must be a value between 0 and 1 (inclusive)");
-
- if ( minAltProportionToTriggerVariant < 0.0 || minAltProportionToTriggerVariant > 1.0 )
- throw new UserException.BadArgumentValue("--minimum_alt_proportion_to_trigger_variant", "must be a value between 0 and 1 (inclusive)");
-
- if ( SampleUtils.getSAMFileSamples(getToolkit().getSAMFileHeader()).size() > 1 && !ALLOW_MULTIPLE_SAMPLES )
- throw new UserException.BadInput("Reduce Reads is not meant to be run for more than 1 sample at a time except for the specific case of tumor/normal pairs in cancer analysis. If that is what you want to do, use the -cancer_mode flag.");
-
- if ( known.isEmpty() )
- knownSnpPositions = null;
- else
- knownSnpPositions = new ObjectAVLTreeSet();
-
- GenomeAnalysisEngine toolkit = getToolkit();
- this.resetReadNameHash(); // prepare the read name hash to keep track of what reads have had their read names compressed
- intervalList = new ObjectAVLTreeSet(); // get the interval list from the engine. If no interval list was provided, the walker will work in WGS mode
-
- if (toolkit.getIntervals() != null)
- intervalList.addAll(toolkit.getIntervals());
-
- final boolean indexOnTheFly = true;
- final SAMFileHeader.SortOrder sortOrder = SAMFileHeader.SortOrder.coordinate;
- if (nwayout) {
- SAMProgramRecord programRecord = NO_PG_TAG ? null : Utils.createProgramRecord(toolkit, this, PROGRAM_RECORD_NAME);
- writerToUse = new BySampleSAMFileWriter(toolkit, PROGRAM_FILENAME_EXTENSION, sortOrder, false, indexOnTheFly, NO_PG_TAG, programRecord, true);
- }
- else {
- writerToUse = out;
- out.setPresorted(false);
- if (!NO_PG_TAG) {
- Utils.setupWriter(out, toolkit, toolkit.getSAMFileHeader(), false, this, PROGRAM_RECORD_NAME);
- }
- }
- }
-
- /** Initializer for {@link #readNameHash}. */
- private void resetReadNameHash() {
- // If the hash grows large, subsequent clear operations can be very expensive, so trim the hash down if it grows beyond its default.
- if (readNameHash == null || readNameHash.size() > READ_NAME_HASH_DEFAULT_SIZE) {
- readNameHash = new Object2LongOpenHashMap(READ_NAME_HASH_DEFAULT_SIZE);
- } else {
- readNameHash.clear();
- }
- }
-
- /**
- * Takes in a read and prepares it for the SlidingWindow machinery by performing the
- * following optional clipping operations:
- * 1. Hard clip adaptor sequences
- * 2. Hard clip low quality tails
- * 3. Hard clip all remaining soft clipped bases
- * 4. Hard clip read to the intervals in the interval list (this step may produce multiple reads)
- *
- * @param ref default map parameter
- * @param read default map parameter
- * @param metaDataTracker default map parameter
- * @return a linked list with all the reads produced by the clipping operations
- */
- @Override
- public ObjectArrayList map(ReferenceContext ref, GATKSAMRecord read, RefMetaDataTracker metaDataTracker) {
- ObjectArrayList mappedReads;
- if (!debugRead.isEmpty() && read.getReadName().contains(debugRead))
- System.out.println("Found debug read!");
-
- if (debugLevel == 1)
- System.out.printf("\nOriginal: %s %s %d %d\n", read, read.getCigar(), read.getAlignmentStart(), read.getAlignmentEnd());
-
- // we write the actual alignment starts to their respective alignment shift tags in the temporary
- // attribute hash so we can determine later if we need to write down the alignment shift to the reduced BAM file
- read.setTemporaryAttribute(GATKSAMRecord.REDUCED_READ_ORIGINAL_ALIGNMENT_START_SHIFT, read.getAlignmentStart());
- read.setTemporaryAttribute(GATKSAMRecord.REDUCED_READ_ORIGINAL_ALIGNMENT_END_SHIFT, read.getAlignmentEnd());
-
- // Check if the read goes beyond the boundaries of the chromosome, and hard clip those boundaries.
- int chromosomeLength = ref.getGenomeLocParser().getContigInfo(read.getReferenceName()).getSequenceLength();
- if (read.getSoftStart() < 0)
- read = ReadClipper.hardClipByReadCoordinates(read, 0, -read.getSoftStart());
- if (read.getSoftEnd() > chromosomeLength)
- read = ReadClipper.hardClipByReadCoordinates(read, chromosomeLength - read.getSoftStart() + 1, read.getReadLength() - 1);
-
- if (!DONT_SIMPLIFY_READS)
- read.simplify(); // Clear all unnecessary attributes
- if (!DONT_CLIP_ADAPTOR_SEQUENCES)
- read = ReadClipper.hardClipAdaptorSequence(read); // Strip away adaptor sequences, if any.
- if (!DONT_CLIP_LOW_QUAL_TAILS)
- read = ReadClipper.hardClipLowQualEnds(read, minTailQuality); // Clip low quality tails
- if (!isWholeGenome()) {
- if (HARD_CLIP_TO_INTERVAL)
- mappedReads = hardClipReadToInterval(read); // Hard clip the remainder of the read to the desired interval
- else {
- mappedReads = new ObjectArrayList();
- mappedReads.add(read);
- }
- }
- else {
- mappedReads = new ObjectArrayList();
- if (!read.isEmpty())
- mappedReads.add(read);
- }
-
- if (!mappedReads.isEmpty() && !DONT_USE_SOFTCLIPPED_BASES) {
- ObjectArrayList tempList = new ObjectArrayList();
- for (GATKSAMRecord mRead : mappedReads) {
- GATKSAMRecord clippedRead = ReadClipper.hardClipLowQualitySoftClips(mRead, minBaseQual);
- if (!clippedRead.isEmpty())
- tempList.add(clippedRead);
- }
- mappedReads = tempList;
- }
-
- if (debugLevel == 1)
- for (GATKSAMRecord mappedRead : mappedReads)
- System.out.printf("MAPPED: %s %d %d\n", mappedRead.getCigar(), mappedRead.getAlignmentStart(), mappedRead.getAlignmentEnd());
-
- // add the SNPs to the list of known positions
- populateKnownSNPs(metaDataTracker);
-
- return mappedReads;
- }
-
- /*
- * Add the positions of known SNPs to the set so that we can keep track of it
- *
- * @param metaDataTracker the ref meta data tracker
- */
- protected void populateKnownSNPs(final RefMetaDataTracker metaDataTracker) {
- for ( final VariantContext vc : metaDataTracker.getValues(known) ) {
- if ( vc.isSNP() )
- knownSnpPositions.add(getToolkit().getGenomeLocParser().createGenomeLoc(vc));
- }
- }
-
- /**
- * Initializes the ReduceReadsStash that keeps track of all reads that are waiting to
- * enter the SlidingWindow machinery. The stash makes sure reads are served in order
- * even though map() may generate reads that are only supposed to enter the machinery
- * in the future.
- *
- * @return the empty stash
- */
- @Override
- public ReduceReadsStash reduceInit() {
- return new ReduceReadsStash(new MultiSampleCompressor(getToolkit().getSAMFileHeader(), contextSize, downsampleCoverage, minMappingQuality, minAltPValueToTriggerVariant, minAltProportionToTriggerVariant, minIndelProportionToTriggerVariant, minBaseQual, downsampleStrategy));
- }
-
- /**
- * Takes the list of reads produced by map(), adds them to the stash (which keeps them sorted) and process
- * all reads that come before the original read (the read that was passed to map) including the original
- * read. This is where we send reads, in order, to the SlidingWindow machinery.
- *
- * @param mappedReads the list of reads sent by map
- * @param stash the stash that keeps the reads in order for processing
- * @return the stash with all reads that have not been processed yet
- */
- public ReduceReadsStash reduce(ObjectArrayList mappedReads, ReduceReadsStash stash) {
- if (debugLevel == 1)
- stash.print();
-
- boolean firstRead = true;
- for (GATKSAMRecord read : mappedReads) {
- boolean originalRead = firstRead && isOriginalRead(mappedReads, read);
-
- if (read.getReadLength() == 0)
- throw new ReviewedStingException("Empty read sent to reduce, this should never happen! " + read.getReadName() + " -- " + read.getCigar() + " -- " + read.getReferenceName() + ":" + read.getAlignmentStart() + "-" + read.getAlignmentEnd());
-
- if (originalRead) {
- ObjectArrayList readsReady = new ObjectArrayList();
- readsReady.addAll(stash.getAllReadsBefore(read));
- readsReady.add(read);
-
- for (GATKSAMRecord readReady : readsReady) {
- if (debugLevel == 1)
- System.out.println("REDUCE: " + readReady.getCigar() + " " + readReady.getAlignmentStart() + " " + readReady.getAlignmentEnd());
-
- for (GATKSAMRecord compressedRead : stash.compress(readReady, knownSnpPositions))
- outputRead(compressedRead);
-
- // We only care about maintaining the link between read pairs if they are in the same variant
- // region. Since an entire variant region's worth of reads is returned in a single call to
- // stash.compress(), the readNameHash can be cleared after the for() loop above.
- // The advantage of clearing the hash is that otherwise it holds all reads that have been encountered,
- // which can use a lot of memory and cause RR to slow to a crawl and/or run out of memory.
- this.resetReadNameHash();
-
- }
- } else
- stash.add(read);
-
- firstRead = false;
- }
-
- // reduce memory requirements by removing old positions
- if ( !mappedReads.isEmpty() )
- clearStaleKnownPositions(mappedReads.get(0));
-
- return stash;
- }
-
- /**
- * Now that now more reads will come, we process all the remaining reads in the stash, in order.
- *
- * @param stash the ReduceReadsStash with all unprocessed reads (from reduce)
- */
- @Override
- public void onTraversalDone(ReduceReadsStash stash) {
-
- // output any remaining reads in the compressor
- for (GATKSAMRecord read : stash.close(knownSnpPositions))
- outputRead(read);
-
- if (nwayout)
- writerToUse.close();
- }
-
- /**
- * Removes known positions that are no longer relevant for use with het compression.
- *
- * @param read the current read, used for checking whether there are stale positions we can remove
- */
- protected void clearStaleKnownPositions(final GATKSAMRecord read) {
- // nothing to clear if not used or empty
- if ( knownSnpPositions == null || knownSnpPositions.isEmpty() )
- return;
-
- // not ready to be cleared until we encounter a read from a different contig
- final int contigIndexOfRead = read.getReferenceIndex();
- if ( knownSnpPositions.first().getContigIndex() == contigIndexOfRead )
- return;
-
- // because we expect most elements to be stale, it's not going to be efficient to remove them one at a time
- final ObjectAVLTreeSet goodLocs = new ObjectAVLTreeSet();
- for ( final GenomeLoc loc : knownSnpPositions ) {
- if ( loc.getContigIndex() == contigIndexOfRead )
- goodLocs.add(loc);
- }
- knownSnpPositions.clear();
- knownSnpPositions.addAll(goodLocs);
- }
-
- /**
- * Hard clips away all parts of the read that doesn't agree with the intervals selected.
- *
- * Note: If read overlaps more than one interval, it will be hard clipped to all
- * the intervals it overlaps with
- *
- * @param read the read to be hard clipped to the interval.
- * @return a shallow copy of the read hard clipped to the interval
- */
- private ObjectArrayList hardClipReadToInterval(GATKSAMRecord read) {
- ObjectArrayList clippedReads = new ObjectArrayList();
-
- GenomeLoc intervalOverlapped = null; // marks the interval to which the original read overlapped (so we can cut all previous intervals from the list)
-
- boolean originalRead = true; // false if this is the right tail of the original read
- boolean overlap; // keeps track of the interval that overlapped the original read
- boolean doneClipping; // triggers an early exit if we are done clipping this read
-
- if (isWholeGenome())
- clippedReads.add(read); // if we don't have intervals (wgs) the read goes in unchanged
-
- for (GenomeLoc interval : intervalList) {
-
- if (read.isEmpty()) // nothing to do with an empty read (could have been fully clipped before)
- break;
-
- GATKSAMRecord clippedRead = null; // this will hold the read clipped to the interval to be added in the end of the switch
-
- switch (ReadUtils.getReadAndIntervalOverlapType(read, interval)) {
- case NO_OVERLAP_RIGHT: // no reads on this interval, check the next interval if this is the original read
- if (!originalRead) // something went wrong if this is the tail of the read
- throw new ReviewedStingException("tail of the read should never NO_OVERLAP_RIGHT the following interval. " + read.getReadName() + " -- " + read.getReferenceName() + ":" + read.getAlignmentStart() + "-" + read.getAlignmentEnd() + " x " + interval.getLocation().toString());
- overlap = false;
- doneClipping = false;
- break;
-
-
- case NO_OVERLAP_HARDCLIPPED_RIGHT: // read used to overlap but got hard clipped and doesn't overlap anymore
- if (originalRead) {
- overlap = true; // effectively, we have found the read's location and now we are going to try and match it's tail (which happens to be the entire read).
- clippedRead = GATKSAMRecord.emptyRead(read);
- } else
- overlap = false;
-
- doneClipping = false;
- break;
-
- case NO_OVERLAP_CONTIG: // read is in a different contig
- if (originalRead) { // the original read can be in a bigger contig, but not on a smaller one.
- if (read.getReferenceIndex() < interval.getContigIndex())
- throw new ReviewedStingException("read is behind interval list. (contig) " + read.getReadName() + " -- " + read.getReferenceName() + ":" + read.getAlignmentStart() + "-" + read.getAlignmentEnd() + " x " + interval.getLocation().toString());
- else {
- overlap = false;
- doneClipping = false;
- }
- } // tail read CANNOT be in a different contig.
- else {
- if (read.getReferenceIndex() < interval.getContigIndex()) {
- overlap = false;
- doneClipping = true;
- } else
- throw new ReviewedStingException("Tail read is in bigger contig than interval traversal. " + read.getReadName() + " -- " + read.getReferenceName() + ":" + read.getAlignmentStart() + "-" + read.getAlignmentEnd() + " x " + interval.getLocation().toString());
-
- }
- break;
-
- case NO_OVERLAP_LEFT:
- if (originalRead) // if this is the first read this should never happen.
- throw new ReviewedStingException("original read cannot be behind the first interval. (position) " + read.getReadName() + " -- " + read.getReferenceName() + ":" + read.getAlignmentStart() + "-" + read.getAlignmentEnd() + " x " + interval.getLocation().toString());
-
- overlap = false;
- doneClipping = true;
- break;
-
- case NO_OVERLAP_HARDCLIPPED_LEFT: // read used to overlap but got hard clipped and doesn't overlap anymore
- overlap = originalRead; // if this is the original read, we should not advance the interval list, the original overlap was here.
- doneClipping = true;
- break;
-
- case OVERLAP_LEFT: // clip the left tail of the read
- clippedRead = ReadClipper.hardClipByReferenceCoordinatesLeftTail(read, interval.getStart() - 1);
-
- overlap = true;
- doneClipping = true;
- break;
-
- case OVERLAP_RIGHT: // clip the right tail of the read and try to match it to the next interval
- clippedRead = ReadClipper.hardClipByReferenceCoordinatesRightTail(read, interval.getStop() + 1);
- read = ReadClipper.hardClipByReferenceCoordinatesLeftTail(read, interval.getStop());
-
- overlap = true;
- doneClipping = false;
- break;
-
- case OVERLAP_LEFT_AND_RIGHT: // clip both left and right ends of the read
- clippedRead = ReadClipper.hardClipBothEndsByReferenceCoordinates(read, interval.getStart() - 1, interval.getStop() + 1);
- read = ReadClipper.hardClipByReferenceCoordinatesLeftTail(read, interval.getStop());
-
- overlap = true;
- doneClipping = false;
- break;
-
- case OVERLAP_CONTAINED: // don't do anything to the read
- clippedRead = read;
-
- overlap = true;
- doneClipping = true;
- break;
-
- default:
- throw new ReviewedStingException("interval overlap returned an unknown / unhandled state. If new state was added to intervalOverlap, it should be handled by hardClipReadToInterval.");
- }
-
- if (overlap && originalRead)
- intervalOverlapped = interval;
-
- if (clippedRead != null) {
- originalRead = false;
-
- if (!clippedRead.isEmpty())
- clippedReads.add(clippedRead); // if the read overlaps the interval entirely within a deletion, it will be entirely clipped off
- }
-
- if (doneClipping)
- break;
- }
-
- if (intervalOverlapped != null)
- intervalList = intervalList.tailSet(intervalOverlapped);
-
- return clippedReads;
- }
-
- /**
- * Compresses the read name and adds it to output BAM file (reduced BAM)
- * after performing some quality control
- *
- * @param read any read
- */
- private void outputRead(GATKSAMRecord read) {
- if (debugLevel == 2) {
- checkForHighMismatch(read);
- checkCigar(read);
- }
-
- if (read.isReducedRead())
- nCompressedReads++;
- else {
- int originalAlignmentStart = (Integer) read.getTemporaryAttribute(GATKSAMRecord.REDUCED_READ_ORIGINAL_ALIGNMENT_START_SHIFT);
- int originalAlignmentEnd = (Integer) read.getTemporaryAttribute(GATKSAMRecord.REDUCED_READ_ORIGINAL_ALIGNMENT_END_SHIFT);
-
- int startShift = originalAlignmentStart - read.getUnclippedStart(); // we annotate the shifts for better compression
- int endShift = read.getUnclippedEnd() - originalAlignmentEnd; // we annotate the shifts for better compression
-
- if (startShift > 0)
- read.setAttribute(GATKSAMRecord.REDUCED_READ_ORIGINAL_ALIGNMENT_START_SHIFT, startShift); // If the read had any soft clips before getting chopped (variant region) annotate it's original alignment (start)
- if (endShift > 0)
- read.setAttribute(GATKSAMRecord.REDUCED_READ_ORIGINAL_ALIGNMENT_END_SHIFT, endShift); // If the read had any soft clips before getting chopped (variant region) annotate it's original alignment (end)
- }
-
- if (debugLevel == 1)
- System.out.println("BAM: " + read.getCigar() + " " + read.getAlignmentStart() + " " + read.getAlignmentEnd());
-
- if (!DONT_COMPRESS_READ_NAMES)
- nextReadNumber = compressReadName(readNameHash, read, nextReadNumber);
-
- writerToUse.addAlignment(read);
- }
-
- /**
- * Quality control procedure that checks if the consensus reads contains too many
- * mismatches with the reference. This should never happen and is a good trigger for
- * errors with the algorithm.
- *
- * @param read any read
- */
- private void checkForHighMismatch(GATKSAMRecord read) {
- final int start = read.getAlignmentStart();
- final int stop = read.getAlignmentEnd();
- final byte[] ref = getToolkit().getReferenceDataSource().getReference().getSubsequenceAt(read.getReferenceName(), start, stop).getBases();
- final int nm = SequenceUtil.countMismatches(read, ref, start - 1);
- final int readLen = read.getReadLength();
- final double nmFraction = nm / (1.0 * readLen);
- if (nmFraction > 0.4 && readLen > 20 && read.getAttribute(GATKSAMRecord.REDUCED_READ_CONSENSUS_TAG) != null && read.getReadName().startsWith("Consensus"))
- throw new ReviewedStingException("BUG: High mismatch fraction found in read " + read.getReadName() + " position: " + read.getReferenceName() + ":" + read.getAlignmentStart() + "-" + read.getAlignmentEnd());
- }
-
- private void checkCigar (GATKSAMRecord read) {
- if (read.getCigar().isValid(null, -1) != null) {
- throw new ReviewedStingException("BUG: cigar string is not valid: " + read.getCigarString());
- }
-
- }
-
-
- /**
- * Compresses the read name using the readNameHash if we have already compressed
- * this read name before.
- *
- * @param hash the hash table containing the read name to compressed read name map
- * @param read any read
- * @param nextReadNumber the number to use in the compressed read name in case this is a new read name
- * @return the next number to use in the compressed read name
- */
- protected static long compressReadName(final Object2LongOpenHashMap hash, final GATKSAMRecord read, final long nextReadNumber) {
- final String name = read.getReadName();
- final StringBuilder compressedName = new StringBuilder();
- long result = nextReadNumber;
- if (read.isReducedRead()) {
- compressedName.append("C");
- }
- final Long readNumber = hash.get(name);
- if (readNumber != null) {
- compressedName.append(readNumber);
- } else {
- hash.put(name, nextReadNumber);
- compressedName.append(nextReadNumber);
- result++;
- }
- read.setReadName(compressedName.toString());
- return result;
- }
-
- /**
- * Returns true if the read is the original read that went through map().
- *
- * This is important to know so we can decide what reads to pull from the stash. Only reads that came before the original read should be pulled.
- *
- * @param list the list
- * @param read the read
- * @return Returns true if the read is the original read that went through map().
- */
- private boolean isOriginalRead(ObjectArrayList list, GATKSAMRecord read) {
- return isWholeGenome() || list.get(0).equals(read);
- }
-
- /**
- * Checks whether or not the intervalList is empty, meaning we're running in WGS mode.
- *
- * @return whether or not we're running in WGS mode.
- */
- private boolean isWholeGenome() {
- return intervalList.isEmpty();
- }
-
-}
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReadsStash.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReadsStash.java
deleted file mode 100644
index 52c5f0903..000000000
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReadsStash.java
+++ /dev/null
@@ -1,160 +0,0 @@
-/*
-* By downloading the PROGRAM you agree to the following terms of use:
-*
-* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
-*
-* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
-*
-* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
-* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
-* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
-*
-* 1. DEFINITIONS
-* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
-*
-* 2. LICENSE
-* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
-* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
-* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
-* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
-*
-* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
-* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
-* Copyright 2012 Broad Institute, Inc.
-* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
-* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
-*
-* 4. INDEMNIFICATION
-* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
-*
-* 5. NO REPRESENTATIONS OR WARRANTIES
-* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
-* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
-*
-* 6. ASSIGNMENT
-* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
-*
-* 7. MISCELLANEOUS
-* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
-* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
-* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
-* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
-* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
-* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
-* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
-*/
-
-package org.broadinstitute.sting.gatk.walkers.compression.reducereads;
-
-import it.unimi.dsi.fastutil.objects.ObjectSortedSet;
-import org.broadinstitute.sting.utils.GenomeLoc;
-import org.broadinstitute.sting.utils.sam.AlignmentStartWithNoTiesComparator;
-import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
-import org.broadinstitute.sting.utils.sam.ReadUtils;
-
-import java.util.LinkedList;
-import java.util.List;
-import java.util.SortedSet;
-import java.util.TreeSet;
-
-/**
- * This class implements a "read stash" that keeps reads always sorted in alignment order. Useful
- * for read walkers that alter the alignment information of the incoming reads, but need to
- * maintain the reads sorted for the reduce step. (e.g. ReduceReads)
- */
-
-public class ReduceReadsStash {
- protected MultiSampleCompressor compressor;
- SortedSet outOfOrderReads;
-
- /**
- * Creates a stash with the default sorting order (read alignment)
- * @param compressor the MultiSampleCompressor object to be used with this stash (for stash.close())
- */
- public ReduceReadsStash(MultiSampleCompressor compressor) {
- this.compressor = compressor;
- this.outOfOrderReads = new TreeSet(new AlignmentStartWithNoTiesComparator());
- }
-
- /**
- * Get all reads before a given read (for processing)
- *
- * @param read the original read
- * @return all reads that have alignment start before the original read.
- */
- public List getAllReadsBefore(GATKSAMRecord read) {
- List result = new LinkedList();
- GATKSAMRecord newHead = null;
-
- for (GATKSAMRecord stashedRead : outOfOrderReads) {
- if (ReadUtils.compareSAMRecords(stashedRead, read) <= 0)
- result.add(stashedRead);
- else {
- newHead = stashedRead;
- break;
- }
- }
-
- if (result.size() > 0) {
- if (result.size() == outOfOrderReads.size())
- outOfOrderReads.clear();
- else
- outOfOrderReads = new TreeSet(outOfOrderReads.tailSet(newHead));
- }
-
- return result;
- }
-
- /**
- * sends the read to the MultiSampleCompressor
- *
- * @param read the read to be compressed
- * @param knownSnpPositions the set of known SNP positions
- * @return any compressed reads that may have resulted from adding this read to the machinery (due to the sliding window)
- */
- public Iterable compress(final GATKSAMRecord read, final ObjectSortedSet knownSnpPositions) {
- return compressor.addAlignment(read, knownSnpPositions);
- }
-
- /**
- * Add a read to the stash
- *
- * @param read any read
- */
- public void add(GATKSAMRecord read) {
- outOfOrderReads.add(read);
- }
-
- /**
- * Close the stash, processing all remaining reads in order
- *
- * @param knownSnpPositions the set of known SNP positions
- * @return a list of all the reads produced by the SlidingWindow machinery)
- */
- public Iterable close(final ObjectSortedSet knownSnpPositions) {
- LinkedList result = new LinkedList();
-
- // compress all the stashed reads (in order)
- for (GATKSAMRecord read : outOfOrderReads)
- for (GATKSAMRecord compressedRead : compressor.addAlignment(read, knownSnpPositions))
- result.add(compressedRead);
-
- // output any remaining reads from the compressor
- for (GATKSAMRecord read : compressor.close(knownSnpPositions))
- result.add(read);
-
- return result;
- }
-
- /**
- * Useful debug functionality, outputs all elements in the stash
- */
- public void print() {
- int i = 1;
- System.out.println("Stash Contents:");
- for (GATKSAMRecord read : outOfOrderReads)
- System.out.println(String.format("%3d: %s %d %d", i++, read.getCigarString(), read.getAlignmentStart(), read.getAlignmentEnd()));
- System.out.println();
- }
-
-}
\ No newline at end of file
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/SingleSampleCompressor.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/SingleSampleCompressor.java
deleted file mode 100644
index 61c34b6a0..000000000
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/SingleSampleCompressor.java
+++ /dev/null
@@ -1,153 +0,0 @@
-/*
-* By downloading the PROGRAM you agree to the following terms of use:
-*
-* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
-*
-* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
-*
-* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
-* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
-* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
-*
-* 1. DEFINITIONS
-* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
-*
-* 2. LICENSE
-* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
-* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
-* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
-* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
-*
-* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
-* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
-* Copyright 2012 Broad Institute, Inc.
-* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
-* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
-*
-* 4. INDEMNIFICATION
-* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
-*
-* 5. NO REPRESENTATIONS OR WARRANTIES
-* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
-* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
-*
-* 6. ASSIGNMENT
-* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
-*
-* 7. MISCELLANEOUS
-* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
-* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
-* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
-* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
-* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
-* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
-* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
-*/
-
-package org.broadinstitute.sting.gatk.walkers.compression.reducereads;
-
-import com.google.java.contract.Ensures;
-import it.unimi.dsi.fastutil.objects.*;
-import org.broadinstitute.sting.utils.GenomeLoc;
-import org.broadinstitute.sting.utils.collections.Pair;
-import org.broadinstitute.sting.utils.sam.AlignmentStartWithNoTiesComparator;
-import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
-
-/**
- *
- * @author carneiro, depristo
- * @version 3.0
- */
-public class SingleSampleCompressor {
- final private int contextSize;
- final private int downsampleCoverage;
- final private int minMappingQuality;
- final private double minAltPValueToTriggerVariant;
- final private double minAltProportionToTriggerVariant;
- final private double minIndelProportionToTriggerVariant;
- final private int minBaseQual;
- final private ReduceReads.DownsampleStrategy downsampleStrategy;
-
- private SlidingWindow slidingWindow;
- private int slidingWindowCounter;
-
- public static Pair, CompressionStash> emptyPair = new Pair,CompressionStash>(new ObjectAVLTreeSet(), new CompressionStash());
-
- public SingleSampleCompressor(final int contextSize,
- final int downsampleCoverage,
- final int minMappingQuality,
- final double minAltPValueToTriggerVariant,
- final double minAltProportionToTriggerVariant,
- final double minIndelProportionToTriggerVariant,
- final int minBaseQual,
- final ReduceReads.DownsampleStrategy downsampleStrategy) {
- this.contextSize = contextSize;
- this.downsampleCoverage = downsampleCoverage;
- this.minMappingQuality = minMappingQuality;
- this.slidingWindowCounter = 0;
- this.minAltPValueToTriggerVariant = minAltPValueToTriggerVariant;
- this.minAltProportionToTriggerVariant = minAltProportionToTriggerVariant;
- this.minIndelProportionToTriggerVariant = minIndelProportionToTriggerVariant;
- this.minBaseQual = minBaseQual;
- this.downsampleStrategy = downsampleStrategy;
- }
-
- /**
- * Add an alignment to the compressor
- *
- * @param read the read to be added
- * @param knownSnpPositions the set of known SNP positions
- * @return any compressed reads that may have resulted from adding this read to the machinery (due to the sliding window)
- */
- public Pair, CompressionStash> addAlignment( final GATKSAMRecord read, final ObjectSortedSet knownSnpPositions ) {
- ObjectSet reads = new ObjectAVLTreeSet(new AlignmentStartWithNoTiesComparator());
- CompressionStash stash = new CompressionStash();
- int readOriginalStart = read.getUnclippedStart();
-
- // create a new window if:
- if ((slidingWindow != null) &&
- ( ( read.getReferenceIndex() != slidingWindow.getContigIndex() ) || // this is a brand new contig
- (readOriginalStart - contextSize > slidingWindow.getStopLocation()))) { // this read is too far away from the end of the current sliding window
-
- // close the current sliding window
- Pair, CompressionStash> readsAndStash = slidingWindow.close(knownSnpPositions);
- reads = readsAndStash.getFirst();
- stash = readsAndStash.getSecond();
- slidingWindow = null; // so we create a new one on the next if
- }
-
- if ( slidingWindow == null) { // this is the first read
- slidingWindow = new SlidingWindow(read.getReferenceName(), read.getReferenceIndex(), contextSize, read.getHeader(), read.getReadGroup(),
- slidingWindowCounter, minAltPValueToTriggerVariant, minAltProportionToTriggerVariant, minIndelProportionToTriggerVariant,
- minBaseQual, minMappingQuality, downsampleCoverage, downsampleStrategy, read.hasBaseIndelQualities());
- slidingWindowCounter++;
- }
-
- stash.addAll(slidingWindow.addRead(read));
- return new Pair, CompressionStash>(reads, stash);
- }
-
- /**
- * Properly closes the compressor.
- *
- * @param knownSnpPositions the set of known SNP positions
- * @return A non-null set/list of all reads generated
- */
- @Ensures("result != null")
- public Pair, CompressionStash> close(final ObjectSortedSet knownSnpPositions) {
- return (slidingWindow != null) ? slidingWindow.close(knownSnpPositions) : emptyPair;
- }
-
- /**
- * Finalizes current variant regions.
- *
- * @param knownSnpPositions the set of known SNP positions
- * @return A non-null set/list of all reads generated
- */
- @Ensures("result != null")
- public ObjectSet closeVariantRegions(final CompressionStash regions, final ObjectSortedSet knownSnpPositions) {
- return slidingWindow == null ? ObjectSets.EMPTY_SET : slidingWindow.closeVariantRegions(regions, knownSnpPositions);
- }
-
-}
-
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/SlidingWindow.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/SlidingWindow.java
deleted file mode 100644
index d5aa8f944..000000000
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/SlidingWindow.java
+++ /dev/null
@@ -1,1110 +0,0 @@
-/*
-* By downloading the PROGRAM you agree to the following terms of use:
-*
-* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
-*
-* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
-*
-* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
-* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
-* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
-*
-* 1. DEFINITIONS
-* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
-*
-* 2. LICENSE
-* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
-* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
-* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
-* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
-*
-* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
-* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
-* Copyright 2012 Broad Institute, Inc.
-* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
-* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
-*
-* 4. INDEMNIFICATION
-* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
-*
-* 5. NO REPRESENTATIONS OR WARRANTIES
-* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
-* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
-*
-* 6. ASSIGNMENT
-* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
-*
-* 7. MISCELLANEOUS
-* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
-* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
-* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
-* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
-* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
-* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
-* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
-*/
-
-package org.broadinstitute.sting.gatk.walkers.compression.reducereads;
-
-import com.google.java.contract.Ensures;
-import com.google.java.contract.Requires;
-import it.unimi.dsi.fastutil.bytes.Byte2IntArrayMap;
-import it.unimi.dsi.fastutil.bytes.Byte2IntMap;
-import it.unimi.dsi.fastutil.objects.*;
-import net.sf.samtools.CigarElement;
-import net.sf.samtools.CigarOperator;
-import net.sf.samtools.SAMFileHeader;
-import org.broadinstitute.sting.gatk.downsampling.ReservoirDownsampler;
-import org.broadinstitute.sting.utils.BaseUtils;
-import org.broadinstitute.sting.utils.GenomeLoc;
-import org.broadinstitute.sting.utils.UnvalidatingGenomeLoc;
-import org.broadinstitute.sting.utils.collections.Pair;
-import org.broadinstitute.sting.utils.sam.AlignmentStartWithNoTiesComparator;
-import org.broadinstitute.sting.utils.sam.GATKSAMReadGroupRecord;
-import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
-import org.broadinstitute.sting.utils.sam.ReadUtils;
-
-import java.util.*;
-
-
-/**
- * Created by IntelliJ IDEA.
- * User: roger
- * Date: 8/3/11
- * Time: 2:24 PM
- */
-public class SlidingWindow {
-
- // Sliding Window data
- final protected PriorityQueue readsInWindow;
- final protected LinkedList windowHeader;
- protected int contextSize; // the largest context size (between mismatches and indels)
- protected String contig;
- protected int contigIndex;
- protected SAMFileHeader samHeader;
- protected GATKSAMReadGroupRecord readGroupAttribute;
- protected int downsampleCoverage;
-
- // Running consensus data
- protected int consensusCounter;
- protected String consensusReadName;
-
- // Filtered Data Consensus data
- protected int filteredDataConsensusCounter;
- protected String filteredDataReadName;
-
- // Additional parameters
- protected double MIN_ALT_PVALUE_TO_TRIGGER_VARIANT; // pvalue has to be greater than this value to trigger variant region due to mismatches
- protected double MIN_ALT_PROPORTION_TO_TRIGGER_VARIANT; // proportion has to be greater than this value to trigger variant region due to mismatches
- protected double MIN_INDEL_BASE_PROPORTION_TO_TRIGGER_VARIANT; // proportion has to be greater than this value to trigger variant region due to deletions
- protected int MIN_BASE_QUAL_TO_COUNT; // qual has to be greater than or equal to this value
- protected int MIN_MAPPING_QUALITY;
-
- protected ReduceReads.DownsampleStrategy downsampleStrategy;
- private boolean hasIndelQualities;
-
- private static CompressionStash emptyRegions = new CompressionStash();
-
- /**
- * The types of synthetic reads
- */
- protected enum ConsensusType {
- POSITIVE_CONSENSUS,
- NEGATIVE_CONSENSUS,
- FILTERED
- }
-
- public int getStopLocation() {
- return getStopLocation(windowHeader);
- }
-
- private int getStopLocation(final LinkedList header) {
- return header.isEmpty() ? -1 : header.peekLast().getLocation();
- }
-
- public String getContig() {
- return contig;
- }
-
- public int getContigIndex() {
- return contigIndex;
- }
-
- public int getStartLocation(final LinkedList header) {
- return header.isEmpty() ? -1 : header.peek().getLocation();
- }
-
- // for testing only
- protected SlidingWindow(final String contig, final int contigIndex, final int startLocation) {
- this.contig = contig;
- this.contigIndex = contigIndex;
-
- contextSize = 10;
-
- this.windowHeader = new LinkedList<>();
- windowHeader.addFirst(new HeaderElement(startLocation));
- this.readsInWindow = new PriorityQueue<>(100, new Comparator() {
- @Override
- public int compare(GATKSAMRecord read1, GATKSAMRecord read2) {
- return read1.getSoftEnd() - read2.getSoftEnd();
- }
- });
- }
-
- public SlidingWindow(final String contig, final int contigIndex, final int contextSize, final SAMFileHeader samHeader,
- final GATKSAMReadGroupRecord readGroupAttribute, final int windowNumber,
- final double minAltPValueToTriggerVariant, final double minAltProportionToTriggerVariant, final double minIndelProportionToTriggerVariant,
- final int minBaseQual, final int minMappingQuality, final int downsampleCoverage,
- final ReduceReads.DownsampleStrategy downsampleStrategy, final boolean hasIndelQualities) {
- this.contextSize = contextSize;
- this.downsampleCoverage = downsampleCoverage;
-
- this.MIN_ALT_PVALUE_TO_TRIGGER_VARIANT = minAltPValueToTriggerVariant;
- this.MIN_ALT_PROPORTION_TO_TRIGGER_VARIANT = minAltProportionToTriggerVariant;
- this.MIN_INDEL_BASE_PROPORTION_TO_TRIGGER_VARIANT = minIndelProportionToTriggerVariant;
- this.MIN_BASE_QUAL_TO_COUNT = minBaseQual;
- this.MIN_MAPPING_QUALITY = minMappingQuality;
-
- this.windowHeader = new LinkedList<>();
- this.readsInWindow = new PriorityQueue<>(1000, new Comparator() {
- @Override
- public int compare(GATKSAMRecord read1, GATKSAMRecord read2) {
- return read1.getSoftEnd() - read2.getSoftEnd();
- }
- });
-
- this.contig = contig;
- this.contigIndex = contigIndex;
- this.samHeader = samHeader;
- this.readGroupAttribute = readGroupAttribute;
-
- this.consensusCounter = 0;
- this.consensusReadName = "Consensus-" + windowNumber + "-";
-
- this.filteredDataConsensusCounter = 0;
- this.filteredDataReadName = "Filtered-" + windowNumber + "-";
-
- this.downsampleStrategy = downsampleStrategy;
- this.hasIndelQualities = hasIndelQualities;
- }
-
- /**
- * Add a read to the sliding window and slides the window accordingly.
- *
- * Reads are assumed to be in order, therefore, when a read is added the sliding window can
- * assume that no more reads will affect read.getUnclippedStart() - contextSizeMismatches. The window
- * slides forward to that position and returns all reads that may have been finalized in the
- * sliding process.
- *
- * @param read the read
- * @return a non-null list of reads (in the CompressionStash) that have been finished by sliding the window.
- */
- @Requires({"read != null"})
- @Ensures("result != null")
- public CompressionStash addRead(GATKSAMRecord read) {
- addToHeader(windowHeader, read); // update the window header counts
- // no need to track low mapping quality reads
- if ( read.getMappingQuality() >= MIN_MAPPING_QUALITY )
- readsInWindow.add(read); // add read to sliding reads
- return slideWindow(read.getUnclippedStart());
- }
-
- /**
- * Returns the next complete (or incomplete if closeLastRegion is true) variant region between 'from' (inclusive) and 'to' (exclusive)
- * but converted to global coordinates.
- *
- * @param from beginning window header index of the search window (inclusive) in local (to the windowHeader) coordinates
- * @param to end window header index of the search window (exclusive) in local (to the windowHeader) coordinates
- * @param variantSite boolean array with true marking variant regions
- * @param closeLastRegion if the last index is variant (so it's an incomplete region), should we close (and return as an interval) the location or ignore it?
- * @return null if nothing is variant, start/stop if there is a complete variant region, start/-1 if there is an incomplete variant region. All coordinates returned are global.
- */
- @Requires({"from >= 0", "from <= to", "to <= variantSite.length"})
- private FinishedGenomeLoc findNextVariantRegion(int from, int to, boolean[] variantSite, boolean closeLastRegion) {
- boolean foundStart = false;
- final int windowHeaderStart = getStartLocation(windowHeader);
- int variantRegionStartIndex = 0;
- for (int i=from; i= 0", "from <= to", "to <= variantSite.length"})
- @Ensures("result != null")
- protected CompressionStash findVariantRegions(int from, int to, boolean[] variantSite, boolean closeLastRegion) {
- final int windowHeaderStart = getStartLocation(windowHeader);
-
- CompressionStash regions = new CompressionStash();
- int index = from;
- while(index < to) {
- // returns results in global coordinates
- FinishedGenomeLoc result = findNextVariantRegion(index, to, variantSite, closeLastRegion);
- if (result == null)
- break;
-
- regions.add(result);
- if (!result.isFinished())
- break;
-
- index = result.getStop() - windowHeaderStart + 1; // go back to local coordinates
- }
- return regions;
- }
-
- /**
- * Determines if the window can be slid given the new incoming read.
- *
- * We check from the start of the window to the (unclipped) start of the new incoming read if there
- * is any variant.
- * If there are variant sites, we check if it's time to close the variant region.
- *
- * @param incomingReadUnclippedStart the incoming read's start position. Must be the unclipped start!
- * @return all reads that have fallen to the left of the sliding window after the slide
- */
- protected CompressionStash slideWindow(final int incomingReadUnclippedStart) {
- final int windowHeaderStartLocation = getStartLocation(windowHeader);
- CompressionStash regions = emptyRegions;
- boolean forceClose = true;
-
- if (incomingReadUnclippedStart - contextSize > windowHeaderStartLocation) {
- markSites(incomingReadUnclippedStart);
- int readStartHeaderIndex = incomingReadUnclippedStart - windowHeaderStartLocation;
- int breakpoint = Math.max(readStartHeaderIndex - contextSize - 1, 0); // this is the limit of what we can close/send to consensus (non-inclusive)
-
- regions = findVariantRegions(0, breakpoint, markedSites.getVariantSiteBitSet(), !forceClose);
- }
-
- while (!readsInWindow.isEmpty() && readsInWindow.peek().getSoftEnd() < windowHeaderStartLocation) {
- readsInWindow.poll();
- }
-
- return regions;
- }
-
-
- protected final class MarkedSites {
-
- private boolean[] siteIsVariant = new boolean[0];
- private int startLocation = 0;
-
- public MarkedSites() {}
-
- public boolean[] getVariantSiteBitSet() { return siteIsVariant; }
-
- protected int getStartLocation() { return startLocation; }
-
- /**
- * Updates the variant site bitset given the new startlocation and size of the region to mark.
- *
- * @param newStartLocation the new start location of the bitset
- * @param sizeOfRegion the new size of the region to be represented
- *
- * @return the end position (newStartLocation + index) of the region marked by this method; the calling method is responsible for the remainder.
- */
- public int updateRegion(final int newStartLocation, final int sizeOfRegion) {
- int lastPositionMarked = sizeOfRegion;
-
- // if this is the first time we set the array and we can't reuse anything, just create a new array from scratch
- if ( newStartLocation >= this.startLocation + siteIsVariant.length || newStartLocation < this.startLocation ) {
- siteIsVariant = new boolean[sizeOfRegion];
- lastPositionMarked = 0;
- }
- // if the dimensions change, copy what we can and continue
- else if ( newStartLocation != this.startLocation || sizeOfRegion != siteIsVariant.length ) {
- final boolean[] tempArray = new boolean[sizeOfRegion];
- final int differenceInStartPositions = newStartLocation - this.startLocation;
- lastPositionMarked = Math.min(siteIsVariant.length - differenceInStartPositions, sizeOfRegion);
- System.arraycopy(siteIsVariant, differenceInStartPositions, tempArray, 0, lastPositionMarked);
- siteIsVariant = null; // explicitly allow garbage collection
- siteIsVariant = tempArray;
- }
-
- this.startLocation = newStartLocation;
-
- return lastPositionMarked + newStartLocation;
- }
- }
-
- private final MarkedSites markedSites = new MarkedSites();
-
- /**
- * returns the MarkedSites object so that it can be tested after adding data to the Sliding Window
- *
- * @return the Marked Sites object used by this Sliding Window
- */
- protected MarkedSites getMarkedSitesForTesting() { return markedSites; }
-
- /**
- * returns an array marked with variant and non-variant regions (it uses markVariantRegion to make the marks)
- *
- * @param stop check the window from start to stop (not-inclusive); given in global coordinates
- */
- protected void markSites(final int stop) {
-
- final int windowHeaderStartLocation = getStartLocation(windowHeader);
- final int sizeOfMarkedRegion = stop - windowHeaderStartLocation + contextSize + 1;
-
- // copy over as many bits as we can from the previous calculation. Note that we can't trust the
- // last (contextSize - 1) worth of bits because we may not have actually looked at variant regions there.
- final int lastPositionMarked = markedSites.updateRegion(windowHeaderStartLocation, sizeOfMarkedRegion) - contextSize - 1;
- final int locationToProcess = Math.max(windowHeaderStartLocation, Math.min(lastPositionMarked, stop - contextSize));
-
- final ListIterator headerElementIterator = windowHeader.listIterator(locationToProcess - windowHeaderStartLocation);
-
- // process a contextSize worth of region from scratch in case there's a variant there
- for (int i = locationToProcess; i < stop; i++) {
- if (headerElementIterator.hasNext()) {
- HeaderElement headerElement = headerElementIterator.next();
-
- if (headerElement.isVariant(MIN_ALT_PVALUE_TO_TRIGGER_VARIANT, MIN_ALT_PROPORTION_TO_TRIGGER_VARIANT, MIN_INDEL_BASE_PROPORTION_TO_TRIGGER_VARIANT))
- markVariantRegion(i - windowHeaderStartLocation);
-
- } else
- break;
- }
- }
-
- /**
- * Marks the sites around the variant site (as true)
- *
- * @param variantSiteLocation the location where a variant site was found
- */
- protected void markVariantRegion(final int variantSiteLocation) {
- int from = (variantSiteLocation < contextSize) ? 0 : variantSiteLocation - contextSize;
- int to = (variantSiteLocation + contextSize + 1 > markedSites.getVariantSiteBitSet().length) ? markedSites.getVariantSiteBitSet().length - 1 : variantSiteLocation + contextSize;
- markRegionAs(from, to, true);
- }
-
- /**
- * Marks the sites around the variant site (as true)
- *
- * @param from the start index (inclusive) to mark
- * @param to the end index (inclusive) to mark
- * @param isVariant mark the region with this boolean value
- */
- private void markRegionAs(final int from, final int to, final boolean isVariant) {
- for (int i = from; i <= to; i++)
- markedSites.getVariantSiteBitSet()[i] = isVariant;
- }
-
- /**
- * Adds bases to the running consensus
- *
- * If adding a sequence with gaps, it will finalize multiple consensus reads and keep the last running consensus
- *
- * @param header the header to use
- * @param start the first header index to add to consensus
- * @param end the first header index NOT TO add to consensus
- * @param consensusType the consensus type to use
- * @return a non-null list of consensus reads generated by this call. Empty list if no consensus was generated.
- */
- @Requires({"start >= 0 && (end >= start || end == 0)"})
- @Ensures("result != null")
- protected ObjectArrayList addToSyntheticReads(final LinkedList header, final int start, final int end, final ConsensusType consensusType) {
- final ObjectArrayList reads = new ObjectArrayList<>();
-
- SyntheticRead consensus = null;
- final ListIterator headerElementIterator = header.listIterator(start);
- boolean wasInConsensus = false;
-
- for ( int currentPosition = start; currentPosition < end; currentPosition++ ) {
-
- if ( ! headerElementIterator.hasNext() )
- throw new IllegalStateException(String.format("Requested to add to synthetic reads a region that contains no header element at index: %d - %d / %d", start, windowHeader.size(), end));
- final HeaderElement headerElement = headerElementIterator.next();
-
- if ( headerElement.hasConsensusData(consensusType) ) {
- wasInConsensus = true;
-
- // add to running consensus
- if ( consensus == null )
- consensus = createNewConsensus(consensusType, headerElement.getLocation());
-
- genericAddBaseToConsensus(consensus, headerElement.getBaseCounts(consensusType));
-
- } else {
-
- // add any outstanding consensus data
- if ( wasInConsensus ) {
- reads.addAll(finalizeAndAdd(consensus, consensusType));
- consensus = null;
- }
-
- wasInConsensus = false;
- }
- }
-
- // add any outstanding consensus data
- reads.addAll(finalizeAndAdd(consensus, consensusType));
-
- return reads;
- }
-
- private SyntheticRead createNewConsensus(final ConsensusType consensusType, final int start) {
- if ( consensusType == ConsensusType.FILTERED )
- return new SyntheticRead(samHeader, readGroupAttribute, contig, contigIndex, filteredDataReadName + filteredDataConsensusCounter++, start, hasIndelQualities, SyntheticRead.StrandType.STRANDLESS);
- return new SyntheticRead(samHeader, readGroupAttribute, contig, contigIndex, consensusReadName + consensusCounter++, start, hasIndelQualities, consensusType == ConsensusType.POSITIVE_CONSENSUS ? SyntheticRead.StrandType.POSITIVE : SyntheticRead.StrandType.NEGATIVE);
- }
-
- /**
- * Finalizes a synthetic read.
- *
- * @param consensus the consensus to finalize
- * @param type the synthetic reads you want to close
- * @return a possibly empty list of GATKSAMRecords generated by finalizing the synthetic reads
- */
- private ObjectArrayList finalizeAndAdd(final SyntheticRead consensus, final ConsensusType type) {
-
- final ObjectArrayList list = new ObjectArrayList<>();
-
- final GATKSAMRecord read;
- if ( type == ConsensusType.FILTERED )
- read = finalizeFilteredDataConsensus(consensus);
- else
- read = finalizeRunningConsensus(consensus);
-
- if ( read != null )
- list.add(read);
-
- return list;
- }
-
- /**
- * Generic accessor to add base and qualities to a synthetic read
- *
- * @param syntheticRead the synthetic read to add to
- * @param baseCounts the base counts object in the header element
- */
- private void genericAddBaseToConsensus(final SyntheticRead syntheticRead, final BaseAndQualsCounts baseCounts) {
- final BaseIndex base = baseCounts.baseIndexWithMostProbability();
- final int count = baseCounts.countOfBase(base);
- final byte qual = baseCounts.averageQualsOfBase(base);
- final byte insQual = baseCounts.averageInsertionQualsOfBase(base);
- final byte delQual = baseCounts.averageDeletionQualsOfBase(base);
- syntheticRead.add(base, count, qual, insQual, delQual, baseCounts.getRMS());
- }
-
- /**
- * Method to compress a variant region and return the associated reduced reads
- *
- * @param start the first window header index in the variant region (inclusive)
- * @param stop the last window header index of the variant region (inclusive)
- * @param knownSnpPositions the set of known SNPs used to determine whether to allow polyploid consensus creation here; can be null (to allow polyploid consensus anywhere)
- * @return a non-null object representing all reads contained in the variant region
- */
- @Requires({"start >= 0 && (stop >= start || stop == 0)"})
- @Ensures("result != null")
- protected CloseVariantRegionResult compressVariantRegion(final int start, final int stop, final ObjectSortedSet knownSnpPositions) {
- final CloseVariantRegionResult allReads = new CloseVariantRegionResult(stop);
-
- // Try to compress into a polyploid consensus
- // Optimization: don't bother if there are no known SNPs here
- final int hetRefPosition = (knownSnpPositions != null && knownSnpPositions.isEmpty()) ? -1 : findSinglePolyploidCompressiblePosition(start, stop);
-
- // Note that using the hetRefPosition protects us from trying to compress variant regions that are created by
- // insertions (which we don't want because we can't confirm that they represent the same allele).
- // Also, we only allow polyploid consensus creation at known sites if provided.
- if ( hetRefPosition != -1 && matchesKnownPosition(windowHeader.get(hetRefPosition).getLocation(), knownSnpPositions) ) {
- // try to create the polyploid consensus
- allReads.reads.addAll(createPolyploidConsensus(hetRefPosition));
- allReads.stopPerformed = hetRefPosition; // we stopped at the het position
- }
- // if we can't create a polyploid consensus here, return all reads that overlap the variant region and remove them
- // from the window header entirely; also remove all reads preceding the variant region (since they will be output
- // as consensus right after compression)
- else {
- final int refStart = windowHeader.get(start).getLocation();
- final int refStop = windowHeader.get(stop).getLocation();
-
- final ObjectList toRemoveFromWindow = new ObjectArrayList<>();
- final ObjectList toEmit = new ObjectArrayList<>();
- for ( final GATKSAMRecord read : readsInWindow ) {
- if ( read.getSoftStart() <= refStop ) {
- if ( read.getAlignmentEnd() >= refStart ) {
- toEmit.add(read);
- removeFromHeader(windowHeader, read);
- }
- toRemoveFromWindow.add(read);
- }
- }
-
- // remove all used reads
- for ( final GATKSAMRecord read : toRemoveFromWindow )
- readsInWindow.remove(read);
-
- // down-sample the unreduced reads if needed
- allReads.reads.addAll(downsampleCoverage > 0 ? downsampleVariantRegion(toEmit) : toEmit);
- }
-
- return allReads;
- }
-
- /**
- * Determines whether the given position match one of the known sites
- *
- * @param targetPosition the position of the het site
- * @param knownSnpPositions the set of known SNPs used to determine whether to allow polyploid consensus creation here; can be null (to allow polyploid consensus anywhere)
- * @return true if the targetPosition matches a known SNP position, false otherwise
- */
- @Requires({"targetPosition >= 1 && knownSnpPositions != null"})
- protected boolean matchesKnownPosition(final int targetPosition, final ObjectSortedSet knownSnpPositions) {
- final GenomeLoc targetLoc = new UnvalidatingGenomeLoc(contig, contigIndex, targetPosition, targetPosition);
- return knownSnpPositions == null || knownSnpPositions.contains(targetLoc);
- }
-
- /*
- * Finds the het variant position located within start and stop (inclusive) if one exists.
- *
- * @param start the first header index in the region to check (inclusive)
- * @param stop the last header index of the region to check (inclusive)
- * @return the window header index of the single het position or -1 if either none or more than one exists
- */
- @Requires("start >= 0 && (stop >= start || stop == 0)")
- protected int findSinglePolyploidCompressiblePosition(final int start, final int stop) {
- int hetRefPosition = -1;
-
- for ( int i = start; i <= stop; i++ ) {
-
- final int nAlleles = windowHeader.get(i).getNumberOfBaseAlleles(MIN_ALT_PVALUE_TO_TRIGGER_VARIANT, MIN_ALT_PROPORTION_TO_TRIGGER_VARIANT);
-
- // we will only work on diploid non-indel cases because we just don't want to handle/test other scenarios
- if ( nAlleles > 2 || nAlleles == -1 )
- return -1;
-
- if ( nAlleles == 2 ) {
-
- // make sure that there is only 1 site in the region that contains more than one allele
- if ( hetRefPosition != -1 )
- return -1;
-
- hetRefPosition = i;
- }
- }
-
- return hetRefPosition;
- }
-
- /*
- * Checks whether there's a position in the header with a significant number of softclips or a variant.
- *
- * @param header the window header to examine
- * @param positionToSkip the global position to skip in the examination (use negative number if you don't want to make use of this argument)
- * @return true if there exists a position with significant softclips, false otherwise
- */
- @Requires("header != null")
- protected boolean hasPositionWithSignificantSoftclipsOrVariant(final List header, final int positionToSkip) {
-
- for ( final HeaderElement headerElement : header ) {
-
- if ( headerElement.getLocation() == positionToSkip )
- continue;
-
- if ( headerElement.hasSignificantSoftclips(MIN_ALT_PVALUE_TO_TRIGGER_VARIANT, MIN_ALT_PROPORTION_TO_TRIGGER_VARIANT) ||
- headerElement.getNumberOfBaseAlleles(MIN_ALT_PVALUE_TO_TRIGGER_VARIANT, MIN_ALT_PROPORTION_TO_TRIGGER_VARIANT) != 1 )
- return true;
- }
-
- return false;
- }
-
- /**
- * Finalizes a variant region, any adjacent synthetic reads.
- *
- * @param start the first window header index in the variant region (inclusive)
- * @param stop the last window header index of the variant region (inclusive)
- * @param knownSnpPositions the set of known SNPs used to determine whether to allow polyploid consensus creation here; can be null (to allow polyploid consensus anywhere)
- * @return a non-null object representing all reads contained in the variant region plus any adjacent synthetic reads
- */
- @Requires({"start >= 0 && (stop >= start || stop == 0)"})
- @Ensures("result != null")
- protected CloseVariantRegionResult closeVariantRegion(final int start, final int stop, final ObjectSortedSet knownSnpPositions) {
- final CloseVariantRegionResult allReads = compressVariantRegion(start, stop, knownSnpPositions);
- allReads.reads.addAll(addAllSyntheticReadTypes(0, allReads.stopPerformed + 1));
- return allReads;
- }
-
- /**
- * Adds reads for all possible strands (positive, negative, filtered) from the global windowHeader object
- *
- * @param start the start position (inclusive)
- * @param end the end position (exclusive)
- * @return non-null but possibly empty array list with reduced reads
- */
- private ObjectArrayList addAllSyntheticReadTypes(final int start, final int end) {
- final ObjectArrayList reads = new ObjectArrayList<>();
- reads.addAll(addToSyntheticReads(windowHeader, start, end, ConsensusType.POSITIVE_CONSENSUS));
- reads.addAll(addToSyntheticReads(windowHeader, start, end, ConsensusType.NEGATIVE_CONSENSUS));
- reads.addAll(addToSyntheticReads(windowHeader, start, end, ConsensusType.FILTERED));
- return reads;
- }
-
- /*
- * @see #closeVariantRegions(CompressionStash, ObjectSortedSet, boolean) with forceCloseFullRegions set to false
- */
- public ObjectSet closeVariantRegions(final CompressionStash regions, final ObjectSortedSet knownSnpPositions) {
- return closeVariantRegions(regions, knownSnpPositions, false);
- }
-
- private static final class CloseVariantRegionResult {
- final private ObjectList reads = new ObjectArrayList<>();
- private int stopPerformed;
-
- public CloseVariantRegionResult(final int stopPerformed) { this.stopPerformed = stopPerformed; }
- }
-
- /*
- * Finalizes the list of regions requested (and any regions preceding them)
- *
- * @param regions the list of regions to finalize
- * @param knownSnpPositions the set of known SNP positions; can be null (to allow polyploid consensus anywhere)
- * @param forceCloseFullRegions if true, requires this method to make sure all regions are fully closed; otherwise, we may decide not to close up to the very end (e.g. during het compression)
- * @return a non-null set of reduced reads representing the finalized regions
- */
- public ObjectSet closeVariantRegions(final CompressionStash regions, final ObjectSortedSet knownSnpPositions, final boolean forceCloseFullRegions) {
- final ObjectAVLTreeSet allReads = new ObjectAVLTreeSet<>(new AlignmentStartWithNoTiesComparator());
- if ( !regions.isEmpty() ) {
-
- int windowHeaderStart = getStartLocation(windowHeader);
- HeaderElement lastCleanedElement = null;
-
- for ( final GenomeLoc region : regions ) {
- if (((FinishedGenomeLoc)region).isFinished() && region.getContig().equals(contig) && region.getStart() >= windowHeaderStart && region.getStop() < windowHeaderStart + windowHeader.size()) {
- final int start = region.getStart() - windowHeaderStart;
- int stop = region.getStop() - windowHeaderStart;
-
- // make sure the bitset is complete given the region (it might not be in multi-sample mode)
- if ( region.getStop() > markedSites.getStartLocation() + markedSites.getVariantSiteBitSet().length - 1 )
- markSites(region.getStop());
-
- CloseVariantRegionResult closeVariantRegionResult = closeVariantRegion(start, stop, knownSnpPositions);
- allReads.addAll(closeVariantRegionResult.reads);
-
- // check whether we didn't close the whole region that was requested
- if ( stop > 0 && closeVariantRegionResult.stopPerformed < stop ) {
- // we should update the variant sites bitset because the context size's worth of bases after the variant position are no longer "variant"
- markRegionAs(closeVariantRegionResult.stopPerformed + 1, stop, false);
-
- // if the calling method said that it didn't care then we are okay so update the stop
- if ( !forceCloseFullRegions ) {
- stop = closeVariantRegionResult.stopPerformed;
- }
- // otherwise, we need to forcibly push the stop that we originally requested
- else {
- while ( closeVariantRegionResult.stopPerformed < stop ) {
- // first clean up used header elements so they don't get reused
- for ( int i = 0; i <= closeVariantRegionResult.stopPerformed; i++ )
- windowHeader.remove();
- stop -= (closeVariantRegionResult.stopPerformed + 1);
-
- closeVariantRegionResult = closeVariantRegion(0, stop, knownSnpPositions);
- allReads.addAll(closeVariantRegionResult.reads);
- }
- }
- }
-
- // We need to clean up the window header elements up until the end of the requested region so that they don't get used for future regions.
- // Note that this cleanup used to happen outside the above for-loop, but that was causing an occasional doubling of the reduced reads
- // (in the case where there are multiple regions to close we'd reuse the reads for each region).
- if ( stop >= 0 ) {
- for ( int i = 0; i < stop; i++ )
- windowHeader.remove();
- lastCleanedElement = windowHeader.remove();
- windowHeaderStart = getStartLocation(windowHeader);
- }
- }
- }
-
- // we need to keep the last element of the last cleaned region in the event that the following element has a read that starts with an insertion.
- if ( lastCleanedElement != null && lastCleanedElement.hasInsertionToTheRight() )
- windowHeader.addFirst(new HeaderElement(lastCleanedElement.getLocation(), lastCleanedElement.numInsertionsToTheRight()));
- }
-
- return allReads;
- }
-
- /**
- * Downsamples a variant region to the downsample coverage of the sliding window.
- *
- * It will use the downsampling strategy defined by the SlidingWindow
- *
- * @param allReads a non-null list of reads to select from (all reads that cover the window)
- * @return a non-null list of reads selected by the downsampler to cover the window to at least the desired coverage
- */
- @Requires({"allReads != null"})
- @Ensures("result != null")
- protected ObjectList downsampleVariantRegion(final ObjectList allReads) {
- int nReads = allReads.size();
- if (nReads == 0)
- return allReads;
-
- if (downsampleCoverage >= nReads)
- return allReads;
-
- ReservoirDownsampler downsampler = new ReservoirDownsampler<>(downsampleCoverage);
- downsampler.submit(allReads);
- return new ObjectArrayList<>(downsampler.consumeFinalizedItems());
- }
-
-
- /**
- * Properly closes a Sliding Window, finalizing all consensus and variant
- * regions that still exist regardless of being able to fulfill the
- * context size requirement in the end.
- *
- * @param knownSnpPositions the set of known SNP positions; can be null (to allow polyploid consensus anywhere)
- * @return A non-null set/list of all reads generated
- */
- @Ensures("result != null")
- public Pair, CompressionStash> close(final ObjectSortedSet knownSnpPositions) {
- // mark variant regions
- ObjectSet finalizedReads = new ObjectAVLTreeSet<>(new AlignmentStartWithNoTiesComparator());
- CompressionStash regions = new CompressionStash();
-
- if (!windowHeader.isEmpty()) {
- markSites(getStopLocation(windowHeader) + 1);
- regions = findVariantRegions(0, windowHeader.size(), markedSites.getVariantSiteBitSet(), true);
- finalizedReads = closeVariantRegions(regions, knownSnpPositions, true);
-
- if (!windowHeader.isEmpty())
- finalizedReads.addAll(addAllSyntheticReadTypes(0, windowHeader.size()));
- }
-
- return new Pair<>(finalizedReads, regions);
- }
-
- /**
- * generates the SAM record for the running consensus read and resets it (to null)
- *
- * @param runningConsensus the consensus to finalize
- * @return the read contained in the running consensus or null
- */
- protected GATKSAMRecord finalizeRunningConsensus(final SyntheticRead runningConsensus) {
- GATKSAMRecord finalizedRead = null;
-
- if ( runningConsensus != null ) {
- if ( runningConsensus.size() > 0 )
- finalizedRead = runningConsensus.close();
- else
- consensusCounter--;
- }
-
- return finalizedRead;
- }
-
- /**
- * generates the SAM record for the filtered data consensus and resets it (to null)
- *
- * @param filteredDataConsensus the consensus to finalize
- * @return the read contained in the running consensus or null
- */
- protected GATKSAMRecord finalizeFilteredDataConsensus(final SyntheticRead filteredDataConsensus) {
- GATKSAMRecord finalizedRead = null;
- if (filteredDataConsensus != null) {
- if (filteredDataConsensus.size() > 0)
- finalizedRead = filteredDataConsensus.close();
- else
- filteredDataConsensusCounter--;
- }
- return finalizedRead;
- }
-
- // define this so that we can use Java generics below
- private final static class HeaderElementList extends LinkedList {}
-
- private final static class SingleStrandConsensusData {
- final HeaderElementList consensus = new HeaderElementList();
- final ObjectList reads = new ObjectArrayList<>();
- }
-
- /**
- * Finalizes a variant region - and any adjacent synthetic reads - for point mutations (indel sites are not
- * supported) with polyploid compression.
- *
- * @param hetRefPosition window header index of the het site; MUST NOT BE AN INDEL SITE!
- * @return a non-null list of all reads contained in the variant region as a polyploid consensus
- */
- @Requires({"start >= 0 && (stop >= start || stop == 0)"})
- @Ensures({"result != null"})
- protected ObjectList createPolyploidConsensus(final int hetRefPosition) {
- // we will create two (positive strand, negative strand) headers for each haplotype
- final SingleStrandConsensusData[] headersPosStrand = new SingleStrandConsensusData[2];
- final SingleStrandConsensusData[] headersNegStrand = new SingleStrandConsensusData[2];
-
- final int globalHetRefPosition = windowHeader.get(hetRefPosition).getLocation();
-
- // initialize the mapping from base (allele) to header
- final Byte2IntMap alleleHeaderMap = new Byte2IntArrayMap(2);
- alleleHeaderMap.defaultReturnValue(-1);
- for ( final BaseIndex allele : windowHeader.get(hetRefPosition).getAlleles(MIN_ALT_PVALUE_TO_TRIGGER_VARIANT, MIN_ALT_PROPORTION_TO_TRIGGER_VARIANT) ) {
- final int currentIndex = alleleHeaderMap.size();
- if ( currentIndex > 1 )
- throw new IllegalStateException("There are more than 2 alleles present when creating a diploid consensus");
-
- alleleHeaderMap.put(allele.b, currentIndex);
- headersPosStrand[currentIndex] = new SingleStrandConsensusData();
- headersNegStrand[currentIndex] = new SingleStrandConsensusData();
- }
-
- // sanity check that we saw 2 alleles
- if ( alleleHeaderMap.size() != 2 )
- throw new IllegalStateException("We expected to see 2 alleles when creating a diploid consensus but saw " + alleleHeaderMap.size());
-
- final ObjectList readsToRemove = new ObjectArrayList<>();
-
- for ( final GATKSAMRecord read : readsInWindow ) {
-
- // if the read falls after the het position, just skip it for now (we'll get to it later)
- if ( read.getSoftStart() > globalHetRefPosition )
- continue;
-
- // remove all other reads from the read cache since we're going to use them here
- readsToRemove.add(read);
-
- // if the read falls before the het position or has low MQ, we don't need to look at it
- if ( read.getSoftEnd() < globalHetRefPosition || read.getMappingQuality() < MIN_MAPPING_QUALITY)
- continue;
-
- // remove all spanning reads from the consensus header since we're going to incorporate them into a consensus here instead
- removeFromHeader(windowHeader, read);
-
- // where on the read is the het position?
- final int readPosOfHet = ReadUtils.getReadCoordinateForReferenceCoordinate(read, globalHetRefPosition, ReadUtils.ClippingTail.LEFT_TAIL);
-
- // this is safe because indels are not supported
- final byte base = read.getReadBases()[readPosOfHet];
-
- // check which allele this read represents
- final int allele = alleleHeaderMap.get(base);
-
- // ignore the read if it represents a base that's not part of the consensus
- if ( allele != -1 ) {
- // add to the appropriate polyploid header
- final SingleStrandConsensusData header = read.getReadNegativeStrandFlag() ? headersNegStrand[allele] : headersPosStrand[allele];
- header.reads.add(read);
- addToHeader(header.consensus, read);
- }
- }
-
- for ( final GATKSAMRecord read : readsToRemove )
- readsInWindow.remove(read);
-
- // create the polyploid synthetic reads if we can
- final ObjectList hetReads = new ObjectArrayList<>();
-
- // sanity check that no new "variant region" exists on just a single consensus strand due to softclips
- // or multi-allelic sites now that we've broken everything out into their component parts. if one does
- // exist then we need to back out the consensus for that strand only.
- for ( final SingleStrandConsensusData header : headersPosStrand ) {
- if ( hasPositionWithSignificantSoftclipsOrVariant(header.consensus, globalHetRefPosition) )
- hetReads.addAll(header.reads);
- else
- finalizeHetConsensus(header.consensus, false, hetReads);
- }
- for ( final SingleStrandConsensusData header : headersNegStrand ) {
- if ( hasPositionWithSignificantSoftclipsOrVariant(header.consensus, globalHetRefPosition) )
- hetReads.addAll(header.reads);
- else
- finalizeHetConsensus(header.consensus, true, hetReads);
- }
-
- return hetReads;
- }
-
- /*
- * Finalizes a particular het consensus for the given header representation
- *
- * @param header the list of header elements representing the header for the consensus
- * @param isNegativeStrand does this header represent reads on the negative strand?
- * @param result list in which to store results
- */
- protected void finalizeHetConsensus(final LinkedList header, final boolean isNegativeStrand, final ObjectList result) {
- if ( header.size() > 0 ) {
- if ( isNegativeStrand )
- result.addAll(addToSyntheticReads(header, 0, header.size(), ConsensusType.NEGATIVE_CONSENSUS));
- else
- result.addAll(addToSyntheticReads(header, 0, header.size(), ConsensusType.POSITIVE_CONSENSUS));
- }
- }
-
- private void addToHeader(LinkedList header, GATKSAMRecord read) {
- updateHeaderCounts(header, read, false);
- }
-
- private void removeFromHeader(LinkedList header, GATKSAMRecord read) {
- updateHeaderCounts(header, read, true);
- }
-
- /**
- * Updates the sliding window's header counts with the incoming read bases, insertions
- * and deletions.
- *
- * @param header the sliding window header to use
- * @param read the incoming read to be added to the sliding window
- * @param removeRead if we are removing the read from the header or adding
- */
- protected void updateHeaderCounts(final LinkedList header, final GATKSAMRecord read, final boolean removeRead) {
- final int readStart = read.getSoftStart();
- final int headerStart = getStartLocation(header);
- int locationIndex = headerStart < 0 ? 0 : readStart - headerStart;
-
- if ( removeRead && locationIndex < 0 )
- throw new IllegalStateException("Provided read is behind the Sliding Window! Read = " + read + ", readStart = " + readStart + ", cigar = " + read.getCigarString() + ", window = " + headerStart + "-" + getStopLocation(header));
-
- // we only need to create new header elements if we are adding the read, not when we're removing it
- if ( !removeRead )
- locationIndex = createNewHeaderElements(header, read, locationIndex);
-
- actuallyUpdateHeaderForRead(header, read, removeRead, locationIndex);
- }
-
- /*
- * Creates new header elements if needed for the given read.
- *
- * @param header the sliding window header to use
- * @param read the incoming read to be added to the sliding window
- * @param startIndex the start location index into the header for this read
- *
- * @return an updated index into the modified header
- */
- @Requires("header != null && read != null")
- protected int createNewHeaderElements(final LinkedList header, final GATKSAMRecord read, final int startIndex) {
-
- int headerStart = getStartLocation(header);
- int locationIndex = startIndex;
-
- // Do we need to add extra elements before the start of the header? This could happen if the previous read was
- // clipped and this alignment starts before the beginning of the window
- final int readStart = read.getSoftStart();
- if ( startIndex < 0 ) {
- for ( int i = 1; i <= -startIndex; i++ )
- header.addFirst(new HeaderElement(headerStart - i));
-
- // update the start location accordingly
- headerStart = readStart;
- locationIndex = 0;
- }
-
- // Do we need to add extra elements to the end of the header?
- final int headerStop = getStopLocation(header);
- final int readEnd = read.getSoftEnd();
- if ( headerStop < readEnd ) {
- final int elementsToAdd = (headerStop < 0) ? readEnd - readStart + 1 : readEnd - headerStop;
- for ( int i = elementsToAdd - 1; i >= 0; i-- )
- header.addLast(new HeaderElement(readEnd - i));
- }
-
- // Special case for leading insertions before the beginning of the sliding read
- if ( (readStart == headerStart || headerStart < 0) && ReadUtils.readStartsWithInsertion(read.getCigar(), false) != null ) {
- // create a new first element to the window header with no bases added
- header.addFirst(new HeaderElement(readStart - 1));
- // this allows the first element (I) to look at locationIndex - 1 when we update the header and do the right thing
- locationIndex = 1;
- }
-
- return locationIndex;
- }
-
- /*
- * Actually updates the sliding window's header counts with the incoming read bases and quals (including insertion and deletion quals).
- *
- * @param header the sliding window header to use
- * @param read the incoming read to be added to the sliding window
- * @param removeRead if we are removing the read from the header or adding
- * @param startIndex the start location index into the header for this read
- */
- @Requires("header != null && read != null && startIndex >= 0")
- protected void actuallyUpdateHeaderForRead(final LinkedList header, final GATKSAMRecord read, final boolean removeRead, final int startIndex) {
-
- final Iterator headerElementIterator = header.listIterator(startIndex);
- final int mappingQuality = read.getMappingQuality();
- final boolean isNegativeStrand = read.getReadNegativeStrandFlag();
-
- // iterator variables
- int locationIndex = startIndex;
- int readBaseIndex = 0;
- HeaderElement headerElement;
-
- for ( final CigarElement cigarElement : read.getCigar().getCigarElements() ) {
- switch ( cigarElement.getOperator() ) {
- case H:
- break;
- case I:
- readBaseIndex += cigarElement.getLength();
-
- // special case, if we don't have the previous header element anymore, don't worry about it.
- if ( locationIndex == 0 )
- break;
-
- // insertions are added to the base to the left (previous element)
- headerElement = header.get(locationIndex - 1);
-
- if ( removeRead )
- headerElement.removeInsertionToTheRight();
- else
- headerElement.addInsertionToTheRight();
-
- break;
- case D:
- // deletions are added to the baseCounts with the read mapping quality as its quality score
- final int nDeletionBases = cigarElement.getLength();
- final byte MQbyte = mappingQuality > Byte.MAX_VALUE ? Byte.MAX_VALUE : (byte)mappingQuality;
- for ( int i = 0; i < nDeletionBases; i++ ) {
- headerElement = headerElementIterator.next();
- if (removeRead)
- headerElement.removeBase(BaseUtils.Base.D.base, MQbyte, MQbyte, MQbyte, mappingQuality, MIN_BASE_QUAL_TO_COUNT, MIN_MAPPING_QUALITY, false, isNegativeStrand);
- else
- headerElement.addBase(BaseUtils.Base.D.base, MQbyte, MQbyte, MQbyte, mappingQuality, MIN_BASE_QUAL_TO_COUNT, MIN_MAPPING_QUALITY, false, isNegativeStrand);
- }
- locationIndex += nDeletionBases;
- break;
- case S:
- case M:
- case P:
- case EQ:
- case X:
- final int nBasesToAdd = cigarElement.getLength();
- final boolean isSoftClip = cigarElement.getOperator() == CigarOperator.S;
- final byte[] readBases = read.getReadBases();
- final byte[] readQuals = read.getBaseQualities();
- final boolean readHasIndelQuals = read.hasBaseIndelQualities();
- final byte[] insertionQuals = readHasIndelQuals ? read.getBaseInsertionQualities() : null;
- final byte[] deletionQuals = readHasIndelQuals ? read.getBaseDeletionQualities() : null;
-
- for ( int i = 0; i < nBasesToAdd; i++ ) {
- headerElement = headerElementIterator.next();
- final byte insertionQuality = readHasIndelQuals ? insertionQuals[readBaseIndex] : -1;
- final byte deletionQuality = readHasIndelQuals ? deletionQuals[readBaseIndex] : -1;
-
- if ( removeRead )
- headerElement.removeBase(readBases[readBaseIndex], readQuals[readBaseIndex], insertionQuality, deletionQuality, mappingQuality, MIN_BASE_QUAL_TO_COUNT, MIN_MAPPING_QUALITY, isSoftClip, isNegativeStrand);
- else
- headerElement.addBase(readBases[readBaseIndex], readQuals[readBaseIndex], insertionQuality, deletionQuality, mappingQuality, MIN_BASE_QUAL_TO_COUNT, MIN_MAPPING_QUALITY, isSoftClip, isNegativeStrand);
-
- readBaseIndex++;
- }
- locationIndex += nBasesToAdd;
- break;
- default:
- break;
- }
- }
- }
-}
-
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/SyntheticRead.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/SyntheticRead.java
deleted file mode 100644
index 9d16ea06f..000000000
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/SyntheticRead.java
+++ /dev/null
@@ -1,369 +0,0 @@
-/*
-* By downloading the PROGRAM you agree to the following terms of use:
-*
-* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
-*
-* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
-*
-* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
-* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
-* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
-*
-* 1. DEFINITIONS
-* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
-*
-* 2. LICENSE
-* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
-* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
-* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
-* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
-*
-* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
-* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
-* Copyright 2012 Broad Institute, Inc.
-* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
-* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
-*
-* 4. INDEMNIFICATION
-* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
-*
-* 5. NO REPRESENTATIONS OR WARRANTIES
-* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
-* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
-*
-* 6. ASSIGNMENT
-* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
-*
-* 7. MISCELLANEOUS
-* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
-* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
-* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
-* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
-* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
-* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
-* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
-*/
-
-package org.broadinstitute.sting.gatk.walkers.compression.reducereads;
-
-import com.google.java.contract.Requires;
-import it.unimi.dsi.fastutil.objects.ObjectArrayList;
-import net.sf.samtools.Cigar;
-import net.sf.samtools.CigarElement;
-import net.sf.samtools.CigarOperator;
-import net.sf.samtools.SAMFileHeader;
-import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
-import org.broadinstitute.sting.utils.recalibration.EventType;
-import org.broadinstitute.sting.utils.sam.GATKSAMReadGroupRecord;
-import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
-
-import java.util.Iterator;
-
-
-/**
- * Running Consensus is a read that is compressed as a sliding window travels over the reads
- * and keeps track of all the bases that are outside of variant regions.
- *
- * Consensus reads have qual fields that correspond to the number of reads that had the base
- * and passed the minimum quality threshold.
- *
- * The mapping quality of a consensus read is the average RMS of the mapping qualities of all reads
- * that compose the consensus
- *
- * @author Mauricio Carneiro
- * @since 8/26/11
- */
-public class SyntheticRead {
-
- /**
- * The types of strandedness for synthetic reads
- */
- public enum StrandType {
- POSITIVE,
- NEGATIVE,
- STRANDLESS
- }
-
- // Rather than storing a separate list for each attribute in SingleBaseInfo, store one list to reduce memory footprint.
- private static class SingleBaseInfo {
- byte baseIndexOrdinal; // enum BaseIndex.ordinal
- int count;
- byte qual;
- byte insertionQual;
- byte deletionQual;
-
- SingleBaseInfo(byte baseIndexOrdinal, int count, byte qual, byte insertionQual, byte deletionQual) {
- this.baseIndexOrdinal = baseIndexOrdinal;
- this.count = count;
- this.qual = qual;
- this.insertionQual = insertionQual;
- this.deletionQual = deletionQual;
- }
- }
-
- // This class is merely sharing of code for convertVariableGivenBases().
- private abstract class SingleBaseInfoIterator implements Iterator {
- final Iterator it;
-
- SingleBaseInfoIterator() {
- this.it = basesCountsQuals.iterator();
- }
-
- public boolean hasNext() {
- return it.hasNext();
- }
-
- public void remove() {
- throw new UnsupportedOperationException();
- }
- }
-
-
- // Map from ordinal to enum value.
- private static final BaseIndex[] BaseIndexByOrdinal = new BaseIndex[BaseIndex.values().length];
- static
- {
- for (final BaseIndex baseIndex : BaseIndex.values()) {
- BaseIndexByOrdinal[baseIndex.ordinal()] = baseIndex;
- }
- }
-
-
- private final ObjectArrayList basesCountsQuals;
- private double mappingQuality;
-
- // Information to produce a GATKSAMRecord
- private SAMFileHeader header;
- private GATKSAMReadGroupRecord readGroupRecord;
- private String contig;
- private int contigIndex;
- private String readName;
- private int refStart;
- private boolean hasIndelQualities = false;
- private StrandType strandType = StrandType.STRANDLESS;
-
- /**
- * Full initialization of the running consensus if you have all the information and are ready to
- * start adding to the running consensus.
- *
- * @param header GATKSAMRecord file header
- * @param readGroupRecord Read Group for the GATKSAMRecord
- * @param contig the read's contig name
- * @param contigIndex the read's contig index
- * @param readName the read's name
- * @param refStart the alignment start (reference based)
- */
- public SyntheticRead(SAMFileHeader header, GATKSAMReadGroupRecord readGroupRecord, String contig, int contigIndex, String readName, int refStart, boolean hasIndelQualities, StrandType strandType) {
- final int initialCapacity = 10000;
- basesCountsQuals = new ObjectArrayList(initialCapacity);
- mappingQuality = 0.0;
-
- this.header = header;
- this.readGroupRecord = readGroupRecord;
- this.contig = contig;
- this.contigIndex = contigIndex;
- this.readName = readName;
- this.refStart = refStart;
- this.hasIndelQualities = hasIndelQualities;
- this.strandType = strandType;
- }
-
- /**
- * Easy access to keep adding to a running consensus that has already been
- * initialized with the correct read name and refStart
- *
- * @param base the base to add
- * @param count number of reads with this base
- */
- @Requires("count <= Byte.MAX_VALUE")
- public void add(BaseIndex base, int count, byte qual, byte insQual, byte delQual, double mappingQuality) {
- basesCountsQuals.add(new SingleBaseInfo(base.getOrdinalByte(), count, qual, insQual, delQual));
- this.mappingQuality += mappingQuality;
- }
-
- public BaseIndex getBase(final int readCoordinate) {
- return BaseIndexByOrdinal[basesCountsQuals.get(readCoordinate).baseIndexOrdinal];
- }
-
- public int getRefStart() {
- return refStart;
- }
-
- /**
- * Creates a GATKSAMRecord of the synthetic read. Will return null if the read is invalid.
- *
- * Invalid reads are :
- * - exclusively composed of deletions
- *
- * @return a GATKSAMRecord or null
- */
- public GATKSAMRecord close () {
- if (isAllDeletions())
- return null;
-
- GATKSAMRecord read = new GATKSAMRecord(header);
- read.setReferenceName(contig);
- read.setReferenceIndex(contigIndex);
- read.setReadPairedFlag(false);
- read.setReadUnmappedFlag(false);
- if ( strandType != StrandType.STRANDLESS ) {
- read.setAttribute(GATKSAMRecord.REDUCED_READ_STRANDED_TAG, '1'); // must come before next line
- read.setReadNegativeStrandFlag(strandType == StrandType.NEGATIVE);
- }
- read.setCigar(buildCigar()); // the alignment start may change while building the cigar (leading deletions)
- read.setAlignmentStart(refStart);
- read.setReadName(readName);
- read.setBaseQualities(convertBaseQualities(), EventType.BASE_SUBSTITUTION);
- read.setReadBases(convertReadBases());
- read.setMappingQuality((int) Math.ceil(mappingQuality / basesCountsQuals.size()));
- read.setReadGroup(readGroupRecord);
- read.setReducedReadCountsTag(convertBaseCounts());
-
- if (hasIndelQualities) {
- read.setBaseQualities(convertInsertionQualities(), EventType.BASE_INSERTION);
- read.setBaseQualities(convertDeletionQualities(), EventType.BASE_DELETION);
- }
-
- return read;
- }
-
- /**
- * Checks if the synthetic read is composed exclusively of deletions
- *
- * @return true if it is, false if it isn't.
- */
- private boolean isAllDeletions() {
- for (SingleBaseInfo b : basesCountsQuals)
- if (b.baseIndexOrdinal != BaseIndex.D.getOrdinalByte())
- return false;
- return true;
- }
-
- public int size () {
- return basesCountsQuals.size();
- }
-
- private byte [] convertBaseQualities() {
- return convertVariableGivenBases(new SingleBaseInfoIterator() {
- public Byte next() {
- return it.next().qual;
- }
- });
- }
-
- private byte [] convertInsertionQualities() {
- return convertVariableGivenBases(new SingleBaseInfoIterator() {
- public Byte next() {
- return it.next().insertionQual;
- }
- });
- }
-
- private byte [] convertDeletionQualities() {
- return convertVariableGivenBases(new SingleBaseInfoIterator() {
- public Byte next() {
- return it.next().deletionQual;
- }
- });
- }
-
- protected int[] convertBaseCounts() {
- int[] variableArray = new int[getReadLengthWithNoDeletions()];
- int i = 0;
- for (final SingleBaseInfo singleBaseInfo : basesCountsQuals) {
- if (singleBaseInfo.baseIndexOrdinal != BaseIndex.D.getOrdinalByte())
- variableArray[i++] = singleBaseInfo.count;
- }
- return variableArray;
- }
-
- private byte [] convertReadBases() {
- byte [] readArray = new byte[getReadLengthWithNoDeletions()];
- int i = 0;
- for (final SingleBaseInfo singleBaseInfo : basesCountsQuals) {
- final BaseIndex baseIndex = BaseIndexByOrdinal[singleBaseInfo.baseIndexOrdinal];
- if (baseIndex != BaseIndex.D)
- readArray[i++] = baseIndex.getByte();
- }
-
- return readArray;
- }
-
- /**
- * Builds the cigar string for the synthetic read
- *
- * Warning: if the synthetic read has leading deletions, it will shift the refStart (alignment start) of the read.
- *
- * @return the cigar string for the synthetic read
- */
- private Cigar buildCigar() {
- ObjectArrayList cigarElements = new ObjectArrayList();
- CigarOperator cigarOperator = null;
- int length = 0;
- for (final SingleBaseInfo singleBaseInfo : basesCountsQuals) {
- final BaseIndex b = BaseIndexByOrdinal[singleBaseInfo.baseIndexOrdinal];
- CigarOperator op;
- switch (b) {
- case D:
- op = CigarOperator.DELETION;
- break;
- case I:
- throw new ReviewedStingException("Trying to create an insertion in a synthetic read. This operation is currently unsupported.");
- default:
- op = CigarOperator.MATCH_OR_MISMATCH;
- break;
- }
- if (cigarOperator == null) {
- if (op == CigarOperator.D) // read cannot start with a deletion
- refStart++; // if it does, we need to move the reference start forward
- else
- cigarOperator = op;
- }
- else if (cigarOperator != op) { // if this is a new operator, we need to close the previous one
- cigarElements.add(new CigarElement(length, cigarOperator)); // close previous operator
- cigarOperator = op;
- length = 0;
- }
-
- if (cigarOperator != null) // only increment the length of the cigar element if we really added it to the read (no leading deletions)
- length++;
- }
- if (length > 0 && cigarOperator != CigarOperator.D) // read cannot end with a deletion
- cigarElements.add(new CigarElement(length, cigarOperator)); // add the last cigar element
-
- return new Cigar(cigarElements);
- }
-
- /**
- * Shared functionality for all conversion utilities
- *
- * @param variableIterator the list to convert
- * @return a converted variable given the bases and skipping deletions
- */
-
- private byte [] convertVariableGivenBases (Iterator variableIterator) {
- byte [] variableArray = new byte[getReadLengthWithNoDeletions()];
- int i = 0;
- for (final SingleBaseInfo singleBaseInfo : basesCountsQuals) {
- byte count = variableIterator.next();
- if (singleBaseInfo.baseIndexOrdinal != BaseIndex.D.getOrdinalByte())
- variableArray[i++] = count;
- }
- return variableArray;
- }
-
- /**
- * Shared functionality for all conversion utilities
- *
- * @return the length of the read with no deletions
- */
- private int getReadLengthWithNoDeletions() {
- int readLength = basesCountsQuals.size();
- for (final SingleBaseInfo singleBaseInfo : basesCountsQuals)
- if (singleBaseInfo.baseIndexOrdinal == BaseIndex.D.getOrdinalByte())
- readLength--;
- return readLength;
- }
-
-
-}
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidSNPGenotypeLikelihoods.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidSNPGenotypeLikelihoods.java
index f3b26f295..77c51f88b 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidSNPGenotypeLikelihoods.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidSNPGenotypeLikelihoods.java
@@ -187,18 +187,6 @@ public class DiploidSNPGenotypeLikelihoods implements Cloneable {
if ( qual == 0 )
return 0;
- if ( elt.getRead().isReducedRead() ) {
- // reduced read representation
- if ( BaseUtils.isRegularBase( obsBase )) {
- int representativeCount = elt.getRepresentativeCount();
- add(obsBase, qual, (byte)0, (byte)0, representativeCount); // fast calculation of n identical likelihoods
- return representativeCount; // we added nObs bases here
- }
-
- // odd bases or deletions => don't use them
- return 0;
- }
-
return add(obsBase, qual, (byte)0, (byte)0, 1);
}
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/ErrorModel.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/ErrorModel.java
index 7ce736b0c..a57502bc0 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/ErrorModel.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/ErrorModel.java
@@ -140,11 +140,10 @@ public class ErrorModel {
Allele refAllele = refSampleVC.getReference();
if ( refSampleVC.isIndel()) {
- final int readCounts[] = new int[refSamplePileup.getNumberOfElements()];
//perReadLikelihoods = new double[readCounts.length][refSampleVC.getAlleles().size()];
final int eventLength = IndelGenotypeLikelihoodsCalculationModel.getEventLength(refSampleVC.getAlleles());
if (!haplotypeMap.isEmpty())
- perReadLikelihoods = pairModel.computeGeneralReadHaplotypeLikelihoods(refSamplePileup,haplotypeMap,refContext, eventLength, perReadAlleleLikelihoodMap, readCounts);
+ perReadLikelihoods = pairModel.computeGeneralReadHaplotypeLikelihoods(refSamplePileup,haplotypeMap,refContext, eventLength, perReadAlleleLikelihoodMap);
}
int idx = 0;
for (PileupElement refPileupElement : refSamplePileup) {
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/GeneralPloidyIndelGenotypeLikelihoods.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/GeneralPloidyIndelGenotypeLikelihoods.java
index 2f2a93fa4..530ba3ef8 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/GeneralPloidyIndelGenotypeLikelihoods.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/GeneralPloidyIndelGenotypeLikelihoods.java
@@ -193,8 +193,7 @@ public class GeneralPloidyIndelGenotypeLikelihoods extends GeneralPloidyGenotype
if (!hasReferenceSampleData) {
- final int readCounts[] = new int[pileup.getNumberOfElements()];
- readHaplotypeLikelihoods = pairModel.computeGeneralReadHaplotypeLikelihoods(pileup, haplotypeMap, refContext, eventLength, perReadAlleleLikelihoodMap, readCounts);
+ readHaplotypeLikelihoods = pairModel.computeGeneralReadHaplotypeLikelihoods(pileup, haplotypeMap, refContext, eventLength, perReadAlleleLikelihoodMap);
n = readHaplotypeLikelihoods.length;
} else {
Allele refAllele = null;
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/GenotypeLikelihoodsCalculationModel.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/GenotypeLikelihoodsCalculationModel.java
index f48ae81cf..95d3fb78b 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/GenotypeLikelihoodsCalculationModel.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/GenotypeLikelihoodsCalculationModel.java
@@ -132,7 +132,7 @@ public abstract class GenotypeLikelihoodsCalculationModel implements Cloneable {
int count = 0;
for ( PileupElement p : pileup ) {
if ( BaseUtils.isRegularBase( p.getBase() ) )
- count += p.getRepresentativeCount();
+ count++;
}
return count;
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java
index 4a3231b3e..ae2ea427b 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java
@@ -253,7 +253,7 @@ public class IndelGenotypeLikelihoodsCalculationModel extends GenotypeLikelihood
int count = 0;
for (PileupElement p : pileup) {
if (p.isDeletion() || BaseUtils.isRegularBase(p.getBase()))
- count += p.getRepresentativeCount();
+ count++;
}
return count;
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java
index aa334f680..c5070a76f 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java
@@ -608,7 +608,7 @@ public class UnifiedGenotyperEngine {
int numDeletions = 0;
for ( final PileupElement p : rawContext.getBasePileup() ) {
if ( p.isDeletion() )
- numDeletions += p.getRepresentativeCount();
+ numDeletions++;
}
if ( ((double) numDeletions) / ((double) rawContext.getBasePileup().depthOfCoverage()) > UAC.MAX_DELETION_FRACTION ) {
return null;
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java
index 276103277..74164415e 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java
@@ -278,6 +278,13 @@ public class HaplotypeCaller extends ActiveRegionWalker, In
@Argument(fullName="numPruningSamples", shortName="numPruningSamples", doc="The number of samples that must pass the minPuning factor in order for the path to be kept", required = false)
protected int numPruningSamples = 1;
+ /**
+ * This mode is currently experimental and should only be used in the RNA-seq calling pipeline.
+ */
+ @Advanced
+ @Argument(fullName="recoverDanglingHeads", shortName="recoverDanglingHeads", doc="Should we enable dangling head recovery in the read threading assembler?", required = false)
+ protected boolean recoverDanglingHeads = false;
+
@Hidden
@Argument(fullName="dontRecoverDanglingTails", shortName="dontRecoverDanglingTails", doc="Should we disable dangling tail recovery in the read threading assembler?", required = false)
protected boolean dontRecoverDanglingTails = false;
@@ -452,7 +459,7 @@ public class HaplotypeCaller extends ActiveRegionWalker, In
@Argument(fullName="debugGraphTransformations", shortName="debugGraphTransformations", doc="If specified, we will write DOT formatted graph files out of the assembler for only this graph size", required = false)
protected boolean debugGraphTransformations = false;
- @Hidden
+ @Advanced
@Argument(fullName="dontUseSoftClippedBases", shortName="dontUseSoftClippedBases", doc="If specified, we will not analyze soft clipped bases in the reads", required = false)
protected boolean dontUseSoftClippedBases = false;
@@ -634,6 +641,7 @@ public class HaplotypeCaller extends ActiveRegionWalker, In
assemblyEngine.setDebugGraphTransformations(debugGraphTransformations);
assemblyEngine.setAllowCyclesInKmerGraphToGeneratePaths(allowCyclesInKmerGraphToGeneratePaths);
assemblyEngine.setRecoverDanglingTails(!dontRecoverDanglingTails);
+ assemblyEngine.setRecoverDanglingHeads(recoverDanglingHeads);
assemblyEngine.setMinBaseQualityToUseInAssembly(MIN_BASE_QUALTY_SCORE);
MIN_TAIL_QUALITY = (byte)(MIN_BASE_QUALTY_SCORE - 1);
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LocalAssemblyEngine.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LocalAssemblyEngine.java
index d0e28d878..d49827405 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LocalAssemblyEngine.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LocalAssemblyEngine.java
@@ -86,6 +86,7 @@ public abstract class LocalAssemblyEngine {
protected boolean allowCyclesInKmerGraphToGeneratePaths = false;
protected boolean debugGraphTransformations = false;
protected boolean recoverDanglingTails = true;
+ protected boolean recoverDanglingHeads = true;
protected byte minBaseQualityToUseInAssembly = DEFAULT_MIN_BASE_QUALITY_TO_USE;
protected int pruneFactor = 2;
@@ -456,4 +457,12 @@ public abstract class LocalAssemblyEngine {
public void setRecoverDanglingTails(boolean recoverDanglingTails) {
this.recoverDanglingTails = recoverDanglingTails;
}
+
+ public boolean isRecoverDanglingHeads() {
+ return recoverDanglingHeads;
+ }
+
+ public void setRecoverDanglingHeads(boolean recoverDanglingHeads) {
+ this.recoverDanglingHeads = recoverDanglingHeads;
+ }
}
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/PairHMMLikelihoodCalculationEngine.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/PairHMMLikelihoodCalculationEngine.java
index c534764a5..9de532dcf 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/PairHMMLikelihoodCalculationEngine.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/PairHMMLikelihoodCalculationEngine.java
@@ -59,7 +59,6 @@ import org.broadinstitute.sting.utils.pairhmm.*;
import org.broadinstitute.sting.utils.recalibration.covariates.RepeatCovariate;
import org.broadinstitute.sting.utils.recalibration.covariates.RepeatLengthCovariate;
import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
-import org.broadinstitute.sting.utils.sam.ReadUtils;
import org.broadinstitute.variant.variantcontext.*;
import java.io.File;
@@ -433,8 +432,7 @@ public class PairHMMLikelihoodCalculationEngine implements LikelihoodCalculation
for( final Map.Entry> entry : stratifiedReadMap.get(sample).getLikelihoodReadMap().entrySet() ) {
// Compute log10(10^x1/2 + 10^x2/2) = log10(10^x1+10^x2)-log10(2)
// First term is approximated by Jacobian log with table lookup.
- haplotypeLikelihood += ReadUtils.getMeanRepresentativeReadCount( entry.getKey() ) *
- ( MathUtils.approximateLog10SumLog10(entry.getValue().get(iii_allele), entry.getValue().get(jjj_allele)) + MathUtils.LOG_ONE_HALF );
+ haplotypeLikelihood += ( MathUtils.approximateLog10SumLog10(entry.getValue().get(iii_allele), entry.getValue().get(jjj_allele)) + MathUtils.LOG_ONE_HALF );
}
}
haplotypeLikelihoodMatrix[iii][jjj] = haplotypeLikelihood;
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/ReadErrorCorrector.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/ReadErrorCorrector.java
index e1471ab33..a48ac9ee0 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/ReadErrorCorrector.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/ReadErrorCorrector.java
@@ -233,10 +233,6 @@ public class ReadErrorCorrector {
*/
@Requires("inputRead != null")
private GATKSAMRecord correctRead(final GATKSAMRecord inputRead) {
- // no support for reduced reads (which shouldn't need to be error-corrected anyway!)
- if (inputRead.isReducedRead())
- return inputRead;
-
// do actual correction
boolean corrected = false;
final byte[] correctedBases = inputRead.getReadBases();
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/ReferenceConfidenceModel.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/ReferenceConfidenceModel.java
index 7f7e65817..5ef310498 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/ReferenceConfidenceModel.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/ReferenceConfidenceModel.java
@@ -297,13 +297,13 @@ public class ReferenceConfidenceModel {
if( hqSoftClips != null && p.isNextToSoftClip() ) {
hqSoftClips.add(AlignmentUtils.calcNumHighQualitySoftClips(p.getRead(), (byte) 28));
}
- result.AD_Ref_Any[1] += p.getRepresentativeCount();
+ result.AD_Ref_Any[1]++;
} else {
- result.AD_Ref_Any[0] += p.getRepresentativeCount();
+ result.AD_Ref_Any[0]++;
}
- result.genotypeLikelihoods[AA] += p.getRepresentativeCount() * QualityUtils.qualToProbLog10(qual);
- result.genotypeLikelihoods[AB] += p.getRepresentativeCount() * MathUtils.approximateLog10SumLog10( QualityUtils.qualToProbLog10(qual) + MathUtils.LOG_ONE_HALF, QualityUtils.qualToErrorProbLog10(qual) + MathUtils.LOG_ONE_THIRD + MathUtils.LOG_ONE_HALF );
- result.genotypeLikelihoods[BB] += p.getRepresentativeCount() * QualityUtils.qualToErrorProbLog10(qual) + MathUtils.LOG_ONE_THIRD;
+ result.genotypeLikelihoods[AA] += QualityUtils.qualToProbLog10(qual);
+ result.genotypeLikelihoods[AB] += MathUtils.approximateLog10SumLog10( QualityUtils.qualToProbLog10(qual) + MathUtils.LOG_ONE_HALF, QualityUtils.qualToErrorProbLog10(qual) + MathUtils.LOG_ONE_THIRD + MathUtils.LOG_ONE_HALF );
+ result.genotypeLikelihoods[BB] += QualityUtils.qualToErrorProbLog10(qual) + MathUtils.LOG_ONE_THIRD;
}
}
@@ -484,7 +484,7 @@ public class ReferenceConfidenceModel {
// todo -- this code really should handle CIGARs directly instead of relying on the above tests
if ( isReadInformativeAboutIndelsOfSize(read.getReadBases(), read.getBaseQualities(), offset, ref, pileupOffsetIntoRef, maxIndelSize) ) {
- nInformative += p.getRepresentativeCount();
+ nInformative++;
if( nInformative > MAX_N_INDEL_INFORMATIVE_READS ) {
return MAX_N_INDEL_INFORMATIVE_READS;
}
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/BaseGraph.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/BaseGraph.java
index edd8dbb16..c9d51b81b 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/BaseGraph.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/BaseGraph.java
@@ -646,24 +646,6 @@ public class BaseGraph extends Default
'}';
}
- /**
- * The base sequence for the given path.
- * Note, this assumes that the path does not start with a source node.
- *
- * @param path the list of vertexes that make up the path
- * @return non-null sequence of bases corresponding to the given path
- */
- @Ensures({"result != null"})
- public byte[] getBasesForPath(final List path) {
- if ( path == null ) throw new IllegalArgumentException("Path cannot be null");
-
- final StringBuffer sb = new StringBuffer();
- for ( final DeBruijnVertex v : path )
- sb.append((char)v.getSuffix());
-
- return sb.toString().getBytes();
- }
-
/**
* Get the set of vertices within distance edges of source, regardless of edge direction
*
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/DeBruijnVertex.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/DeBruijnVertex.java
index ec2ccff20..cf95f6a5a 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/DeBruijnVertex.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/graphs/DeBruijnVertex.java
@@ -65,14 +65,6 @@ public class DeBruijnVertex extends BaseVertex {
super(sequence);
}
- /**
- * For testing purposes only
- * @param sequence
- */
- protected DeBruijnVertex( final String sequence ) {
- this(sequence.getBytes());
- }
-
/**
* Get the kmer size for this DeBruijnVertex
* @return integer >= 1
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/DanglingChainMergingGraph.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/DanglingChainMergingGraph.java
new file mode 100644
index 000000000..c696c50ae
--- /dev/null
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/DanglingChainMergingGraph.java
@@ -0,0 +1,526 @@
+/*
+* By downloading the PROGRAM you agree to the following terms of use:
+*
+* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
+*
+* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
+*
+* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
+* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
+* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
+*
+* 1. DEFINITIONS
+* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
+*
+* 2. LICENSE
+* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
+* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
+* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
+* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
+*
+* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
+* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
+* Copyright 2012 Broad Institute, Inc.
+* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
+* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
+*
+* 4. INDEMNIFICATION
+* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
+*
+* 5. NO REPRESENTATIONS OR WARRANTIES
+* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
+* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
+*
+* 6. ASSIGNMENT
+* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
+*
+* 7. MISCELLANEOUS
+* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
+* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
+* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
+* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
+* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
+* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
+* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
+*/
+
+package org.broadinstitute.sting.gatk.walkers.haplotypecaller.readthreading;
+
+import com.google.java.contract.Ensures;
+import net.sf.samtools.Cigar;
+import net.sf.samtools.CigarElement;
+import net.sf.samtools.CigarOperator;
+import org.broadinstitute.sting.gatk.walkers.haplotypecaller.graphs.*;
+import org.broadinstitute.sting.utils.sam.AlignmentUtils;
+import org.broadinstitute.sting.utils.smithwaterman.*;
+import org.jgrapht.EdgeFactory;
+
+import java.util.*;
+
+public abstract class DanglingChainMergingGraph extends BaseGraph {
+
+ private static final int MAX_CIGAR_COMPLEXITY = 3;
+ private static final int MIN_DANGLING_TAIL_LENGTH = 5; // SNP + 3 stabilizing nodes + the LCA
+ private static final int MAXIMUM_MISMATCHES_IN_DANGLING_HEAD_MERGE = 1;
+
+ protected boolean alreadyBuilt;
+
+ /**
+ * Create a new ReadThreadingAssembler using kmerSize for matching
+ * @param kmerSize must be >= 1
+ */
+ protected DanglingChainMergingGraph(final int kmerSize, final EdgeFactory edgeFactory) {
+ super(kmerSize, edgeFactory);
+ }
+
+ /**
+ * Edge factory that encapsulates the numPruningSamples assembly parameter
+ */
+ protected static class MyEdgeFactory implements EdgeFactory {
+ final int numPruningSamples;
+
+ public MyEdgeFactory(int numPruningSamples) {
+ this.numPruningSamples = numPruningSamples;
+ }
+
+ @Override
+ public MultiSampleEdge createEdge(final MultiDeBruijnVertex sourceVertex, final MultiDeBruijnVertex targetVertex) {
+ return new MultiSampleEdge(false, 1, numPruningSamples);
+ }
+
+ public MultiSampleEdge createEdge(final boolean isRef, final int multiplicity) {
+ return new MultiSampleEdge(isRef, multiplicity, numPruningSamples);
+ }
+
+ }
+
+ /**
+ * Class to keep track of the important dangling chain merging data
+ */
+ protected static final class DanglingChainMergeHelper {
+ final List danglingPath, referencePath;
+ final byte[] danglingPathString, referencePathString;
+ final Cigar cigar;
+
+ public DanglingChainMergeHelper(final List danglingPath,
+ final List referencePath,
+ final byte[] danglingPathString,
+ final byte[] referencePathString,
+ final Cigar cigar) {
+ this.danglingPath = danglingPath;
+ this.referencePath = referencePath;
+ this.danglingPathString = danglingPathString;
+ this.referencePathString = referencePathString;
+ this.cigar = cigar;
+ }
+ }
+
+ /**
+ * Try to recover dangling tails
+ *
+ * @param pruneFactor the prune factor to use in ignoring chain pieces
+ */
+ public void recoverDanglingTails(final int pruneFactor) {
+ if ( ! alreadyBuilt ) throw new IllegalStateException("recoverDanglingTails requires the graph be already built");
+
+ int attempted = 0;
+ int nRecovered = 0;
+ for ( final MultiDeBruijnVertex v : vertexSet() ) {
+ if ( outDegreeOf(v) == 0 && ! isRefSink(v) ) {
+ attempted++;
+ nRecovered += recoverDanglingTail(v, pruneFactor);
+ }
+ }
+
+ logger.debug("Recovered " + nRecovered + " of " + attempted + " dangling tails");
+ }
+
+ /**
+ * Try to recover dangling heads
+ *
+ * @param pruneFactor the prune factor to use in ignoring chain pieces
+ */
+ public void recoverDanglingHeads(final int pruneFactor) {
+ if ( ! alreadyBuilt ) throw new IllegalStateException("recoverDanglingHeads requires the graph be already built");
+
+ // we need to build a list of dangling heads because that process can modify the graph (and otherwise generate
+ // a ConcurrentModificationException if we do it while iterating over the vertexes)
+ final List danglingHeads = new ArrayList<>();
+
+ int attempted = 0;
+ int nRecovered = 0;
+ for ( final MultiDeBruijnVertex v : vertexSet() ) {
+ if ( inDegreeOf(v) == 0 && ! isRefSource(v) )
+ danglingHeads.add(v);
+ }
+
+ // now we can try to recover the dangling heads
+ for ( final MultiDeBruijnVertex v : danglingHeads ) {
+ attempted++;
+ nRecovered += recoverDanglingHead(v, pruneFactor);
+ }
+
+ logger.debug("Recovered " + nRecovered + " of " + attempted + " dangling heads");
+ }
+
+ /**
+ * Attempt to attach vertex with out-degree == 0 to the graph
+ *
+ * @param vertex the vertex to recover
+ * @param pruneFactor the prune factor to use in ignoring chain pieces
+ * @return 1 if we successfully recovered the vertex and 0 otherwise
+ */
+ protected int recoverDanglingTail(final MultiDeBruijnVertex vertex, final int pruneFactor) {
+ if ( outDegreeOf(vertex) != 0 ) throw new IllegalStateException("Attempting to recover a dangling tail for " + vertex + " but it has out-degree > 0");
+
+ // generate the CIGAR string from Smith-Waterman between the dangling tail and reference paths
+ final DanglingChainMergeHelper danglingTailMergeResult = generateCigarAgainstDownwardsReferencePath(vertex, pruneFactor);
+
+ // if the CIGAR is too complex (or couldn't be computed) then we do not allow the merge into the reference path
+ if ( danglingTailMergeResult == null || ! cigarIsOkayToMerge(danglingTailMergeResult.cigar, false, true) )
+ return 0;
+
+ // merge
+ return mergeDanglingTail(danglingTailMergeResult);
+ }
+
+ /**
+ * Attempt to attach vertex with in-degree == 0, or a vertex on its path, to the graph
+ *
+ * @param vertex the vertex to recover
+ * @param pruneFactor the prune factor to use in ignoring chain pieces
+ * @return 1 if we successfully recovered a vertex and 0 otherwise
+ */
+ protected int recoverDanglingHead(final MultiDeBruijnVertex vertex, final int pruneFactor) {
+ if ( inDegreeOf(vertex) != 0 ) throw new IllegalStateException("Attempting to recover a dangling head for " + vertex + " but it has in-degree > 0");
+
+ // generate the CIGAR string from Smith-Waterman between the dangling tail and reference paths
+ final DanglingChainMergeHelper danglingHeadMergeResult = generateCigarAgainstUpwardsReferencePath(vertex, pruneFactor);
+
+ // if the CIGAR is too complex (or couldn't be computed) then we do not allow the merge into the reference path
+ if ( danglingHeadMergeResult == null || ! cigarIsOkayToMerge(danglingHeadMergeResult.cigar, true, false) )
+ return 0;
+
+ // merge
+ return mergeDanglingHead(danglingHeadMergeResult);
+ }
+
+ /**
+ * Determine whether the provided cigar is okay to merge into the reference path
+ *
+ * @param cigar the cigar to analyze
+ * @param requireFirstElementM if true, require that the first cigar element be an M operator in order for it to be okay
+ * @param requireLastElementM if true, require that the last cigar element be an M operator in order for it to be okay
+ * @return true if it's okay to merge, false otherwise
+ */
+ protected boolean cigarIsOkayToMerge(final Cigar cigar, final boolean requireFirstElementM, final boolean requireLastElementM) {
+
+ final List elements = cigar.getCigarElements();
+ final int numElements = elements.size();
+
+ // don't allow more than a couple of different ops
+ if ( numElements == 0 || numElements > MAX_CIGAR_COMPLEXITY )
+ return false;
+
+ // the first element must be an M
+ if ( requireFirstElementM && elements.get(0).getOperator() != CigarOperator.M )
+ return false;
+
+ // the last element must be an M
+ if ( requireLastElementM && elements.get(numElements - 1).getOperator() != CigarOperator.M )
+ return false;
+
+ // TODO -- do we want to check whether the Ms mismatch too much also?
+
+ return true;
+ }
+
+ /**
+ * Actually merge the dangling tail if possible
+ *
+ * @param danglingTailMergeResult the result from generating a Cigar for the dangling tail against the reference
+ * @return 1 if merge was successful, 0 otherwise
+ */
+ protected int mergeDanglingTail(final DanglingChainMergeHelper danglingTailMergeResult) {
+
+ final List elements = danglingTailMergeResult.cigar.getCigarElements();
+ final CigarElement lastElement = elements.get(elements.size() - 1);
+ if ( lastElement.getOperator() != CigarOperator.M )
+ throw new IllegalArgumentException("The last Cigar element must be an M");
+
+ final int lastRefIndex = danglingTailMergeResult.cigar.getReferenceLength() - 1;
+ final int matchingSuffix = Math.min(GraphUtils.longestSuffixMatch(danglingTailMergeResult.referencePathString, danglingTailMergeResult.danglingPathString, lastRefIndex), lastElement.getLength());
+ if ( matchingSuffix == 0 )
+ return 0;
+
+ final int altIndexToMerge = Math.max(danglingTailMergeResult.cigar.getReadLength() - matchingSuffix - 1, 0);
+
+ // there is an important edge condition that we need to handle here: Smith-Waterman correctly calculates that there is a
+ // deletion, that deletion is left-aligned such that the LCA node is part of that deletion, and the rest of the dangling
+ // tail is a perfect match to the suffix of the reference path. In this case we need to push the reference index to merge
+ // down one position so that we don't incorrectly cut a base off of the deletion.
+ final boolean firstElementIsDeletion = elements.get(0).getOperator() == CigarOperator.D;
+ final boolean mustHandleLeadingDeletionCase = firstElementIsDeletion && (elements.get(0).getLength() + matchingSuffix == lastRefIndex + 1);
+ final int refIndexToMerge = lastRefIndex - matchingSuffix + 1 + (mustHandleLeadingDeletionCase ? 1 : 0);
+
+ // another edge condition occurs here: if Smith-Waterman places the whole tail into an insertion then it will try to
+ // merge back to the LCA, which results in a cycle in the graph. So we do not want to merge in such a case.
+ if ( refIndexToMerge == 0 )
+ return 0;
+
+ // it's safe to merge now
+ addEdge(danglingTailMergeResult.danglingPath.get(altIndexToMerge), danglingTailMergeResult.referencePath.get(refIndexToMerge), ((MyEdgeFactory)getEdgeFactory()).createEdge(false, 1));
+
+ return 1;
+ }
+
+ /**
+ * Actually merge the dangling head if possible
+ *
+ * @param danglingHeadMergeResult the result from generating a Cigar for the dangling head against the reference
+ * @return 1 if merge was successful, 0 otherwise
+ */
+ protected int mergeDanglingHead(final DanglingChainMergeHelper danglingHeadMergeResult) {
+
+ final List elements = danglingHeadMergeResult.cigar.getCigarElements();
+ final CigarElement firstElement = elements.get(0);
+ if ( firstElement.getOperator() != CigarOperator.M )
+ throw new IllegalArgumentException("The first Cigar element must be an M");
+
+ final int indexesToMerge = bestPrefixMatch(danglingHeadMergeResult.referencePathString, danglingHeadMergeResult.danglingPathString, firstElement.getLength());
+ if ( indexesToMerge <= 0 )
+ return 0;
+
+ // we can't push back the reference path
+ if ( indexesToMerge >= danglingHeadMergeResult.referencePath.size() - 1 )
+ return 0;
+
+ // but we can manipulate the dangling path if we need to
+ if ( indexesToMerge >= danglingHeadMergeResult.danglingPath.size() &&
+ ! extendDanglingPathAgainstReference(danglingHeadMergeResult, indexesToMerge - danglingHeadMergeResult.danglingPath.size() + 2) )
+ return 0;
+
+ addEdge(danglingHeadMergeResult.referencePath.get(indexesToMerge+1), danglingHeadMergeResult.danglingPath.get(indexesToMerge), ((MyEdgeFactory)getEdgeFactory()).createEdge(false, 1));
+
+ return 1;
+ }
+
+ /**
+ * Generates the CIGAR string from the Smith-Waterman alignment of the dangling path (where the
+ * provided vertex is the sink) and the reference path.
+ *
+ * @param vertex the sink of the dangling chain
+ * @param pruneFactor the prune factor to use in ignoring chain pieces
+ * @return a SmithWaterman object which can be null if no proper alignment could be generated
+ */
+ protected DanglingChainMergeHelper generateCigarAgainstDownwardsReferencePath(final MultiDeBruijnVertex vertex, final int pruneFactor) {
+
+ // find the lowest common ancestor path between vertex and the reference sink if available
+ final List altPath = findPathUpwardsToLowestCommonAncestorOfReference(vertex, pruneFactor);
+ if ( altPath == null || isRefSource(altPath.get(0)) || altPath.size() < MIN_DANGLING_TAIL_LENGTH )
+ return null;
+
+ // now get the reference path from the LCA
+ final List refPath = getReferencePath(altPath.get(0), TraversalDirection.downwards);
+
+ // create the Smith-Waterman strings to use
+ final byte[] refBases = getBasesForPath(refPath, false);
+ final byte[] altBases = getBasesForPath(altPath, false);
+
+ // run Smith-Waterman to determine the best alignment (and remove trailing deletions since they aren't interesting)
+ final SmithWaterman alignment = new SWPairwiseAlignment(refBases, altBases, SWParameterSet.STANDARD_NGS, SWPairwiseAlignment.OVERHANG_STRATEGY.LEADING_INDEL);
+ return new DanglingChainMergeHelper(altPath, refPath, altBases, refBases, AlignmentUtils.removeTrailingDeletions(alignment.getCigar()));
+ }
+
+ /**
+ * Generates the CIGAR string from the Smith-Waterman alignment of the dangling path (where the
+ * provided vertex is the source) and the reference path.
+ *
+ * @param vertex the source of the dangling head
+ * @param pruneFactor the prune factor to use in ignoring chain pieces
+ * @return a SmithWaterman object which can be null if no proper alignment could be generated
+ */
+ protected DanglingChainMergeHelper generateCigarAgainstUpwardsReferencePath(final MultiDeBruijnVertex vertex, final int pruneFactor) {
+
+ // find the highest common descendant path between vertex and the reference source if available
+ final List altPath = findPathDownwardsToHighestCommonDescendantOfReference(vertex, pruneFactor);
+ if ( altPath == null || isRefSink(altPath.get(0)) )
+ return null;
+
+ // now get the reference path from the LCA
+ final List refPath = getReferencePath(altPath.get(0), TraversalDirection.upwards);
+
+ // create the Smith-Waterman strings to use
+ final byte[] refBases = getBasesForPath(refPath, true);
+ final byte[] altBases = getBasesForPath(altPath, true);
+
+ // run Smith-Waterman to determine the best alignment (and remove trailing deletions since they aren't interesting)
+ final SmithWaterman alignment = new SWPairwiseAlignment(refBases, altBases, SWParameterSet.STANDARD_NGS, SWPairwiseAlignment.OVERHANG_STRATEGY.LEADING_INDEL);
+ return new DanglingChainMergeHelper(altPath, refPath, altBases, refBases, AlignmentUtils.removeTrailingDeletions(alignment.getCigar()));
+ }
+
+ /**
+ * Finds the path upwards in the graph from this vertex to the reference sequence, including the lowest common ancestor vertex.
+ * Note that nodes are excluded if their pruning weight is less than the pruning factor.
+ *
+ * @param vertex the original vertex
+ * @param pruneFactor the prune factor to use in ignoring chain pieces
+ * @return the path if it can be determined or null if this vertex either doesn't merge onto the reference path or
+ * has an ancestor with multiple incoming edges before hitting the reference path
+ */
+ protected List findPathUpwardsToLowestCommonAncestorOfReference(final MultiDeBruijnVertex vertex, final int pruneFactor) {
+ final LinkedList path = new LinkedList<>();
+
+ MultiDeBruijnVertex v = vertex;
+ while ( ! isReferenceNode(v) && inDegreeOf(v) == 1 ) {
+ final MultiSampleEdge edge = incomingEdgeOf(v);
+ // if it has too low a weight, don't use it (or previous vertexes) for the path
+ if ( edge.getPruningMultiplicity() < pruneFactor )
+ path.clear();
+ // otherwise it is safe to use
+ else
+ path.addFirst(v);
+ v = getEdgeSource(edge);
+ }
+ path.addFirst(v);
+
+ return isReferenceNode(v) ? path : null;
+ }
+
+ /**
+ * Finds the path downwards in the graph from this vertex to the reference sequence, including the highest common descendant vertex.
+ * However note that the path is reversed so that this vertex ends up at the end of the path.
+ * Also note that nodes are excluded if their pruning weight is less than the pruning factor.
+ *
+ * @param vertex the original vertex
+ * @param pruneFactor the prune factor to use in ignoring chain pieces
+ * @return the path if it can be determined or null if this vertex either doesn't merge onto the reference path or
+ * has a descendant with multiple outgoing edges before hitting the reference path
+ */
+ protected List findPathDownwardsToHighestCommonDescendantOfReference(final MultiDeBruijnVertex vertex, final int pruneFactor) {
+ final LinkedList path = new LinkedList<>();
+
+ MultiDeBruijnVertex v = vertex;
+ while ( ! isReferenceNode(v) && outDegreeOf(v) == 1 ) {
+ final MultiSampleEdge edge = outgoingEdgeOf(v);
+ // if it has too low a weight, don't use it (or previous vertexes) for the path
+ if ( edge.getPruningMultiplicity() < pruneFactor )
+ path.clear();
+ // otherwise it is safe to use
+ else
+ path.addFirst(v);
+ v = getEdgeTarget(edge);
+ }
+ path.addFirst(v);
+
+ return isReferenceNode(v) ? path : null;
+ }
+
+ private enum TraversalDirection {
+ downwards,
+ upwards
+ }
+
+ /**
+ * Finds the path in the graph from this vertex to the reference sink, including this vertex
+ *
+ * @param start the reference vertex to start from
+ * @param direction describes which direction to move in the graph (i.e. down to the reference sink or up to the source)
+ * @return the path (non-null, non-empty)
+ */
+ protected List getReferencePath(final MultiDeBruijnVertex start, final TraversalDirection direction) {
+ if ( ! isReferenceNode(start) ) throw new IllegalArgumentException("Cannot construct the reference path from a vertex that is not on that path");
+
+ final List path = new ArrayList<>();
+
+ MultiDeBruijnVertex v = start;
+ while ( v != null ) {
+ path.add(v);
+ v = (direction == TraversalDirection.downwards ? getNextReferenceVertex(v) : getPrevReferenceVertex(v));
+ }
+
+ return path;
+ }
+
+ /**
+ * The base sequence for the given path.
+ *
+ * @param path the list of vertexes that make up the path
+ * @param reverseIfSource if true and if we encounter a source node, then reverse the character sequence for that node
+ * @return non-null sequence of bases corresponding to the given path
+ */
+ @Ensures({"result != null"})
+ public byte[] getBasesForPath(final List path, final boolean reverseIfSource) {
+ if ( path == null ) throw new IllegalArgumentException("Path cannot be null");
+
+ final StringBuilder sb = new StringBuilder();
+ for ( final MultiDeBruijnVertex v : path ) {
+ if ( isSource(v) ) {
+ final String seq = v.getSequenceString();
+ sb.append(reverseIfSource ? new StringBuilder(seq).reverse().toString() : seq);
+ } else {
+ sb.append((char)v.getSuffix());
+ }
+ }
+
+ return sb.toString().getBytes();
+ }
+
+ /**
+ * Finds the index of the best extent of the prefix match between the provided paths, for dangling head merging.
+ * Assumes that path1.length >= maxIndex and path2.length >= maxIndex.
+ *
+ * @param path1 the first path
+ * @param path2 the second path
+ * @param maxIndex the maximum index to traverse (not inclusive)
+ * @return the index of the ideal prefix match or -1 if it cannot find one, must be less than maxIndex
+ */
+ protected static int bestPrefixMatch(final byte[] path1, final byte[] path2, final int maxIndex) {
+ int mismatches = 0;
+ int index = 0;
+ int lastGoodIndex = -1;
+ while ( index < maxIndex ) {
+ if ( path1[index] != path2[index] ) {
+ if ( ++mismatches > MAXIMUM_MISMATCHES_IN_DANGLING_HEAD_MERGE )
+ return lastGoodIndex;
+ lastGoodIndex = index;
+ }
+ index++;
+ }
+ // if we got here then we hit the max index
+ return lastGoodIndex;
+ }
+
+ protected boolean extendDanglingPathAgainstReference(final DanglingChainMergeHelper danglingHeadMergeResult, final int numNodesToExtend) {
+
+ final int indexOfLastDanglingNode = danglingHeadMergeResult.danglingPath.size() - 1;
+ final int indexOfRefNodeToUse = indexOfLastDanglingNode + numNodesToExtend;
+ if ( indexOfRefNodeToUse >= danglingHeadMergeResult.referencePath.size() )
+ return false;
+
+ final MultiDeBruijnVertex danglingSource = danglingHeadMergeResult.danglingPath.remove(indexOfLastDanglingNode);
+ final StringBuilder sb = new StringBuilder();
+ final byte[] refSourceSequence = danglingHeadMergeResult.referencePath.get(indexOfRefNodeToUse).getSequence();
+ for ( int i = 0; i < numNodesToExtend; i++ )
+ sb.append((char)refSourceSequence[i]);
+ sb.append(danglingSource.getSequenceString());
+ final byte[] sequenceToExtend = sb.toString().getBytes();
+
+ // clean up the source and edge
+ final MultiSampleEdge sourceEdge = outgoingEdgeOf(danglingSource);
+ MultiDeBruijnVertex prevV = getEdgeTarget(sourceEdge);
+ removeEdge(danglingSource, prevV);
+
+ // extend the path
+ for ( int i = numNodesToExtend; i > 0; i-- ) {
+ final MultiDeBruijnVertex newV = new MultiDeBruijnVertex(Arrays.copyOfRange(sequenceToExtend, i, i+kmerSize));
+ addVertex(newV);
+ final MultiSampleEdge newE = addEdge(newV, prevV);
+ newE.setMultiplicity(sourceEdge.getMultiplicity());
+ danglingHeadMergeResult.danglingPath.add(newV);
+ prevV = newV;
+ }
+
+ return true;
+ }
+}
\ No newline at end of file
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/HaplotypeGraph.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/HaplotypeGraph.java
index 55ff2f978..150cdc826 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/HaplotypeGraph.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/HaplotypeGraph.java
@@ -239,14 +239,14 @@ public class HaplotypeGraph extends ReadThreadingGraph {
super(kmerSize);
referenceHaplotype = findReferenceHaplotypeOrFail(haplotypes);
this.haplotypes = new LinkedHashSet<>(haplotypes);
- addSequence("anonymous", referenceHaplotype.getBases(), null, true);
+ addSequence("anonymous", referenceHaplotype.getBases(), true);
for (final Haplotype h : haplotypes) {
if (h.isReference())
continue;
if (h.length() < kmerSize) {
Utils.warnUser(logger, "haplotype shorter than kmerSize " + h.length() + " < " + kmerSize + " will be dropped");
} else
- addSequence("anonymous", h.getBases(), null, false);
+ addSequence("anonymous", h.getBases(), false);
}
buildGraphIfNecessary();
@@ -391,10 +391,9 @@ public class HaplotypeGraph extends ReadThreadingGraph {
graphWriter.println("}");
}
-
@Override
- public Pair findStart(final SequenceForKmers seqForKmers) {
- return getOrCreateKmerVertex(seqForKmers.sequence, 0, true);
+ protected int findStart(final SequenceForKmers seqForKmers) {
+ return 0;
}
/**
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/ReadThreadingAssembler.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/ReadThreadingAssembler.java
index f33a4883f..a932f8a96 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/ReadThreadingAssembler.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/ReadThreadingAssembler.java
@@ -151,14 +151,12 @@ public class ReadThreadingAssembler extends LocalAssemblyEngine {
final ReadThreadingGraph rtgraph = new ReadThreadingGraph(kmerSize, debugGraphTransformations, minBaseQualityToUseInAssembly, numPruningSamples);
// add the reference sequence to the graph
- rtgraph.addSequence("ref", refHaplotype.getBases(), null, true);
+ rtgraph.addSequence("ref", refHaplotype.getBases(), true);
// add the artificial GGA haplotypes to the graph
int hapCount = 0;
for ( final Haplotype h : activeAlleleHaplotypes ) {
- final int[] counts = new int[h.length()];
- Arrays.fill(counts, GGA_MODE_ARTIFICIAL_COUNTS);
- rtgraph.addSequence("activeAllele" + hapCount++, h.getBases(), counts, false);
+ rtgraph.addSequence("activeAllele" + hapCount++, h.getBases(), GGA_MODE_ARTIFICIAL_COUNTS, false);
}
// Next pull kmers out of every read and throw them on the graph
@@ -188,9 +186,10 @@ public class ReadThreadingAssembler extends LocalAssemblyEngine {
// tails that we'll ultimately just trim away anyway, as the dangling tail edges have weight of 1
rtgraph.pruneLowWeightChains(pruneFactor);
- // look at all chains in the graph that terminate in a non-ref node (dangling sinks) and see if
+ // look at all chains in the graph that terminate in a non-ref node (dangling sources and sinks) and see if
// we can recover them by merging some N bases from the chain back into the reference
if ( recoverDanglingTails ) rtgraph.recoverDanglingTails(pruneFactor);
+ if ( recoverDanglingHeads ) rtgraph.recoverDanglingHeads(pruneFactor);
// remove all heading and trailing paths
if ( removePathsNotConnectedToRef ) rtgraph.removePathsNotConnectedToRef();
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/ReadThreadingGraph.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/ReadThreadingGraph.java
index dc057294e..a7989ac2c 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/ReadThreadingGraph.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/haplotypecaller/readthreading/ReadThreadingGraph.java
@@ -46,21 +46,12 @@
package org.broadinstitute.sting.gatk.walkers.haplotypecaller.readthreading;
-import net.sf.samtools.Cigar;
-import net.sf.samtools.CigarElement;
-import net.sf.samtools.CigarOperator;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.gatk.walkers.haplotypecaller.KMerCounter;
import org.broadinstitute.sting.gatk.walkers.haplotypecaller.Kmer;
import org.broadinstitute.sting.gatk.walkers.haplotypecaller.graphs.*;
import org.broadinstitute.sting.utils.BaseUtils;
-import org.broadinstitute.sting.utils.collections.Pair;
-import org.broadinstitute.sting.utils.sam.AlignmentUtils;
import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
-import org.broadinstitute.sting.utils.smithwaterman.SWPairwiseAlignment;
-import org.broadinstitute.sting.utils.smithwaterman.SWParameterSet;
-import org.broadinstitute.sting.utils.smithwaterman.SmithWaterman;
-import org.jgrapht.EdgeFactory;
import org.jgrapht.alg.CycleDetector;
import java.io.File;
@@ -68,28 +59,7 @@ import java.util.*;
import java.util.regex.Matcher;
import java.util.regex.Pattern;
-public class ReadThreadingGraph extends BaseGraph implements KmerSearchableGraph {
-
- /**
- * Edge factory that encapsulates the numPruningSamples assembly parameter
- */
- private static class MyEdgeFactory implements EdgeFactory {
- final int numPruningSamples;
-
- public MyEdgeFactory(int numPruningSamples) {
- this.numPruningSamples = numPruningSamples;
- }
-
- @Override
- public MultiSampleEdge createEdge(final MultiDeBruijnVertex sourceVertex, final MultiDeBruijnVertex targetVertex) {
- return new MultiSampleEdge(false, 1, numPruningSamples);
- }
-
- public MultiSampleEdge createEdge(final boolean isRef, final int multiplicity) {
- return new MultiSampleEdge(isRef, multiplicity, numPruningSamples);
- }
-
- }
+public class ReadThreadingGraph extends DanglingChainMergingGraph implements KmerSearchableGraph {
private final static Logger logger = Logger.getLogger(ReadThreadingGraph.class);
@@ -97,9 +67,6 @@ public class ReadThreadingGraph extends BaseGraph start + kmerSize in seqForKmers
- */
- private int getCountGivenKmerStart(final SequenceForKmers seqForKmers, final int kmerStart) {
- return seqForKmers.getCount(kmerStart + kmerSize - 1);
+ sampleSequences.add(new SequenceForKmers(seqName, sequence, start, stop, count, isRef));
}
/**
@@ -254,12 +216,11 @@ public class ReadThreadingGraph extends BaseGraph startingInfo = findStart(seqForKmers);
- if ( startingInfo == null )
+ final int uniqueStartPos = findStart(seqForKmers);
+ if ( uniqueStartPos == -1 )
return;
- final MultiDeBruijnVertex startingVertex = startingInfo.getFirst();
- final int uniqueStartPos = startingInfo.getSecond();
+ final MultiDeBruijnVertex startingVertex = getOrCreateKmerVertex(seqForKmers.sequence, uniqueStartPos);
// increase the counts of all edges incoming into the starting vertex supported by going back in sequence
if ( increaseCountsBackwards )
@@ -276,185 +237,28 @@ public class ReadThreadingGraph extends BaseGraph danglingPath, referencePath;
- final byte[] danglingPathString, referencePathString;
- final Cigar cigar;
-
- public DanglingTailMergeResult(final List danglingPath,
- final List referencePath,
- final byte[] danglingPathString,
- final byte[] referencePathString,
- final Cigar cigar) {
- this.danglingPath = danglingPath;
- this.referencePath = referencePath;
- this.danglingPathString = danglingPathString;
- this.referencePathString = referencePathString;
- this.cigar = cigar;
- }
- }
-
- /**
- * Attempt to attach vertex with out-degree == 0 to the graph
+ * Find vertex and its position in seqForKmers where we should start assembling seqForKmers
*
- * @param vertex the vertex to recover
- * @param pruneFactor the prune factor to use in ignoring chain pieces
- * @return 1 if we successfully recovered the vertex and 0 otherwise
+ * @param seqForKmers the sequence we want to thread into the graph
+ * @return the position of the starting vertex in seqForKmer, or -1 if it cannot find one
*/
- protected int recoverDanglingChain(final MultiDeBruijnVertex vertex, final int pruneFactor) {
- if ( outDegreeOf(vertex) != 0 ) throw new IllegalStateException("Attempting to recover a dangling tail for " + vertex + " but it has out-degree > 0");
-
- // generate the CIGAR string from Smith-Waterman between the dangling tail and reference paths
- final DanglingTailMergeResult danglingTailMergeResult = generateCigarAgainstReferencePath(vertex, pruneFactor);
-
- // if the CIGAR is too complex (or couldn't be computed) then we do not allow the merge into the reference path
- if ( danglingTailMergeResult == null || ! cigarIsOkayToMerge(danglingTailMergeResult.cigar) )
+ protected int findStart(final SequenceForKmers seqForKmers) {
+ if ( seqForKmers.isRef )
return 0;
- // merge
- return mergeDanglingTail(danglingTailMergeResult);
- }
-
- /**
- * Determine whether the provided cigar is okay to merge into the reference path
- *
- * @param cigar the cigar to analyze
- * @return true if it's okay to merge, false otherwise
- */
- protected boolean cigarIsOkayToMerge(final Cigar cigar) {
-
- final List elements = cigar.getCigarElements();
- final int numElements = elements.size();
-
- // don't allow more than a couple of different ops
- if ( numElements > MAX_CIGAR_COMPLEXITY )
- return false;
-
- // the last element must be an M
- if ( elements.get(numElements - 1).getOperator() != CigarOperator.M )
- return false;
-
- // TODO -- do we want to check whether the Ms mismatch too much also?
-
- return true;
- }
-
- /**
- * Actually merge the dangling tail if possible
- *
- * @param danglingTailMergeResult the result from generating a Cigar for the dangling tail against the reference
- * @return 1 if merge was successful, 0 otherwise
- */
- protected int mergeDanglingTail(final DanglingTailMergeResult danglingTailMergeResult) {
-
- final List elements = danglingTailMergeResult.cigar.getCigarElements();
- final CigarElement lastElement = elements.get(elements.size() - 1);
- if ( lastElement.getOperator() != CigarOperator.M )
- throw new IllegalArgumentException("The last Cigar element must be an M");
-
- final int lastRefIndex = danglingTailMergeResult.cigar.getReferenceLength() - 1;
- final int matchingSuffix = Math.min(GraphUtils.longestSuffixMatch(danglingTailMergeResult.referencePathString, danglingTailMergeResult.danglingPathString, lastRefIndex), lastElement.getLength());
- if ( matchingSuffix == 0 )
- return 0;
-
- final int altIndexToMerge = Math.max(danglingTailMergeResult.cigar.getReadLength() - matchingSuffix - 1, 0);
-
- // there is an important edge condition that we need to handle here: Smith-Waterman correctly calculates that there is a
- // deletion, that deletion is left-aligned such that the LCA node is part of that deletion, and the rest of the dangling
- // tail is a perfect match to the suffix of the reference path. In this case we need to push the reference index to merge
- // down one position so that we don't incorrectly cut a base off of the deletion.
- final boolean firstElementIsDeletion = elements.get(0).getOperator() == CigarOperator.D;
- final boolean mustHandleLeadingDeletionCase = firstElementIsDeletion && (elements.get(0).getLength() + matchingSuffix == lastRefIndex + 1);
- final int refIndexToMerge = lastRefIndex - matchingSuffix + 1 + (mustHandleLeadingDeletionCase ? 1 : 0);
-
- addEdge(danglingTailMergeResult.danglingPath.get(altIndexToMerge), danglingTailMergeResult.referencePath.get(refIndexToMerge), ((MyEdgeFactory)getEdgeFactory()).createEdge(false, 1));
-
- return 1;
- }
-
- /**
- * Generates the CIGAR string from the Smith-Waterman alignment of the dangling path (where the
- * provided vertex is the sink) and the reference path.
- *
- * @param vertex the sink of the dangling tail
- * @param pruneFactor the prune factor to use in ignoring chain pieces
- * @return a SmithWaterman object which can be null if no proper alignment could be generated
- */
- protected DanglingTailMergeResult generateCigarAgainstReferencePath(final MultiDeBruijnVertex vertex, final int pruneFactor) {
-
- // find the lowest common ancestor path between vertex and the reference sink if available
- final List altPath = findPathToLowestCommonAncestorOfReference(vertex, pruneFactor);
- if ( altPath == null || isRefSource(altPath.get(0)) || altPath.size() < MIN_DANGLING_TAIL_LENGTH )
- return null;
-
- // now get the reference path from the LCA
- final List refPath = getReferencePath(altPath.get(0));
-
- // create the Smith-Waterman strings to use
- final byte[] refBases = getBasesForPath(refPath);
- final byte[] altBases = getBasesForPath(altPath);
-
- // run Smith-Waterman to determine the best alignment (and remove trailing deletions since they aren't interesting)
- final SmithWaterman alignment = new SWPairwiseAlignment(refBases, altBases, SWParameterSet.STANDARD_NGS, SWPairwiseAlignment.OVERHANG_STRATEGY.LEADING_INDEL);
- return new DanglingTailMergeResult(altPath, refPath, altBases, refBases, AlignmentUtils.removeTrailingDeletions(alignment.getCigar()));
- }
-
- /**
- * Finds the path upwards in the graph from this vertex to the reference sequence, including the lowest common ancestor vertex.
- * Note that nodes are excluded if their pruning weight is less than the pruning factor.
- *
- * @param vertex the original vertex
- * @param pruneFactor the prune factor to use in ignoring chain pieces
- * @return the path if it can be determined or null if this vertex either doesn't merge onto the reference path or
- * has an ancestor with multiple incoming edges before hitting the reference path
- */
- protected List findPathToLowestCommonAncestorOfReference(final MultiDeBruijnVertex vertex, final int pruneFactor) {
- final LinkedList path = new LinkedList<>();
-
- MultiDeBruijnVertex v = vertex;
- while ( ! isReferenceNode(v) && inDegreeOf(v) == 1 ) {
- final MultiSampleEdge edge = incomingEdgeOf(v);
- // if it has too low a weight, don't use it (or previous vertexes) for the path
- if ( edge.getPruningMultiplicity() < pruneFactor )
- path.clear();
- // otherwise it is safe to use
- else
- path.addFirst(v);
- v = getEdgeSource(edge);
- }
- path.addFirst(v);
-
- return isReferenceNode(v) ? path : null;
- }
-
- /**
- * Finds the path downwards in the graph from this vertex to the reference sink, including this vertex
- *
- * @param start the reference vertex to start from
- * @return the path (non-null, non-empty)
- */
- protected List getReferencePath(final MultiDeBruijnVertex start) {
- if ( ! isReferenceNode(start) ) throw new IllegalArgumentException("Cannot construct the reference path from a vertex that is not on that path");
-
- final List path = new ArrayList<>();
-
- MultiDeBruijnVertex v = start;
- while ( v != null ) {
- path.add(v);
- v = getNextReferenceVertex(v);
+ for ( int i = seqForKmers.start; i < seqForKmers.stop - kmerSize; i++ ) {
+ final Kmer kmer1 = new Kmer(seqForKmers.sequence, i, kmerSize);
+ if ( !nonUniqueKmers.contains(kmer1) )
+ return i;
}
- return path;
+ return -1;
}
/**
@@ -532,26 +336,6 @@ public class ReadThreadingGraph extends BaseGraph uniqueKmers.size();
}
- /**
- * Try to recover dangling tails
- *
- * @param pruneFactor the prune factor to use in ignoring chain pieces
- */
- public void recoverDanglingTails(final int pruneFactor) {
- if ( ! alreadyBuilt ) throw new IllegalStateException("recoverDanglingTails requires the graph be already built");
-
- int attempted = 0;
- int nRecovered = 0;
- for ( final MultiDeBruijnVertex v : vertexSet() ) {
- if ( outDegreeOf(v) == 0 && ! isRefSink(v) ) {
- attempted++;
- nRecovered += recoverDanglingChain(v, pruneFactor);
- }
- }
-
- if ( debugGraphTransformations ) logger.info("Recovered " + nRecovered + " of " + attempted + " dangling tails");
- }
-
/** structure that keeps track of the non-unique kmers for a given kmer size */
private static class NonUniqueResult {
final Set nonUniques;
@@ -574,7 +358,7 @@ public class ReadThreadingGraph extends BaseGraph withNonUniques = getAllPendingSequences();
- final Set nonUniqueKmers = new HashSet();
+ final Set nonUniqueKmers = new HashSet<>();
// go through the sequences and determine which kmers aren't unique within each read
int kmerSize = minKmerSize;
@@ -612,7 +396,7 @@ public class ReadThreadingGraph extends BaseGraph getAllPendingSequences() {
- final LinkedList result = new LinkedList();
+ final LinkedList result = new LinkedList<>();
for ( final List oneSampleWorth : pending.values() ) result.addAll(oneSampleWorth);
return result;
}
@@ -648,7 +432,7 @@ public class ReadThreadingGraph extends BaseGraph vertexMap = new HashMap();
+ final Map vertexMap = new HashMap<>();
// create all of the equivalent seq graph vertices
@@ -683,58 +467,22 @@ public class ReadThreadingGraph extends BaseGraph %s via %s at %d with suffix %s vs. %s",
// prev, vertex, edge, offset, (char)suffix, (char)seqBase));
if ( suffix == seqBase && (increaseCountsThroughBranches || inDegreeOf(vertex) == 1) ) {
- edge.incMultiplicity(seqForKmers.getCount(offset));
+ edge.incMultiplicity(seqForKmers.count);
increaseCountsInMatchedKmers(seqForKmers, prev, originalKmer, offset-1);
}
}
}
- /**
- * Find vertex and its position in seqForKmers where we should start assembling seqForKmers
- *
- * @param seqForKmers the sequence we want to thread into the graph
- * @return a pair of the starting vertex and its position in seqForKmer
- */
- protected Pair findStart(final SequenceForKmers seqForKmers) {
- final int uniqueStartPos = seqForKmers.isRef ? 0 : findUniqueStartPosition(seqForKmers.sequence, seqForKmers.start, seqForKmers.stop);
-
- if ( uniqueStartPos == -1 )
- return null;
-
- return getOrCreateKmerVertex(seqForKmers.sequence, uniqueStartPos, true);
- }
-
- /**
- * Find a starting point in sequence that begins a unique kmer among all kmers in the graph
- * @param sequence the sequence of bases
- * @param start the first base to use in sequence
- * @param stop the last base to use in sequence
- * @return the index into sequence that begins a unique kmer of size kmerSize, or -1 if none could be found
- */
- private int findUniqueStartPosition(final byte[] sequence, final int start, final int stop) {
- for ( int i = start; i < stop - kmerSize; i++ ) {
- final Kmer kmer1 = new Kmer(sequence, i, kmerSize);
- if ( uniqueKmers.containsKey(kmer1) )
- return i;
- }
- return -1;
- }
-
/**
* Get the vertex for the kmer in sequence starting at start
* @param sequence the sequence
* @param start the position of the kmer start
- * @param allowRefSource if true, we will allow matches to the kmer that represents the reference starting kmer
* @return a non-null vertex
*/
- protected Pair getOrCreateKmerVertex(final byte[] sequence, final int start, final boolean allowRefSource) {
+ private MultiDeBruijnVertex getOrCreateKmerVertex(final byte[] sequence, final int start) {
final Kmer kmer = new Kmer(sequence, start, kmerSize);
- final MultiDeBruijnVertex vertex = getUniqueKmerVertex(kmer, allowRefSource);
- if ( vertex != null ) {
- return new Pair<>(vertex, start);
- } else {
- return new Pair<>(createVertex(kmer), start);
- }
+ final MultiDeBruijnVertex vertex = getUniqueKmerVertex(kmer, true);
+ return ( vertex != null ) ? vertex : createVertex(kmer);
}
/**
@@ -780,7 +528,7 @@ public class ReadThreadingGraph extends BaseGraph 1 for reduced reads)
+ * @param count the number of observations of this kmer in graph (can be > 1 for GGA)
* @param isRef is this the reference sequence?
* @return a non-null vertex connecting prevVertex to in the graph based on sequence
*/
@@ -819,7 +567,6 @@ public class ReadThreadingGraph extends BaseGraph= kmerSize ) {
// if the sequence is long enough to get some value out of, add it to the graph
final String name = read.getReadName() + "_" + start + "_" + end;
- addSequence(name, read.getReadGroup().getSample(), read.getReadBases(), start, end, reducedReadCounts, false);
+ addSequence(name, read.getReadGroup().getSample(), read.getReadBases(), start, end, 1, false);
}
lastGood = -1; // reset the last good base
@@ -885,11 +632,11 @@ public class ReadThreadingGraph extends BaseGraph
* Note: only used for testing.
- * Checkout {@link HaplotypeGraphUnitTest} for examples.
+ * Checkout {@see HaplotypeGraphUnitTest} for examples.
*
* @param s the string representation of the graph {@code null}.
*/
@@ -920,7 +667,7 @@ public class ReadThreadingGraph extends BaseGraph
* Note: this is done just for testing purposes.
- * Checkout {@link HaplotypeGraphUnitTest} for examples.
+ * Checkout {@see HaplotypeGraphUnitTest} for examples.
*
* @param str the string representation.
*/
@@ -941,9 +688,9 @@ public class ReadThreadingGraph extends BaseGraph 1 because sequence may be a reduced read and therefore count as N observations
- *
- * @param i the offset into sequence for the start of the kmer
- * @return a count >= 1 that indicates the number of observations of kmer starting at i in this sequence.
- */
- public int getCount(final int i) {
- if ( i < 0 || i > sequence.length ) throw new ArrayIndexOutOfBoundsException("i must be >= 0 and <= " + sequence.length + " but got " + i);
- return counts == null ? 1 : counts[i];
}
}
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModel.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModel.java
index 318779cd2..aa8b46312 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModel.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModel.java
@@ -261,10 +261,9 @@ public class PairHMMIndelErrorModel {
final double downsamplingFraction) {
final int numHaplotypes = haplotypeMap.size();
- final int readCounts[] = new int[pileup.getNumberOfElements()];
- final double[][] readLikelihoods = computeGeneralReadHaplotypeLikelihoods(pileup, haplotypeMap, ref, eventLength, perReadAlleleLikelihoodMap, readCounts);
+ final double[][] readLikelihoods = computeGeneralReadHaplotypeLikelihoods(pileup, haplotypeMap, ref, eventLength, perReadAlleleLikelihoodMap);
perReadAlleleLikelihoodMap.performPerAlleleDownsampling(downsamplingFraction);
- return getDiploidHaplotypeLikelihoods(numHaplotypes, readCounts, readLikelihoods);
+ return getDiploidHaplotypeLikelihoods(numHaplotypes, readLikelihoods);
}
@@ -295,16 +294,13 @@ public class PairHMMIndelErrorModel {
final LinkedHashMap haplotypeMap,
final ReferenceContext ref,
final int eventLength,
- final PerReadAlleleLikelihoodMap perReadAlleleLikelihoodMap,
- final int[] readCounts) {
+ final PerReadAlleleLikelihoodMap perReadAlleleLikelihoodMap) {
final double readLikelihoods[][] = new double[pileup.getNumberOfElements()][haplotypeMap.size()];
final LinkedList readList = new LinkedList<>();
final Map readGCPArrayMap = new LinkedHashMap<>();
int readIdx=0;
for (PileupElement p: pileup) {
- // > 1 when the read is a consensus read representing multiple independent observations
- readCounts[readIdx] = p.getRepresentativeCount();
// check if we've already computed likelihoods for this pileup element (i.e. for this read at this location)
if (perReadAlleleLikelihoodMap.containsPileupElement(p)) {
@@ -499,7 +495,7 @@ public class PairHMMIndelErrorModel {
// return b1.length;
// }
- private static double[] getDiploidHaplotypeLikelihoods(final int numHaplotypes, final int readCounts[], final double readLikelihoods[][]) {
+ private static double[] getDiploidHaplotypeLikelihoods(final int numHaplotypes, final double readLikelihoods[][]) {
final double[][] haplotypeLikehoodMatrix = new double[numHaplotypes][numHaplotypes];
// todo: MAD 09/26/11 -- I'm almost certain this calculation can be simplified to just a single loop without the intermediate NxN matrix
@@ -515,8 +511,7 @@ public class PairHMMIndelErrorModel {
continue;
final double li = readLikelihoods[readIdx][i];
final double lj = readLikelihoods[readIdx][j];
- final int readCount = readCounts[readIdx];
- haplotypeLikehoodMatrix[i][j] += readCount * (MathUtils.approximateLog10SumLog10(li, lj) + MathUtils.LOG_ONE_HALF);
+ haplotypeLikehoodMatrix[i][j] += MathUtils.approximateLog10SumLog10(li, lj) + MathUtils.LOG_ONE_HALF;
}
}
}
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/qc/AssessReducedCoverage.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/qc/AssessReducedCoverage.java
deleted file mode 100644
index 2f8295008..000000000
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/qc/AssessReducedCoverage.java
+++ /dev/null
@@ -1,175 +0,0 @@
-/*
-* By downloading the PROGRAM you agree to the following terms of use:
-*
-* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
-*
-* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
-*
-* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
-* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
-* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
-*
-* 1. DEFINITIONS
-* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
-*
-* 2. LICENSE
-* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
-* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
-* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
-* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
-*
-* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
-* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
-* Copyright 2012 Broad Institute, Inc.
-* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
-* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
-*
-* 4. INDEMNIFICATION
-* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
-*
-* 5. NO REPRESENTATIONS OR WARRANTIES
-* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
-* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
-*
-* 6. ASSIGNMENT
-* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
-*
-* 7. MISCELLANEOUS
-* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
-* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
-* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
-* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
-* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
-* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
-* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
-*/
-
-package org.broadinstitute.sting.gatk.walkers.qc;
-
-import org.broadinstitute.sting.commandline.Argument;
-import org.broadinstitute.sting.commandline.Hidden;
-import org.broadinstitute.sting.commandline.Output;
-import org.broadinstitute.sting.gatk.CommandLineGATK;
-import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
-import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
-import org.broadinstitute.sting.gatk.filters.*;
-import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
-import org.broadinstitute.sting.gatk.walkers.LocusWalker;
-import org.broadinstitute.sting.gatk.walkers.ReadFilters;
-import org.broadinstitute.sting.gatk.walkers.TreeReducible;
-import org.broadinstitute.sting.utils.GenomeLoc;
-import org.broadinstitute.sting.utils.help.DocumentedGATKFeature;
-import org.broadinstitute.sting.utils.pileup.PileupElement;
-import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
-
-import java.io.PrintStream;
-import java.util.HashSet;
-import java.util.Set;
-
-/**
- * Emits intervals present in either the original or reduced bam but not the other.
- *
- * Input
- *
- * The original and reduced BAM files.
- *
- *
- * Output
- *
- * A list of intervals present in one bam but not the other.
- *
- *
- * Examples
- *
- * java -Xmx2g -jar GenomeAnalysisTK.jar \
- * -I:original original.bam \
- * -I:reduced reduced.bam \
- * -R ref.fasta \
- * -T AssessReducedCoverage \
- * -o output.intervals
- *
- *
- * @author ebanks
- */
-@DocumentedGATKFeature( groupName = "Quality Control and Simple Analysis Tools", extraDocs = {CommandLineGATK.class} )
-@ReadFilters({UnmappedReadFilter.class, NotPrimaryAlignmentFilter.class, DuplicateReadFilter.class, FailsVendorQualityCheckFilter.class, BadCigarFilter.class})
-@Hidden
-public class AssessReducedCoverage extends LocusWalker implements TreeReducible {
-
- private static final String original = "original";
- private static final String reduced = "reduced";
-
- @Output
- protected PrintStream out;
-
- @Override
- public boolean includeReadsWithDeletionAtLoci() { return true; }
-
- @Argument(fullName = "output_reduced_only_coverage", shortName = "output_reduced_only_coverage", doc = "Output an interval if the reduced bam has coverage where the original does not", required = false)
- public boolean OUTPUT_REDUCED_ONLY_INTERVALS = false;
-
- public void initialize() {}
-
- public GenomeLoc map(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
-
- if ( tracker == null )
- return null;
-
- final Set tags = getAllTags(context.getBasePileup());
- return (tags.contains(original) && !tags.contains(reduced)) ||
- (OUTPUT_REDUCED_ONLY_INTERVALS && tags.contains(reduced) && !tags.contains(original)) ? ref.getLocus() : null;
- }
-
- private Set getAllTags(final ReadBackedPileup pileup) {
-
- final Set tags = new HashSet(10);
-
- for ( final PileupElement p : pileup ) {
- if ( (int)p.getQual() > 2 && p.getMappingQual() > 0 && !p.isDeletion() )
- tags.addAll(getToolkit().getReaderIDForRead(p.getRead()).getTags().getPositionalTags());
- }
-
- return tags;
- }
-
- public void onTraversalDone(GenomeLoc sum) {
- if ( sum != null )
- out.println(sum);
- }
-
- public GenomeLoc reduceInit() {
- return null;
- }
-
- public GenomeLoc treeReduce(GenomeLoc lhs, GenomeLoc rhs) {
- if ( lhs == null )
- return rhs;
-
- if ( rhs == null )
- return lhs;
-
- // if contiguous, just merge them
- if ( lhs.contiguousP(rhs) )
- return getToolkit().getGenomeLocParser().createGenomeLoc(lhs.getContig(), lhs.getStart(), rhs.getStop());
-
- // otherwise, print the lhs and start over with the rhs
- out.println(lhs);
- return rhs;
- }
-
- public GenomeLoc reduce(GenomeLoc value, GenomeLoc sum) {
- if ( value == null )
- return sum;
-
- if ( sum == null )
- return value;
-
- // if contiguous, just merge them
- if ( sum.contiguousP(value) )
- return getToolkit().getGenomeLocParser().createGenomeLoc(sum.getContig(), sum.getStart(), value.getStop());
-
- // otherwise, print the sum and start over with the value
- out.println(sum);
- return value;
- }
-}
\ No newline at end of file
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/rnaseq/OverhangFixingManager.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/rnaseq/OverhangFixingManager.java
new file mode 100644
index 000000000..581a9e426
--- /dev/null
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/rnaseq/OverhangFixingManager.java
@@ -0,0 +1,380 @@
+/*
+* By downloading the PROGRAM you agree to the following terms of use:
+*
+* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
+*
+* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
+*
+* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
+* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
+* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
+*
+* 1. DEFINITIONS
+* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
+*
+* 2. LICENSE
+* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
+* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
+* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
+* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
+*
+* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
+* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
+* Copyright 2012 Broad Institute, Inc.
+* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
+* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
+*
+* 4. INDEMNIFICATION
+* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
+*
+* 5. NO REPRESENTATIONS OR WARRANTIES
+* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
+* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
+*
+* 6. ASSIGNMENT
+* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
+*
+* 7. MISCELLANEOUS
+* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
+* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
+* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
+* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
+* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
+* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
+* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
+*/
+
+package org.broadinstitute.sting.gatk.walkers.rnaseq;
+
+import net.sf.samtools.SAMFileWriter;
+import net.sf.samtools.SAMRecordCoordinateComparator;
+import org.apache.log4j.Logger;
+import org.broadinstitute.sting.utils.GenomeLoc;
+import org.broadinstitute.sting.utils.GenomeLocParser;
+import org.broadinstitute.sting.utils.clipping.ReadClipper;
+import org.broadinstitute.sting.utils.fasta.CachingIndexedFastaSequenceFile;
+import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
+
+import java.util.*;
+
+/**
+ * The class manages reads and splices and tries to apply overhang clipping when appropriate.
+ * Important note: although for efficiency the manager does try to send reads to the underlying writer in coordinate
+ * sorted order, it does NOT guarantee that it will do so in every case! So unless there's a good reason not to,
+ * methods that instantiate this manager should pass in a writer that does not assume the reads are pre-sorted.
+ */
+public class OverhangFixingManager {
+
+ protected static final Logger logger = Logger.getLogger(OverhangFixingManager.class);
+ private static final boolean DEBUG = false;
+
+ // how many reads should we store in memory before flushing the queue?
+ private final int MAX_RECORDS_IN_MEMORY;
+
+ // how many mismatches do we tolerate in the overhangs?
+ private final int MAX_MISMATCHES_IN_OVERHANG;
+
+ // how many bases do we tolerate in the overhang before deciding not to clip?
+ private final int MAX_BASES_IN_OVERHANG;
+
+ // should we not bother fixing overhangs?
+ private final boolean doNotFixOverhangs;
+
+ // where we ultimately write out our records
+ private final SAMFileWriter writer;
+
+ // fasta reference reader to check overhanging edges in the exome reference sequence
+ private final CachingIndexedFastaSequenceFile referenceReader;
+
+ // the genome loc parser
+ private final GenomeLocParser genomeLocParser;
+
+ // the read cache
+ private final static int initialCapacity = 5000;
+ private PriorityQueue waitingReads = new PriorityQueue<>(initialCapacity, new SplitReadComparator());
+
+ // the set of current splices to use
+ private final Set splices = new TreeSet<>(new SpliceComparator());
+
+ protected static final int MAX_SPLICES_TO_KEEP = 1000;
+
+
+ /**
+ *
+ * @param writer actual writer
+ * @param genomeLocParser the GenomeLocParser object
+ * @param referenceReader the reference reader
+ * @param maxRecordsInMemory max records to keep in memory
+ * @param maxMismatchesInOverhangs max number of mismatches permitted in the overhangs before requiring clipping
+ * @param maxBasesInOverhangs max number of bases permitted in the overhangs before deciding not to clip
+ * @param doNotFixOverhangs if true, don't clip overhangs at all
+ */
+ public OverhangFixingManager(final SAMFileWriter writer,
+ final GenomeLocParser genomeLocParser,
+ final CachingIndexedFastaSequenceFile referenceReader,
+ final int maxRecordsInMemory,
+ final int maxMismatchesInOverhangs,
+ final int maxBasesInOverhangs,
+ final boolean doNotFixOverhangs) {
+ this.writer = writer;
+ this.genomeLocParser = genomeLocParser;
+ this.referenceReader = referenceReader;
+ this.MAX_RECORDS_IN_MEMORY = maxRecordsInMemory;
+ this.MAX_MISMATCHES_IN_OVERHANG = maxMismatchesInOverhangs;
+ this.MAX_BASES_IN_OVERHANG = maxBasesInOverhangs;
+ this.doNotFixOverhangs = doNotFixOverhangs;
+ }
+
+ public final int getNReadsInQueue() { return waitingReads.size(); }
+
+ /**
+ * For testing purposes only
+ *
+ * @return the list of reads currently in the queue
+ */
+ protected List getReadsInQueueForTesting() {
+ return new ArrayList<>(waitingReads);
+ }
+
+ /**
+ * For testing purposes only
+ *
+ * @return the list of splices currently in the queue
+ */
+ protected List getSplicesForTesting() {
+ return new ArrayList<>(splices);
+ }
+
+ /**
+ * Add a new observed split to the list to use
+ *
+ * @param contig the contig
+ * @param start the start of the split, inclusive
+ * @param end the end of the split, inclusive
+ */
+ public void addSplicePosition(final String contig, final int start, final int end) {
+ if ( doNotFixOverhangs )
+ return;
+
+ // is this a new splice? if not, we are done
+ final Splice splice = new Splice(contig, start, end);
+ if ( splices.contains(splice) )
+ return;
+
+ // initialize it with the reference context
+ // we don't want to do this until we know for sure that it's a new splice position
+ splice.initialize(referenceReader);
+
+ // clear the set of old split positions seen if we hit a new contig
+ final boolean sameContig = splices.isEmpty() || splices.iterator().next().loc.getContig().equals(contig);
+ if ( !sameContig )
+ splices.clear();
+
+ // run this position against the existing reads
+ for ( final SplitRead read : waitingReads )
+ fixSplit(read, splice);
+
+ splices.add(splice);
+
+ if ( splices.size() > MAX_SPLICES_TO_KEEP )
+ cleanSplices();
+ }
+
+ /**
+ * Add a read to the manager
+ *
+ * @param read the read to add
+ */
+ public void addRead(final GATKSAMRecord read) {
+ if ( read == null ) throw new IllegalArgumentException("read added to manager is null, which is not allowed");
+
+ // if the new read is on a different contig or we have too many reads, then we need to flush the queue and clear the map
+ final boolean tooManyReads = getNReadsInQueue() >= MAX_RECORDS_IN_MEMORY;
+ final boolean encounteredNewContig = getNReadsInQueue() > 0 && !waitingReads.peek().read.getReferenceIndex().equals(read.getReferenceIndex());
+
+ if ( tooManyReads || encounteredNewContig ) {
+ if ( DEBUG ) logger.warn("Flushing queue on " + (tooManyReads ? "too many reads" : ("move to new contig: " + read.getReferenceName() + " from " + waitingReads.peek().read.getReferenceName())) + " at " + read.getAlignmentStart());
+
+ final int targetQueueSize = encounteredNewContig ? 0 : MAX_RECORDS_IN_MEMORY / 2;
+
+ // write the required number of waiting reads to disk
+ while ( getNReadsInQueue() > targetQueueSize )
+ writer.addAlignment(waitingReads.poll().read);
+ }
+
+ final SplitRead splitRead = new SplitRead(read);
+
+ // fix overhangs, as needed
+ for ( final Splice splice : splices)
+ fixSplit(splitRead, splice);
+
+ // add the new read to the queue
+ waitingReads.add(splitRead);
+ }
+
+ /**
+ * Clean up the list of splices
+ */
+ private void cleanSplices() {
+ final int targetQueueSize = splices.size() / 2;
+ final Iterator iter = splices.iterator();
+ for ( int i = 0; i < targetQueueSize; i++ ) {
+ iter.next();
+ iter.remove();
+ }
+ }
+
+ /**
+ * Try to fix the given read using the given split
+ *
+ * @param read the read to fix
+ * @param splice the split (bad region to clip out)
+ */
+ private void fixSplit(final SplitRead read, final Splice splice) {
+ // if the read doesn't even overlap the split position then we can just exit
+ if ( !splice.loc.overlapsP(read.loc) )
+ return;
+
+ if ( isLeftOverhang(read.loc, splice.loc) ) {
+ final int overhang = splice.loc.getStop() - read.loc.getStart() + 1;
+ if ( overhangingBasesMismatch(read.read.getReadBases(), 0, splice.reference, splice.reference.length - overhang, overhang) ) {
+ final GATKSAMRecord clippedRead = ReadClipper.hardClipByReadCoordinates(read.read, 0, overhang - 1);
+ read.setRead(clippedRead);
+ }
+ }
+ else if ( isRightOverhang(read.loc, splice.loc) ) {
+ final int overhang = read.loc.getStop() - splice.loc.getStart() + 1;
+ if ( overhangingBasesMismatch(read.read.getReadBases(), read.read.getReadLength() - overhang, splice.reference, 0, overhang) ) {
+ final GATKSAMRecord clippedRead = ReadClipper.hardClipByReadCoordinates(read.read, read.read.getReadLength() - overhang, read.read.getReadLength() - 1);
+ read.setRead(clippedRead);
+ }
+ }
+ }
+
+ /**
+ * Is this a proper overhang on the left side of the read?
+ *
+ * @param readLoc the read's loc
+ * @param spliceLoc the split's loc
+ * @return true if it's a left side overhang
+ */
+ protected static boolean isLeftOverhang(final GenomeLoc readLoc, final GenomeLoc spliceLoc) {
+ return readLoc.getStart() <= spliceLoc.getStop() && readLoc.getStart() > spliceLoc.getStart() && readLoc.getStop() > spliceLoc.getStop();
+ }
+
+ /**
+ * Is this a proper overhang on the right side of the read?
+ *
+ * @param readLoc the read's loc
+ * @param spliceLoc the split's loc
+ * @return true if it's a right side overhang
+ */
+ protected static boolean isRightOverhang(final GenomeLoc readLoc, final GenomeLoc spliceLoc) {
+ return readLoc.getStop() >= spliceLoc.getStart() && readLoc.getStop() < spliceLoc.getStop() && readLoc.getStart() < spliceLoc.getStart();
+ }
+
+ /**
+ * Are there too many mismatches to the reference among the overhanging bases?
+ *
+ * @param read the read bases
+ * @param readStartIndex where to start on the read
+ * @param reference the reference bases
+ * @param referenceStartIndex where to start on the reference
+ * @param spanToTest how many bases to test
+ * @return true if too many overhanging bases mismatch, false otherwise
+ */
+ protected boolean overhangingBasesMismatch(final byte[] read,
+ final int readStartIndex,
+ final byte[] reference,
+ final int referenceStartIndex,
+ final int spanToTest) {
+ // don't process too small a span, too large a span, or a span that is most of a read
+ if ( spanToTest < 1 || spanToTest > MAX_BASES_IN_OVERHANG || spanToTest > read.length / 2 )
+ return false;
+
+ int numMismatchesSeen = 0;
+ for ( int i = 0; i < spanToTest; i++ ) {
+ if ( read[readStartIndex + i] != reference[referenceStartIndex + i] ) {
+ if ( ++numMismatchesSeen > MAX_MISMATCHES_IN_OVERHANG )
+ return true;
+ }
+ }
+
+ // we can still mismatch overall if at least half of the bases mismatch
+ return numMismatchesSeen >= ((spanToTest+1)/2);
+ }
+
+ /**
+ * Close out the manager stream by clearing the read cache
+ */
+ public void close() {
+ // write out all of the remaining reads
+ while ( ! waitingReads.isEmpty() )
+ writer.addAlignment(waitingReads.poll().read);
+ }
+
+ // class to represent the reads with their soft-clip-included GenomeLocs
+ protected final class SplitRead {
+
+ public GATKSAMRecord read;
+ public GenomeLoc loc;
+
+ public SplitRead(final GATKSAMRecord read) {
+ setRead(read);
+ }
+
+ public void setRead(final GATKSAMRecord read) {
+ if ( !read.isEmpty() ) {
+ this.read = read;
+ loc = genomeLocParser.createGenomeLoc(read.getReferenceName(), read.getSoftStart(), read.getSoftEnd());
+ }
+ }
+ }
+
+ // class to represent the comparator for the split reads
+ private final class SplitReadComparator implements Comparator {
+
+ private final SAMRecordCoordinateComparator readComparator;
+
+ public SplitReadComparator() {
+ readComparator = new SAMRecordCoordinateComparator();
+ }
+
+ public int compare(final SplitRead read1, final SplitRead read2) {
+ return readComparator.compare(read1.read, read2.read);
+ }
+ }
+
+ // class to represent the split positions
+ protected final class Splice {
+
+ public final GenomeLoc loc;
+ public byte[] reference;
+
+ public Splice(final String contig, final int start, final int end) {
+ loc = genomeLocParser.createGenomeLoc(contig, start, end);
+ }
+
+ public void initialize(final CachingIndexedFastaSequenceFile referenceReader) {
+ reference = referenceReader.getSubsequenceAt(loc.getContig(), loc.getStart(), loc.getStop()).getBases();
+ }
+
+ @Override
+ public boolean equals(final Object other) {
+ return other != null && (other instanceof Splice) && this.loc.equals(((Splice)other).loc);
+ }
+
+ @Override
+ public int hashCode() {
+ return loc.hashCode();
+ }
+ }
+
+ // class to represent the comparator for the split reads
+ private final class SpliceComparator implements Comparator {
+
+ public int compare(final Splice position1, final Splice position2) {
+ return position1.loc.compareTo(position2.loc);
+ }
+ }
+}
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/rnaseq/SplitNCigarReads.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/rnaseq/SplitNCigarReads.java
index 223745431..6b9fca312 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/rnaseq/SplitNCigarReads.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/gatk/walkers/rnaseq/SplitNCigarReads.java
@@ -47,30 +47,28 @@
package org.broadinstitute.sting.gatk.walkers.rnaseq;
import net.sf.samtools.*;
+import org.broadinstitute.sting.commandline.Advanced;
import org.broadinstitute.sting.commandline.Argument;
+import org.broadinstitute.sting.commandline.Hidden;
import org.broadinstitute.sting.commandline.Output;
import org.broadinstitute.sting.gatk.CommandLineGATK;
import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.io.StingSAMFileWriter;
-import org.broadinstitute.sting.gatk.iterators.ReadTransformer;
-import org.broadinstitute.sting.gatk.iterators.ReadTransformersMode;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
-import org.broadinstitute.sting.gatk.walkers.BAQMode;
import org.broadinstitute.sting.gatk.walkers.DataSource;
import org.broadinstitute.sting.gatk.walkers.ReadWalker;
import org.broadinstitute.sting.gatk.walkers.Requires;
import org.broadinstitute.sting.utils.Utils;
-import org.broadinstitute.sting.utils.baq.BAQ;
import org.broadinstitute.sting.utils.clipping.ReadClipper;
import org.broadinstitute.sting.utils.exceptions.UserException;
+import org.broadinstitute.sting.utils.fasta.CachingIndexedFastaSequenceFile;
import org.broadinstitute.sting.utils.help.DocumentedGATKFeature;
import org.broadinstitute.sting.utils.help.HelpConstants;
import org.broadinstitute.sting.utils.sam.CigarUtils;
import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
-import java.io.PrintStream;
-import java.util.*;
+import java.io.FileNotFoundException;
/**
*
@@ -88,119 +86,138 @@ import java.util.*;
@DocumentedGATKFeature( groupName = HelpConstants.DOCS_CAT_DATA, extraDocs = {CommandLineGATK.class} )
@Requires({DataSource.READS, DataSource.REFERENCE})
-public class SplitNCigarReads extends ReadWalker, SAMFileWriter> {
+public class SplitNCigarReads extends ReadWalker {
+
+ // The name that will go in the @PG tag
+ public static final String PROGRAM_RECORD_NAME = "GATK SplitNCigarReads";
+
@Output(doc="Write output to this BAM filename instead of STDOUT")
- StingSAMFileWriter out;
-
- @Argument(required = false)
- PrintStream splitPositionsOutput = System.out;
-
- @Argument(fullName="outputAsBED", shortName="bed", doc="Output as BED file", required=false)
- boolean outputAsBED = false;
-
- @Argument(fullName="printSplitPositions", shortName="splitPosition", doc="print the split positions", required=false)
- boolean printSplitPositions = false;
-
- public static final String PROGRAM_RECORD_NAME = "GATK SplitNCigarReads"; // The name that will go in the @PG tag
- // public static SplitPositions splitPositions = null;
- public static String results = "";
+ protected StingSAMFileWriter writer;
/**
- * The initialize function.
+ * For expert users only! To minimize memory consumption you can lower this number, but then the tool may skip
+ * overhang fixing in regions with too much coverage. Just make sure to give Java enough memory! 4Gb should be
+ * enough with the default value.
*/
+ @Advanced
+ @Argument(fullName="maxReadsInMemory", shortName="maxInMemory", doc="max reads allowed to be kept in memory at a time by the BAM writer", required=false)
+ protected int MAX_RECORDS_IN_MEMORY = 150000;
+
+ /**
+ * If there are more than this many mismatches within the overhang regions, the whole overhang will get hard-clipped out.
+ * It is still possible in some cases that the overhang could get clipped if the number of mismatches do not exceed this
+ * value, e.g. if most of the overhang mismatches.
+ */
+ @Advanced
+ @Argument(fullName="maxMismatchesInOverhang", shortName="maxMismatches", doc="max number of mismatches allowed in the overhang", required=false)
+ protected int MAX_MISMATCHES_IN_OVERHANG = 1;
+
+ /**
+ * If there are more than this many bases in the overhang, we won't try to hard-clip them out
+ */
+ @Advanced
+ @Argument(fullName="maxBasesInOverhang", shortName="maxOverhang", doc="max number of bases allowed in the overhang", required=false)
+ protected int MAX_BASES_TO_CLIP = 40;
+
+ @Argument(fullName="doNotFixOverhangs", shortName="doNotFixOverhangs", doc="do not have the walker hard-clip overhanging sections of the reads", required=false)
+ protected boolean doNotFixOverhangs = false;
+
+ @Hidden
+ @Argument(fullName = "no_pg_tag", shortName = "npt", doc = "Necessary for integration tests", required = false)
+ protected boolean NO_PG_TAG = false;
+
+ /**
+ * This stores all of the already-split reads and manages any processing (e.g. clipping overhangs) that happens to them.
+ * It will emit reads to the underlying writer as needed so we don't need to worry about any of that in this class.
+ */
+ protected OverhangFixingManager overhangManager;
+
+
+ @Override
public void initialize() {
final GenomeAnalysisEngine toolkit = getToolkit();
- final boolean preSorted = false;
- if (getToolkit() != null) {
- Utils.setupWriter(out, toolkit, toolkit.getSAMFileHeader(), preSorted, this, PROGRAM_RECORD_NAME);
+ if ( !NO_PG_TAG ) {
+ // we don't want to assume that reads will be written in order by the manager because in deep, deep pileups it won't work
+ Utils.setupWriter(writer, toolkit, toolkit.getSAMFileHeader(), false, this, PROGRAM_RECORD_NAME);
+ }
+
+ try {
+ final CachingIndexedFastaSequenceFile referenceReader = new CachingIndexedFastaSequenceFile(toolkit.getArguments().referenceFile);
+ overhangManager = new OverhangFixingManager(writer, toolkit.getGenomeLocParser(), referenceReader, MAX_RECORDS_IN_MEMORY, MAX_MISMATCHES_IN_OVERHANG, MAX_BASES_TO_CLIP, doNotFixOverhangs);
+ }
+ catch (FileNotFoundException ex) {
+ throw new UserException.CouldNotReadInputFile(toolkit.getArguments().referenceFile, ex);
}
- out.setPresorted(preSorted);
- // splitPositions = new SplitPositions();
}
- /**
- * The reads map function.
- *
- * @param ref the reference bases that correspond to our read, if a reference was provided
- * @param read the read itself, as a GATKSAMRecord
- *
- * @return a list of split read if there are N's in the cigar string, or the read itself.
- */
- public List map(final ReferenceContext ref,final GATKSAMRecord read,final RefMetaDataTracker metaDataTracker) {
- return splitNCigarRead(read);
+ @Override
+ public GATKSAMRecord map(final ReferenceContext ref, final GATKSAMRecord read, final RefMetaDataTracker metaDataTracker) {
+ return read;
+ }
+
+ @Override
+ public OverhangFixingManager reduceInit() {
+ return overhangManager;
+ }
+
+ @Override
+ public OverhangFixingManager reduce(final GATKSAMRecord read, final OverhangFixingManager manager) {
+ splitNCigarRead(read, manager);
+ return manager;
+ }
+
+ @Override
+ public void onTraversalDone(final OverhangFixingManager manager) {
+ manager.close();
}
/**
* Goes through the cigar string of the read and create new reads for each consecutive non-N elements (while hard clipping the rest of the read).
* For example: for a read with cigar '1H2M2D1M2N1M2I1N1M2S' 3 new reads will be created with cigar strings: 1H2M2D1M, 1M2I and 1M2S
*
- * @param read the read to split, as a GATKSAMRecord
- * @return a list of split read if there are N's in the cigar string, or the read itself.
+ * @param read the read to split
+ * @param manager the output manager
*/
+ public static void splitNCigarRead(final GATKSAMRecord read, final OverhangFixingManager manager) {
+ final int numCigarElements = read.getCigar().numCigarElements();
- public static List splitNCigarRead(final GATKSAMRecord read){
- final List splitReads = new ArrayList<>();
int firstCigarIndex = 0;
- for (int i = 0; i < read.getCigar().numCigarElements(); i ++){
+ for ( int i = 0; i < numCigarElements; i++ ) {
final CigarElement cigarElement = read.getCigar().getCigarElement(i);
- if(cigarElement.getOperator() == CigarOperator.N){
- final boolean addToSplitPositions = true;
- splitReads.add(splitReadBasedOnCigar(read,firstCigarIndex,i, addToSplitPositions));
+ if (cigarElement.getOperator() == CigarOperator.N) {
+ manager.addRead(splitReadBasedOnCigar(read, firstCigarIndex, i, manager));
firstCigarIndex = i+1;
}
}
- //if there are no N's in the read
- if (firstCigarIndex == 0)
- splitReads.add(read);
+ // if there are no N's in the read
+ if (firstCigarIndex == 0) {
+ manager.addRead(read);
+ }
//add the last section of the read: from the last N to the the end of the read
// (it will be done for all the usual cigar string that does not end with N)
- else if(firstCigarIndex < read.getCigar().numCigarElements()){
- final boolean addToSplitPositions = false;
- splitReads.add(splitReadBasedOnCigar(read,firstCigarIndex,read.getCigar().numCigarElements(), addToSplitPositions));
+ else if (firstCigarIndex < numCigarElements) {
+ manager.addRead(splitReadBasedOnCigar(read, firstCigarIndex, numCigarElements, null));
}
- return splitReads;
- }
-
-
- /**
- * reduceInit is called once before any calls to the map function. We use it here to setup the splitPositionsOutput
- * bam file, if it was specified on the command line
- *
- * @return SAMFileWriter, set to the BAM splitPositionsOutput file if the command line option was set, null otherwise
- */
- public SAMFileWriter reduceInit() {
- return out;
}
/**
- * given a read and a splitPositionsOutput location, reduce by emitting the read
+ * Pull out an individual split position for a read
*
- * @param reads the split reads itself
- * @param output the splitPositionsOutput source
- * @return the SAMFileWriter, so that the next reduce can emit to the same source
+ * @param read the read being split
+ * @param cigarStartIndex the index of the first cigar element to keep
+ * @param cigarEndIndex the index of the last cigar element to keep
+ * @param forSplitPositions the manager for keeping track of split positions; can be null
+ * @return a non-null read representing the section of the original read being split out
*/
- public SAMFileWriter reduce(final List reads,final SAMFileWriter output ) {
- for (final GATKSAMRecord read: reads)
- output.addAlignment(read);
- return output;
- }
-
- public void onTraversalDone(SAMFileWriter readResult) {
- super.onTraversalDone(readResult);
- if(printSplitPositions)
- splitPositionsOutput.println(results);
- // splitPositionsOutput.println(splitPositions);
-
- }
- private static GATKSAMRecord splitReadBasedOnCigar(final GATKSAMRecord read, final int cigarStartIndex, final int cigarEndIndex, final boolean addToSplitPositions){
+ private static GATKSAMRecord splitReadBasedOnCigar(final GATKSAMRecord read, final int cigarStartIndex, final int cigarEndIndex, final OverhangFixingManager forSplitPositions) {
int cigarFirstIndex = cigarStartIndex;
int cigarSecondIndex = cigarEndIndex;
- //in case a section of the read is end or start with D (for example the first section in 1M1D1N1M is 1M1D), we should trim this cigar element
- // it can be if, but it was kept as while to make sure the code can work with Cigar string that were not "cleaned"
+ //in case a section of the read ends or starts with D (for example the first section in 1M1D1N1M is 1M1D), we should trim this cigar element
+ // it can be 'if', but it was kept as 'while' to make sure the code can work with Cigar strings that were not "cleaned"
while(read.getCigar().getCigarElement(cigarFirstIndex).getOperator().equals(CigarOperator.D))
cigarFirstIndex++;
while(read.getCigar().getCigarElement(cigarSecondIndex-1).getOperator().equals(CigarOperator.D))
@@ -213,54 +230,13 @@ public class SplitNCigarReads extends ReadWalker, SAMFileWri
final int startRefIndex = read.getOriginalAlignmentStart() + CigarUtils.countRefBasesBasedOnCigar(read,0,cigarFirstIndex); //goes through the prefix of the cigar (up to cigarStartIndex) and move the reference index.
final int stopRefIndex = startRefIndex + CigarUtils.countRefBasesBasedOnCigar(read,cigarFirstIndex,cigarSecondIndex)-1; //goes through a consecutive non-N section of the cigar (up to cigarEndIndex) and move the reference index.
- if(addToSplitPositions){
- final int splitPosition = startRefIndex + CigarUtils.countRefBasesBasedOnCigar(read,cigarFirstIndex,cigarEndIndex); //we use cigarEndIndex instead of cigarSecondIndex so we won't take into account the D's at the end.
+ if ( forSplitPositions != null ) {
final String contig = read.getReferenceName();
-// results += String.format("%s:%d-%d\n", contig, splitPosition, splitPosition );
-// splitPositions.addSplitPosition(contig,splitPosition);
+ final int splitStart = startRefIndex + CigarUtils.countRefBasesBasedOnCigar(read,cigarFirstIndex,cigarEndIndex); //we use cigarEndIndex instead of cigarSecondIndex so we won't take into account the D's at the end.
+ final int splitEnd = splitStart + read.getCigar().getCigarElement(cigarEndIndex).getLength() - 1;
+ forSplitPositions.addSplicePosition(contig, splitStart, splitEnd);
}
return ReadClipper.hardClipToRegionIncludingClippedBases(read, startRefIndex, stopRefIndex);
-
}
-
- private class SplitPosition {
- public final String contig;
- public final int start;
- public final int end;
-
- public SplitPosition(final String c, final int position) {
- contig = c;
- start = position;
- end = position;
- }
- }
-
-
- private class SplitPositions {
- private final HashSet splitPositions;
-
- public SplitPositions() {
- splitPositions = new HashSet<>();
- }
-
- public void addSplitPosition(final String contig, final int position) {
- final SplitPosition newSplitPosition = new SplitPosition(contig, position);
- splitPositions.add(newSplitPosition);
- }
-
- public String toString() {
- String result = ""; // = "Contig\tstart\tstop\n";
-
- for (SplitPosition position: splitPositions) {
- if (outputAsBED)
- result += String.format("%s\t%d\t%d\n", position.contig, position.start-1, position.end );
- else
- result += String.format("%s:%d-%d\n", position.contig, position.start, position.end );
- }
- return result;
- }
- }
-
-
}
diff --git a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/utils/recalibration/RecalUtils.java b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/utils/recalibration/RecalUtils.java
index 56f7e8257..325237d05 100644
--- a/protected/gatk-protected/src/main/java/org/broadinstitute/sting/utils/recalibration/RecalUtils.java
+++ b/protected/gatk-protected/src/main/java/org/broadinstitute/sting/utils/recalibration/RecalUtils.java
@@ -48,12 +48,10 @@ package org.broadinstitute.sting.utils.recalibration;
import com.google.java.contract.Ensures;
import com.google.java.contract.Requires;
-import net.sf.samtools.SAMFileHeader;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.gatk.report.GATKReport;
import org.broadinstitute.sting.gatk.report.GATKReportTable;
import org.broadinstitute.sting.gatk.walkers.bqsr.RecalibrationArgumentCollection;
-import org.broadinstitute.sting.gatk.walkers.compression.reducereads.ReduceReads;
import org.broadinstitute.sting.utils.classloader.JVMUtils;
import org.broadinstitute.sting.utils.recalibration.covariates.*;
import org.broadinstitute.sting.utils.BaseUtils;
@@ -1063,20 +1061,4 @@ public class RecalUtils {
private static RecalDatum createDatumObject(final byte reportedQual, final double isError) {
return new RecalDatum(1, isError, reportedQual);
}
-
- /**
- * Checks for invalid BAMs that are being used with BQSR and fails with a UserException if it finds one
- *
- * @param headers sam file headers being passed into the GATK engine
- * @param allowBqsrOnReducedBams should we allow BQSR on reduced bams?
- */
- public static void checkForInvalidRecalBams(final List headers, final boolean allowBqsrOnReducedBams) {
- // for now, the only check we make is against reduced bams
- if ( !allowBqsrOnReducedBams ) {
- for ( final SAMFileHeader header : headers ) {
- if ( header.getProgramRecord(ReduceReads.PROGRAM_RECORD_NAME) != null )
- throw new UserException.BadInput("base quality score recalibration should absolutely not be run on reduced BAM files! Please run ReduceReads only after BQSR has been performed");
- }
- }
- }
}
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/annotator/RankSumUnitTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/annotator/RankSumUnitTest.java
index fec83e1a8..b1c280748 100644
--- a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/annotator/RankSumUnitTest.java
+++ b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/annotator/RankSumUnitTest.java
@@ -46,8 +46,6 @@
package org.broadinstitute.sting.gatk.walkers.annotator;
-import org.broadinstitute.sting.gatk.walkers.compression.reducereads.*;
-import org.broadinstitute.sting.gatk.walkers.compression.reducereads.BaseCounts;
import org.broadinstitute.sting.utils.MannWhitneyU;
import org.testng.Assert;
import org.testng.annotations.BeforeClass;
@@ -122,13 +120,15 @@ public class RankSumUnitTest {
final List dist2 = new ArrayList<>(distribution2);
if ( numToReduceIn2 > 0 ) {
- final org.broadinstitute.sting.gatk.walkers.compression.reducereads.BaseCounts counts = new BaseCounts();
+ int counts = 0;
+ int quals = 0;
+
for ( int i = 0; i < numToReduceIn2; i++ ) {
- final int value = dist2.remove(0);
- counts.incr(BaseIndex.A, (byte)value, 0, false);
+ counts++;
+ quals += dist2.remove(0);
}
- final int qual = (int)counts.averageQualsOfBase(BaseIndex.A);
+ final int qual = quals / counts;
for ( int i = 0; i < numToReduceIn2; i++ )
dist2.add(qual);
}
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/annotator/VariantAnnotatorIntegrationTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/annotator/VariantAnnotatorIntegrationTest.java
index 287cd45d0..0ed7eb2e8 100644
--- a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/annotator/VariantAnnotatorIntegrationTest.java
+++ b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/annotator/VariantAnnotatorIntegrationTest.java
@@ -358,7 +358,7 @@ public class VariantAnnotatorIntegrationTest extends WalkerTest {
final File outputVCFNoQD = executeTest("testQualByDepth calling without QD", specNoQD).getFirst().get(0);
final String baseAnn = String.format("-T VariantAnnotator -R %s -V %s", REF, outputVCFNoQD.getAbsolutePath()) + " --no_cmdline_in_header -o %s -L 20:10130000-10134800 -A QualByDepth";
- final WalkerTestSpec specAnn = new WalkerTestSpec(baseAnn, 1, Arrays.asList("4ccdbebcfd02be87ae5b4ad94666f011"));
+ final WalkerTestSpec specAnn = new WalkerTestSpec(baseAnn, 1, Arrays.asList("388200d107fb47326df78a971a52698f"));
specAnn.disableShadowBCF();
final File outputVCFAnn = executeTest("testQualByDepth re-annotation of QD", specAnn).getFirst().get(0);
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/BaseCountsUnitTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/BaseCountsUnitTest.java
deleted file mode 100644
index f988471a0..000000000
--- a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/BaseCountsUnitTest.java
+++ /dev/null
@@ -1,201 +0,0 @@
-/*
-* By downloading the PROGRAM you agree to the following terms of use:
-*
-* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
-*
-* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
-*
-* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
-* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
-* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
-*
-* 1. DEFINITIONS
-* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
-*
-* 2. LICENSE
-* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
-* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
-* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
-* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
-*
-* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
-* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
-* Copyright 2012 Broad Institute, Inc.
-* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
-* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
-*
-* 4. INDEMNIFICATION
-* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
-*
-* 5. NO REPRESENTATIONS OR WARRANTIES
-* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
-* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
-*
-* 6. ASSIGNMENT
-* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
-*
-* 7. MISCELLANEOUS
-* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
-* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
-* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
-* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
-* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
-* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
-* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
-*/
-
-package org.broadinstitute.sting.gatk.walkers.compression.reducereads;
-
-
-import org.broadinstitute.sting.BaseTest;
-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.List;
-
-/**
- * Basic unit test for BaseCounts in reduced reads
- */
-public class BaseCountsUnitTest extends BaseTest {
-
- private class BaseCountsTest {
- public String bases;
- public byte mostCountBase;
- public int mostCommonCount;
-
- private BaseCountsTest(String bases, char mostCountBase, int mostCommonCount) {
- this.mostCommonCount = mostCommonCount;
- this.mostCountBase = (byte)mostCountBase;
- this.bases = bases;
- }
- }
-
- @DataProvider(name = "counting")
- public Object[][] createCountingData() {
- List params = new ArrayList();
-
- params.add(new BaseCountsTest("A", 'A', 1 ));
- params.add(new BaseCountsTest("AA", 'A', 2 ));
- params.add(new BaseCountsTest("AC", 'A', 1 ));
- params.add(new BaseCountsTest("AAC", 'A', 2 ));
- params.add(new BaseCountsTest("AAA", 'A', 3 ));
- params.add(new BaseCountsTest("AAAN", 'A', 3 ));
- params.add(new BaseCountsTest("AAANNNN", 'N', 4 ));
- params.add(new BaseCountsTest("AACTG", 'A', 2 ));
- params.add(new BaseCountsTest("D", 'D', 1 ));
- params.add(new BaseCountsTest("DDAAD", 'D', 3));
- params.add(new BaseCountsTest("", (char)BaseCounts.MAX_BASE_WITH_NO_COUNTS, 0 ));
- params.add(new BaseCountsTest("AAIIIAI", 'I', 4 ));
-
- List params2 = new ArrayList();
- for ( BaseCountsTest x : params ) params2.add(new Object[]{x});
- return params2.toArray(new Object[][]{});
- }
-
- @Test(dataProvider = "counting", enabled = true)
- public void testCounting(BaseCountsTest params) {
- BaseCounts counts = new BaseCounts();
-
- for ( byte base : params.bases.getBytes() )
- counts.incr(base);
-
- String name = String.format("Test-%s", params.bases);
- Assert.assertEquals(counts.totalCount(), params.bases.length(), name);
- Assert.assertEquals(counts.countOfBase(counts.baseIndexWithMostCounts()), params.mostCommonCount, name);
- Assert.assertEquals((char)counts.baseWithMostCounts(), (char)params.mostCountBase, name);
-
- // test the static creation
- final int[] countsArray = new int[] { counts.countOfBase(BaseIndex.A), counts.countOfBase(BaseIndex.C),
- counts.countOfBase(BaseIndex.G), counts.countOfBase(BaseIndex.T)};
- final BaseCounts countsFromArray = BaseCounts.createWithCounts(countsArray);
- Assert.assertEquals(counts.countOfBase(BaseIndex.A), countsFromArray.countOfBase(BaseIndex.A));
- Assert.assertEquals(counts.countOfBase(BaseIndex.C), countsFromArray.countOfBase(BaseIndex.C));
- Assert.assertEquals(counts.countOfBase(BaseIndex.G), countsFromArray.countOfBase(BaseIndex.G));
- Assert.assertEquals(counts.countOfBase(BaseIndex.T), countsFromArray.countOfBase(BaseIndex.T));
- Assert.assertEquals(ACGTcounts(counts), countsFromArray.totalCount());
-
- // test addition
- counts.add(countsFromArray);
- Assert.assertEquals(counts.countOfBase(BaseIndex.A), 2 * countsFromArray.countOfBase(BaseIndex.A));
- Assert.assertEquals(counts.countOfBase(BaseIndex.C), 2 * countsFromArray.countOfBase(BaseIndex.C));
- Assert.assertEquals(counts.countOfBase(BaseIndex.G), 2 * countsFromArray.countOfBase(BaseIndex.G));
- Assert.assertEquals(counts.countOfBase(BaseIndex.T), 2 * countsFromArray.countOfBase(BaseIndex.T));
- Assert.assertEquals(ACGTcounts(counts), 2 * countsFromArray.totalCount());
-
- // test subtraction
- counts.sub(countsFromArray);
- Assert.assertEquals(counts.countOfBase(BaseIndex.A), countsFromArray.countOfBase(BaseIndex.A));
- Assert.assertEquals(counts.countOfBase(BaseIndex.C), countsFromArray.countOfBase(BaseIndex.C));
- Assert.assertEquals(counts.countOfBase(BaseIndex.G), countsFromArray.countOfBase(BaseIndex.G));
- Assert.assertEquals(counts.countOfBase(BaseIndex.T), countsFromArray.countOfBase(BaseIndex.T));
- Assert.assertEquals(ACGTcounts(counts), countsFromArray.totalCount());
-
- // test decrementing
- if ( counts.countOfBase(BaseIndex.A) > 0 ) {
- counts.decr((byte)'A');
- Assert.assertEquals(counts.countOfBase(BaseIndex.A), countsFromArray.countOfBase(BaseIndex.A) - 1);
- }
- }
-
- private static int ACGTcounts(final BaseCounts baseCounts) {
- return baseCounts.totalCountWithoutIndels() - baseCounts.countOfBase(BaseIndex.N);
- }
-
-
- //////////////////////////////////
- // TEST FOR QUALS IN BASECOUNTS //
- //////////////////////////////////
-
- private class BaseCountsQualsTest {
- public final List quals;
-
- private BaseCountsQualsTest(final List quals) {
- this.quals = quals;
- }
- }
-
- @DataProvider(name = "quals")
- public Object[][] createQualsData() {
- List tests = new ArrayList();
-
- final int[] quals = new int[]{ 0, 5, 10, 15, 20, 30, 40, 50 };
-
- for ( final int qual1 : quals ) {
- for ( final int qual2 : quals ) {
- for ( final int qual3 : quals ) {
- tests.add(new Object[]{new BaseCountsQualsTest(Arrays.asList(qual1, qual2, qual3))});
- }
- }
- }
-
- return tests.toArray(new Object[][]{});
- }
-
- @Test(dataProvider = "quals", enabled = true)
- public void testQuals(BaseCountsQualsTest test) {
- BaseCounts counts = new BaseCounts();
-
- for ( int qual : test.quals )
- counts.incr(BaseIndex.A, (byte)qual, 20, false);
-
- final int actualSum = (int)counts.getSumQuals((byte)'A');
- final int expectedSum = qualSum(test.quals);
- Assert.assertEquals(actualSum, expectedSum);
-
- final int actualAverage = (int)counts.averageQuals((byte)'A');
- Assert.assertEquals(actualAverage, expectedSum / test.quals.size());
-
- // test both proportion methods
- Assert.assertEquals(counts.baseCountProportion(BaseIndex.A), counts.baseCountProportion((byte)'A'));
- }
-
- private static int qualSum(final List quals) {
- int sum = 0;
- for ( final int qual : quals )
- sum += qual;
- return sum;
- }
-}
\ No newline at end of file
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/HeaderElementUnitTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/HeaderElementUnitTest.java
deleted file mode 100644
index 4f5b7477c..000000000
--- a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/HeaderElementUnitTest.java
+++ /dev/null
@@ -1,214 +0,0 @@
-/*
-* By downloading the PROGRAM you agree to the following terms of use:
-*
-* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
-*
-* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
-*
-* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
-* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
-* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
-*
-* 1. DEFINITIONS
-* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
-*
-* 2. LICENSE
-* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
-* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
-* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
-* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
-*
-* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
-* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
-* Copyright 2012 Broad Institute, Inc.
-* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
-* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
-*
-* 4. INDEMNIFICATION
-* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
-*
-* 5. NO REPRESENTATIONS OR WARRANTIES
-* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
-* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
-*
-* 6. ASSIGNMENT
-* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
-*
-* 7. MISCELLANEOUS
-* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
-* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
-* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
-* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
-* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
-* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
-* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
-*/
-
-package org.broadinstitute.sting.gatk.walkers.compression.reducereads;
-
-
-import org.broadinstitute.sting.BaseTest;
-import org.broadinstitute.sting.utils.MathUtils;
-import org.testng.Assert;
-import org.testng.annotations.DataProvider;
-import org.testng.annotations.Test;
-
-import java.util.ArrayList;
-import java.util.List;
-
-public class HeaderElementUnitTest extends BaseTest {
-
- private class HETest {
- public byte base, baseQual, insQual, delQual;
- public int MQ;
- public boolean isClip;
-
- private HETest(final byte base, final byte baseQual, final byte insQual, final byte delQual, final int MQ, final boolean isClip) {
- this.base = base;
- this.baseQual = baseQual;
- this.insQual = insQual;
- this.delQual = delQual;
- this.MQ = MQ;
- this.isClip = isClip;
- }
- }
-
- private static final byte byteA = (byte)'A';
- private static final byte byte10 = (byte)10;
- private static final byte byte20 = (byte)20;
- private static final int minBaseQual = 20;
- private static final int minMappingQual = 20;
-
- @DataProvider(name = "data")
- public Object[][] createData() {
- List tests = new ArrayList();
-
- tests.add(new Object[]{new HETest(byteA, byte20, byte20, byte20, 20, false)});
- tests.add(new Object[]{new HETest(byteA, byte10, byte20, byte20, 20, false)});
- tests.add(new Object[]{new HETest(byteA, byte20, byte20, byte20, 10, false)});
- tests.add(new Object[]{new HETest(byteA, byte20, byte20, byte20, 20, true)});
-
- return tests.toArray(new Object[][]{});
- }
-
- @Test(dataProvider = "data", enabled = true)
- public void testHE(HETest test) {
-
- HeaderElement headerElement = new HeaderElement(1000, 0);
-
- // first test that if we add and then remove it, we have no data
- headerElement.addBase(test.base, test.baseQual, test.insQual, test.delQual, test.MQ, minBaseQual, minMappingQual, test.isClip, false);
- headerElement.addInsertionToTheRight();
- headerElement.removeBase(test.base, test.baseQual, test.insQual, test.delQual, test.MQ, minBaseQual, minMappingQual, test.isClip, false);
- headerElement.removeInsertionToTheRight();
- testHeaderIsEmpty(headerElement);
-
- // now, test that the data was added as expected
- for ( int i = 0; i < 10; i++ )
- headerElement.addBase(test.base, test.baseQual, test.insQual, test.delQual, test.MQ, minBaseQual, minMappingQual, test.isClip, false);
- testHeaderData(headerElement, test);
-
- // test the insertion adding functionality
- for ( int i = 0; i < 10; i++ )
- headerElement.addInsertionToTheRight();
- Assert.assertEquals(headerElement.numInsertionsToTheRight(), 10);
- }
-
- private void testHeaderIsEmpty(final HeaderElement headerElement) {
- Assert.assertFalse(headerElement.hasConsensusData(SlidingWindow.ConsensusType.POSITIVE_CONSENSUS));
- Assert.assertFalse(headerElement.hasConsensusData(SlidingWindow.ConsensusType.FILTERED));
- Assert.assertFalse(headerElement.hasInsertionToTheRight());
- Assert.assertTrue(headerElement.isEmpty());
- }
-
- private void testHeaderData(final HeaderElement headerElement, final HETest test) {
- Assert.assertEquals(headerElement.isVariantFromSoftClips(), test.isClip);
- Assert.assertFalse(headerElement.isEmpty());
- Assert.assertFalse(headerElement.hasInsertionToTheRight());
- Assert.assertEquals(headerElement.hasConsensusData(SlidingWindow.ConsensusType.POSITIVE_CONSENSUS), test.MQ >= minMappingQual);
- Assert.assertEquals(headerElement.hasConsensusData(SlidingWindow.ConsensusType.FILTERED), test.MQ < minMappingQual);
- Assert.assertEquals(headerElement.getBaseCounts(headerElement.hasConsensusData(SlidingWindow.ConsensusType.POSITIVE_CONSENSUS) ? SlidingWindow.ConsensusType.POSITIVE_CONSENSUS : SlidingWindow.ConsensusType.FILTERED).getRMS(), (double)test.MQ);
- Assert.assertFalse(headerElement.isVariantFromMismatches(0.05, 0.05));
- Assert.assertEquals(headerElement.isVariant(0.05, 0.05, 0.05), test.isClip);
- }
-
-
- private class AllelesTest {
- public final int[] counts;
- public final double pvalue;
-
- private AllelesTest(final int[] counts, final double pvalue) {
- this.counts = counts;
- this.pvalue = pvalue;
- }
- }
-
- @DataProvider(name = "alleles")
- public Object[][] createAllelesData() {
- List tests = new ArrayList<>();
-
- final int[] counts = new int[]{ 0, 5, 10, 15, 20 };
- final double [] pvalues = new double[]{ 0.0, 0.01, 0.05, 0.20, 1.0 };
-
- for ( final int countA : counts ) {
- for ( final int countC : counts ) {
- for ( final int countG : counts ) {
- for ( final int countT : counts ) {
- for ( final int countD : counts ) {
- for ( final double pvalue : pvalues ) {
- tests.add(new Object[]{new AllelesTest(new int[]{countA, countC, countG, countT, countD}, pvalue)});
- }
- }
- }
- }
- }
- }
-
- return tests.toArray(new Object[][]{});
- }
-
- @Test(dataProvider = "alleles", enabled = true)
- public void testAlleles(AllelesTest test) {
-
- HeaderElement headerElement = new HeaderElement(1000, 0);
- for ( int i = 0; i < test.counts.length; i++ ) {
- final BaseIndex base = BaseIndex.values()[i];
- for ( int j = 0; j < test.counts[i]; j++ )
- headerElement.addBase(base.b, byte20, byte10, byte10, byte20, minBaseQual, minMappingQual, false, false);
- }
-
- final int nAllelesSeen = headerElement.getNumberOfBaseAlleles(test.pvalue, test.pvalue);
- final int nAllelesExpected = calculateExpectedAlleles(test.counts, test.pvalue);
-
- Assert.assertEquals(nAllelesSeen, nAllelesExpected);
- }
-
- private static int calculateExpectedAlleles(final int[] counts, final double targetPvalue) {
- int total = 0;
- for ( final int count : counts ) {
- total += count;
- }
-
- int result = 0;
- for ( int index = 0; index < counts.length; index++ ) {
- final int count = counts[index];
- if ( count == 0 )
- continue;
-
- final boolean isSignificant;
- if ( count <= HeaderElement.MIN_COUNT_FOR_USING_PVALUE ) {
- isSignificant = MathUtils.binomialCumulativeProbability(total, 0, count) > targetPvalue;
- } else {
- isSignificant = (count >= targetPvalue * total);
- }
-
- if ( isSignificant ) {
- if ( index == BaseIndex.D.index )
- return -1;
- result++;
- }
- }
-
- return result;
- }
-}
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReadsIntegrationTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReadsIntegrationTest.java
deleted file mode 100644
index 067f36d58..000000000
--- a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReadsIntegrationTest.java
+++ /dev/null
@@ -1,347 +0,0 @@
-/*
-* By downloading the PROGRAM you agree to the following terms of use:
-*
-* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
-*
-* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
-*
-* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
-* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
-* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
-*
-* 1. DEFINITIONS
-* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
-*
-* 2. LICENSE
-* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
-* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
-* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
-* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
-*
-* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
-* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
-* Copyright 2012 Broad Institute, Inc.
-* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
-* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
-*
-* 4. INDEMNIFICATION
-* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
-*
-* 5. NO REPRESENTATIONS OR WARRANTIES
-* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
-* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
-*
-* 6. ASSIGNMENT
-* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
-*
-* 7. MISCELLANEOUS
-* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
-* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
-* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
-* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
-* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
-* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
-* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
-*/
-
-package org.broadinstitute.sting.gatk.walkers.compression.reducereads;
-
-import net.sf.samtools.SAMFileReader;
-import net.sf.samtools.SAMRecord;
-import org.broadinstitute.sting.WalkerTest;
-import org.broadinstitute.sting.utils.collections.Pair;
-import org.broadinstitute.sting.utils.exceptions.UserException;
-import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
-import org.broadinstitute.sting.utils.sam.GATKSamRecordFactory;
-import org.testng.Assert;
-import org.testng.annotations.Test;
-
-import java.io.File;
-import java.util.Arrays;
-import java.util.Collections;
-import java.util.List;
-
-public class ReduceReadsIntegrationTest extends WalkerTest {
- final static String REF = b37KGReference;
- final static String DBSNP = b37dbSNP132;
- final String BAM = validationDataLocation + "NA12878.HiSeq.b37.chr20.10_11mb.bam";
- final String DELETION_BAM = validationDataLocation + "filtered_deletion_for_reduce_reads.bam";
- final String STASH_BAM = validationDataLocation + "ReduceReadsStashBug.bam";
- final String STASH_L = " -L 14:73718184-73718284 -L 14:73718294-73718330 -L 14:73718360-73718556";
- final String DIVIDEBYZERO_BAM = validationDataLocation + "ReduceReadsDivideByZeroBug.bam";
- final String DIVIDEBYZERO_L = " -L " + validationDataLocation + "ReduceReadsDivideByZeroBug.intervals";
- final String L = " -L 20:10,100,000-10,120,000 ";
- final String COREDUCTION_BAM_A = validationDataLocation + "coreduction.test.A.bam";
- final String COREDUCTION_BAM_B = validationDataLocation + "coreduction.test.B.bam";
- final String COREDUCTION_L = " -L 1:1,853,860-1,854,354 -L 1:1,884,131-1,892,057";
- final String OFFCONTIG_BAM = privateTestDir + "readOffb37contigMT.bam";
- final String HIGH_COVERAGE_BAM = privateTestDir + "NA20313.highCoverageRegion.bam";
- final String HIGH_COVERAGE_L = " -L 1:1650830-1650870";
- final String BOTH_ENDS_OF_PAIR_IN_VARIANT_REGION_BAM = privateTestDir + "bothEndsOfPairInVariantRegion.bam";
- final String INSERTIONS_AT_EDGE_OF_CONSENSUS_BAM = privateTestDir + "rr-too-many-insertions.bam";
-
- final static String emptyFileMd5 = "d41d8cd98f00b204e9800998ecf8427e";
-
- protected Pair, List> executeTest(final String name, final WalkerTestSpec spec) {
- return executeTest(name, spec, emptyFileMd5);
- }
-
- protected Pair, List> executeTest(final String name, final WalkerTestSpec spec, final String qualsTestMD5) {
- final Pair, List> result = super.executeTest(name, spec);
-
- // perform some Reduce Reads specific testing now
- if ( result != null ) {
-
- // generate a new command-line based on the old one
- spec.disableImplicitArgs();
- final String[] originalArgs = spec.getArgsWithImplicitArgs().split(" ");
-
- final StringBuilder reducedInputs = new StringBuilder();
- for ( final File file : result.getFirst() ) {
- reducedInputs.append(" -I:reduced ");
- reducedInputs.append(file.getAbsolutePath());
- }
-
- // the coverage test is a less stricter version of the quals test so we can safely ignore it for now
- //final String coverageCommand = createCommandLine("AssessReducedCoverage", originalArgs);
- //super.executeTest(name + " : COVERAGE_TEST", new WalkerTestSpec(coverageCommand + reducedInputs.toString(), Arrays.asList(emptyFileMd5)));
-
- // run the quals test
- final String qualsCommand = createCommandLine("AssessReducedQuals", originalArgs);
- super.executeTest(name + " : QUALS_TEST", new WalkerTestSpec(qualsCommand + reducedInputs.toString(), Arrays.asList(qualsTestMD5)));
- }
-
- return result;
- }
-
- /*
- * Generate a new command-line based on the old one
- *
- * @param walkerName the new walker name to use
- * @param originalArgs the original arguments used for the test
- * @return the new command line
- */
- private String createCommandLine(final String walkerName, final String[] originalArgs) {
-
- final StringBuilder newArgs = new StringBuilder();
-
- for ( int i = 0; i < originalArgs.length; i++ ) {
- final String arg = originalArgs[i];
-
- if ( arg.equals("-T") ) {
- newArgs.append("-T ");
- newArgs.append(walkerName);
- } else if ( arg.startsWith("-I") ) {
- newArgs.append("-I:original ");
- newArgs.append(originalArgs[++i]);
- } else if ( arg.equals("-R") || arg.equals("-L") ) {
- newArgs.append(arg);
- newArgs.append(" ");
- newArgs.append(originalArgs[++i]);
- }
-
- // always add a trailing space
- newArgs.append(" ");
- }
-
- newArgs.append("-o %s");
-
- return newArgs.toString();
- }
-
- protected Pair, List> executeTestWithoutAdditionalRRTests(final String name, final WalkerTestSpec spec) {
- return super.executeTest(name, spec);
- }
-
- private void RRTest(final String testName, final String args, final String md5, final boolean useKnowns) {
- this.RRTest(testName, args, md5, useKnowns, emptyFileMd5);
- }
-
- private void RRTest(final String testName, final String args, final String md5, final boolean useKnowns, final String qualsTestMD5) {
- String base = String.format("-T ReduceReads -npt -R %s -I %s ", REF, BAM) + " -o %s" + (useKnowns ? " -known " + DBSNP : "") + " ";
- WalkerTestSpec spec = new WalkerTestSpec(base + args, Arrays.asList("bam"), Arrays.asList(md5));
- executeTest(testName, spec, qualsTestMD5);
- }
-
- @Test(enabled = true)
- public void testDefaultCompression() {
- RRTest("testDefaultCompression ", L, "0e503f7b79ace4c89d74f0943a0de1c0", false);
- }
-
- @Test(enabled = true)
- public void testDefaultCompressionWithKnowns() {
- RRTest("testDefaultCompressionWithKnowns ", L, "6db7ce2733d006f8bd61c42a40d23728", true);
- }
-
- private final String intervals = "-L 20:10,100,000-10,100,500 -L 20:10,200,000-10,200,500 -L 20:10,300,000-10,300,500 -L 20:10,400,000-10,500,000 -L 20:10,500,050-10,500,060 -L 20:10,600,000-10,600,015 -L 20:10,700,000-10,700,110";
-
- @Test(enabled = true)
- public void testMultipleIntervals() {
- RRTest("testMultipleIntervals ", intervals, "207f2c6d3db956e19412a45a231ca367", false, "043b2838c27d8f9580379b54c18ff40a");
- }
-
- @Test(enabled = true)
- public void testMultipleIntervalsWithKnowns() {
- RRTest("testMultipleIntervalsWithKnowns ", intervals, "f3b11a8a7673b301e27137936fafc6b6", true, "043b2838c27d8f9580379b54c18ff40a");
- }
-
- @Test(enabled = true)
- public void testHighCompression() {
- RRTest("testHighCompression ", " -cs 10 -min_pvalue 0.3 -minvar 0.3 -mindel 0.3 " + L, "dcc3716b3665aa1c2dbe6b22d6534aef", false);
- }
-
- @Test(enabled = true)
- public void testHighCompressionWithKnowns() {
- RRTest("testHighCompressionWithKnowns ", " -cs 10 -min_pvalue 0.3 -minvar 0.3 -mindel 0.3 " + L, "97ae655bf0e483ea227b1aac67ced024", true);
- }
-
- @Test(enabled = true)
- public void testLowCompression() {
- RRTest("testLowCompression ", " -cs 30 -min_pvalue 0.001 -minvar 0.01 -mindel 0.01 -minmap 5 -minqual 5 " + L, "a1377eb922e0b09a03a280b691b0b3ff", false);
- }
-
- @Test(enabled = true)
- public void testLowCompressionWithKnowns() {
- RRTest("testLowCompressionWithKnowns ", " -cs 30 -min_pvalue 0.001 -minvar 0.01 -mindel 0.01 -minmap 5 -minqual 5 " + L, "bd7c5b0b210694f364ca6a41f5b89870", true);
- }
-
- @Test(enabled = true)
- public void testBadPvalueInput() {
- final String cmd = String.format("-T ReduceReads -npt -R %s -I %s ", REF, BAM) + "-o %s -min_pvalue -0.01";
- WalkerTestSpec spec = new WalkerTestSpec(cmd, 1, UserException.BadArgumentValue.class);
- executeTest("testBadPvalueInput", spec);
- }
-
- @Test(enabled = true)
- public void testIndelCompression() {
- final String md5 = "9c9305eda5e4e7f22246ec8a4b242c97";
- RRTest("testIndelCompression ", " -cs 50 -L 20:10,100,500-10,100,600 ", md5, false);
- RRTest("testIndelCompressionWithKnowns ", " -cs 50 -L 20:10,100,500-10,100,600 ", md5, true);
- }
-
- @Test(enabled = true)
- public void testFilteredDeletionCompression() {
- String base = String.format("-T ReduceReads -npt -R %s -I %s ", REF, DELETION_BAM) + " -o %s ";
- executeTest("testFilteredDeletionCompression", new WalkerTestSpec(base, Arrays.asList("bam"), Arrays.asList("1bda512143be1016dfaca1f7020b6398")), "4f916da29d91852077f0a2fdbdd2c7f6");
- }
-
- @Test(enabled = true)
- public void testCoReduction() {
- String base = String.format("-T ReduceReads %s --cancer_mode -npt -R %s -I %s -I %s", COREDUCTION_L, REF, COREDUCTION_BAM_A, COREDUCTION_BAM_B) + " -o %s ";
- executeTest("testCoReduction", new WalkerTestSpec(base, Arrays.asList("bam"), Arrays.asList("58c2bae5a339af2ea3c22a46ce8faa68")));
- }
-
- @Test(enabled = true)
- public void testCoReductionWithKnowns() {
- String base = String.format("-T ReduceReads %s --cancer_mode -npt -R %s -I %s -I %s -known %s", COREDUCTION_L, REF, COREDUCTION_BAM_A, COREDUCTION_BAM_B, DBSNP) + " -o %s ";
- executeTest("testCoReductionWithKnowns", new WalkerTestSpec(base, Arrays.asList("bam"), Arrays.asList("5c251932b49d99a810581e3a6f762878")));
- }
-
- @Test(enabled = true)
- public void testInsertionsAtEdgeOfConsensus() {
- String base = String.format("-T ReduceReads -npt -R %s -I %s ", REF, INSERTIONS_AT_EDGE_OF_CONSENSUS_BAM) + " -o %s ";
- executeTest("testInsertionsAtEdgeOfConsensus", new WalkerTestSpec(base, Arrays.asList("bam"), Arrays.asList("c10653a8c21fb32b5cf580d3704b0edd")));
- }
-
- /**
- * Bug reported by Adam where a read that got clipped before actually belongs 2 intervals ahead
- * and a subsequent tail leaves only this read in the stash. The next read to come in is in fact
- * before (alignment start) than this read, so the TreeSet breaks with a Key out of Range error
- * that was freaking hard to catch.
- *
- * This bam is simplified to replicate the exact bug with the three provided intervals.
- */
- @Test(enabled = true)
- public void testAddingReadAfterTailingTheStash() {
- String base = String.format("-T ReduceReads %s -npt -R %s -I %s", STASH_L, REF, STASH_BAM) + " -o %s ";
- executeTest("testAddingReadAfterTailingTheStash", new WalkerTestSpec(base, Arrays.asList("bam"), Arrays.asList("fddbec29d0945afbbb34b42994614c15")), "3eab32c215ba68e75efd5ab7e9f7a2e7");
- }
-
- /**
- * Divide by zero bug reported by GdA and users in the forum. Happens when the downsampler goes over a region where all reads get
- * filtered out.
- */
- @Test(enabled = true)
- public void testDivideByZero() {
- String base = String.format("-T ReduceReads %s -npt -R %s -I %s", DIVIDEBYZERO_L, REF, DIVIDEBYZERO_BAM) + " -o %s ";
- // we expect to lose coverage due to the downsampling so don't run the systematic tests
- executeTestWithoutAdditionalRRTests("testDivideByZero", new WalkerTestSpec(base, Arrays.asList("bam"), Arrays.asList("7dfe2647992ce1154db340fc742d523a")));
- }
-
- /**
- * Bug happens when reads are soft-clipped off the contig (usually in the MT). This test guarantees no changes to the upstream code will
- * break the current hard-clipping routine that protects reduce reads from such reads.
- */
- @Test(enabled = true)
- public void testReadOffContig() {
- String base = String.format("-T ReduceReads -npt -R %s -I %s ", REF, OFFCONTIG_BAM) + " -o %s ";
- executeTest("testReadOffContig", new WalkerTestSpec(base, Arrays.asList("bam"), Arrays.asList("595e5812c37189930cae93e45765def4")));
- }
-
- /**
- * Confirm that if both ends of pair are in same variant region, compressed names of both ends of pair are the same.
- */
- @Test(enabled = true)
- public void testPairedReadsInVariantRegion() {
- String base = String.format("-T ReduceReads -npt -R %s -I %s ", hg19Reference, BOTH_ENDS_OF_PAIR_IN_VARIANT_REGION_BAM) +
- " -o %s --downsample_coverage 250 -dcov 50 ";
- executeTest("testPairedReadsInVariantRegion", new WalkerTestSpec(base, Arrays.asList("bam"), Arrays.asList("b005727119eee27995705959a637085e")), "2af063d1bd3c322b03405dbb3ecf59a9");
- }
-
- /**
- * Confirm that this bam does not fail when multi-sample mode is enabled. The provided example is tricky and used to cause
- * us to exception out in the code.
- */
- @Test(enabled = true)
- public void testMultiSampleDoesNotFailWithFlag() {
- String cmd = "-T ReduceReads --cancer_mode -npt -R " + b37KGReference + " -I " + privateTestDir + "rr_multisample.bam -o /dev/null";
- executeTestWithoutAdditionalRRTests("testMultiSampleDoesNotFailWithFlag", new WalkerTestSpec(cmd, 0, Collections.emptyList()));
- }
-
- /**
- * Confirm that this bam fails when multi-sample mode is not enabled
- */
- @Test(enabled = true)
- public void testMultiSampleFailsWithoutFlag() {
- String cmd = "-T ReduceReads -npt -R " + b37KGReference + " -I " + privateTestDir + "rr_multisample.bam -o /dev/null";
- executeTestWithoutAdditionalRRTests("testMultiSampleDoesNotFailWithFlag", new WalkerTestSpec(cmd, 0, UserException.BadInput.class));
- }
-
- /**
- * Confirm that compression is not capping coverage counts to max byte
- */
- @Test(enabled = true)
- public void testCompressionWorksForHighDepth() {
- final String base = String.format("-T ReduceReads -npt -R %s -I %s %s", b37KGReference, HIGH_COVERAGE_BAM, HIGH_COVERAGE_L) + " -o %s";
- final File outputBam = executeTestWithoutAdditionalRRTests("testCompressionWorksForHighDepth",
- new WalkerTestSpec(base, 1, Arrays.asList(""))).first.get(0); // No MD5s; we only want to check the coverage
-
- boolean sawHighCoveragePosition = false;
- final SAMFileReader reader = new SAMFileReader(outputBam);
- reader.setSAMRecordFactory(new GATKSamRecordFactory());
-
- for ( final SAMRecord rawRead : reader ) {
- final GATKSAMRecord read = (GATKSAMRecord)rawRead;
- read.setAttribute(GATKSAMRecord.REDUCED_READ_CONSENSUS_TAG, rawRead.getByteArrayAttribute(GATKSAMRecord.REDUCED_READ_CONSENSUS_TAG));
-
- if ( ! read.isReducedRead() )
- continue;
-
- final int[] decodedCounts = read.getReducedReadCounts();
- for ( final int count : decodedCounts ) {
- if ( count > Byte.MAX_VALUE ) {
- sawHighCoveragePosition = true;
- break;
- }
- }
-
- if ( sawHighCoveragePosition )
- break;
- }
-
- reader.close();
-
- Assert.assertTrue(sawHighCoveragePosition, "No positions were found with coverage over max byte (127); the coverage is incorrectly being capped somewhere!");
- }
-}
-
diff --git a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReadsUnitTest.java b/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReadsUnitTest.java
deleted file mode 100644
index 6032affa7..000000000
--- a/protected/gatk-protected/src/test/java/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReadsUnitTest.java
+++ /dev/null
@@ -1,214 +0,0 @@
-/*
-* By downloading the PROGRAM you agree to the following terms of use:
-*
-* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
-*
-* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
-*
-* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
-* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
-* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
-*
-* 1. DEFINITIONS
-* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
-*
-* 2. LICENSE
-* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
-* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
-* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
-* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
-*
-* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
-* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
-* Copyright 2012 Broad Institute, Inc.
-* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
-* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
-*
-* 4. INDEMNIFICATION
-* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
-*
-* 5. NO REPRESENTATIONS OR WARRANTIES
-* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
-* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
-*
-* 6. ASSIGNMENT
-* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
-*
-* 7. MISCELLANEOUS
-* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
-* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
-* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
-* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
-* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
-* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
-* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
-*/
-
-package org.broadinstitute.sting.gatk.walkers.compression.reducereads;
-
-import it.unimi.dsi.fastutil.objects.*;
-import net.sf.samtools.SAMFileHeader;
-import org.broad.tribble.Feature;
-import org.broadinstitute.sting.BaseTest;
-import org.broadinstitute.sting.commandline.RodBinding;
-import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
-import org.broadinstitute.sting.gatk.refdata.RODRecordListImpl;
-import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
-import org.broadinstitute.sting.gatk.refdata.utils.GATKFeature;
-import org.broadinstitute.sting.gatk.refdata.utils.RODRecordList;
-import org.broadinstitute.sting.utils.GenomeLoc;
-import org.broadinstitute.sting.utils.GenomeLocParser;
-import org.broadinstitute.sting.utils.sam.ArtificialSAMUtils;
-import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
-import org.broadinstitute.variant.variantcontext.Allele;
-import org.broadinstitute.variant.variantcontext.VariantContext;
-import org.broadinstitute.variant.variantcontext.VariantContextBuilder;
-import org.testng.Assert;
-import org.testng.annotations.BeforeClass;
-import org.testng.annotations.DataProvider;
-import org.testng.annotations.Test;
-
-import java.util.ArrayList;
-import java.util.Arrays;
-import java.util.List;
-import java.util.Random;
-
-
-public class ReduceReadsUnitTest extends BaseTest {
-
- Random random = new Random(987743);
- Object2LongOpenHashMap hash = new Object2LongOpenHashMap();
- long nextNumber = 0L;
-
- /**
- * Combinatorial unit test data provider example.
- *
- * Creates data for testMyData test function, containing two arguments, start and size at each value
- *
- * @return Object[][] for testng DataProvider
- */
- @DataProvider(name = "ReadNameProvider")
- public Object[][] readNameProvider() {
- final int readNameLength = 4;
- final int nReads = 100000;
- final int charVariety = 20;
- ObjectArrayList tests = new ObjectArrayList();
- ObjectOpenHashSet