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<X> rather than new Quad<X,X,X,X>
@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
This commit is contained in:
chartl
parent
6c7a300664
commit
2e237a12e9
|
|
@ -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<Integer, Integer>,
|
|||
Pair<List<SAMRecord>,List<Integer>> readsFilteredByQuality = filterByQuality(context.getReads(),context.getOffsets(), this.getMinQualityScore());
|
||||
filteredContext = new AlignmentContext(context.getLocation(),readsFilteredByQuality.getFirst(),readsFilteredByQuality.getSecond());
|
||||
}
|
||||
Pair<Pair<List<SAMRecord>,List<SAMRecord>>,Pair<List<Integer>,List<Integer>>> readsByDirection = PoolUtils.splitReadsByReadDirection(filteredContext.getReads(),filteredContext.getOffsets());
|
||||
ReadOffsetQuad splitReads = PoolUtils.splitReadsByReadDirection(filteredContext.getReads(), filteredContext.getOffsets());
|
||||
Pair<Pair<List<SAMRecord>,List<SAMRecord>>,Pair<List<Integer>,List<Integer>>> readsByDirection = new Pair(new Pair(splitReads.getFirstReads(),splitReads.getSecondReads()),new Pair(splitReads.getFirstOffsets(),splitReads.getSecondOffsets()));
|
||||
if ( ! suppress_printing && ! outputUnthreshCvg ) {
|
||||
Pair<double[],byte[]> 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<Pair<Integer, Integer>,
|
|||
powers.getFirst()[0], powers.getFirst()[1], powers.getFirst()[2]);
|
||||
} else if (! suppress_printing && outputUnthreshCvg) {
|
||||
Pair<double[],byte[]> powers = calculatePower(readsByDirection, useBootstrap, filteredContext);
|
||||
Pair<Pair<List<SAMRecord>,List<SAMRecord>>,Pair<List<Integer>,List<Integer>>> ufReadsByDirection = PoolUtils.splitReadsByReadDirection(context.getReads(),context.getOffsets());
|
||||
ReadOffsetQuad readsByDir = PoolUtils.splitReadsByReadDirection(context.getReads(),context.getOffsets());
|
||||
Pair<Pair<List<SAMRecord>,List<SAMRecord>>,Pair<List<Integer>,List<Integer>>> 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(),
|
||||
|
|
|
|||
|
|
@ -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<Integer>, SQuad<Long>> {
|
||||
|
||||
@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<Long> reduceInit() {
|
||||
return new SQuad<Long>(0l,0l,0l,0l);
|
||||
}
|
||||
|
||||
/* REDUCE
|
||||
* Add up filtered and unfiltered coverage in the forward
|
||||
* and reverse directions.
|
||||
*/
|
||||
public SQuad<Long> reduce( SQuad<Integer> newCvg, SQuad<Long> oldCvg) {
|
||||
return new SQuad<Long>(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<Integer> 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<Byte> medianQScoresAboveThresh = getMedianQScores(readsByDirectionAboveThresh);
|
||||
SQuad<Double> powerAboveThresh = calculatePower(readsByDirectionAboveThresh, medianQScoresAboveThresh);
|
||||
SQuad<Byte> medianQScoresUnthresh;
|
||||
SQuad<Double> 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<Integer>(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<Double> pow, SQuad<Double> rawPow, SQuad<Byte> medQ, SQuad<Byte> 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<SAMRecord> reads, List<Integer> 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<Byte> getMedianQScores(ReadOffsetQuad reads) {
|
||||
return new SQuad<Byte>(ListUtils.getQScoreMedian(reads.getFirstReads(), reads.getFirstOffsets()),
|
||||
ListUtils.getQScoreMedian(reads.getSecondReads(), reads.getSecondOffsets()),
|
||||
ListUtils.getQScoreMedian(reads.getReadsCombined(), reads.getOffsetsCombined()),
|
||||
(byte) 0);
|
||||
}
|
||||
|
||||
private SQuad<Double> calculatePower(ReadOffsetQuad readQuad, SQuad<Byte> medianQs) {
|
||||
SQuad<Double> power;
|
||||
if( bootstrapIterations > 0 ) {
|
||||
power = bootstrapPowerByDirection(readQuad);
|
||||
} else {
|
||||
power = theoreticalPowerByDirection(readQuad,medianQs);
|
||||
}
|
||||
return power;
|
||||
}
|
||||
|
||||
private SQuad<Double> theoreticalPowerByDirection(ReadOffsetQuad readQuad, SQuad<Byte> medianQs) {
|
||||
return new SQuad<Double>( theoreticalPower(readQuad.numReadsFirst(), medianQs.getFirst()),
|
||||
theoreticalPower(readQuad.numReadsSecond(), medianQs.getSecond()),
|
||||
theoreticalPower(readQuad.numReads(), medianQs.getThird()),
|
||||
0.0 );
|
||||
}
|
||||
|
||||
private SQuad<Double> bootstrapPowerByDirection(ReadOffsetQuad readQuad) {
|
||||
return new SQuad<Double> ( 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<SAMRecord> reads, List<Integer> 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);
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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<Pair<List<SAMRecord>, List<SAMRecord>>, Pair<List<Integer>, List<Integer>>> splitReadsByReadDirection(List<SAMRecord> reads, List<Integer> offsets) {
|
||||
public static ReadOffsetQuad splitReadsByReadDirection(List<SAMRecord> reads, List<Integer> offsets) {
|
||||
ArrayList<SAMRecord> forwardReads;
|
||||
ArrayList<SAMRecord> reverseReads;
|
||||
ArrayList<Integer> 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<SAMRecord>[], List<Integer>[]> splitReadsByBase(List<SAMRecord> reads, List<Integer> offsets) {
|
||||
|
|
@ -117,7 +117,11 @@ public class PoolUtils {
|
|||
}
|
||||
}
|
||||
|
||||
return new Pair(threshReads, threshOffsets);
|
||||
return new Pair<List<SAMRecord>,List<Integer>>(threshReads, threshOffsets);
|
||||
}
|
||||
|
||||
public static Pair<List<SAMRecord>,List<Integer>> thresholdReadsByQuality(Pair<List<SAMRecord>,List<Integer>> 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<Byte> qSNP = getReadBaseQualities(readsBySupport.getFirstReads(), readsBySupport.getFirstOffsets());
|
||||
List<Byte> qRef = getReadBaseQualities(readsBySupport.getSecondReads(), readsBySupport.getSecondOffsets());
|
||||
Pair<Double,Double> logsumSNP = qListToSumLogProbabilities(true, qSNP, 1/alleleFreq);
|
||||
Pair<Double,Double> logsumRef = qListToSumLogProbabilities(false, qRef, 1/alleleFreq);
|
||||
return 0.0 - logsumSNP.getFirst() - logsumRef.getFirst() + logsumSNP.getSecond() + logsumRef.getSecond();
|
||||
}
|
||||
|
||||
public static double calculateLogLikelihoodOfSample(Pair<Pair<List<SAMRecord>,List<SAMRecord>>,Pair<List<Integer>,List<Integer>>> snpReadsRefReads, int nIndivids) {
|
||||
List<Byte> qListSnps = getReadBaseQualities(snpReadsRefReads.getFirst().getFirst(),snpReadsRefReads.getSecond().getFirst());
|
||||
List<Byte> qListRefs = getReadBaseQualities(snpReadsRefReads.getFirst().getSecond(),snpReadsRefReads.getSecond().getSecond());
|
||||
|
|
@ -194,4 +206,23 @@ public class PoolUtils {
|
|||
return new Pair<Double,Double>(logProbObserveXAndSNPTrue,logProbObserveXAndRefTrue);
|
||||
}
|
||||
|
||||
public static ReadOffsetQuad coinTossPartition(List<SAMRecord> reads, List<Integer> offsets, double alleleFreq) {
|
||||
ArrayList<SAMRecord> snpReads = new ArrayList<SAMRecord>();
|
||||
ArrayList<Integer> snpOffsets = new ArrayList<Integer>();
|
||||
ArrayList<SAMRecord> refReads = new ArrayList<SAMRecord>();
|
||||
ArrayList<Integer> refOffsets = new ArrayList<Integer>();
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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<W,X,Y,Z> {
|
||||
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<W,X> a, Pair<Y,Z> 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<W,X> getFirstPair() { return new Pair<W,X>(first,second); }
|
||||
public Pair<Y,Z> getSecondPair() { return new Pair<Y,Z>(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;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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<SAMRecord>,List<Integer>,List<SAMRecord>, List<Integer>> {
|
||||
/*
|
||||
* 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<SAMRecord> a, List<Integer> b, List<SAMRecord> c, List<Integer> d) {
|
||||
super(a,b,c,d);
|
||||
}
|
||||
|
||||
// another constructor that IntelliJ wants
|
||||
public ReadOffsetQuad(Pair<List<SAMRecord>,List<Integer>> a, Pair<List<SAMRecord>,List<Integer>> 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<SAMRecord> getFirstReads() {
|
||||
return this.first;
|
||||
}
|
||||
|
||||
public List<SAMRecord> getSecondReads() {
|
||||
return this.third;
|
||||
}
|
||||
|
||||
public List<SAMRecord> getReadsCombined() {
|
||||
ArrayList<SAMRecord> combined = new ArrayList<SAMRecord>(first);
|
||||
combined.addAll(third);
|
||||
return combined;
|
||||
}
|
||||
|
||||
public List<Integer> getOffsetsCombined() {
|
||||
ArrayList<Integer> combined = new ArrayList<Integer>(second);
|
||||
combined.addAll(fourth);
|
||||
return combined;
|
||||
}
|
||||
|
||||
public List<Integer> getFirstOffsets() {
|
||||
return this.second;
|
||||
}
|
||||
|
||||
public List<Integer> getSecondOffsets() {
|
||||
return this.fourth;
|
||||
}
|
||||
|
||||
public Pair<List<SAMRecord>,List<Integer>> getFirstReadOffsetPair() {
|
||||
return new Pair<List<SAMRecord>,List<Integer>>(first, second);
|
||||
}
|
||||
|
||||
public Pair<List<SAMRecord>,List<Integer>> getSecondReadOffsetPair() {
|
||||
return new Pair<List<SAMRecord>,List<Integer>>(third,fourth);
|
||||
}
|
||||
|
||||
public AlignmentContext getFirstPairAsAlignmentContext(GenomeLoc loc) {
|
||||
return new AlignmentContext(loc, first, second);
|
||||
}
|
||||
|
||||
public AlignmentContext getSecondPairAsAlignmentContext(GenomeLoc loc) {
|
||||
return new AlignmentContext(loc, third, fourth);
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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<X> extends Quad<X,X,X,X> {
|
||||
/* 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
|
||||
|
||||
}
|
||||
Loading…
Reference in New Issue