From 2e237a12e9abc7b93fd2f57e94ada7d76b9b8ebf Mon Sep 17 00:00:00 2001 From: chartl Date: Sun, 13 Sep 2009 01:00:04 +0000 Subject: [PATCH] This commit has a bunch to do with cleaning up the CoverageAndPowerWalker code: implementing some new printing options, but mostly altering the code so it's much more readable and understandable, and much less hacky-looking. ADDED: @Quad: This is just like Pair, except with four fields. In the original CoverageAndPowerWalker I often used a pair of pairs to hold things, which made the code nigh unreadable. @SQuad: An extension of Quad for when you want to store objects of the same type. Let's you simply declare new SQuad rather than new Quad @ReadOffsetQuad: An extension of Quad specifically for holding two lists of reads and two lists of offsets Supports construction from AlignmentContexts and conversion to AlignmentContexts (given a GenomeLoc). There are methods that make it very clear what the code is doing (getSecondRead() rather than the cryptic getThird() ) @PowerAndCoverageWalker: The new version of CoverageAndPowerWalker. If the tests all go well, then I'll remove the old version. New to this version is the ability to give an output file directly to the walker, so that locus information prints to the file, while the final reduce prints to standard out. Bootstrap iterations are now a command line argument rather than a final int; and users can instruct the walker to print out the coverage/power statistics for both the original reads, and those reads whose quality score exceeds a user-defined threshold. CHANGES: @PoolUtils: Altered methods to accept as argumetns, and return, Quad objects. Added a random partition method for bootstrapping. @CoverageAndPowerWalker: Altered methods to work with the new PoolUtils methods. git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1602 348d0f76-0448-11de-a6fe-93d51630548a --- .../poolseq/CoverageAndPowerWalker.java | 7 +- .../poolseq/PowerAndCoverageWalker.java | 261 ++++++++++++++++++ .../sting/playground/utils/PoolUtils.java | 37 ++- .../sting/playground/utils/Quad.java | 90 ++++++ .../playground/utils/ReadOffsetQuad.java | 100 +++++++ .../sting/playground/utils/SQuad.java | 20 ++ 6 files changed, 510 insertions(+), 5 deletions(-) create mode 100644 java/src/org/broadinstitute/sting/playground/gatk/walkers/poolseq/PowerAndCoverageWalker.java create mode 100644 java/src/org/broadinstitute/sting/playground/utils/Quad.java create mode 100644 java/src/org/broadinstitute/sting/playground/utils/ReadOffsetQuad.java create mode 100644 java/src/org/broadinstitute/sting/playground/utils/SQuad.java diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/poolseq/CoverageAndPowerWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/poolseq/CoverageAndPowerWalker.java index 13f810804..9d35f0297 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/poolseq/CoverageAndPowerWalker.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/poolseq/CoverageAndPowerWalker.java @@ -12,6 +12,7 @@ import org.broadinstitute.sting.gatk.walkers.By; import org.broadinstitute.sting.gatk.contexts.ReferenceContext; import org.broadinstitute.sting.gatk.contexts.AlignmentContext; import org.broadinstitute.sting.playground.utils.PoolUtils; +import org.broadinstitute.sting.playground.utils.ReadOffsetQuad; import net.sf.samtools.SAMRecord; import java.util.*; @@ -71,7 +72,8 @@ public class CoverageAndPowerWalker extends LocusWalker, Pair,List> readsFilteredByQuality = filterByQuality(context.getReads(),context.getOffsets(), this.getMinQualityScore()); filteredContext = new AlignmentContext(context.getLocation(),readsFilteredByQuality.getFirst(),readsFilteredByQuality.getSecond()); } - Pair,List>,Pair,List>> readsByDirection = PoolUtils.splitReadsByReadDirection(filteredContext.getReads(),filteredContext.getOffsets()); + ReadOffsetQuad splitReads = PoolUtils.splitReadsByReadDirection(filteredContext.getReads(), filteredContext.getOffsets()); + Pair,List>,Pair,List>> readsByDirection = new Pair(new Pair(splitReads.getFirstReads(),splitReads.getSecondReads()),new Pair(splitReads.getFirstOffsets(),splitReads.getSecondOffsets())); if ( ! suppress_printing && ! outputUnthreshCvg ) { Pair powers = calculatePower(readsByDirection, useBootstrap, filteredContext); out.printf("%s %d %d %d %d %d %d %f %f %f%n", filteredContext.getLocation(), readsByDirection.getFirst().getFirst().size(), readsByDirection.getFirst().getSecond().size(), @@ -79,7 +81,8 @@ public class CoverageAndPowerWalker extends LocusWalker, powers.getFirst()[0], powers.getFirst()[1], powers.getFirst()[2]); } else if (! suppress_printing && outputUnthreshCvg) { Pair powers = calculatePower(readsByDirection, useBootstrap, filteredContext); - Pair,List>,Pair,List>> ufReadsByDirection = PoolUtils.splitReadsByReadDirection(context.getReads(),context.getOffsets()); + ReadOffsetQuad readsByDir = PoolUtils.splitReadsByReadDirection(context.getReads(),context.getOffsets()); + Pair,List>,Pair,List>> ufReadsByDirection = new Pair(new Pair(readsByDir.getFirstReads(), readsByDir.getSecondReads()), new Pair(readsByDir.getFirstOffsets(), readsByDir.getSecondOffsets())); out.printf("%s %d %d %d %d %d %d %d %d %d %f %f %f%n", filteredContext.getLocation(), ufReadsByDirection.getFirst().getFirst().size(), ufReadsByDirection.getFirst().getSecond().size(), context.getReads().size(), readsByDirection.getFirst().getFirst().size(), readsByDirection.getFirst().getSecond().size(), diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/poolseq/PowerAndCoverageWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/poolseq/PowerAndCoverageWalker.java new file mode 100644 index 000000000..67b5f8b2a --- /dev/null +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/poolseq/PowerAndCoverageWalker.java @@ -0,0 +1,261 @@ +package org.broadinstitute.sting.playground.gatk.walkers.poolseq; + +import org.broadinstitute.sting.playground.utils.SQuad; +import org.broadinstitute.sting.playground.utils.ReadOffsetQuad; +import org.broadinstitute.sting.playground.utils.PoolUtils; +import org.broadinstitute.sting.gatk.walkers.LocusWalker; +import org.broadinstitute.sting.gatk.walkers.By; +import org.broadinstitute.sting.gatk.walkers.DataSource; +import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker; +import org.broadinstitute.sting.gatk.contexts.ReferenceContext; +import org.broadinstitute.sting.gatk.contexts.AlignmentContext; +import org.broadinstitute.sting.utils.cmdLine.Argument; +import org.broadinstitute.sting.utils.*; +import net.sf.samtools.SAMRecord; + +import java.io.PrintStream; +import java.io.FileNotFoundException; +import java.util.List; + +/** + * Created by IntelliJ IDEA. + * User: Ghost + * Date: Sep 12, 2009 + * Time: 12:20:42 PM + * To change this template use File | Settings | File Templates. + */ +@By(DataSource.REFERENCE) +public class PowerAndCoverageWalker extends LocusWalker, SQuad> { + + @Argument(fullName="suppressLocusPrinting", doc="Walker does not output locus data to any file. Only total coverage information is given at the end", required=false) + public boolean suppressPrinting = false; + + @Argument(fullName="outputFile", shortName="of", doc="Locus information is printed to this file instead of standard out", required = false) + public String outputFile = null; + + @Argument(fullName="bootStrapIterations", shortName="bs", doc="Use a bootstrap method with this many iterations", required=false) + public int bootstrapIterations = 0; + + @Argument(fullName="lodThreshold", shortName="lod", doc="Threshold for log likelihood ratio to be called a SNP. Defaults to 3.0", required = false) + public double lodThresh = 3.0; + + @Argument(fullName="minimumQScore", shortName="qm", doc="Use bases whose phred (Q) score meets or exceeds this number. Defaults to 0", required = false) + public byte minQ = 0; + + @Argument(fullName="aboveQScoreOutputOnly", doc="Output only the information for reads that exceed the q-score threshold", required=false) + public boolean aboveQScoreOutputOnly = false; + + @Argument(fullName="outputRawStatistics", shortName="rp", doc="Walker outputs the median quality score and power for the un-thresholded reads as well as those that are thresholded", required=false) + public boolean outputRawStatistics = false; + + @Argument(fullName="poolSize", shortName="pf", doc="Number of individuals in the pool", required = true) + public int numIndividuals = 0; + + protected PrintStream outputWriter = null; + + public void initialize() { + if(numIndividuals <= 0) { + throw new StingException("Pool size must be greater than 1. You input "+numIndividuals); + } + if( outputFile != null ) { + try { + outputWriter = new PrintStream(outputFile); + } catch (FileNotFoundException e) { + String errMsg = "The filepath input as an argument to -of, "+outputFile+" was a directory or could not be found or created."; + throw new StingException(errMsg, e); + } + outputWriter.printf("%s%n", generateHeaderString()); + } else if( ! suppressPrinting ) { + out.printf("%s%n", generateHeaderString()); + } + } + + public SQuad reduceInit() { + return new SQuad(0l,0l,0l,0l); + } + + /* REDUCE + * Add up filtered and unfiltered coverage in the forward + * and reverse directions. + */ + public SQuad reduce( SQuad newCvg, SQuad oldCvg) { + return new SQuad(newCvg.getFirst() + oldCvg.getFirst(), newCvg.getSecond() + oldCvg.getSecond(), + newCvg.getThird() + oldCvg.getThird(), newCvg.getFourth() + oldCvg.getFourth()); + } + /* MAP + * Get the coverage in each direction above and below threshold and calculate the estimated power + * if printing is not suppressed print this information to the file + * pass coverage information on to Reduce + */ + public SQuad map( RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) { + ReadOffsetQuad readsByDirection = splitReadsByDirection(context.getReads(),context.getOffsets()); + ReadOffsetQuad readsByDirectionAboveThresh; + + if(minQ > 0) { + readsByDirectionAboveThresh = thresholdReadsByQuality(readsByDirection); + } else { + readsByDirectionAboveThresh = readsByDirection; + } + + + SQuad medianQScoresAboveThresh = getMedianQScores(readsByDirectionAboveThresh); + SQuad powerAboveThresh = calculatePower(readsByDirectionAboveThresh, medianQScoresAboveThresh); + SQuad medianQScoresUnthresh; + SQuad powerUnthresh; + + if(! outputRawStatistics) { + medianQScoresUnthresh = null; + powerUnthresh = null; + } else { + medianQScoresUnthresh = getMedianQScores(readsByDirection); + powerUnthresh= calculatePower(readsByDirection, medianQScoresUnthresh); + } + + if ( ! suppressPrinting && outputWriter == null) { + out.printf("%s%n",outputString(readsByDirection,readsByDirectionAboveThresh,powerAboveThresh, powerUnthresh, + medianQScoresAboveThresh, medianQScoresUnthresh, context.getLocation())); + } else if (! suppressPrinting && outputWriter != null) { + outputWriter.printf("%s%n",outputString(readsByDirection,readsByDirectionAboveThresh,powerAboveThresh, powerUnthresh, + medianQScoresAboveThresh, medianQScoresUnthresh, context.getLocation())); + } + + return new SQuad(readsByDirection.numReadsFirst(), readsByDirection.numReadsSecond(), + readsByDirectionAboveThresh.numReadsFirst(), readsByDirectionAboveThresh.numReadsSecond()); + + } + + // private helper methods + + /* + * Make different headers for the output file depending on printing-related command line arguments + */ + private String generateHeaderString() { + String header; + if (minQ <= 0) { + header = "Chrom:Pos CoverageFwd CoverageRev CoverageCom MedianQFwd MedianQRev MedianQCom PowerFwd PowerRev PowerCom"; + } else if (! aboveQScoreOutputOnly && ! outputRawStatistics){ + String m = Byte.toString(minQ); + header = "Chrom:Pos CoverageFwd CoverageRev CoverageCom CoverageAboveQ"+m+"Fwd CoverageAboveQ"+m+"Rev CoverageAboveQ"+m+"Com MedianQFwd MedianQRev MedianQCom PowerFwd PowerRev PowerCom"; + } else if (! outputRawStatistics){ + header = "Chrom:Pos CoverageFwd CoverageRev CoverageCom MedianQFwd MedianQRev MedianQCom PowerFwd PowerRev PowerCom"; + } else { + String m = Byte.toString(minQ); + header = "Chrom:Pos CoverageFwd CoverageRev CoverageCom CoverageAboveQ"+m+"Fwd CoverageAboveQ"+m+"Rev "+ + "CoverageAboveQ"+m+"Com RawMedianQFwd RawMedianQRev RawMedianQCom MedianQAboveQ"+m+"Fwd "+ + "MedianQAboveQ"+m+"Rev MedianQAboveQ"+m+"Com RawPowerFwd RawPowerRev RawPowerCom "+ + "PowerAboveQ"+m+"Fwd PowerAboveQ"+m+"Rev PowerAboveQ"+m+"Com"; + } + + return header; + } + /* + * Make different output strings depending on printing-related command-line options + */ + private String outputString(ReadOffsetQuad reads, ReadOffsetQuad threshReads, SQuad pow, SQuad rawPow, SQuad medQ, SQuad rawQ, GenomeLoc loc) { + String outputString; + if(minQ <= 0) { + outputString = String.format("%s %d %d %d %d %d %d %f %f %f", loc.toString(), reads.numReadsFirst(), reads.numReadsSecond(), reads.numReads(), + medQ.getFirst(), medQ.getSecond(), medQ.getThird(), pow.getFirst(), pow.getSecond(), pow.getThird()); + } else if (! aboveQScoreOutputOnly && ! outputRawStatistics ){ + outputString = String.format("%s %d %d %d %d %d %d %d %d %d %f %f %f", loc.toString(), reads.numReadsFirst(), reads.numReadsSecond(), reads.numReads(), + threshReads.numReadsFirst(), threshReads.numReadsSecond(), threshReads.numReads(), medQ.getFirst(), + medQ.getSecond(), medQ.getThird(), pow.getFirst(), pow.getSecond(), pow.getThird()); + } else if (! outputRawStatistics ){ + outputString = String.format("%s %d %d %d %d %d %d %f %f %f", loc.toString(), threshReads.numReadsFirst(), threshReads.numReadsSecond(), + threshReads.numReads(), medQ.getFirst(), medQ.getSecond(), medQ.getThird(), pow.getFirst(), + pow.getSecond(), pow.getThird()); + } else { + outputString = String.format("%s %d %d %d %d %d %d %d %d %d %d %d %d %f %f %f %f %f %f", loc.toString(), reads.numReadsFirst(), + reads.numReadsSecond(), reads.numReads(), threshReads.numReadsFirst(), threshReads.numReadsSecond(), + threshReads.numReads(), rawQ.getFirst(), rawQ.getSecond(), rawQ.getThird(), medQ.getFirst(), + medQ.getSecond(), medQ.getThird(), rawPow.getFirst(), rawPow.getSecond(), rawPow.getThird(), + pow.getFirst(), pow.getSecond(), pow.getThird()); + } + return outputString; + } + + private ReadOffsetQuad splitReadsByDirection(List reads, List offsets) { + return PoolUtils.splitReadsByReadDirection(reads,offsets); + } + + private ReadOffsetQuad thresholdReadsByQuality(ReadOffsetQuad readQuad) { + return new ReadOffsetQuad(PoolUtils.thresholdReadsByQuality(readQuad.getFirstReadOffsetPair(), minQ), + PoolUtils.thresholdReadsByQuality(readQuad.getSecondReadOffsetPair(), minQ)); + } + + private SQuad getMedianQScores(ReadOffsetQuad reads) { + return new SQuad(ListUtils.getQScoreMedian(reads.getFirstReads(), reads.getFirstOffsets()), + ListUtils.getQScoreMedian(reads.getSecondReads(), reads.getSecondOffsets()), + ListUtils.getQScoreMedian(reads.getReadsCombined(), reads.getOffsetsCombined()), + (byte) 0); + } + + private SQuad calculatePower(ReadOffsetQuad readQuad, SQuad medianQs) { + SQuad power; + if( bootstrapIterations > 0 ) { + power = bootstrapPowerByDirection(readQuad); + } else { + power = theoreticalPowerByDirection(readQuad,medianQs); + } + return power; + } + + private SQuad theoreticalPowerByDirection(ReadOffsetQuad readQuad, SQuad medianQs) { + return new SQuad( theoreticalPower(readQuad.numReadsFirst(), medianQs.getFirst()), + theoreticalPower(readQuad.numReadsSecond(), medianQs.getSecond()), + theoreticalPower(readQuad.numReads(), medianQs.getThird()), + 0.0 ); + } + + private SQuad bootstrapPowerByDirection(ReadOffsetQuad readQuad) { + return new SQuad ( bootstrapPower(readQuad.getFirstReads(), readQuad.getFirstOffsets()), + bootstrapPower(readQuad.getSecondReads(), readQuad.getSecondOffsets()), + bootstrapPower(readQuad.getReadsCombined(), readQuad.getOffsetsCombined()), + 0.0 ); + } + + private double theoreticalPower(int depth, byte q) { + double power; + if( depth <= 0 ) { + power = 0.0; + } else { + double p_error = QualityUtils.qualToErrorProb(q); + double snpProp = getSNPProportion(); + double kterm_num = Math.log10( snpProp * (1 - p_error) + (1 - snpProp) * (p_error/3) ); + double kterm_denom = Math.log10( p_error/3 ); + double dkterm_num = Math.log10( snpProp * (p_error/3) + (1 - snpProp) * (1 - p_error) ); + double dkterm_denom = Math.log10( 1 - p_error); + + int kaccept = (int) Math.ceil( ( lodThresh - ( (double) depth ) * ( dkterm_num - dkterm_denom ) ) / + ( kterm_num - kterm_denom- dkterm_num + dkterm_denom ) ); + + // we will reject the null hypothesis if we see kaccept or more SNPs, the power is the probability that this occurs + // we can optimize this by checking to see which sum is smaller + + if ( depth - kaccept < kaccept ) {// kaccept > depth/2 - calculate power as P[hits between kaccept and depth] + power = MathUtils.cumBinomialProbLog(kaccept, depth, depth, snpProp); + } else { // kaccept < depth/2 - calculate power as 1-P[hits between 0 and kaccept] + power = 1-MathUtils.cumBinomialProbLog(0,kaccept,depth,snpProp); + } + } + + return power; + } + + private double bootstrapPower(List reads, List offsets) { + double power = 0.0; + double snpProp = getSNPProportion(); + for ( int boot = 0; boot < bootstrapIterations; boot ++) { + ReadOffsetQuad partition = PoolUtils.coinTossPartition(reads, offsets, getSNPProportion()); + if ( PoolUtils.calculateLogLikelihoodOfSample(partition, snpProp) >= lodThresh ) { + power += 1.0/bootstrapIterations; + } + } + return power; + } + + private double getSNPProportion() { + return 1/(2.0*numIndividuals); + } + +} diff --git a/java/src/org/broadinstitute/sting/playground/utils/PoolUtils.java b/java/src/org/broadinstitute/sting/playground/utils/PoolUtils.java index 2614e18c8..baa8a2fb1 100755 --- a/java/src/org/broadinstitute/sting/playground/utils/PoolUtils.java +++ b/java/src/org/broadinstitute/sting/playground/utils/PoolUtils.java @@ -26,7 +26,7 @@ public class PoolUtils { public static final int BASE_T_OFFSET = 3; public static final int BASE_INDEXED_ARRAY_SIZE = 4; - public static Pair, List>, Pair, List>> splitReadsByReadDirection(List reads, List offsets) { + public static ReadOffsetQuad splitReadsByReadDirection(List reads, List offsets) { ArrayList forwardReads; ArrayList reverseReads; ArrayList forwardOffsets; @@ -54,7 +54,7 @@ public class PoolUtils { } } - return new Pair(new Pair(forwardReads, reverseReads), new Pair(forwardOffsets, reverseOffsets)); + return new ReadOffsetQuad(forwardReads,forwardOffsets,reverseReads,reverseOffsets); } public static Pair[], List[]> splitReadsByBase(List reads, List offsets) { @@ -117,7 +117,11 @@ public class PoolUtils { } } - return new Pair(threshReads, threshOffsets); + return new Pair,List>(threshReads, threshOffsets); + } + + public static Pair,List> thresholdReadsByQuality(Pair,List> readPair, byte qThresh) { + return thresholdReadsByQuality(readPair.getFirst(),readPair.getSecond(),qThresh); } public static int getBaseOffset(char base) { @@ -165,6 +169,14 @@ public class PoolUtils { return qualities; } + public static double calculateLogLikelihoodOfSample(ReadOffsetQuad readsBySupport, double alleleFreq) { + List qSNP = getReadBaseQualities(readsBySupport.getFirstReads(), readsBySupport.getFirstOffsets()); + List qRef = getReadBaseQualities(readsBySupport.getSecondReads(), readsBySupport.getSecondOffsets()); + Pair logsumSNP = qListToSumLogProbabilities(true, qSNP, 1/alleleFreq); + Pair logsumRef = qListToSumLogProbabilities(false, qRef, 1/alleleFreq); + return 0.0 - logsumSNP.getFirst() - logsumRef.getFirst() + logsumSNP.getSecond() + logsumRef.getSecond(); + } + public static double calculateLogLikelihoodOfSample(Pair,List>,Pair,List>> snpReadsRefReads, int nIndivids) { List qListSnps = getReadBaseQualities(snpReadsRefReads.getFirst().getFirst(),snpReadsRefReads.getSecond().getFirst()); List qListRefs = getReadBaseQualities(snpReadsRefReads.getFirst().getSecond(),snpReadsRefReads.getSecond().getSecond()); @@ -194,4 +206,23 @@ public class PoolUtils { return new Pair(logProbObserveXAndSNPTrue,logProbObserveXAndRefTrue); } + public static ReadOffsetQuad coinTossPartition(List reads, List offsets, double alleleFreq) { + ArrayList snpReads = new ArrayList(); + ArrayList snpOffsets = new ArrayList(); + ArrayList refReads = new ArrayList(); + ArrayList refOffsets = new ArrayList(); + + for( int i = 0; i < reads.size(); i ++ ) { + if ( Math.random() < alleleFreq ) { + snpReads.add(reads.get(i)); + snpOffsets.add(offsets.get(i)); + } else { + refReads.add(reads.get(i)); + refOffsets.add(offsets.get(i)); + } + } + + return new ReadOffsetQuad(snpReads,snpOffsets,refReads,refOffsets); + } + } diff --git a/java/src/org/broadinstitute/sting/playground/utils/Quad.java b/java/src/org/broadinstitute/sting/playground/utils/Quad.java new file mode 100644 index 000000000..458090ed8 --- /dev/null +++ b/java/src/org/broadinstitute/sting/playground/utils/Quad.java @@ -0,0 +1,90 @@ +package org.broadinstitute.sting.playground.utils; + +import org.broadinstitute.sting.utils.Pair; + +/** + * Created by IntelliJ IDEA. + * User: chartl + * Date: Sep 11, 2009 + * Time: 5:12:29 PM + * To change this template use File | Settings | File Templates. + */ +public class Quad { + public W first; + public X second; + public Y third; + public Z fourth; + + public Quad() { + first = null; + second = null; + third = null; + fourth = null; + } + + public Quad(W w, X x, Y y, Z z) { + first = w; + second = x; + third = y; + fourth = z; + } + + public Quad(Pair a, Pair b) { + first = a.getFirst(); + second = a.getSecond(); + third = b.getFirst(); + fourth = b.getSecond(); + } + + public boolean equals(Object o) { + if(o == null) { + return false; + } else if (! (o instanceof Quad) ) { + return false; + } + + Quad other = (Quad) o; + + return ( equalToNotNull(this.first,other.first) && equalToNotNull(this.second,other.second) + && equalToNotNull(this.third,other.third) && equalToNotNull(this.fourth,other.fourth)); + } + + public int hashCode() { + return getHash(first) ^ getHash(second) ^ getHash(third) ^ getHash(fourth); + } + + public W getFirst() { return first; } + public X getSecond() { return second; } + public Y getThird() { return third; } + public Z getFourth() { return fourth; } + + public Pair getFirstPair() { return new Pair(first,second); } + public Pair getSecondPair() { return new Pair(third,fourth); } + + public void setFirst(Object o) { first = (W) o; } + public void setSecond(Object o) { second = (X) o; } + public void setThird(Object o) { third = (Y) o; } + public void setFourth(Object o) { fourth = (Z) o; } + + private int getHash(Object o) { + int hash = 0; + if(o != null) { + hash = o.hashCode(); + } + return hash; + } + + private boolean equalToNotNull(Object a, Object b) { + boolean areEqual = false; + if ( a != null && b != null ) { + areEqual = a.equals(b); + } else if (a == null && b == null ) { + areEqual = true; + // todo -- make sure we don't want to check for instanceOf here... + // todo -- maybe this statement should be eliminated + } + + return areEqual; + } + +} diff --git a/java/src/org/broadinstitute/sting/playground/utils/ReadOffsetQuad.java b/java/src/org/broadinstitute/sting/playground/utils/ReadOffsetQuad.java new file mode 100644 index 000000000..cb003e7cf --- /dev/null +++ b/java/src/org/broadinstitute/sting/playground/utils/ReadOffsetQuad.java @@ -0,0 +1,100 @@ +package org.broadinstitute.sting.playground.utils; + +import org.broadinstitute.sting.gatk.contexts.AlignmentContext; +import org.broadinstitute.sting.utils.GenomeLoc; +import org.broadinstitute.sting.utils.Pair; +import net.sf.samtools.SAMRecord; +import java.util.List; +import java.util.ArrayList; + +/** + * Created by IntelliJ IDEA. + * User: chartl + * Date: Sep 12, 2009 + * Time: 1:21:38 AM + * To change this template use File | Settings | File Templates. + */ +public class ReadOffsetQuad extends Quad,List,List, List> { + /* + * ReadOffsetQuad separates the user from specifying the types of objects required to + * store two sets of read/offset pairs in a quad. + * + * Implements methods that return read/offset pairs as AlignmentContexts + * and allows ReadOffsetQuad to be constructed from two AlignmentContexts + * + */ + + // constructor that IntelliJ wants + public ReadOffsetQuad(List a, List b, List c, List d) { + super(a,b,c,d); + } + + // another constructor that IntelliJ wants + public ReadOffsetQuad(Pair,List> a, Pair,List> b) { + super(a,b); + } + + public ReadOffsetQuad(AlignmentContext a, AlignmentContext b) { + first = a.getReads(); + second = a.getOffsets(); + third = b.getReads(); + fourth = b.getOffsets(); + } + + public int numReads() { + return first.size() + third.size(); + } + + public int numReadsFirst() { + return first.size(); + } + + public int numReadsSecond() { + return third.size(); + } + + public List getFirstReads() { + return this.first; + } + + public List getSecondReads() { + return this.third; + } + + public List getReadsCombined() { + ArrayList combined = new ArrayList(first); + combined.addAll(third); + return combined; + } + + public List getOffsetsCombined() { + ArrayList combined = new ArrayList(second); + combined.addAll(fourth); + return combined; + } + + public List getFirstOffsets() { + return this.second; + } + + public List getSecondOffsets() { + return this.fourth; + } + + public Pair,List> getFirstReadOffsetPair() { + return new Pair,List>(first, second); + } + + public Pair,List> getSecondReadOffsetPair() { + return new Pair,List>(third,fourth); + } + + public AlignmentContext getFirstPairAsAlignmentContext(GenomeLoc loc) { + return new AlignmentContext(loc, first, second); + } + + public AlignmentContext getSecondPairAsAlignmentContext(GenomeLoc loc) { + return new AlignmentContext(loc, third, fourth); + } + +} diff --git a/java/src/org/broadinstitute/sting/playground/utils/SQuad.java b/java/src/org/broadinstitute/sting/playground/utils/SQuad.java new file mode 100644 index 000000000..110955ade --- /dev/null +++ b/java/src/org/broadinstitute/sting/playground/utils/SQuad.java @@ -0,0 +1,20 @@ +package org.broadinstitute.sting.playground.utils; + +/** + * Created by IntelliJ IDEA. + * User: Ghost + * Date: Sep 12, 2009 + * Time: 1:04:40 PM + * To change this template use File | Settings | File Templates. + */ +public class SQuad extends Quad { + /* SQuad - single object type Quad + * or "Simple Quad". Makes it less code-intensive + * for user to use quads to hold objects. + */ + + public SQuad(X a, X b, X c, X d) { super(a,b,c,d); } + // don't know why this method even neds to be written + // but IntelliJ wants it to be there + +}