diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/poolseq/ArtificialPoolWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/poolseq/ArtificialPoolWalker.java index 7298a8673..3668c74c8 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/poolseq/ArtificialPoolWalker.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/poolseq/ArtificialPoolWalker.java @@ -1,356 +1,154 @@ -/* - * Copyright (c) 2009 The Broad Institute - * Permission is hereby granted, free of charge, to any person - * obtaining a copy of this software and associated documentation - * files (the "Software"), to deal in the Software without - * restriction, including without limitation the rights to use, - * copy, modify, merge, publish, distribute, sublicense, and/or sell - * copies of the Software, and to permit persons to whom the - * Software is furnished to do so, subject to the following - * conditions: - * - * The above copyright notice and this permission notice shall be - * included in all copies or substantial portions of the Software. - * - * TERMS OF USE: WE ARE NOT LIABLE FOR ANYTHING - */ - - package org.broadinstitute.sting.playground.gatk.walkers.poolseq; -import org.broadinstitute.sting.gatk.contexts.AlignmentContext; -import org.broadinstitute.sting.gatk.contexts.ReferenceContext; -import org.broadinstitute.sting.gatk.walkers.genotyper.SingleSampleGenotyper; import org.broadinstitute.sting.gatk.walkers.LocusWalker; +import org.broadinstitute.sting.gatk.walkers.genotyper.SingleSampleGenotyper; import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker; +import org.broadinstitute.sting.gatk.contexts.ReferenceContext; +import org.broadinstitute.sting.gatk.contexts.AlignmentContext; +import net.sf.samtools.SAMFileWriter; +import net.sf.samtools.SAMRecord; +import org.broadinstitute.sting.playground.utils.ArtificialPoolContext; import org.broadinstitute.sting.utils.cmdLine.Argument; import org.broadinstitute.sting.utils.StingException; -import org.broadinstitute.sting.utils.Pair; -import org.broadinstitute.sting.utils.genotype.GenotypeCall; -import java.util.*; import java.io.FileOutputStream; import java.io.PrintWriter; import java.io.FileNotFoundException; +import java.io.IOException; +import java.util.Set; +import java.util.List; +import java.util.ListIterator; +import java.util.LinkedList; -import net.sf.samtools.SAMFileWriter; -import net.sf.samtools.SAMRecord; - -// Display the depth of coverage at a given locus and output the power - -public class ArtificialPoolWalker extends LocusWalker[], SAMFileWriter> { - @Argument(fullName="suppressLocusPrinting",doc="Suppress printing",required=false) - public boolean suppressPrinting = false; - @Argument(fullName="reductionProportion", shortName="rp", doc="Proportion of total coverage reflected in pool file",required=false) - public double red_prop = -1; - @Argument(fullName="auxOutputFile", shortName="af", doc="Filepath for Auxiliary pool information (true genotypes, etc.)", required=true) - public String auxFilepath = null; - @Argument(fullName="outputBamFile",shortName="of",doc="Output to this file rather than to standard output") - SAMFileWriter outputBamFile = null; - - private List> readGroupSets; - private PrintWriter auxWrite; - private Pair[] local_genotypes; - private LinkedList[] living_reads; - private SingleSampleGenotyper ssg; - private int npeople; - //TODO: LOCAL CLASS FOR ALL THIS - //@param local_genotypes - holds the genotype (A A/ A C/ etc) for each individual. Updates at each locus. - //@param auxWrite - the writer to the auxiliary file - //@param readGroupSets : holds the readgroups (identifiers for individuals from each read) - //@param ssg - the instance of the single sample genotyper - //@param living_reads - holds on to the selected reads from each person until the location passes the read endpoint - // Updates at each locus. - //@param npeople - number of people represented in this pool (size of readGroupSets) +/** + * Created by IntelliJ IDEA. + * User: chartl + * Date: Aug 26, 2009 + * Time: 11:28:26 AM + * To change this template use File | Settings | File Templates. + */ +public class ArtificialPoolWalker extends LocusWalker { + @Argument(fullName = "AuxOutputFile", shortName = "af", doc = "Auxiliary file for genotyp & coverage output", required = true) + String auxFilePath = null; + @Argument(fullName = "OutputBamFile", shortName = "of", doc = "Output to this file rather than standard output", required = false) + SAMFileWriter outputBamFile = null; public void initialize() { - // initialize the output writer + } + + public ArtificialPoolContext reduceInit() { // try to initialize the file writer + ArtificialPoolContext apContext = new ArtificialPoolContext(); + apContext.setSingleSampleGenotyper(new SingleSampleGenotyper()); + apContext.setReadGroupSets(getToolkit().getMergedReadGroupsByReaders()); + apContext.setAuxWriter(initializeAuxFileWriter(apContext.getTotalNumberOfPeople())); + apContext.setSAMFileWriter(outputBamFile); + apContext.initializeSSG(); + + return apContext; + } + + public ArtificialPoolContext map(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) { + return new ArtificialPoolContext(tracker,ref,context); + } + + public ArtificialPoolContext reduce(ArtificialPoolContext mapCon, ArtificialPoolContext redCon){ + /* ArtificialPoolContext updatedContext = ArtificialPoolContext.mapReduceMerge(mapCon,redCon); + List[] newReads = updatedContext.splitReadsByGroup(updatedContext.getNewReads()); + long[] newCvg = updateRunningCoverage(updatedContext.getRunningCoverage(), getCoverageByGroup(newReads)); + updatedContext.setRunningCoverage(newCvg); + List[] sampledReads = ArtificialPoolContext.sampleReads(newReads,runningCoverageToDouble(newCvg)); + printToFiles(sampledReads,updatedContext);*/ + AlignmentContext context = redCon.getAlignmentContext(); + SAMFileWriter writer = redCon.getSAMFileWriter(); + for(SAMRecord read : context.getReads()) { + writer.addAlignment(read); + } + PrintWriter auxWrite = redCon.getWriterToAuxiliaryFile(); + auxWrite.print("This is a test."); + ArtificialPoolContext updatedContext = redCon; + return updatedContext; + } + + // Helper methods follow + + private PrintWriter initializeAuxFileWriter(int nFiles) { + PrintWriter auxFileWriter; try { - auxWrite = new PrintWriter(new FileOutputStream(auxFilepath)); - } - catch(FileNotFoundException e) { - String ermsg = "auxiliary file for Artificial Pool Walker was either a folder or could not be created"; - throw new StingException(ermsg, e); + auxFileWriter = new PrintWriter(new FileOutputStream(auxFilePath)); + auxFileWriter.println(createAuxFileHeader(nFiles)); + } catch(FileNotFoundException e) { + String errmsg = "The filepath you entered "+auxFilePath+" could not be opened. Please double-check that the input is correct."; + throw new StingException(errmsg, e); + } catch(IOException e) { + String errmsg = "The file you entered "+auxFilePath+" could not be written to. Please check your permissions to write to this file."; + throw new StingException(errmsg,e); } - // create the file header and write to file - - auxWrite.printf("%s%n",createFileHeader()); - - // obtain the readgroup sets - - readGroupSets = getToolkit().getMergedReadGroupsByReaders(); - npeople = readGroupSets.size(); - - // if a reduction proportion was unspecified, set it to 1/num_files - - if(red_prop <= 0) { - red_prop = 1.0/npeople; - } - - // initialize the local genotype array - - local_genotypes = new Pair[npeople]; - - // initilize the single sample genotyper - - ssg = new SingleSampleGenotyper(); - ssg.initialize(); - - // initialize the living reads array - - living_reads = new LinkedList[npeople]; - - for(int iter=0; iter < readGroupSets.size(); iter++) { - living_reads[iter] = new LinkedList(); - } + return auxFileWriter; } - public List[] map(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) { - - updateLiving(); - // each time we move to the next locus, remove from the coverage count those reads that ended - auxWrite.printf("%s:%s",context.getContig(),context.getPosition()); // TODO: PUT IN REDUCE - - return getNewReadsAndGenotypesByGroup(tracker, ref, context); - } - - public SAMFileWriter reduceInit() { - - return outputBamFile; - } - - public SAMFileWriter reduce(List[] readsByReadGroup, SAMFileWriter outFile) { - // want to generate a sub sample of the new reads - int total_new_coverage = 0; - - for(int pers = 0; pers < npeople; pers++) { - total_new_coverage += readsByReadGroup[pers].size(); - } - - int sought_coverage = (int) Math.ceil(red_prop * (double)total_new_coverage); - List[] randomReadsByGroup = drawReadsRandomlyFromReadsByGroup(readsByReadGroup,sought_coverage); - printToFileAndAuxFile(randomReadsByGroup,sought_coverage,outFile); - return outFile; - } - - - - // ArtificalPoolWalker helper methods begin below - - private List[] drawReadsRandomlyFromReadsByGroup(List[] readsByReadGroup, int numToDraw){ - - int[] simulatedCoverage = simulateRandomUniformDraws(numToDraw, readsByReadGroup); - - return removeReadsToCoverage(readsByReadGroup,simulatedCoverage); - } - - private List[] removeReadsToCoverage(List[] readsByReadGroup, int[] simulatedCoverage) { - for(int group = 0; group < npeople; group++) { - int cvgAtGroup = readsByReadGroup[group].size(); - for(int drawDown = simulatedCoverage[group]; drawDown > 0; drawDown --) { - readsByReadGroup[group].remove((int)Math.floor((double) cvgAtGroup*Math.random())); - cvgAtGroup--; - } - } - - return readsByReadGroup; - } - - private int[] simulateRandomUniformDraws(int numToDraw, List[] readsByReadGroup) { - int[] drawsByPerson = new int[npeople]; - - for(int pers = 0; pers < npeople; pers++) { //array instantiation - drawsByPerson[pers] = 0; - } - - for(int draw = 0; draw < numToDraw; draw++) { - int personToDrawFrom = (int) Math.floor((double) npeople * Math.random()); - drawsByPerson[personToDrawFrom]++; - } - - return redistributeSimulatedCoverage(drawsByPerson,readsByReadGroup,0, 0, null); - } - - private int[] redistributeSimulatedCoverage(int[] drawsByPerson, List[] readsByReadGroup, int excess, int nOffLimits, boolean[] offLimits) { - int [] result; - if(offLimits == null) { // first time calling - offLimits = new boolean[npeople]; - for(int i = 0; i < npeople; i++) { - offLimits[i] = false; - } - result = redistributeSimulatedCoverage(drawsByPerson,readsByReadGroup,excess,0, offLimits); - } else if (excess == 0) { // no overt excess, need to check individual persons - // get excess coverage - int newExcess = 0; - - for(int pers = 0; pers < npeople; pers++) { - if(!offLimits[pers] && drawsByPerson[pers] > readsByReadGroup[pers].size()) { - newExcess += drawsByPerson[pers] - readsByReadGroup[pers].size(); - drawsByPerson[pers] = readsByReadGroup[pers].size(); - offLimits[pers] = true; - nOffLimits++; - } else { - // do nothing - } - } - - if(newExcess == 0) { // we are done! - result = drawsByPerson; - } else { // there is now overt excess - result = redistributeSimulatedCoverage(drawsByPerson,readsByReadGroup,newExcess,nOffLimits,offLimits); - } - } else { // there are overt excess reads to distribute - int nRemaining = npeople - nOffLimits; - int[] drawsToAdd = new int[nRemaining]; - - for(int j = 0; j < nRemaining; j++) { - drawsToAdd[j] = 0; - } - - for(int draw = excess; draw > 0; draw--) { - drawsToAdd[(int) Math.floor((double)nRemaining * Math.random())]++; - } - - int notOffLimitsCounter = 0; - - for(int pers = 0; pers < npeople; pers++) { - if(! offLimits[pers]) { - drawsByPerson[pers] += drawsToAdd[notOffLimitsCounter]; - notOffLimitsCounter++; - } else { - // do nothing - } - } - result = redistributeSimulatedCoverage(drawsByPerson,readsByReadGroup,0,nOffLimits,offLimits); - } - return result; - } - - private int[] createNewPersonIndex(int peopleWithReadsLeft, int[] indexLeft, int indexToRemoveFrom) { - int[] newPersonIndex = new int[peopleWithReadsLeft-1]; - for(int i = 0; i < peopleWithReadsLeft; i++) { - if(i < indexToRemoveFrom) { - newPersonIndex[i] = indexLeft[i]; - } else { - newPersonIndex[i] = indexLeft[i+1]; - } - } - return newPersonIndex; - } - - - private List[] getNewReadsAndGenotypesByGroup(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) - { - // at the same time we will: - // 1) split reads up by read group - // 2) get new reads - // 3) get genotypes for each individual - // then we will: - // return the new reads (a subsample), stratified by read group - - List allReads = context.getReads(); - int n_reads_total = allReads.size(); - List allOffsets = context.getOffsets(); - ListIterator readIterator = allReads.listIterator(); - ListIterator offsetIterator = allOffsets.listIterator(); - LinkedList[] reads_by_group = new LinkedList[npeople]; - LinkedList[] offsets_by_group = new LinkedList[npeople]; - LinkedList[] new_reads_by_group = new LinkedList[npeople]; - - for(int j = 0; j < npeople; j++) { - reads_by_group[j] = new LinkedList(); - offsets_by_group[j] = new LinkedList(); - new_reads_by_group[j] = new LinkedList(); - } - - while(readIterator.hasNext()) { - int iter_person = 0; - SAMRecord read = (SAMRecord) readIterator.next(); - int offset = (Integer) offsetIterator.next(); - while(!readGroupSets.get(iter_person).contains((String) read.getAttribute("RG"))) { - iter_person++; - // for debug purposes only: - if(iter_person > npeople) { - throw new StingException("Read not found in any read group in ArtificialPoolWalker. Life sucks."); - } - } - // here, iter_person is the person from whom the read came so: - reads_by_group[iter_person].add(read); - offsets_by_group[iter_person].add(offset); - - // check to see if the read is new - if(offset == 0) { - new_reads_by_group[iter_person].add(read); - // add this to the new reads set to be randomly selected and printed - } - } - // now we obtain the genotypes - for(int iter=0; iter[] readsByGroup) { + int[] coverage = new int[readsByGroup.length]; + for(int iterator = 0; iterator < readsByGroup.length; iterator ++) { + coverage[iterator] = readsByGroup[iterator].size(); + } - int readLengthLeft; - ListIterator updater; + return coverage; + } - for(int j = 0; j < npeople; j++) { - updater = living_reads[j].listIterator(); + private long[] updateRunningCoverage(long[] cvgUpToNow, int[] newCvgByGroup) { + long[] newCvg = new long[cvgUpToNow.length]; + for(int iter = 0; iter < cvgUpToNow.length; iter++) { + newCvg[iter] = cvgUpToNow[iter] + newCvgByGroup[iter]; + } - while(updater.hasNext()) { - readLengthLeft = (Integer) updater.next(); - if(readLengthLeft <= 1) { - updater.remove(); - } else { - updater.set(readLengthLeft-1); - } + return newCvg; + } - } // end while(updater.hasNext()) - } // end for(int j=0; j < npeople; j++) - } + private double[] runningCoverageToDouble(long[] cvg) { + double[] avgProp = new double[cvg.length]; + long sum = 0; + for(long elem : cvg) { + sum += elem; + } + for(int iter = 0; iter < cvg.length; iter++) { + avgProp[iter] = cvg[iter]/sum; + } - private void printToFileAndAuxFile(List[] downSampledReadsByGroup, int cvg, SAMFileWriter outFile) - { - // NOTE: the chromosome and position were already printed in the auxFile during map. - for(int pers = 0; pers < npeople; pers++) { - List personReads = downSampledReadsByGroup[pers]; - for(SAMRecord read : personReads) { - outFile.addAlignment(read); - living_reads[pers].add(read.getReadLength()); + return avgProp; + } + + private void printToFiles(List[] sampledNewReads, ArtificialPoolContext context) { + SAMFileWriter samWrite = context.getSAMFileWriter(); + String sp = " "; + PrintWriter auxWrite = context.getWriterToAuxiliaryFile(); + int readGroupInt = 0; + for(List readGroup : sampledNewReads) { + for(SAMRecord read : readGroup) { + samWrite.addAlignment(read); } + auxWrite.print(context.getAlignmentContext().getLocation().toString() + sp); + auxWrite.print(context.genotypeAndConfidenceToString(readGroupInt,sp)); + readGroupInt++; } - String auxformstring = " %s %f %d"; - // now print to the aux file - int totalcvg = 0; - for(int pers = 0; pers < npeople; pers++) { - auxWrite.printf(auxformstring, local_genotypes[pers].first, local_genotypes[pers].second,living_reads[pers].size()); - totalcvg += living_reads[pers].size(); - } - auxWrite.printf(" %d%n", totalcvg); } + } + 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 1e489e77e..0dab39fdd 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 @@ -35,8 +35,11 @@ public class CoverageAndPowerWalker extends LocusWalker, @Argument(fullName="bootStrap", shortName="bs", doc="Use a bootstrap method", required=false) public boolean useBootstrap = false; - @Argument(fullName="lodThreshold", shortName="lt", doc="Threshold for LOD score for calls") - public double threshold = 3.0; + @Argument(fullName="lodThreshold", shortName="lt", doc="Threshold for LOD score for calls", required=false) + public double lodThreshold = 3.0; + + @Argument(fullName="minQScore", shortName="qm", doc="Threshold for the minimum quality score, filter out bases below",required=false) + public byte minQScore = 0; private static final int BOOTSTRAP_ITERATIONS = 300; @@ -58,11 +61,18 @@ public class CoverageAndPowerWalker extends LocusWalker, } public Pair map(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) { - Pair,List>,Pair,List>> readsByDirection = PoolUtils.splitReadsByReadDirection(context.getReads(),context.getOffsets()); + AlignmentContext filteredContext; + if (this.getMinQualityScore() <= 0) { + filteredContext = context; + } else { + 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()); if ( ! suppress_printing) { - Pair powers = calculatePower(readsByDirection, useBootstrap, context); - out.printf("%s: %d %d %d %d %d %d %f %f %f%n", context.getLocation(), readsByDirection.getFirst().getFirst().size(), readsByDirection.getFirst().getSecond().size(), - context.getReads().size(), powers.getSecond()[0], powers.getSecond()[1], powers.getSecond()[2], + 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(), + filteredContext.getReads().size(), powers.getSecond()[0], powers.getSecond()[1], powers.getSecond()[2], powers.getFirst()[0], powers.getFirst()[1], powers.getFirst()[2]); } return new Pair(readsByDirection.getFirst().getFirst().size(),readsByDirection.getFirst().getSecond().size()); @@ -73,9 +83,9 @@ public class CoverageAndPowerWalker extends LocusWalker, public Pair calculatePower(Pair,List>,Pair,List>> dirReads, boolean bootStrap, AlignmentContext context) { Pair powerAndMedianQ; if( bootStrap ) { - powerAndMedianQ = calculateDirectionalPowersByBootstrap(dirReads,threshold,context); + powerAndMedianQ = calculateDirectionalPowersByBootstrap(dirReads,lodThreshold,context); } else { - powerAndMedianQ = calculateDirectionalPowersTheoretical(dirReads,threshold,context); + powerAndMedianQ = calculateDirectionalPowersTheoretical(dirReads,lodThreshold,context); } return powerAndMedianQ; @@ -162,7 +172,15 @@ public class CoverageAndPowerWalker extends LocusWalker, } public String createHeaderString() { - return "Chrom:Pos CvgF CvgR CvgC QmedF QmedR QmedC PowF PowR PowC"; + return "Chrom:Pos ForwardCoverage ReverseCoverage CombinedCoverage ForwardMedianQ ReverseMedianQ CombinedMedianQ PowerForward PowerReverse PowerCombined"; + } + + public byte getMinQualityScore() { + return minQScore; + } + + public Pair,List> filterByQuality(List reads, List offsets, byte qMin) { + return PoolUtils.thresholdReadsByQuality(reads,offsets,qMin); } // class methods diff --git a/java/src/org/broadinstitute/sting/playground/utils/PoolUtils.java b/java/src/org/broadinstitute/sting/playground/utils/PoolUtils.java index b87f89b7f..9319b7a7c 100755 --- a/java/src/org/broadinstitute/sting/playground/utils/PoolUtils.java +++ b/java/src/org/broadinstitute/sting/playground/utils/PoolUtils.java @@ -18,6 +18,12 @@ public class PoolUtils { private PoolUtils () {} + public static final int BASE_A_OFFSET = 0; + public static final int BASE_C_OFFSET = 1; + public static final int BASE_G_OFFSET = 2; + 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) { ArrayList forwardReads; ArrayList reverseReads; @@ -48,5 +54,108 @@ public class PoolUtils { return new Pair(new Pair(forwardReads,reverseReads), new Pair(forwardOffsets,reverseOffsets)); } - + + public static Pair[], List[]> splitReadsByBase(List reads, List offsets) { + ArrayList[] readsByBase; + ArrayList[] offsetsByBase; + if ( reads == null ) { + readsByBase = null; + offsetsByBase = null; + } else { + readsByBase = new ArrayList[4]; + offsetsByBase = new ArrayList[4]; + for(int readNum = 0; readNum < reads.size(); readNum++) { + switch (reads.get(readNum).getReadBases()[offsets.get(readNum)]) { + case 'A': + case 'a': readsByBase[0].add(reads.get(readNum)); + offsetsByBase[0].add(offsets.get(readNum)); + break; + case 'C': + case 'c': readsByBase[1].add(reads.get(readNum)); + offsetsByBase[1].add(offsets.get(readNum)); + break; + case 'G': + case 'g': readsByBase[2].add(reads.get(readNum)); + offsetsByBase[2].add(offsets.get(readNum)); + break; + case 'T': + case 't': readsByBase[3].add(reads.get(readNum)); + offsetsByBase[3].add(offsets.get(readNum)); + break; + default: break; // TODO: INDEL AWARENESS + } + } + } + return new Pair(readsByBase,offsetsByBase); + } + + public static Pair, List> thresholdReadsByQuality(List reads, List offsets, byte qThresh) { + List threshReads; + List threshOffsets; + if(reads == null) { + threshReads=null; + threshOffsets = null; + } else if (qThresh <= 0) { + threshReads = reads; + threshOffsets = offsets; + } else { + threshReads = new ArrayList(); + threshOffsets = new ArrayList(); + + for ( int readNo = 0; readNo < reads.size(); readNo ++) { + if ( reads.get(readNo).getBaseQualities()[offsets.get(readNo)] >= qThresh) { + threshReads.add(reads.get(readNo)); + threshOffsets.add(offsets.get(readNo)); + } // else do nothing + } + } + + return new Pair(threshReads,threshOffsets); + } + + public static int getBaseOffset(char base) { + switch(base) { + case 'A': + case 'a': + return getBaseAOffset(); + case 'C': + case 'c': + return getBaseCOffset(); + case 'G': + case 'g': + return getBaseGOffset(); + case 'T': + case 't': + return getBaseTOffset(); + default: + return -1; + } + //TODO: indel offsets + } + + public static int getBaseAOffset() { + return BASE_A_OFFSET; + } + + public static int getBaseCOffset() { + return BASE_C_OFFSET; + } + + public static int getBaseGOffset() { + return BASE_G_OFFSET; + } + + public static int getBaseTOffset() { + return BASE_T_OFFSET; + } + + public static List getListOfBaseQualities(List reads,List offsets) { + //TODO: this is a terrible method name. Change it to something better. + List qualities = new ArrayList(reads.size()); + for (int readNo = 0; readNo < reads.size(); readNo ++) { + qualities.add(reads.get(readNo).getBaseQualities()[offsets.get(readNo)]); + } + + return qualities; + } }