diff --git a/ivy.xml b/ivy.xml
index b7ca65406..1802c1627 100644
--- a/ivy.xml
+++ b/ivy.xml
@@ -61,6 +61,7 @@
+
@@ -82,7 +83,7 @@
-
+
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/downsampling/AlleleBiasedDownsamplingUtils.java b/protected/java/src/org/broadinstitute/sting/gatk/downsampling/AlleleBiasedDownsamplingUtils.java
index a7bb58d0c..02821ab50 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/downsampling/AlleleBiasedDownsamplingUtils.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/downsampling/AlleleBiasedDownsamplingUtils.java
@@ -47,7 +47,6 @@
package org.broadinstitute.sting.gatk.downsampling;
import net.sf.samtools.SAMReadGroupRecord;
-import net.sf.samtools.SAMRecord;
import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.pileup.*;
import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
@@ -78,30 +77,15 @@ public class AlleleBiasedDownsamplingUtils {
for ( int i = 0; i < 4; i++ )
alleleStratifiedElements[i] = new ArrayList();
- // keep all of the reduced reads
- final ArrayList reducedReadPileups = new ArrayList();
-
// start by stratifying the reads by the alleles they represent at this position
- for( final PileupElement pe : pileup ) {
+ for ( final PileupElement pe : pileup ) {
// we do not want to remove a reduced read
- if ( pe.getRead().isReducedRead() )
- reducedReadPileups.add(pe);
-
- final int baseIndex = BaseUtils.simpleBaseToBaseIndex(pe.getBase());
- if ( baseIndex != -1 )
- alleleStratifiedElements[baseIndex].add(pe);
- }
-
- // Unfortunately, we need to maintain the original pileup ordering of reads or FragmentUtils will complain later.
- int numReadsToRemove = (int)(pileup.getNumberOfElements() * downsamplingFraction); // floor
- final TreeSet elementsToKeep = new TreeSet(new Comparator() {
- @Override
- public int compare(PileupElement element1, PileupElement element2) {
- final int difference = element1.getRead().getAlignmentStart() - element2.getRead().getAlignmentStart();
- return difference != 0 ? difference : element1.getRead().getReadName().compareTo(element2.getRead().getReadName());
+ if ( !pe.getRead().isReducedRead() ) {
+ final int baseIndex = BaseUtils.simpleBaseToBaseIndex(pe.getBase());
+ if ( baseIndex != -1 )
+ alleleStratifiedElements[baseIndex].add(pe);
}
- });
- elementsToKeep.addAll(reducedReadPileups);
+ }
// make a listing of allele counts
final int[] alleleCounts = new int[4];
@@ -109,22 +93,30 @@ public class AlleleBiasedDownsamplingUtils {
alleleCounts[i] = alleleStratifiedElements[i].size();
// do smart down-sampling
+ int numReadsToRemove = (int)(pileup.getNumberOfElements() * downsamplingFraction); // floor
final int[] targetAlleleCounts = runSmartDownsampling(alleleCounts, numReadsToRemove);
+ final HashSet readsToRemove = new HashSet(numReadsToRemove);
for ( int i = 0; i < 4; i++ ) {
final ArrayList alleleList = alleleStratifiedElements[i];
- // if we don't need to remove any reads, keep them all
- if ( alleleList.size() <= targetAlleleCounts[i] )
- elementsToKeep.addAll(alleleList);
- else
- elementsToKeep.addAll(downsampleElements(alleleList, alleleList.size() - targetAlleleCounts[i], log));
+ // if we don't need to remove any reads, then don't
+ if ( alleleList.size() > targetAlleleCounts[i] )
+ readsToRemove.addAll(downsampleElements(alleleList, alleleList.size() - targetAlleleCounts[i], log));
}
// clean up pointers so memory can be garbage collected if needed
for ( int i = 0; i < 4; i++ )
alleleStratifiedElements[i].clear();
- return new ReadBackedPileupImpl(pileup.getLocation(), new ArrayList(elementsToKeep));
+ // we need to keep the reads sorted because the FragmentUtils code will expect them in coordinate order and will fail otherwise
+ final List readsToKeep = new ArrayList(pileup.getNumberOfElements() - numReadsToRemove);
+ for ( final PileupElement pe : pileup ) {
+ if ( !readsToRemove.contains(pe) ) {
+ readsToKeep.add(pe);
+ }
+ }
+
+ return new ReadBackedPileupImpl(pileup.getLocation(), new ArrayList(readsToKeep));
}
private static int scoreAlleleCounts(final int[] alleleCounts) {
@@ -188,37 +180,43 @@ public class AlleleBiasedDownsamplingUtils {
}
/**
- * Performs allele biased down-sampling on a pileup and computes the list of elements to keep
+ * Performs allele biased down-sampling on a pileup and computes the list of elements to remove
*
* @param elements original list of records
* @param numElementsToRemove the number of records to remove
* @param log logging output
- * @return the list of pileup elements TO KEEP
+ * @return the list of pileup elements TO REMOVE
*/
- private static List downsampleElements(final ArrayList elements, final int numElementsToRemove, final PrintStream log) {
- if ( numElementsToRemove == 0 )
- return elements;
+ private static List downsampleElements(final List elements, final int numElementsToRemove, final PrintStream log) {
+ ArrayList elementsToRemove = new ArrayList(numElementsToRemove);
+ // are there no elements to remove?
+ if ( numElementsToRemove == 0 )
+ return elementsToRemove;
+
+ // should we remove all of the elements?
final int pileupSize = elements.size();
if ( numElementsToRemove == pileupSize ) {
logAllElements(elements, log);
- return new ArrayList(0);
+ elementsToRemove.addAll(elements);
+ return elementsToRemove;
}
+ // create a bitset describing which elements to remove
final BitSet itemsToRemove = new BitSet(pileupSize);
for ( Integer selectedIndex : MathUtils.sampleIndicesWithoutReplacement(pileupSize, numElementsToRemove) ) {
itemsToRemove.set(selectedIndex);
}
- ArrayList elementsToKeep = new ArrayList(pileupSize - numElementsToRemove);
for ( int i = 0; i < pileupSize; i++ ) {
- if ( itemsToRemove.get(i) )
- logRead(elements.get(i).getRead(), log);
- else
- elementsToKeep.add(elements.get(i));
+ if ( itemsToRemove.get(i) ) {
+ final T element = elements.get(i);
+ logElement(element, log);
+ elementsToRemove.add(element);
+ }
}
- return elementsToKeep;
+ return elementsToRemove;
}
/**
@@ -252,65 +250,30 @@ public class AlleleBiasedDownsamplingUtils {
final List alleleBin = alleleReadMap.get(alleles.get(i));
if ( alleleBin.size() > targetAlleleCounts[i] ) {
- readsToRemove.addAll(downsampleReads(alleleBin, alleleBin.size() - targetAlleleCounts[i], log));
+ readsToRemove.addAll(downsampleElements(alleleBin, alleleBin.size() - targetAlleleCounts[i], log));
}
}
return readsToRemove;
}
- /**
- * Performs allele biased down-sampling on a pileup and computes the list of elements to remove
- *
- * @param reads original list of records
- * @param numElementsToRemove the number of records to remove
- * @param log logging output
- * @return the list of pileup elements TO REMOVE
- */
- private static List downsampleReads(final List reads, final int numElementsToRemove, final PrintStream log) {
- final ArrayList readsToRemove = new ArrayList(numElementsToRemove);
-
- if ( numElementsToRemove == 0 )
- return readsToRemove;
-
- final int pileupSize = reads.size();
- if ( numElementsToRemove == pileupSize ) {
- logAllReads(reads, log);
- return reads;
- }
-
- final BitSet itemsToRemove = new BitSet(pileupSize);
- for ( Integer selectedIndex : MathUtils.sampleIndicesWithoutReplacement(pileupSize, numElementsToRemove) ) {
- itemsToRemove.set(selectedIndex);
- }
-
- for ( int i = 0; i < pileupSize; i++ ) {
- if ( itemsToRemove.get(i) ) {
- final GATKSAMRecord read = reads.get(i);
- readsToRemove.add(read);
- logRead(read, log);
+ private static void logAllElements(final List elements, final PrintStream log) {
+ if ( log != null ) {
+ for ( final T obj : elements ) {
+ logElement(obj, log);
}
}
-
- return readsToRemove;
}
- private static void logAllElements(final List elements, final PrintStream log) {
+ private static void logElement(final T obj, final PrintStream log) {
if ( log != null ) {
- for ( final PileupElement p : elements )
- logRead(p.getRead(), log);
- }
- }
- private static void logAllReads(final List reads, final PrintStream log) {
- if ( log != null ) {
- for ( final GATKSAMRecord read : reads )
- logRead(read, log);
- }
- }
+ final GATKSAMRecord read;
+ if ( obj instanceof PileupElement )
+ read = ((PileupElement)obj).getRead();
+ else
+ read = (GATKSAMRecord)obj;
- private static void logRead(final SAMRecord read, final PrintStream log) {
- if ( log != null ) {
final SAMReadGroupRecord readGroup = read.getReadGroup();
log.println(String.format("%s\t%s\t%s\t%s", read.getReadName(), readGroup.getSample(), readGroup.getLibrary(), readGroup.getPlatformUnit()));
}
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/DepthOfCoverage.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/DepthOfCoverage.java
index aeec36c18..4adb2ca71 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/DepthOfCoverage.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/DepthOfCoverage.java
@@ -99,7 +99,7 @@ public class DepthOfCoverage extends InfoFieldAnnotation implements StandardAnno
for (PerReadAlleleLikelihoodMap maps : perReadAlleleLikelihoodMap.values() ) {
for (Map.Entry> el : maps.getLikelihoodReadMap().entrySet()) {
final GATKSAMRecord read = el.getKey();
- depth += (read.isReducedRead() ? read.getReducedCount(ReadUtils.getReadCoordinateForReferenceCoordinate(read, vc.getStart(), ReadUtils.ClippingTail.RIGHT_TAIL)) : 1);
+ depth += (read.isReducedRead() ? read.getReducedCount(ReadUtils.getReadCoordinateForReferenceCoordinateUpToEndOfRead(read, vc.getStart(), ReadUtils.ClippingTail.RIGHT_TAIL)) : 1);
}
}
}
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/DepthPerAlleleBySample.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/DepthPerAlleleBySample.java
index a194fe323..5acea12f6 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/DepthPerAlleleBySample.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/DepthPerAlleleBySample.java
@@ -144,7 +144,7 @@ public class DepthPerAlleleBySample extends GenotypeAnnotation implements Standa
continue; // read is non-informative
if (!vc.getAlleles().contains(a))
continue; // sanity check - shouldn't be needed
- alleleCounts.put(a, alleleCounts.get(a) + (read.isReducedRead() ? read.getReducedCount(ReadUtils.getReadCoordinateForReferenceCoordinate(read, vc.getStart(), ReadUtils.ClippingTail.RIGHT_TAIL)) : 1));
+ alleleCounts.put(a, alleleCounts.get(a) + (read.isReducedRead() ? read.getReducedCount(ReadUtils.getReadCoordinateForReferenceCoordinateUpToEndOfRead(read, vc.getStart(), ReadUtils.ClippingTail.RIGHT_TAIL)) : 1));
}
final int[] counts = new int[alleleCounts.size()];
counts[0] = alleleCounts.get(vc.getReference());
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/FisherStrand.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/FisherStrand.java
index 167e5df63..ff3d7940f 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/FisherStrand.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/FisherStrand.java
@@ -116,8 +116,8 @@ public class FisherStrand extends InfoFieldAnnotation implements StandardAnnotat
else if (table1 == null)
return annotationForOneTable(pValueForContingencyTable(table2));
else { // take the one with the best (i.e., least significant pvalue)
- double pvalue1 = Math.max(pValueForContingencyTable(table1), MIN_PVALUE);
- double pvalue2 = Math.max(pValueForContingencyTable(table2), MIN_PVALUE);
+ double pvalue1 = pValueForContingencyTable(table1);
+ double pvalue2 = pValueForContingencyTable(table2);
return annotationForOneTable(Math.max(pvalue1, pvalue2));
}
}
@@ -129,7 +129,7 @@ public class FisherStrand extends InfoFieldAnnotation implements StandardAnnotat
* @return a hash map from FS -> phred-scaled pValue
*/
private Map annotationForOneTable(final double pValue) {
- final Object value = String.format("%.3f", QualityUtils.phredScaleErrorRate(pValue));
+ final Object value = String.format("%.3f", QualityUtils.phredScaleErrorRate(Math.max(pValue, MIN_PVALUE))); // prevent INFINITYs
return Collections.singletonMap(FS, value);
// Map map = new HashMap();
// map.put(FS, String.format("%.3f", QualityUtils.phredScaleErrorRate(pValue)));
@@ -265,24 +265,16 @@ public class FisherStrand extends InfoFieldAnnotation implements StandardAnnotat
for (PerReadAlleleLikelihoodMap maps : stratifiedPerReadAlleleLikelihoodMap.values() ) {
for (Map.Entry> el : maps.getLikelihoodReadMap().entrySet()) {
- final boolean matchesRef = PerReadAlleleLikelihoodMap.getMostLikelyAllele(el.getValue()).equals(ref,true);
- final boolean matchesAlt = PerReadAlleleLikelihoodMap.getMostLikelyAllele(el.getValue()).equals(alt,true);
-
- if ( !matchesRef && !matchesAlt )
- continue;
-
- boolean isFW = el.getKey().getReadNegativeStrandFlag();
-
- int row = matchesRef ? 0 : 1;
- int column = isFW ? 0 : 1;
-
+ final Allele mostLikelyAllele = PerReadAlleleLikelihoodMap.getMostLikelyAllele(el.getValue());
final GATKSAMRecord read = el.getKey();
- table[row][column] += (read.isReducedRead() ? read.getReducedCount(ReadUtils.getReadCoordinateForReferenceCoordinate(read, vc.getStart(), ReadUtils.ClippingTail.RIGHT_TAIL)) : 1);
+ final int representativeCount = read.isReducedRead() ? read.getReducedCount(ReadUtils.getReadCoordinateForReferenceCoordinateUpToEndOfRead(read, vc.getStart(), ReadUtils.ClippingTail.RIGHT_TAIL)) : 1;
+ updateTable(table, mostLikelyAllele, read, ref, alt, representativeCount);
}
}
return table;
}
+
/**
Allocate and fill a 2x2 strand contingency table. In the end, it'll look something like this:
* fw rc
@@ -299,31 +291,36 @@ public class FisherStrand extends InfoFieldAnnotation implements StandardAnnotat
for ( Map.Entry sample : stratifiedContexts.entrySet() ) {
for (PileupElement p : sample.getValue().getBasePileup()) {
- // ignore reduced reads because they are always on the forward strand!
- // TODO -- when het compression is enabled in RR, we somehow need to allow those reads through into the Fisher test
- if ( p.getRead().isReducedRead() )
- continue;
-
if ( ! RankSumTest.isUsableBase(p, false) ) // ignore deletions
continue;
if ( p.getQual() < minQScoreToConsider || p.getMappingQual() < minQScoreToConsider )
continue;
- final Allele base = Allele.create(p.getBase(), false);
- final boolean isFW = !p.getRead().getReadNegativeStrandFlag();
-
- final boolean matchesRef = ref.equals(base, true);
- final boolean matchesAlt = alt.equals(base, true);
- if ( matchesRef || matchesAlt ) {
- int row = matchesRef ? 0 : 1;
- int column = isFW ? 0 : 1;
-
- table[row][column] += p.getRepresentativeCount();
- }
+ updateTable(table, Allele.create(p.getBase(), false), p.getRead(), ref, alt, p.getRepresentativeCount());
}
}
return table;
}
+
+ private static void updateTable(final int[][] table, final Allele allele, final GATKSAMRecord read, final Allele ref, final Allele alt, final int representativeCount) {
+ // ignore reduced reads because they are always on the forward strand!
+ // TODO -- when het compression is enabled in RR, we somehow need to allow those reads through into the Fisher test
+ if ( read.isReducedRead() )
+ return;
+
+ final boolean matchesRef = allele.equals(ref, true);
+ final boolean matchesAlt = allele.equals(alt, true);
+
+ if ( matchesRef || matchesAlt ) {
+
+ final boolean isFW = !read.getReadNegativeStrandFlag();
+
+ int row = matchesRef ? 0 : 1;
+ int column = isFW ? 0 : 1;
+
+ table[row][column] += representativeCount;
+ }
+ }
}
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/GCContent.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/GCContent.java
index 3bb3d7d5a..2b3290595 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/GCContent.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/GCContent.java
@@ -95,9 +95,9 @@ public class GCContent extends InfoFieldAnnotation implements ExperimentalAnnota
for ( byte base : ref.getBases() ) {
int baseIndex = BaseUtils.simpleBaseToBaseIndex(base);
- if ( baseIndex == BaseUtils.gIndex || baseIndex == BaseUtils.cIndex )
+ if ( baseIndex == BaseUtils.Base.G.ordinal() || baseIndex == BaseUtils.Base.C.ordinal() )
gc++;
- else if ( baseIndex == BaseUtils.aIndex || baseIndex == BaseUtils.tIndex )
+ else if ( baseIndex == BaseUtils.Base.A.ordinal() || baseIndex == BaseUtils.Base.T.ordinal() )
at++;
else
; // ignore
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/HaplotypeScore.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/HaplotypeScore.java
index fe4075117..3acba48ae 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/HaplotypeScore.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/HaplotypeScore.java
@@ -56,7 +56,6 @@ import org.broadinstitute.sting.gatk.walkers.annotator.interfaces.InfoFieldAnnot
import org.broadinstitute.sting.gatk.walkers.annotator.interfaces.StandardAnnotation;
import org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap;
import org.broadinstitute.variant.utils.BaseUtils;
-import org.broadinstitute.sting.utils.Haplotype;
import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.QualityUtils;
import org.broadinstitute.variant.vcf.VCFHeaderLineType;
@@ -236,8 +235,8 @@ public class HaplotypeScore extends InfoFieldAnnotation implements StandardAnnot
final byte[] haplotypeBases = new byte[contextSize];
Arrays.fill(haplotypeBases, (byte) REGEXP_WILDCARD);
- final double[] baseQualities = new double[contextSize];
- Arrays.fill(baseQualities, 0.0);
+ final byte[] baseQualities = new byte[contextSize];
+ Arrays.fill(baseQualities, (byte)0);
byte[] readBases = read.getReadBases();
readBases = AlignmentUtils.readToAlignmentByteArray(read.getCigar(), readBases); // Adjust the read bases based on the Cigar string
@@ -267,7 +266,7 @@ public class HaplotypeScore extends InfoFieldAnnotation implements StandardAnnot
readQuals[baseOffset] = (byte) 0;
} // quals less than 5 are used as codes and don't have actual probabilistic meaning behind them
haplotypeBases[i] = readBases[baseOffset];
- baseQualities[i] = (double) readQuals[baseOffset];
+ baseQualities[i] = readQuals[baseOffset];
}
return new Haplotype(haplotypeBases, baseQualities);
@@ -286,10 +285,10 @@ public class HaplotypeScore extends InfoFieldAnnotation implements StandardAnnot
final int length = a.length;
final byte[] consensusChars = new byte[length];
- final double[] consensusQuals = new double[length];
+ final int[] consensusQuals = new int[length];
- final double[] qualsA = haplotypeA.getQuals();
- final double[] qualsB = haplotypeB.getQuals();
+ final int[] qualsA = haplotypeA.getQuals();
+ final int[] qualsB = haplotypeB.getQuals();
for (int i = 0; i < length; i++) {
chA = a[i];
@@ -300,7 +299,7 @@ public class HaplotypeScore extends InfoFieldAnnotation implements StandardAnnot
if ((chA == wc) && (chB == wc)) {
consensusChars[i] = wc;
- consensusQuals[i] = 0.0;
+ consensusQuals[i] = 0;
} else if ((chA == wc)) {
consensusChars[i] = chB;
consensusQuals[i] = qualsB[i];
@@ -433,7 +432,6 @@ public class HaplotypeScore extends InfoFieldAnnotation implements StandardAnnot
}
-
public List getKeyNames() {
return Arrays.asList("HaplotypeScore");
}
@@ -441,4 +439,46 @@ public class HaplotypeScore extends InfoFieldAnnotation implements StandardAnnot
public List getDescriptions() {
return Arrays.asList(new VCFInfoHeaderLine("HaplotypeScore", 1, VCFHeaderLineType.Float, "Consistency of the site with at most two segregating haplotypes"));
}
+
+ private static class Haplotype {
+ private final byte[] bases;
+ private final int[] quals;
+ private int qualitySum = -1;
+
+ public Haplotype( final byte[] bases, final int[] quals ) {
+ this.bases = bases;
+ this.quals = quals;
+ }
+
+ public Haplotype( final byte[] bases, final int qual ) {
+ this.bases = bases;
+ quals = new int[bases.length];
+ Arrays.fill(quals, qual);
+ }
+
+ public Haplotype( final byte[] bases, final byte[] quals ) {
+ this.bases = bases;
+ this.quals = new int[quals.length];
+ for ( int i = 0 ; i < quals.length; i++ )
+ this.quals[i] = (int)quals[i];
+ }
+
+ public double getQualitySum() {
+ if ( qualitySum == -1 ) {
+ qualitySum = 0;
+ for ( final int qual : quals ) {
+ qualitySum += qual;
+ }
+ }
+ return qualitySum;
+ }
+
+ public int[] getQuals() {
+ return quals.clone();
+ }
+
+ public byte[] getBases() {
+ return bases.clone();
+ }
+ }
}
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/RankSumTest.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/RankSumTest.java
index 959a26fba..ec107512a 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/RankSumTest.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/RankSumTest.java
@@ -169,8 +169,7 @@ public abstract class RankSumTest extends InfoFieldAnnotation implements ActiveR
* @return true if this base is part of a meaningful read for comparison, false otherwise
*/
public static boolean isUsableBase(final PileupElement p, final boolean allowDeletions) {
- return !(p.isInsertionAtBeginningOfRead() ||
- (! allowDeletions && p.isDeletion()) ||
+ return !((! allowDeletions && p.isDeletion()) ||
p.getMappingQual() == 0 ||
p.getMappingQual() == QualityUtils.MAPPING_QUALITY_UNAVAILABLE ||
((int) p.getQual()) < QualityUtils.MIN_USABLE_Q_SCORE); // need the unBAQed quality score here
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/VariantAnnotatorEngine.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/VariantAnnotatorEngine.java
index 99dadea54..f03a25c04 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/VariantAnnotatorEngine.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/annotator/VariantAnnotatorEngine.java
@@ -52,6 +52,7 @@ import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.walkers.annotator.interfaces.*;
+import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap;
import org.broadinstitute.variant.vcf.*;
import org.broadinstitute.sting.utils.exceptions.UserException;
@@ -214,10 +215,10 @@ public class VariantAnnotatorEngine {
Map infoAnnotations = new LinkedHashMap(vc.getAttributes());
// annotate db occurrences
- vc = annotateDBs(tracker, ref, vc, infoAnnotations);
+ vc = annotateDBs(tracker, ref.getLocus(), vc, infoAnnotations);
// annotate expressions where available
- annotateExpressions(tracker, ref, infoAnnotations);
+ annotateExpressions(tracker, ref.getLocus(), infoAnnotations);
// go through all the requested info annotationTypes
for ( InfoFieldAnnotation annotationType : requestedInfoAnnotations ) {
@@ -254,10 +255,22 @@ public class VariantAnnotatorEngine {
return builder.genotypes(annotateGenotypes(null, null, null, vc, perReadAlleleLikelihoodMap)).make();
}
- private VariantContext annotateDBs(RefMetaDataTracker tracker, ReferenceContext ref, VariantContext vc, Map infoAnnotations) {
+ public VariantContext annotateDBs(final RefMetaDataTracker tracker, final GenomeLoc loc, VariantContext vc) {
+ final Map newInfoAnnotations = new HashMap(0);
+ vc = annotateDBs(tracker, loc, vc, newInfoAnnotations);
+
+ if ( !newInfoAnnotations.isEmpty() ) {
+ final VariantContextBuilder builder = new VariantContextBuilder(vc).attributes(newInfoAnnotations);
+ vc = builder.make();
+ }
+
+ return vc;
+ }
+
+ private VariantContext annotateDBs(final RefMetaDataTracker tracker, final GenomeLoc loc, VariantContext vc, final Map infoAnnotations) {
for ( Map.Entry, String> dbSet : dbAnnotations.entrySet() ) {
if ( dbSet.getValue().equals(VCFConstants.DBSNP_KEY) ) {
- final String rsID = VCFUtils.rsIDOfFirstRealVariant(tracker.getValues(dbSet.getKey(), ref.getLocus()), vc.getType());
+ final String rsID = VCFUtils.rsIDOfFirstRealVariant(tracker.getValues(dbSet.getKey(), loc), vc.getType());
// add the ID if appropriate
if ( rsID != null ) {
@@ -273,7 +286,7 @@ public class VariantAnnotatorEngine {
}
} else {
boolean overlapsComp = false;
- for ( VariantContext comp : tracker.getValues(dbSet.getKey(), ref.getLocus()) ) {
+ for ( VariantContext comp : tracker.getValues(dbSet.getKey(), loc) ) {
if ( !comp.isFiltered() && ( !requireStrictAlleleMatch || comp.getAlleles().equals(vc.getAlleles()) ) ) {
overlapsComp = true;
break;
@@ -287,9 +300,9 @@ public class VariantAnnotatorEngine {
return vc;
}
- private void annotateExpressions(RefMetaDataTracker tracker, ReferenceContext ref, Map infoAnnotations) {
+ private void annotateExpressions(final RefMetaDataTracker tracker, final GenomeLoc loc, final Map infoAnnotations) {
for ( VAExpression expression : requestedExpressions ) {
- Collection VCs = tracker.getValues(expression.binding, ref.getLocus());
+ Collection VCs = tracker.getValues(expression.binding, loc);
if ( VCs.size() == 0 )
continue;
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/coverage/DepthOfCoverage.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/coverage/DepthOfCoverage.java
index 1e4c55e0d..b10daab58 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/coverage/DepthOfCoverage.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/coverage/DepthOfCoverage.java
@@ -938,7 +938,7 @@ public class DepthOfCoverage extends LocusWalker contexts,
AlignmentContextUtils.ReadOrientation contextType) {
final Map consensusIndelStrings = countConsensusAlleles(ref, contexts, contextType);
-// logger.info("Alleles at " + ref.getLocus());
-// for ( Map.Entry elt : consensusIndelStrings.entrySet() ) {
-// logger.info(" " + elt.getValue() + " => " + elt.getKey());
-// }
return consensusCountsToAlleles(ref, consensusIndelStrings);
}
@@ -138,14 +134,9 @@ public class ConsensusAlleleCounter {
final int nReadsOverall = indelPileup.getNumberOfElements();
if ( nIndelReads == 0 || (nIndelReads / (1.0 * nReadsOverall)) < minFractionInOneSample ) {
-// if ( nIndelReads > 0 )
-// logger.info("Skipping sample " + sample.getKey() + " with nIndelReads " + nIndelReads + " nReads " + nReadsOverall);
continue;
-// } else {
-// logger.info("### Keeping sample " + sample.getKey() + " with nIndelReads " + nIndelReads + " nReads " + nReadsOverall);
}
-
for (PileupElement p : indelPileup) {
final GATKSAMRecord read = ReadClipper.hardClipAdaptorSequence(p.getRead());
if (read == null)
@@ -154,17 +145,10 @@ public class ConsensusAlleleCounter {
continue;
}
-/* if (DEBUG && p.isIndel()) {
- System.out.format("Read: %s, cigar: %s, aln start: %d, aln end: %d, p.len:%d, Type:%s, EventBases:%s\n",
- read.getReadName(),read.getCigar().toString(),read.getAlignmentStart(),read.getAlignmentEnd(),
- p.getEventLength(),p.getType().toString(), p.getEventBases());
- }
- */
- String indelString = p.getEventBases();
-
if ( p.isBeforeInsertion() ) {
- // edge case: ignore a deletion immediately preceding an insertion as p.getEventBases() returns null [EB]
- if ( indelString == null )
+ final String insertionBases = p.getBasesOfImmediatelyFollowingInsertion();
+ // edge case: ignore a deletion immediately preceding an insertion as p.getBasesOfImmediatelyFollowingInsertion() returns null [EB]
+ if ( insertionBases == null )
continue;
boolean foundKey = false;
@@ -182,20 +166,20 @@ public class ConsensusAlleleCounter {
String s = cList.get(k).getFirst();
int cnt = cList.get(k).getSecond();
// case 1: current insertion is prefix of indel in hash map
- if (s.startsWith(indelString)) {
+ if (s.startsWith(insertionBases)) {
cList.set(k,new Pair(s,cnt+1));
foundKey = true;
}
- else if (indelString.startsWith(s)) {
+ else if (insertionBases.startsWith(s)) {
// case 2: indel stored in hash table is prefix of current insertion
// In this case, new bases are new key.
foundKey = true;
- cList.set(k,new Pair(indelString,cnt+1));
+ cList.set(k,new Pair(insertionBases,cnt+1));
}
}
if (!foundKey)
// none of the above: event bases not supported by previous table, so add new key
- cList.add(new Pair(indelString,1));
+ cList.add(new Pair(insertionBases,1));
}
else if (read.getAlignmentStart() == loc.getStart()+1) {
@@ -203,28 +187,28 @@ public class ConsensusAlleleCounter {
for (int k=0; k < cList.size(); k++) {
String s = cList.get(k).getFirst();
int cnt = cList.get(k).getSecond();
- if (s.endsWith(indelString)) {
+ if (s.endsWith(insertionBases)) {
// case 1: current insertion (indelString) is suffix of indel in hash map (s)
cList.set(k,new Pair(s,cnt+1));
foundKey = true;
}
- else if (indelString.endsWith(s)) {
+ else if (insertionBases.endsWith(s)) {
// case 2: indel stored in hash table is prefix of current insertion
// In this case, new bases are new key.
foundKey = true;
- cList.set(k,new Pair(indelString,cnt+1));
+ cList.set(k,new Pair(insertionBases,cnt+1));
}
}
if (!foundKey)
// none of the above: event bases not supported by previous table, so add new key
- cList.add(new Pair(indelString,1));
+ cList.add(new Pair(insertionBases,1));
}
else {
// normal case: insertion somewhere in the middle of a read: add count to arrayList
- int cnt = consensusIndelStrings.containsKey(indelString)? consensusIndelStrings.get(indelString):0;
- cList.add(new Pair(indelString,cnt+1));
+ int cnt = consensusIndelStrings.containsKey(insertionBases)? consensusIndelStrings.get(insertionBases):0;
+ cList.add(new Pair(insertionBases,cnt+1));
}
// copy back arrayList into hashMap
@@ -234,11 +218,10 @@ public class ConsensusAlleleCounter {
}
}
- else if ( p.isBeforeDeletedBase() ) {
- indelString = String.format("D%d",p.getEventLength());
- int cnt = consensusIndelStrings.containsKey(indelString)? consensusIndelStrings.get(indelString):0;
- consensusIndelStrings.put(indelString,cnt+1);
-
+ else if ( p.isBeforeDeletionStart() ) {
+ final String deletionString = String.format("D%d",p.getLengthOfImmediatelyFollowingIndel());
+ int cnt = consensusIndelStrings.containsKey(deletionString)? consensusIndelStrings.get(deletionString):0;
+ consensusIndelStrings.put(deletionString,cnt+1);
}
}
}
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidSNPGenotypeLikelihoods.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidSNPGenotypeLikelihoods.java
index fc7573f21..8c8de2bad 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidSNPGenotypeLikelihoods.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidSNPGenotypeLikelihoods.java
@@ -425,7 +425,7 @@ public class DiploidSNPGenotypeLikelihoods implements Cloneable {
byte qual = p.getQual();
if ( qual > SAMUtils.MAX_PHRED_SCORE )
- throw new UserException.MalformedBAM(p.getRead(), String.format("the maximum allowed quality score is %d, but a quality of %d was observed in read %s. Perhaps your BAM incorrectly encodes the quality scores in Sanger format; see http://en.wikipedia.org/wiki/FASTQ_format for more details", SAMUtils.MAX_PHRED_SCORE, qual, p.getRead().getReadName()));
+ throw new UserException.MisencodedBAM(p.getRead(), "we encountered an extremely high quality score (" + (int)qual + ")");
if ( capBaseQualsAtMappingQual )
qual = (byte)Math.min((int)qual, p.getMappingQual());
if ( (int)qual < minBaseQual )
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ErrorModel.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ErrorModel.java
index 12af7839a..1b004d889 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ErrorModel.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ErrorModel.java
@@ -214,7 +214,7 @@ public class ErrorModel {
if (DEBUG)
System.out.format("PE: base:%s isNextToDel:%b isNextToIns:%b eventBases:%s eventLength:%d Allele:%s RefAllele:%s\n",
pileupElement.getBase(), pileupElement.isBeforeDeletionStart(),
- pileupElement.isBeforeInsertion(),pileupElement.getEventBases(),pileupElement.getEventLength(), allele.toString(), refAllele.toString());
+ pileupElement.isBeforeInsertion(),pileupElement.getBasesOfImmediatelyFollowingInsertion(),pileupElement.getLengthOfImmediatelyFollowingIndel(), allele.toString(), refAllele.toString());
//pileupElement.
// if test allele is ref, any base mismatch, or any insertion/deletion at start of pileup count as mismatch
@@ -238,11 +238,11 @@ public class ErrorModel {
// for non-ref alleles,
byte[] alleleBases = allele.getBases();
int eventLength = alleleBases.length - refAllele.getBases().length;
- if (eventLength < 0 && pileupElement.isBeforeDeletionStart() && pileupElement.getEventLength() == -eventLength)
+ if (eventLength < 0 && pileupElement.isBeforeDeletionStart() && pileupElement.getLengthOfImmediatelyFollowingIndel() == -eventLength)
return true;
if (eventLength > 0 && pileupElement.isBeforeInsertion() &&
- Arrays.equals(pileupElement.getEventBases().getBytes(),Arrays.copyOfRange(alleleBases,1,alleleBases.length))) // allele contains ref byte, but pileupElement's event bases doesn't
+ Arrays.equals(pileupElement.getBasesOfImmediatelyFollowingInsertion().getBytes(),Arrays.copyOfRange(alleleBases,1,alleleBases.length))) // allele contains ref byte, but pileupElement's event bases doesn't
return true;
return false;
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GeneralPloidyIndelGenotypeLikelihoods.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GeneralPloidyIndelGenotypeLikelihoods.java
index 7bbe470f8..c957bb9db 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GeneralPloidyIndelGenotypeLikelihoods.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GeneralPloidyIndelGenotypeLikelihoods.java
@@ -210,7 +210,7 @@ public class GeneralPloidyIndelGenotypeLikelihoods extends GeneralPloidyGenotype
// count number of elements in pileup
for (PileupElement elt : pileup) {
if (VERBOSE)
- System.out.format("base:%s isNextToDel:%b isNextToIns:%b eventBases:%s eventLength:%d\n",elt.getBase(), elt.isBeforeDeletionStart(),elt.isBeforeInsertion(),elt.getEventBases(),elt.getEventLength());
+ System.out.format("base:%s isNextToDel:%b isNextToIns:%b eventBases:%s eventLength:%d\n",elt.getBase(), elt.isBeforeDeletionStart(),elt.isBeforeInsertion(),elt.getBasesOfImmediatelyFollowingInsertion(),elt.getLengthOfImmediatelyFollowingIndel());
int idx =0;
for (Allele allele : alleles) {
int cnt = numSeenBases.get(idx);
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GeneralPloidySNPGenotypeLikelihoods.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GeneralPloidySNPGenotypeLikelihoods.java
index edae18a16..aa117eb3b 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GeneralPloidySNPGenotypeLikelihoods.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GeneralPloidySNPGenotypeLikelihoods.java
@@ -323,22 +323,12 @@ public class GeneralPloidySNPGenotypeLikelihoods extends GeneralPloidyGenotypeLi
public ReadBackedPileup createBAQedPileup( final ReadBackedPileup pileup ) {
final List BAQedElements = new ArrayList();
for( final PileupElement PE : pileup ) {
- final PileupElement newPE = new BAQedPileupElement( PE );
+ final PileupElement newPE = new SNPGenotypeLikelihoodsCalculationModel.BAQedPileupElement( PE );
BAQedElements.add( newPE );
}
return new ReadBackedPileupImpl( pileup.getLocation(), BAQedElements );
}
- public class BAQedPileupElement extends PileupElement {
- public BAQedPileupElement( final PileupElement PE ) {
- super(PE.getRead(), PE.getOffset(), PE.isDeletion(), PE.isBeforeDeletedBase(), PE.isAfterDeletedBase(), PE.isBeforeInsertion(), PE.isAfterInsertion(), PE.isNextToSoftClip());
- }
-
- @Override
- public byte getQual( final int offset ) { return BAQ.calcBAQFromTag(getRead(), offset, true); }
- }
-
-
/**
* Helper function that returns the phred-scaled base quality score we should use for calculating
* likelihoods for a pileup element. May return 0 to indicate that the observation is bad, and may
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java
index 86000f236..84c109c9d 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java
@@ -252,7 +252,7 @@ public class IndelGenotypeLikelihoodsCalculationModel extends GenotypeLikelihood
protected int getFilteredDepth(ReadBackedPileup pileup) {
int count = 0;
for (PileupElement p : pileup) {
- if (p.isDeletion() || p.isInsertionAtBeginningOfRead() || BaseUtils.isRegularBase(p.getBase()))
+ if (p.isDeletion() || BaseUtils.isRegularBase(p.getBase()))
count += p.getRepresentativeCount();
}
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/SNPGenotypeLikelihoodsCalculationModel.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/SNPGenotypeLikelihoodsCalculationModel.java
index c1b790559..7dc3e8ee3 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/SNPGenotypeLikelihoodsCalculationModel.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/SNPGenotypeLikelihoodsCalculationModel.java
@@ -237,11 +237,16 @@ public class SNPGenotypeLikelihoodsCalculationModel extends GenotypeLikelihoodsC
public static class BAQedPileupElement extends PileupElement {
public BAQedPileupElement( final PileupElement PE ) {
- super(PE.getRead(), PE.getOffset(), PE.isDeletion(), PE.isBeforeDeletedBase(), PE.isAfterDeletedBase(), PE.isBeforeInsertion(), PE.isAfterInsertion(), PE.isNextToSoftClip());
+ super(PE);
}
@Override
- public byte getQual( final int offset ) { return BAQ.calcBAQFromTag(getRead(), offset, true); }
+ public byte getQual() {
+ if ( isDeletion() )
+ return super.getQual();
+ else
+ return BAQ.calcBAQFromTag(getRead(), offset, true);
+ }
}
private static class SampleGenotypeData {
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java
index 96f327631..26f2560b7 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java
@@ -47,6 +47,7 @@
package org.broadinstitute.sting.gatk.walkers.haplotypecaller;
import com.google.java.contract.Ensures;
+import net.sf.samtools.*;
import org.broadinstitute.sting.commandline.*;
import org.broadinstitute.sting.gatk.CommandLineGATK;
import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
@@ -55,7 +56,9 @@ import org.broadinstitute.sting.gatk.arguments.StandardCallerArgumentCollection;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.AlignmentContextUtils;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
+import org.broadinstitute.sting.gatk.downsampling.DownsampleType;
import org.broadinstitute.sting.gatk.filters.BadMateFilter;
+import org.broadinstitute.sting.gatk.io.StingSAMFileWriter;
import org.broadinstitute.sting.gatk.iterators.ReadTransformer;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.walkers.*;
@@ -66,6 +69,7 @@ import org.broadinstitute.sting.gatk.walkers.genotyper.UnifiedArgumentCollection
import org.broadinstitute.sting.gatk.walkers.genotyper.UnifiedGenotyperEngine;
import org.broadinstitute.sting.gatk.walkers.genotyper.VariantCallContext;
import org.broadinstitute.sting.utils.*;
+import org.broadinstitute.sting.utils.activeregion.ActiveRegion;
import org.broadinstitute.sting.utils.activeregion.ActiveRegionReadState;
import org.broadinstitute.sting.utils.activeregion.ActivityProfileResult;
import org.broadinstitute.sting.utils.clipping.ReadClipper;
@@ -129,6 +133,7 @@ import java.util.*;
@PartitionBy(PartitionType.LOCUS)
@BAQMode(ApplicationTime = ReadTransformer.ApplicationTime.FORBIDDEN)
@ActiveRegionExtension(extension=65, maxRegion=300)
+//@Downsample(by= DownsampleType.BY_SAMPLE, toCoverage=5)
public class HaplotypeCaller extends ActiveRegionWalker implements AnnotatorCompatible {
/**
@@ -140,6 +145,17 @@ public class HaplotypeCaller extends ActiveRegionWalker implem
@Output(fullName="graphOutput", shortName="graph", doc="File to which debug assembly graph information should be written", required = false)
protected PrintStream graphWriter = null;
+ /**
+ * The assembled haplotypes will be written as BAM to this file if requested. Really for debugging purposes only. Note that the output here
+ * does not include uninformative reads so that not every input read is emitted to the bam.
+ */
+ @Hidden
+ @Output(fullName="bamOutput", shortName="bam", doc="File to which assembled haplotypes should be written", required = false)
+ protected StingSAMFileWriter bamWriter = null;
+ private SAMFileHeader bamHeader = null;
+ private long uniqueNameCounter = 1;
+ private final static String readGroupId = "ArtificialHaplotype";
+
/**
* The PairHMM implementation to use for genotype likelihood calculations. The various implementations balance a tradeoff of accuracy and runtime.
*/
@@ -168,6 +184,16 @@ public class HaplotypeCaller extends ActiveRegionWalker implem
@Argument(fullName="downsampleRegion", shortName="dr", doc="coverage, per-sample, to downsample each active region to", required = false)
protected int DOWNSAMPLE_PER_SAMPLE_PER_REGION = 1000;
+ /**
+ * If this flag is provided, the haplotype caller will include unmapped reads in the assembly and calling
+ * when these reads occur in the region being analyzed. Typically, for paired end analyses, one pair of the
+ * read can map, but if its pair is too divergent then it may be unmapped and placed next to its mate, taking
+ * the mates contig and alignment start. If this flag is provided the haplotype caller will see such reads,
+ * and may make use of them in assembly and calling, where possible.
+ */
+ @Argument(fullName="includeUmappedReads", shortName="unmapped", doc="If provided, unmapped reads with chromosomal coordinates (i.e., those placed to their maps) will be included in the assembly and calling", required = false)
+ protected boolean includeUnmappedReads = false;
+
@Argument(fullName="useAllelesTrigger", shortName="allelesTrigger", doc = "If specified, use additional trigger on variants found in an external alleles file", required=false)
protected boolean USE_ALLELES_TRIGGER = false;
@@ -175,6 +201,10 @@ public class HaplotypeCaller extends ActiveRegionWalker implem
@Argument(fullName="useFilteredReadsForAnnotations", shortName="useFilteredReadsForAnnotations", doc = "If specified, use the contamination-filtered read maps for the purposes of annotating variants", required=false)
protected boolean USE_FILTERED_READ_MAP_FOR_ANNOTATIONS = false;
+ @Hidden
+ @Argument(fullName="justDetermineActiveRegions", shortName="justDetermineActiveRegions", doc = "If specified, the HC won't actually do any assembly or calling, it'll just run the upfront active region determination code. Useful for benchmarking and scalability testing", required=false)
+ protected boolean justDetermineActiveRegions = false;
+
/**
* rsIDs from this file are used to populate the ID column of the output. Also, the DB INFO flag will be set when appropriate.
* dbSNP is not used in any way for the calculations themselves.
@@ -236,6 +266,8 @@ public class HaplotypeCaller extends ActiveRegionWalker implem
// the genotyping engine
private GenotypingEngine genotypingEngine = null;
+ private VariantAnnotatorEngine annotationEngine = null;
+
// fasta reference reader to supplement the edges of the reference sequence
private CachingIndexedFastaSequenceFile referenceReader;
@@ -280,7 +312,7 @@ public class HaplotypeCaller extends ActiveRegionWalker implem
UG_engine_simple_genotyper = new UnifiedGenotyperEngine(getToolkit(), simpleUAC, logger, null, null, samples, GATKVariantContextUtils.DEFAULT_PLOIDY);
// initialize the output VCF header
- final VariantAnnotatorEngine annotationEngine = new VariantAnnotatorEngine(Arrays.asList(annotationClassesToUse), annotationsToUse, annotationsToExclude, this, getToolkit());
+ annotationEngine = new VariantAnnotatorEngine(Arrays.asList(annotationClassesToUse), annotationsToUse, annotationsToExclude, this, getToolkit());
Set headerInfo = new HashSet();
@@ -314,6 +346,9 @@ public class HaplotypeCaller extends ActiveRegionWalker implem
assemblyEngine = new SimpleDeBruijnAssembler( DEBUG, graphWriter, minKmer );
likelihoodCalculationEngine = new LikelihoodCalculationEngine( (byte)gcpHMM, DEBUG, pairHMM );
genotypingEngine = new GenotypingEngine( DEBUG, annotationEngine, USE_FILTERED_READ_MAP_FOR_ANNOTATIONS );
+
+ if ( bamWriter != null )
+ setupBamWriter();
}
//---------------------------------------------------------------------------------------------------------------
@@ -329,11 +364,20 @@ public class HaplotypeCaller extends ActiveRegionWalker implem
// enable non primary and extended reads in the active region
@Override
public EnumSet desiredReadStates() {
- return EnumSet.of(
- ActiveRegionReadState.PRIMARY,
- ActiveRegionReadState.NONPRIMARY,
- ActiveRegionReadState.EXTENDED
- );
+ if ( includeUnmappedReads ) {
+ throw new UserException.BadArgumentValue("includeUmappedReads", "is not yet functional");
+// return EnumSet.of(
+// ActiveRegionReadState.PRIMARY,
+// ActiveRegionReadState.NONPRIMARY,
+// ActiveRegionReadState.EXTENDED,
+// ActiveRegionReadState.UNMAPPED
+// );
+ } else
+ return EnumSet.of(
+ ActiveRegionReadState.PRIMARY,
+ ActiveRegionReadState.NONPRIMARY,
+ ActiveRegionReadState.EXTENDED
+ );
}
@Override
@@ -371,7 +415,7 @@ public class HaplotypeCaller extends ActiveRegionWalker implem
final byte qual = p.getQual();
if( p.isDeletion() || qual > (byte) 18) {
int AA = 0; final int AB = 1; int BB = 2;
- if( p.getBase() != ref.getBase() || p.isDeletion() || p.isBeforeDeletedBase() || p.isAfterDeletedBase() || p.isBeforeInsertion() || p.isAfterInsertion() || p.isNextToSoftClip() ) {
+ if( p.getBase() != ref.getBase() || p.isDeletion() || p.isBeforeDeletionStart() || p.isAfterDeletionEnd() || p.isBeforeInsertion() || p.isAfterInsertion() || p.isNextToSoftClip() ) {
AA = 2;
BB = 0;
if( p.isNextToSoftClip() ) {
@@ -402,7 +446,10 @@ public class HaplotypeCaller extends ActiveRegionWalker implem
//---------------------------------------------------------------------------------------------------------------
@Override
- public Integer map( final org.broadinstitute.sting.utils.activeregion.ActiveRegion activeRegion, final RefMetaDataTracker metaDataTracker ) {
+ public Integer map( final ActiveRegion activeRegion, final RefMetaDataTracker metaDataTracker ) {
+ if ( justDetermineActiveRegions )
+ // we're benchmarking ART and/or the active region determination code in the HC, just leave without doing any work
+ return 1;
final ArrayList activeAllelesToGenotype = new ArrayList();
@@ -420,8 +467,7 @@ public class HaplotypeCaller extends ActiveRegionWalker implem
if( UG_engine.getUAC().GenotypingMode == GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES && activeAllelesToGenotype.isEmpty() ) { return 0; } // No alleles found in this region so nothing to do!
finalizeActiveRegion( activeRegion ); // merge overlapping fragments, clip adapter and low qual tails
- final Haplotype referenceHaplotype = new Haplotype(activeRegion.getActiveRegionReference(referenceReader)); // Create the reference haplotype which is the bases from the reference that make up the active region
- referenceHaplotype.setIsReference(true);
+ final Haplotype referenceHaplotype = new Haplotype(activeRegion.getActiveRegionReference(referenceReader), true); // Create the reference haplotype which is the bases from the reference that make up the active region
final byte[] fullReferenceWithPadding = activeRegion.getFullReference(referenceReader, REFERENCE_PADDING);
//int PRUNE_FACTOR = Math.max(MIN_PRUNE_FACTOR, determinePruneFactorFromCoverage( activeRegion ));
final ArrayList haplotypes = assemblyEngine.runLocalAssembly( activeRegion, referenceHaplotype, fullReferenceWithPadding, getPaddedLoc(activeRegion), MIN_PRUNE_FACTOR, activeAllelesToGenotype );
@@ -452,9 +498,31 @@ public class HaplotypeCaller extends ActiveRegionWalker implem
activeRegion.getLocation(),
getToolkit().getGenomeLocParser(),
activeAllelesToGenotype ) ) {
+ annotationEngine.annotateDBs(metaDataTracker, getToolkit().getGenomeLocParser().createGenomeLoc(call), call);
vcfWriter.add( call );
}
+ if ( bamWriter != null ) {
+ // write the haplotypes to the bam
+ final GenomeLoc paddedRefLoc = getPaddedLoc(activeRegion);
+ for ( Haplotype haplotype : haplotypes )
+ writeHaplotype(haplotype, paddedRefLoc, bestHaplotypes.contains(haplotype));
+
+ // we need to remap the Alleles back to the Haplotypes; inefficient but unfortunately this is a requirement currently
+ final Map alleleToHaplotypeMap = new HashMap(haplotypes.size());
+ for ( final Haplotype haplotype : haplotypes )
+ alleleToHaplotypeMap.put(Allele.create(haplotype.getBases()), haplotype);
+
+ // next, output the interesting reads for each sample aligned against the appropriate haplotype
+ for ( final PerReadAlleleLikelihoodMap readAlleleLikelihoodMap : stratifiedReadMap.values() ) {
+ for ( Map.Entry> entry : readAlleleLikelihoodMap.getLikelihoodReadMap().entrySet() ) {
+ final Allele bestAllele = PerReadAlleleLikelihoodMap.getMostLikelyAllele(entry.getValue());
+ if ( bestAllele != Allele.NO_CALL )
+ writeReadAgainstHaplotype(entry.getKey(), alleleToHaplotypeMap.get(bestAllele), paddedRefLoc.getStart());
+ }
+ }
+ }
+
if( DEBUG ) { System.out.println("----------------------------------------------------------------------------------"); }
return 1; // One active region was processed during this map call
@@ -548,6 +616,95 @@ public class HaplotypeCaller extends ActiveRegionWalker implem
return returnMap;
}
+ private void setupBamWriter() {
+ // prepare the bam header
+ bamHeader = new SAMFileHeader();
+ bamHeader.setSequenceDictionary(getToolkit().getSAMFileHeader().getSequenceDictionary());
+ bamHeader.setSortOrder(SAMFileHeader.SortOrder.coordinate);
+
+ // include the original read groups plus a new artificial one for the haplotypes
+ final List readGroups = new ArrayList(getToolkit().getSAMFileHeader().getReadGroups());
+ final SAMReadGroupRecord rg = new SAMReadGroupRecord(readGroupId);
+ rg.setSample("HC");
+ rg.setSequencingCenter("BI");
+ readGroups.add(rg);
+ bamHeader.setReadGroups(readGroups);
+
+ bamWriter.setPresorted(false);
+ bamWriter.writeHeader(bamHeader);
+ }
+
+ private void writeHaplotype(final Haplotype haplotype, final GenomeLoc paddedRefLoc, final boolean isAmongBestHaplotypes) {
+ final GATKSAMRecord record = new GATKSAMRecord(bamHeader);
+ record.setReadBases(haplotype.getBases());
+ record.setAlignmentStart(paddedRefLoc.getStart() + haplotype.getAlignmentStartHapwrtRef());
+ record.setBaseQualities(Utils.dupBytes((byte) '!', haplotype.getBases().length));
+ record.setCigar(haplotype.getCigar());
+ record.setMappingQuality(isAmongBestHaplotypes ? 60 : 0);
+ record.setReadName("HC" + uniqueNameCounter++);
+ record.setReadUnmappedFlag(false);
+ record.setReferenceIndex(paddedRefLoc.getContigIndex());
+ record.setAttribute(SAMTag.RG.toString(), readGroupId);
+ record.setFlags(16);
+ bamWriter.addAlignment(record);
+ }
+
+ private void writeReadAgainstHaplotype(final GATKSAMRecord read, final Haplotype haplotype, final int referenceStart) {
+
+ final SWPairwiseAlignment swPairwiseAlignment = new SWPairwiseAlignment(haplotype.getBases(), read.getReadBases(), 5.0, -10.0, -22.0, -1.2);
+ final int readStartOnHaplotype = swPairwiseAlignment.getAlignmentStart2wrt1();
+ final int readStartOnReference = referenceStart + haplotype.getAlignmentStartHapwrtRef() + readStartOnHaplotype;
+ read.setAlignmentStart(readStartOnReference);
+
+ final Cigar cigar = generateReadCigarFromHaplotype(read, readStartOnHaplotype, haplotype.getCigar());
+ read.setCigar(cigar);
+
+ bamWriter.addAlignment(read);
+ }
+
+ private Cigar generateReadCigarFromHaplotype(final GATKSAMRecord read, final int readStartOnHaplotype, final Cigar haplotypeCigar) {
+
+ int currentReadPos = 0;
+ int currentHapPos = 0;
+ final List readCigarElements = new ArrayList();
+
+ for ( final CigarElement cigarElement : haplotypeCigar.getCigarElements() ) {
+
+ if ( cigarElement.getOperator() == CigarOperator.D ) {
+ if ( currentReadPos > 0 )
+ readCigarElements.add(cigarElement);
+ } else if ( cigarElement.getOperator() == CigarOperator.M || cigarElement.getOperator() == CigarOperator.I ) {
+
+ final int elementLength = cigarElement.getLength();
+ final int nextReadPos = currentReadPos + elementLength;
+ final int nextHapPos = currentHapPos + elementLength;
+
+ // do we want this element?
+ if ( currentReadPos > 0 ) {
+ // do we want the entire element?
+ if ( nextReadPos < read.getReadLength() ) {
+ readCigarElements.add(cigarElement);
+ currentReadPos = nextReadPos;
+ }
+ // otherwise, we can finish up and return the cigar
+ else {
+ readCigarElements.add(new CigarElement(read.getReadLength() - currentReadPos, cigarElement.getOperator()));
+ return new Cigar(readCigarElements);
+ }
+ }
+ // do we want part of the element to start?
+ else if ( currentReadPos == 0 && nextHapPos > readStartOnHaplotype ) {
+ currentReadPos = Math.min(nextHapPos - readStartOnHaplotype, read.getReadLength());
+ readCigarElements.add(new CigarElement(currentReadPos, cigarElement.getOperator()));
+ }
+
+ currentHapPos = nextHapPos;
+ }
+ }
+
+ return new Cigar(readCigarElements);
+ }
+
/*
private int determinePruneFactorFromCoverage( final ActiveRegion activeRegion ) {
final ArrayList readLengthDistribution = new ArrayList();
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java
index 8b844817d..57e071189 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java
@@ -124,9 +124,14 @@ public class LikelihoodCalculationEngine {
}
private PerReadAlleleLikelihoodMap computeReadLikelihoods( final ArrayList haplotypes, final ArrayList reads) {
+ // first, a little set up to get copies of the Haplotypes that are Alleles (more efficient than creating them each time)
+ final int numHaplotypes = haplotypes.size();
+ final Map alleleVersions = new HashMap(numHaplotypes);
+ for ( final Haplotype haplotype : haplotypes ) {
+ alleleVersions.put(haplotype, Allele.create(haplotype.getBases()));
+ }
final PerReadAlleleLikelihoodMap perReadAlleleLikelihoodMap = new PerReadAlleleLikelihoodMap();
- final int numHaplotypes = haplotypes.size();
for( final GATKSAMRecord read : reads ) {
final byte[] overallGCP = new byte[read.getReadLength()];
Arrays.fill( overallGCP, constantGCP ); // Is there a way to derive empirical estimates for this from the data?
@@ -138,20 +143,17 @@ public class LikelihoodCalculationEngine {
readQuals[kkk] = ( readQuals[kkk] > (byte) read.getMappingQuality() ? (byte) read.getMappingQuality() : readQuals[kkk] ); // cap base quality by mapping quality
//readQuals[kkk] = ( readQuals[kkk] > readInsQuals[kkk] ? readInsQuals[kkk] : readQuals[kkk] ); // cap base quality by base insertion quality, needs to be evaluated
//readQuals[kkk] = ( readQuals[kkk] > readDelQuals[kkk] ? readDelQuals[kkk] : readQuals[kkk] ); // cap base quality by base deletion quality, needs to be evaluated
+ // TODO -- why is Q18 hard-coded here???
readQuals[kkk] = ( readQuals[kkk] < (byte) 18 ? QualityUtils.MIN_USABLE_Q_SCORE : readQuals[kkk] );
}
for( int jjj = 0; jjj < numHaplotypes; jjj++ ) {
final Haplotype haplotype = haplotypes.get(jjj);
- // TODO -- need to test against a reference/position with non-standard bases
- //if ( !Allele.acceptableAlleleBases(haplotype.getBases(), false) )
- // continue;
-
final int haplotypeStart = ( previousHaplotypeSeen == null ? 0 : computeFirstDifferingPosition(haplotype.getBases(), previousHaplotypeSeen.getBases()) );
previousHaplotypeSeen = haplotype;
- perReadAlleleLikelihoodMap.add(read, Allele.create(haplotype.getBases()),
+ perReadAlleleLikelihoodMap.add(read, alleleVersions.get(haplotype),
pairHMM.computeReadLikelihoodGivenHaplotypeLog10(haplotype.getBases(), read.getReadBases(),
readQuals, readInsQuals, readDelQuals, overallGCP, haplotypeStart, jjj == 0));
}
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/SimpleDeBruijnAssembler.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/SimpleDeBruijnAssembler.java
index e1a94eee7..e16994fa4 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/SimpleDeBruijnAssembler.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/SimpleDeBruijnAssembler.java
@@ -338,7 +338,7 @@ public class SimpleDeBruijnAssembler extends LocalAssemblyEngine {
for( final DefaultDirectedGraph graph : graphs ) {
for ( final KBestPaths.Path path : KBestPaths.getKBestPaths(graph, NUM_BEST_PATHS_PER_KMER_GRAPH) ) {
- final Haplotype h = new Haplotype( path.getBases( graph ), path.getScore() );
+ final Haplotype h = new Haplotype( path.getBases( graph ) );
if( addHaplotype( h, fullReferenceWithPadding, returnHaplotypes, activeRegionStart, activeRegionStop, false ) ) {
// for GGA mode, add the desired allele into the haplotype if it isn't already present
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/validation/ValidationAmplicons.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/validation/ValidationAmplicons.java
index 9ba74b98a..8547c0993 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/validation/ValidationAmplicons.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/validation/ValidationAmplicons.java
@@ -376,7 +376,7 @@ public class ValidationAmplicons extends RodWalker {
if ( lowerCaseSNPs ) {
sequence.append(Character.toLowerCase((char) ref.getBase()));
} else {
- sequence.append((char) BaseUtils.N);
+ sequence.append((char) BaseUtils.Base.N.base);
}
rawSequence.append(Character.toUpperCase((char) ref.getBase()));
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/ConcordanceMetrics.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/ConcordanceMetrics.java
index bb76006bf..8a87c9957 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/ConcordanceMetrics.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/ConcordanceMetrics.java
@@ -272,22 +272,7 @@ public class ConcordanceMetrics {
@Requires({"evalVC != null","truthVC != null"})
private SiteConcordanceType getMatchType(VariantContext evalVC, VariantContext truthVC) {
- if ( evalVC.isMonomorphicInSamples() )
- return SiteConcordanceType.TRUTH_ONLY;
- if ( truthVC.isMonomorphicInSamples() )
- return SiteConcordanceType.EVAL_ONLY;
-
- boolean evalSusbsetTruth = VariantContextUtils.allelesAreSubset(evalVC,truthVC);
- boolean truthSubsetEval = VariantContextUtils.allelesAreSubset(truthVC,evalVC);
-
- if ( evalSusbsetTruth && truthSubsetEval )
- return SiteConcordanceType.ALLELES_MATCH;
- else if ( evalSusbsetTruth )
- return SiteConcordanceType.EVAL_SUBSET_TRUTH;
- else if ( truthSubsetEval )
- return SiteConcordanceType.EVAL_SUPERSET_TRUTH;
-
- return SiteConcordanceType.ALLELES_DO_NOT_MATCH;
+ return SiteConcordanceType.getConcordanceType(evalVC,truthVC);
}
public int[] getSiteConcordance() {
@@ -305,6 +290,27 @@ public class ConcordanceMetrics {
EVAL_SUBSET_TRUTH,
ALLELES_DO_NOT_MATCH,
EVAL_ONLY,
- TRUTH_ONLY
+ TRUTH_ONLY;
+
+ public static SiteConcordanceType getConcordanceType(VariantContext eval, VariantContext truth) {
+ if ( eval.isMonomorphicInSamples() )
+ return TRUTH_ONLY;
+ if ( truth.isMonomorphicInSamples() )
+ return EVAL_ONLY;
+
+ boolean evalSubsetTruth = VariantContextUtils.allelesAreSubset(eval,truth);
+ boolean truthSubsetEval = VariantContextUtils.allelesAreSubset(eval,truth);
+
+ if ( evalSubsetTruth && truthSubsetEval )
+ return ALLELES_MATCH;
+
+ if ( evalSubsetTruth )
+ return EVAL_SUBSET_TRUTH;
+
+ if ( truthSubsetEval )
+ return EVAL_SUPERSET_TRUTH;
+
+ return ALLELES_DO_NOT_MATCH;
+ }
}
}
diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/GenotypeConcordance.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/GenotypeConcordance.java
index 64eb4deef..e8965dfc8 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/GenotypeConcordance.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/GenotypeConcordance.java
@@ -64,7 +64,7 @@ import java.util.*;
/**
* A simple walker for performing genotype concordance calculations between two callsets
*/
-public class GenotypeConcordance extends RodWalker,ConcordanceMetrics> {
+public class GenotypeConcordance extends RodWalker>,ConcordanceMetrics> {
@Input(fullName="eval",shortName="eval",doc="The variants and genotypes to evaluate",required=true)
RodBinding evalBinding;
@@ -75,21 +75,41 @@ public class GenotypeConcordance extends RodWalker genotypeFilterExpressionsEval = new ArrayList();
+
+ @Argument(shortName="gfc", fullName="genotypeFilterExpressionComp", doc="One or more criteria to use to set COMP genotypes to no-call. "+
+ "These genotype-level filters are only applied to the COMP rod.", required=false)
+ public ArrayList genotypeFilterExpressionsComp = new ArrayList();
+
+ @Argument(shortName="moltenize",fullName="moltenize",doc="Molten rather than tabular output")
+ public boolean moltenize = false;
+
@Output
PrintStream out;
- List evalSamples;
- List compSamples;
+ private List evalSamples;
+ private List compSamples;
+ private List evalJexls = null;
+ private List compJexls = null;
- // todo -- integration test coverage
- // todo -- deal with occurrences like:
- // Eval: 20 4000 A C
- // Eval: 20 4000 A AC
- // Comp: 20 4000 A C
- // currently this results in a warning and skipping
- // todo -- extend to multiple eval, multiple comp
- // todo -- table with "proportion of overlapping sites" (not just eval/comp margins)
+ // todo -- table with "proportion of overlapping sites" (not just eval/comp margins) [e.g. drop no-calls]
+ // (this will break all the integration tests of course, due to new formatting)
+ public void initialize() {
+ evalJexls = initializeJexl(genotypeFilterExpressionsEval);
+ compJexls = initializeJexl(genotypeFilterExpressionsComp);
+ }
+
+ private List initializeJexl(ArrayList genotypeFilterExpressions) {
+ ArrayList dummyNames = new ArrayList(genotypeFilterExpressions.size());
+ int expCount = 1;
+ for ( String exp : genotypeFilterExpressions ) {
+ dummyNames.add(String.format("gfe%d",expCount++));
+ }
+ return VariantContextUtils.initializeMatchExps(dummyNames, genotypeFilterExpressions);
+ }
public ConcordanceMetrics reduceInit() {
Map headerMap = GATKVCFUtils.getVCFHeadersFromRods(getToolkit(), Arrays.asList(evalBinding,compBinding));
@@ -101,120 +121,320 @@ public class GenotypeConcordance extends RodWalker map(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
- Pair evalCompPair = null;
+ public List> map(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
+ List> evalCompPair = new ArrayList>(3);
if ( tracker != null && (
tracker.getValues(evalBinding,ref.getLocus()).size() > 0 ||
- tracker.getValues(compBinding,ref.getLocus()).size() > 0 ) ) {
+ tracker.getValues(compBinding,ref.getLocus()).size() > 0 ) ) {
List eval = tracker.getValues(evalBinding,ref.getLocus());
List comp = tracker.getValues(compBinding,ref.getLocus());
if ( eval.size() > 1 || comp.size() > 1 ) {
- logger.warn("Eval or Comp Rod at position "+ref.getLocus().toString()+" has multiple records. Site will be skipped.");
- return evalCompPair;
+ if ( noDuplicateTypes(eval) && noDuplicateTypes(comp) ) {
+ logger.info("Eval or Comp Rod at position " + ref.getLocus().toString() + " has multiple records. Resolving.");
+ evalCompPair = resolveMultipleRecords(eval,comp);
+ } else {
+ logger.warn("Eval or Comp Rod at position "+ref.getLocus().toString()+" has multiple records of the same type. This locus will be skipped.");
+ }
+ } else {
+ // if a rod is missing, explicitly create a variant context with 'missing' genotypes. Slow, but correct.
+ // note that if there is no eval rod there must be a comp rod, and also the reverse
+ VariantContext evalContext = eval.size() == 1 ? eval.get(0) : createEmptyContext(comp.get(0),evalSamples);
+ VariantContext compContext = comp.size() == 1 ? comp.get(0) : createEmptyContext(eval.get(0),compSamples);
+ evalContext = filterGenotypes(evalContext,ignoreFilters,evalJexls);
+ compContext = filterGenotypes(compContext,ignoreFilters,compJexls);
+ evalCompPair.add(new Pair(evalContext,compContext));
}
- // if a rod is missing, explicitly create a variant context with 'missing' genotypes. Slow, but correct.
- // note that if there is no eval rod there must be a comp rod, and also the reverse
- VariantContext evalContext = eval.size() == 1 ? eval.get(0) : createEmptyContext(ref,comp.get(0),evalSamples);
- VariantContext compContext = comp.size() == 1 ? comp.get(0) : createEmptyContext(ref,eval.get(0),compSamples);
- evalContext = filterGenotypes(evalContext,ignoreFilters);
- compContext = filterGenotypes(compContext,ignoreFilters);
- evalCompPair = new Pair(evalContext,compContext);
}
return evalCompPair;
}
- public ConcordanceMetrics reduce(Pair evalComp, ConcordanceMetrics metrics) {
- if ( evalComp != null )
+ private boolean noDuplicateTypes(List vcList) {
+ HashSet types = new HashSet(vcList.size());
+ for ( VariantContext vc : vcList ) {
+ VariantContext.Type type = vc.getType();
+ if ( types.contains(type) )
+ return false;
+ types.add(type);
+ }
+
+ return true;
+ }
+
+ /**
+ * The point of this method is to match up pairs of evals and comps by their type (or alternate alleles for mixed).
+ * Basically multiple records could exist for a site such as:
+ * Eval: 20 4000 A C
+ * Eval: 20 4000 A AC
+ * Comp: 20 4000 A C
+ * So for each eval, loop through the comps. If the eval alleles (non-emptily) intersect the comp alleles, pair them up and remove
+ * that comp records. Continue until we're out of evals or comps. This is n^2, but should rarely actually happen.
+ *
+ * The remaining unpaired records get paird with an empty contexts. So in the example above we'd get a list of:
+ * 1 - (20,4000,A/C | 20,4000,A/C)
+ * 2 - (20,4000,A/AC | Empty )
+ * @param evalList - list of eval variant contexts
+ * @param compList - list of comp variant contexts
+ * @return resolved pairs of the input lists
+ */
+ private List> resolveMultipleRecords(List evalList, List compList) {
+ List> resolvedPairs = new ArrayList>(evalList.size()+compList.size()); // oversized but w/e
+ List pairedEval = new ArrayList(evalList.size());
+ for ( VariantContext eval : evalList ) {
+ VariantContext.Type evalType = eval.getType();
+ Set evalAlleles = new HashSet(eval.getAlternateAlleles());
+ VariantContext pairedComp = null;
+ for ( VariantContext comp : compList ) {
+ if ( evalType.equals(comp.getType()) ) {
+ pairedComp = comp;
+ break;
+ } else if ( eval.isMixed() || comp.isMixed() ) {
+ for ( Allele compAllele : comp.getAlternateAlleles() ) {
+ if ( evalAlleles.contains(compAllele) ) {
+ pairedComp = comp;
+ break;
+ }
+ }
+ }
+ }
+ if ( pairedComp != null ) {
+ compList.remove(pairedComp);
+ resolvedPairs.add(new Pair(filterGenotypes(eval,ignoreFilters,evalJexls),filterGenotypes(pairedComp,ignoreFilters,compJexls)));
+ pairedEval.add(eval);
+ if ( compList.size() < 1 )
+ break;
+ }
+ }
+ evalList.removeAll(pairedEval);
+ for ( VariantContext unpairedEval : evalList ) {
+ resolvedPairs.add(new Pair(filterGenotypes(unpairedEval,ignoreFilters,evalJexls),createEmptyContext(unpairedEval,compSamples)));
+ }
+
+ for ( VariantContext unpairedComp : compList ) {
+ resolvedPairs.add(new Pair(createEmptyContext(unpairedComp,evalSamples),filterGenotypes(unpairedComp,ignoreFilters,compJexls)));
+ }
+
+ return resolvedPairs;
+ }
+
+ public ConcordanceMetrics reduce(List> evalCompList, ConcordanceMetrics metrics) {
+ for ( Pair evalComp : evalCompList)
metrics.update(evalComp.getFirst(),evalComp.getSecond());
return metrics;
}
+ private static double repairNaN(double d) {
+ if ( Double.isNaN(d) ) {
+ return 0.0;
+ }
+ return d;
+ }
+
public void onTraversalDone(ConcordanceMetrics metrics) {
+ // todo -- this is over 200 lines of code just to format the output and could use some serious cleanup
GATKReport report = new GATKReport();
GATKReportTable concordanceCounts = new GATKReportTable("GenotypeConcordance_Counts","Per-sample concordance tables: comparison counts",2+GenotypeType.values().length*GenotypeType.values().length);
GATKReportTable concordanceEvalProportions = new GATKReportTable("GenotypeConcordance_EvalProportions", "Per-sample concordance tables: proportions of genotypes called in eval",2+GenotypeType.values().length*GenotypeType.values().length);
GATKReportTable concordanceCompProportions = new GATKReportTable("GenotypeConcordance_CompProportions", "Per-sample concordance tables: proportions of genotypes called in comp",2+GenotypeType.values().length*GenotypeType.values().length);
GATKReportTable concordanceSummary = new GATKReportTable("GenotypeConcordance_Summary","Per-sample summary statistics: NRS and NRD",2);
GATKReportTable siteConcordance = new GATKReportTable("SiteConcordance_Summary","Site-level summary statistics",ConcordanceMetrics.SiteConcordanceType.values().length);
- concordanceCompProportions.addColumn("Sample","%s");
- concordanceCounts.addColumn("Sample","%s");
- concordanceEvalProportions.addColumn("Sample","%s");
- concordanceSummary.addColumn("Sample","%s");
- for ( GenotypeType evalType : GenotypeType.values() ) {
- for ( GenotypeType compType : GenotypeType.values() ) {
- String colKey = String.format("%s_%s", evalType.toString(), compType.toString());
- concordanceCounts.addColumn(colKey,"%d");
- if ( evalType == GenotypeType.HET || evalType == GenotypeType.HOM_REF || evalType == GenotypeType.HOM_VAR)
- concordanceEvalProportions.addColumn(colKey,"%.3f");
- if ( compType == GenotypeType.HET || compType == GenotypeType.HOM_VAR || compType == GenotypeType.HOM_REF )
- concordanceCompProportions.addColumn(colKey,"%.3f");
- }
- }
- concordanceEvalProportions.addColumn("Mismatching_Alleles","%.3f");
- concordanceCompProportions.addColumn("Mismatching_Alleles","%.3f");
- concordanceCounts.addColumn("Mismatching_Alleles","%d");
- concordanceSummary.addColumn("Non-Reference Sensitivity","%.3f");
- concordanceSummary.addColumn("Non-Reference Discrepancy","%.3f");
- for (ConcordanceMetrics.SiteConcordanceType type : ConcordanceMetrics.SiteConcordanceType.values() ) {
- siteConcordance.addColumn(type.toString(),"%d");
- }
+ if ( moltenize ) {
+ concordanceCompProportions.addColumn("Sample","%s");
+ concordanceCounts.addColumn("Sample","%s");
+ concordanceEvalProportions.addColumn("Sample","%s");
+ concordanceSummary.addColumn("Sample","%s");
- for ( Map.Entry entry : metrics.getPerSampleGenotypeConcordance().entrySet() ) {
- ConcordanceMetrics.GenotypeConcordanceTable table = entry.getValue();
- concordanceEvalProportions.set(entry.getKey(),"Sample",entry.getKey());
- concordanceCompProportions.set(entry.getKey(),"Sample",entry.getKey());
- concordanceCounts.set(entry.getKey(),"Sample",entry.getKey());
+ concordanceCompProportions.addColumn("Eval_Genotype","%s");
+ concordanceCounts.addColumn("Eval_Genotype","%s");
+ concordanceEvalProportions.addColumn("Eval_Genotype","%s");
+ concordanceSummary.addColumn("Non-Reference_Discrepancy","%.3f");
+
+ concordanceCompProportions.addColumn("Comp_Genotype","%s");
+ concordanceCounts.addColumn("Comp_Genotype","%s");
+ concordanceEvalProportions.addColumn("Comp_Genotype","%s");
+ concordanceSummary.addColumn("Non-Reference_Sensitivity","%.3f");
+
+ concordanceCompProportions.addColumn("Proportion","%.3f");
+ concordanceCounts.addColumn("Count","%d");
+ concordanceEvalProportions.addColumn("Proportion","%.3f");
+
+ for ( Map.Entry entry : metrics.getPerSampleGenotypeConcordance().entrySet() ) {
+ ConcordanceMetrics.GenotypeConcordanceTable table = entry.getValue();
+ for ( GenotypeType evalType : GenotypeType.values() ) {
+ for ( GenotypeType compType : GenotypeType.values() ) {
+ String rowKey = String.format("%s_%s_%s",entry.getKey(),evalType.toString(),compType.toString());
+ concordanceCounts.set(rowKey,"Sample",entry.getKey());
+ concordanceCounts.set(rowKey,"Eval_Genotype",evalType.toString());
+ concordanceCounts.set(rowKey,"Comp_Genotype",evalType.toString());
+ int count = table.get(evalType, compType);
+ concordanceCounts.set(rowKey,"Count",count);
+ if ( evalType == GenotypeType.HET || evalType == GenotypeType.HOM_REF || evalType == GenotypeType.HOM_VAR) {
+ concordanceEvalProportions.set(rowKey,"Sample",entry.getKey());
+ concordanceEvalProportions.set(rowKey,"Eval_Genotype",evalType.toString());
+ concordanceEvalProportions.set(rowKey,"Comp_Genotype",evalType.toString());
+ concordanceEvalProportions.set(rowKey,"Proportion",repairNaN(( (double) count)/table.getnEvalGenotypes(evalType)));
+ }
+ if ( compType == GenotypeType.HET || compType == GenotypeType.HOM_VAR || compType == GenotypeType.HOM_REF ) {
+ concordanceCompProportions.set(rowKey,"Sample",entry.getKey());
+ concordanceCompProportions.set(rowKey,"Eval_Genotype",evalType.toString());
+ concordanceCompProportions.set(rowKey,"Comp_Genotype",evalType.toString());
+ concordanceCompProportions.set(rowKey,"Proportion",repairNaN(( (double) count)/table.getnCompGenotypes(compType)));
+ }
+ }
+ }
+ String mismatchKey = String.format("%s_%s",entry.getKey(),"Mismatching");
+ concordanceCounts.set(mismatchKey,"Sample",entry.getKey());
+ concordanceCounts.set(mismatchKey,"Eval_Genotype","Mismatching_Alleles");
+ concordanceCounts.set(mismatchKey,"Comp_Genotype","Mismatching_Alleles");
+ concordanceEvalProportions.set(mismatchKey,"Sample",entry.getKey());
+ concordanceEvalProportions.set(mismatchKey,"Eval_Genotype","Mismatching_Alleles");
+ concordanceEvalProportions.set(mismatchKey,"Comp_Genotype","Mismatching_Alleles");
+ concordanceCompProportions.set(mismatchKey,"Sample",entry.getKey());
+ concordanceCompProportions.set(mismatchKey,"Eval_Genotype","Mismatching_Alleles");
+ concordanceCompProportions.set(mismatchKey,"Comp_Genotype","Mismatching_Alleles");
+ concordanceEvalProportions.set(mismatchKey,"Proportion", repairNaN(( (double) table.getnMismatchingAlt() )/table.getnCalledEvalGenotypes()));
+ concordanceCompProportions.set(mismatchKey,"Proportion", repairNaN(( (double) table.getnMismatchingAlt() )/table.getnCalledCompGenotypes()));
+ concordanceCounts.set(mismatchKey,"Count",table.getnMismatchingAlt());
+ }
+
+ String sampleKey = "ALL";
+ ConcordanceMetrics.GenotypeConcordanceTable table = metrics.getOverallGenotypeConcordance();
+ for ( GenotypeType evalType : GenotypeType.values() ) {
+ for ( GenotypeType compType : GenotypeType.values() ) {
+ String rowKey = String.format("%s_%s_%s",sampleKey,evalType.toString(),compType.toString());
+ concordanceCounts.set(rowKey,"Sample",sampleKey);
+ concordanceCounts.set(rowKey,"Eval_Genotype",evalType.toString());
+ concordanceCounts.set(rowKey,"Comp_Genotype",evalType.toString());
+ int count = table.get(evalType, compType);
+ concordanceCounts.set(rowKey,"Count",count);
+ if ( evalType == GenotypeType.HET || evalType == GenotypeType.HOM_REF || evalType == GenotypeType.HOM_VAR) {
+ concordanceEvalProportions.set(rowKey,"Sample",sampleKey);
+ concordanceEvalProportions.set(rowKey,"Eval_Genotype",evalType.toString());
+ concordanceEvalProportions.set(rowKey,"Comp_Genotype",evalType.toString());
+ concordanceEvalProportions.set(rowKey,"Proportion",repairNaN(( (double) count)/table.getnEvalGenotypes(evalType)));
+ }
+ if ( compType == GenotypeType.HET || compType == GenotypeType.HOM_VAR || compType == GenotypeType.HOM_REF ) {
+ concordanceCompProportions.set(rowKey,"Sample",sampleKey);
+ concordanceCompProportions.set(rowKey,"Eval_Genotype",evalType.toString());
+ concordanceCompProportions.set(rowKey,"Comp_Genotype",evalType.toString());
+ concordanceCompProportions.set(rowKey,"Proportion",repairNaN(( (double) count)/table.getnCompGenotypes(compType)));
+ }
+ }
+ }
+ String rowKey = String.format("%s_%s",sampleKey,"Mismatching");
+ concordanceCounts.set(rowKey,"Sample",sampleKey);
+ concordanceCounts.set(rowKey,"Eval_Genotype","Mismatching_Alleles");
+ concordanceCounts.set(rowKey,"Comp_Genotype","Mismatching_Alleles");
+ concordanceEvalProportions.set(rowKey,"Sample",sampleKey);
+ concordanceEvalProportions.set(rowKey,"Eval_Genotype","Mismatching_Alleles");
+ concordanceEvalProportions.set(rowKey,"Comp_Genotype","Mismatching_Alleles");
+ concordanceCompProportions.set(rowKey,"Sample",sampleKey);
+ concordanceCompProportions.set(rowKey,"Eval_Genotype","Mismatching_Alleles");
+ concordanceCompProportions.set(rowKey,"Comp_Genotype","Mismatching_Alleles");
+ concordanceEvalProportions.set(rowKey,"Proportion", repairNaN(( (double) table.getnMismatchingAlt() )/table.getnCalledEvalGenotypes()));
+ concordanceCompProportions.set(rowKey,"Proportion", repairNaN(( (double) table.getnMismatchingAlt() )/table.getnCalledCompGenotypes()));
+ concordanceCounts.set(rowKey,"Count",table.getnMismatchingAlt());
+
+ for ( Map.Entry nrsEntry : metrics.getPerSampleNRS().entrySet() ) {
+ concordanceSummary.set(nrsEntry.getKey(),"Sample",nrsEntry.getKey());
+ concordanceSummary.set(nrsEntry.getKey(),"Non-Reference_Sensitivity",nrsEntry.getValue());
+ }
+ for ( Map.Entry nrdEntry : metrics.getPerSampleNRD().entrySet() ) {
+ concordanceSummary.set(nrdEntry.getKey(),"Non-Reference_Discrepancy",nrdEntry.getValue());
+ }
+ concordanceSummary.set("ALL_NRS_NRD","Sample","ALL");
+ concordanceSummary.set("ALL_NRS_NRD","Non-Reference_Sensitivity",metrics.getOverallNRS());
+ concordanceSummary.set("ALL_NRS_NRD","Non-Reference_Discrepancy",metrics.getOverallNRD());
+
+
+ for (ConcordanceMetrics.SiteConcordanceType type : ConcordanceMetrics.SiteConcordanceType.values() ) {
+ siteConcordance.addColumn(type.toString(),"%d");
+ }
+
+ for (ConcordanceMetrics.SiteConcordanceType type : ConcordanceMetrics.SiteConcordanceType.values() ) {
+ siteConcordance.set("Comparison",type.toString(),metrics.getOverallSiteConcordance().get(type));
+ }
+
+ } else {
+ concordanceCompProportions.addColumn("Sample","%s");
+ concordanceCounts.addColumn("Sample","%s");
+ concordanceEvalProportions.addColumn("Sample","%s");
+ concordanceSummary.addColumn("Sample","%s");
+ for ( GenotypeType evalType : GenotypeType.values() ) {
+ for ( GenotypeType compType : GenotypeType.values() ) {
+ String colKey = String.format("%s_%s", evalType.toString(), compType.toString());
+ concordanceCounts.addColumn(colKey,"%d");
+ if ( evalType == GenotypeType.HET || evalType == GenotypeType.HOM_REF || evalType == GenotypeType.HOM_VAR)
+ concordanceEvalProportions.addColumn(colKey,"%.3f");
+ if ( compType == GenotypeType.HET || compType == GenotypeType.HOM_VAR || compType == GenotypeType.HOM_REF )
+ concordanceCompProportions.addColumn(colKey,"%.3f");
+ }
+ }
+ concordanceEvalProportions.addColumn("Mismatching_Alleles","%.3f");
+ concordanceCompProportions.addColumn("Mismatching_Alleles","%.3f");
+ concordanceCounts.addColumn("Mismatching_Alleles","%d");
+ concordanceSummary.addColumn("Non-Reference Sensitivity","%.3f");
+ concordanceSummary.addColumn("Non-Reference Discrepancy","%.3f");
+ for (ConcordanceMetrics.SiteConcordanceType type : ConcordanceMetrics.SiteConcordanceType.values() ) {
+ siteConcordance.addColumn(type.toString(),"%d");
+ }
+
+ for ( Map.Entry entry : metrics.getPerSampleGenotypeConcordance().entrySet() ) {
+ ConcordanceMetrics.GenotypeConcordanceTable table = entry.getValue();
+ concordanceEvalProportions.set(entry.getKey(),"Sample",entry.getKey());
+ concordanceCompProportions.set(entry.getKey(),"Sample",entry.getKey());
+ concordanceCounts.set(entry.getKey(),"Sample",entry.getKey());
+ for ( GenotypeType evalType : GenotypeType.values() ) {
+ for ( GenotypeType compType : GenotypeType.values() ) {
+ String colKey = String.format("%s_%s",evalType.toString(),compType.toString());
+ int count = table.get(evalType, compType);
+ concordanceCounts.set(entry.getKey(),colKey,count);
+ if ( evalType == GenotypeType.HET || evalType == GenotypeType.HOM_REF || evalType == GenotypeType.HOM_VAR)
+ concordanceEvalProportions.set(entry.getKey(),colKey,repairNaN(( (double) count)/table.getnEvalGenotypes(evalType)));
+ if ( compType == GenotypeType.HET || compType == GenotypeType.HOM_VAR || compType == GenotypeType.HOM_REF )
+ concordanceCompProportions.set(entry.getKey(),colKey,repairNaN(( (double) count)/table.getnCompGenotypes(compType)));
+ }
+ }
+ concordanceEvalProportions.set(entry.getKey(),"Mismatching_Alleles", repairNaN(( (double) table.getnMismatchingAlt() )/table.getnCalledEvalGenotypes()));
+ concordanceCompProportions.set(entry.getKey(),"Mismatching_Alleles", repairNaN(( (double) table.getnMismatchingAlt() )/table.getnCalledCompGenotypes()));
+ concordanceCounts.set(entry.getKey(),"Mismatching_Alleles",table.getnMismatchingAlt());
+ }
+
+ String rowKey = "ALL";
+ concordanceCompProportions.set(rowKey,"Sample",rowKey);
+ concordanceEvalProportions.set(rowKey,"Sample",rowKey);
+ concordanceCounts.set(rowKey,"Sample",rowKey);
+ ConcordanceMetrics.GenotypeConcordanceTable table = metrics.getOverallGenotypeConcordance();
for ( GenotypeType evalType : GenotypeType.values() ) {
for ( GenotypeType compType : GenotypeType.values() ) {
String colKey = String.format("%s_%s",evalType.toString(),compType.toString());
- int count = table.get(evalType, compType);
- concordanceCounts.set(entry.getKey(),colKey,count);
+ int count = table.get(evalType,compType);
+ concordanceCounts.set(rowKey,colKey,count);
if ( evalType == GenotypeType.HET || evalType == GenotypeType.HOM_REF || evalType == GenotypeType.HOM_VAR)
- concordanceEvalProportions.set(entry.getKey(),colKey,( (double) count)/table.getnEvalGenotypes(evalType));
+ concordanceEvalProportions.set(rowKey,colKey,repairNaN(( (double) count)/table.getnEvalGenotypes(evalType)));
if ( compType == GenotypeType.HET || compType == GenotypeType.HOM_VAR || compType == GenotypeType.HOM_REF )
- concordanceCompProportions.set(entry.getKey(),colKey,( (double) count)/table.getnCompGenotypes(compType));
+ concordanceCompProportions.set(rowKey,colKey,repairNaN(( (double) count)/table.getnCompGenotypes(compType)));
}
}
- concordanceEvalProportions.set(entry.getKey(),"Mismatching_Alleles", ( (double) table.getnMismatchingAlt() )/table.getnCalledEvalGenotypes());
- concordanceCompProportions.set(entry.getKey(),"Mismatching_Alleles", ( (double) table.getnMismatchingAlt() )/table.getnCalledCompGenotypes());
- concordanceCounts.set(entry.getKey(),"Mismatching_Alleles",table.getnMismatchingAlt());
- }
+ concordanceEvalProportions.set(rowKey,"Mismatching_Alleles", repairNaN(( (double) table.getnMismatchingAlt() )/table.getnCalledEvalGenotypes()));
+ concordanceCompProportions.set(rowKey,"Mismatching_Alleles", repairNaN(( (double) table.getnMismatchingAlt() )/table.getnCalledCompGenotypes()));
+ concordanceCounts.set(rowKey,"Mismatching_Alleles",table.getnMismatchingAlt());
- String rowKey = "ALL";
- concordanceCompProportions.set(rowKey,"Sample",rowKey);
- concordanceEvalProportions.set(rowKey,"Sample",rowKey);
- concordanceCounts.set(rowKey,"Sample",rowKey);
- ConcordanceMetrics.GenotypeConcordanceTable table = metrics.getOverallGenotypeConcordance();
- for ( GenotypeType evalType : GenotypeType.values() ) {
- for ( GenotypeType compType : GenotypeType.values() ) {
- String colKey = String.format("%s_%s",evalType.toString(),compType.toString());
- int count = table.get(evalType,compType);
- concordanceCounts.set(rowKey,colKey,count);
- if ( evalType == GenotypeType.HET || evalType == GenotypeType.HOM_REF || evalType == GenotypeType.HOM_VAR)
- concordanceEvalProportions.set(rowKey,colKey,( (double) count)/table.getnEvalGenotypes(evalType));
- if ( compType == GenotypeType.HET || compType == GenotypeType.HOM_VAR || compType == GenotypeType.HOM_REF )
- concordanceCompProportions.set(rowKey,colKey,( (double) count)/table.getnCompGenotypes(compType));
+ for ( Map.Entry nrsEntry : metrics.getPerSampleNRS().entrySet() ) {
+ concordanceSummary.set(nrsEntry.getKey(),"Sample",nrsEntry.getKey());
+ concordanceSummary.set(nrsEntry.getKey(),"Non-Reference Sensitivity",nrsEntry.getValue());
}
- }
- concordanceEvalProportions.set(rowKey,"Mismatching_Alleles", ( (double) table.getnMismatchingAlt() )/table.getnCalledEvalGenotypes());
- concordanceCompProportions.set(rowKey,"Mismatching_Alleles", ( (double) table.getnMismatchingAlt() )/table.getnCalledCompGenotypes());
- concordanceCounts.set(rowKey,"Mismatching_Alleles",table.getnMismatchingAlt());
+ for ( Map.Entry nrdEntry : metrics.getPerSampleNRD().entrySet() ) {
+ concordanceSummary.set(nrdEntry.getKey(),"Non-Reference Discrepancy",nrdEntry.getValue());
+ }
+ concordanceSummary.set("ALL","Sample","ALL");
+ concordanceSummary.set("ALL","Non-Reference Sensitivity",metrics.getOverallNRS());
+ concordanceSummary.set("ALL","Non-Reference Discrepancy",metrics.getOverallNRD());
- for ( Map.Entry nrsEntry : metrics.getPerSampleNRS().entrySet() ) {
- concordanceSummary.set(nrsEntry.getKey(),"Sample",nrsEntry.getKey());
- concordanceSummary.set(nrsEntry.getKey(),"Non-Reference Sensitivity",nrsEntry.getValue());
- }
- for ( Map.Entry nrdEntry : metrics.getPerSampleNRD().entrySet() ) {
- concordanceSummary.set(nrdEntry.getKey(),"Non-Reference Discrepancy",nrdEntry.getValue());
- }
- concordanceSummary.set("ALL","Sample","ALL");
- concordanceSummary.set("ALL","Non-Reference Sensitivity",metrics.getOverallNRS());
- concordanceSummary.set("ALL","Non-Reference Discrepancy",metrics.getOverallNRD());
-
- for (ConcordanceMetrics.SiteConcordanceType type : ConcordanceMetrics.SiteConcordanceType.values() ) {
- siteConcordance.set("Comparison",type.toString(),metrics.getOverallSiteConcordance().get(type));
+ for (ConcordanceMetrics.SiteConcordanceType type : ConcordanceMetrics.SiteConcordanceType.values() ) {
+ siteConcordance.set("Comparison",type.toString(),metrics.getOverallSiteConcordance().get(type));
+ }
}
report.addTable(concordanceCompProportions);
@@ -226,7 +446,7 @@ public class GenotypeConcordance extends RodWalker samples) {
+ public VariantContext createEmptyContext(VariantContext other, List samples) {
VariantContextBuilder builder = new VariantContextBuilder();
// set the alleles to be the same
builder.alleles(other.getAlleles());
@@ -239,13 +459,29 @@ public class GenotypeConcordance extends RodWalker exps) {
if ( ! context.isFiltered() || ignoreSiteFilter ) {
- // todo -- add genotype-level jexl filtering here
- return context;
+ List filteredGenotypes = new ArrayList(context.getNSamples());
+ for ( Genotype g : context.getGenotypes() ) {
+ Map matchMap = VariantContextUtils.match(context, g, exps);
+ boolean filtered = false;
+ for ( Boolean b : matchMap.values() ) {
+ if ( b ) {
+ filtered = true;
+ break;
+ }
+ }
+ if ( filtered ) {
+ filteredGenotypes.add(GenotypeBuilder.create(g.getSampleName(),Arrays.asList(Allele.NO_CALL,Allele.NO_CALL),g.getExtendedAttributes()));
+ } else {
+ filteredGenotypes.add(g);
+ }
+ }
+ VariantContextBuilder builder = new VariantContextBuilder(context);
+ builder.genotypes(filteredGenotypes);
+ return builder.make();
}
+
VariantContextBuilder builder = new VariantContextBuilder();
builder.alleles(Arrays.asList(context.getReference()));
builder.loc(context.getChr(),context.getStart(),context.getEnd());
@@ -256,4 +492,4 @@ public class GenotypeConcordance extends RodWalker implements TreeRe
headerLines = new LinkedHashSet(getSelectedHeaders(headerLines));
// Optionally add in the intervals.
- if (includeIntervals && getToolkit().getArguments().intervals != null) {
- for (IntervalBinding intervalBinding : getToolkit().getArguments().intervals) {
- String source = intervalBinding.getSource();
- if (source == null)
- continue;
- File file = new File(source);
- if (file.exists()) {
- headerLines.add(new VCFHeaderLine(VCFHeader.INTERVALS_KEY, FilenameUtils.getBaseName(file.getName())));
- } else {
- headerLines.add(new VCFHeaderLine(VCFHeader.INTERVALS_KEY, source));
+ if (includeIntervals) {
+ IntervalArgumentCollection intervalArguments = getToolkit().getArguments().intervalArguments;
+ if (intervalArguments.intervals != null) {
+ for (IntervalBinding intervalBinding : intervalArguments.intervals) {
+ String source = intervalBinding.getSource();
+ if (source == null)
+ continue;
+ File file = new File(source);
+ if (file.exists()) {
+ headerLines.add(new VCFHeaderLine(VCFHeader.INTERVALS_KEY, FilenameUtils.getBaseName(file.getName())));
+ } else {
+ headerLines.add(new VCFHeaderLine(VCFHeader.INTERVALS_KEY, source));
+ }
}
}
+
+ if (intervalArguments.excludeIntervals != null) {
+ for (IntervalBinding intervalBinding : intervalArguments.excludeIntervals) {
+ String source = intervalBinding.getSource();
+ if (source == null)
+ continue;
+ File file = new File(source);
+ if (file.exists()) {
+ headerLines.add(new VCFHeaderLine(VCFHeader.EXCLUDE_INTERVALS_KEY, FilenameUtils.getBaseName(file.getName())));
+ } else {
+ headerLines.add(new VCFHeaderLine(VCFHeader.EXCLUDE_INTERVALS_KEY, source));
+ }
+ }
+ }
+
+ if (intervalArguments.intervalMerging != IntervalMergingRule.ALL) {
+ headerLines.add(new VCFHeaderLine(VCFHeader.INTERVAL_MERGING_KEY, String.valueOf(intervalArguments.intervalMerging)));
+ }
+
+ if (intervalArguments.intervalSetRule != IntervalSetRule.UNION) {
+ headerLines.add(new VCFHeaderLine(VCFHeader.INTERVAL_SET_RULE_KEY, String.valueOf(intervalArguments.intervalSetRule)));
+ }
+
+ if (intervalArguments.intervalPadding != 0) {
+ headerLines.add(new VCFHeaderLine(VCFHeader.INTERVAL_PADDING_KEY, String.valueOf(intervalArguments.intervalPadding)));
+ }
}
TreeSet vcfSamples = new TreeSet(SampleUtils.getSampleList(vcfRods, VariantContextUtils.GenotypeMergeType.REQUIRE_UNIQUE));
diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/ArtificialReadPileupTestProvider.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/ArtificialReadPileupTestProvider.java
index 80ef7293f..047d69c5f 100644
--- a/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/ArtificialReadPileupTestProvider.java
+++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/ArtificialReadPileupTestProvider.java
@@ -51,6 +51,7 @@ import net.sf.samtools.SAMReadGroupRecord;
import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
+import org.broadinstitute.sting.utils.locusiterator.LocusIteratorByState;
import org.broadinstitute.variant.utils.BaseUtils;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.GenomeLocParser;
@@ -214,8 +215,7 @@ public class ArtificialReadPileupTestProvider {
read.setReadNegativeStrandFlag(false);
read.setReadGroup(sampleRG(sample));
-
- pileupElements.add(new PileupElement(read,readOffset,false,isBeforeDeletion, false, isBeforeInsertion,false,false,altBases,Math.abs(eventLength)));
+ pileupElements.add(LocusIteratorByState.createPileupForReadAndOffset(read, readOffset));
}
return pileupElements;
diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java
index 527e5c5e1..45a42d018 100644
--- a/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java
+++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java
@@ -108,7 +108,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
public void testMultipleSNPAlleles() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
"-T UnifiedGenotyper -R " + b37KGReference + " --no_cmdline_in_header -glm BOTH --dbsnp " + b37dbSNP129 + " -I " + privateTestDir + "multiallelic.snps.bam -o %s -L " + privateTestDir + "multiallelic.snps.intervals", 1,
- Arrays.asList("b41b95aaa2c453c9b75b3b29a9c2718e"));
+ Arrays.asList("35479a79e1ce7c15493bd77e58cadcaa"));
executeTest("test Multiple SNP alleles", spec);
}
@@ -124,7 +124,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
public void testReverseTrim() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
"-T UnifiedGenotyper -R " + b37KGReference + " --no_cmdline_in_header -glm INDEL -I " + validationDataLocation + "CEUTrio.HiSeq.b37.chr20.10_11mb.bam -o %s -L 20:10289124 -L 20:10090289", 1,
- Arrays.asList("44e9f6cf11b4efecb454cd3de8de9877"));
+ Arrays.asList("1e61de694b51d7c0f26da5179ee6bb0c"));
executeTest("test reverse trim", spec);
}
@@ -238,12 +238,12 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
@Test
public void testOutputParameterAllConfident() {
- testOutputParameters("--output_mode EMIT_ALL_CONFIDENT_SITES", "9dbc9389db39cf9697e93e0bf529314f");
+ testOutputParameters("--output_mode EMIT_ALL_CONFIDENT_SITES", "5649f72de04e1391e0f2bb86843d3d72");
}
@Test
public void testOutputParameterAllSites() {
- testOutputParameters("--output_mode EMIT_ALL_SITES", "8b26088a035e579c4afd3b46737291e4");
+ testOutputParameters("--output_mode EMIT_ALL_SITES", "cb151bb9e90680b12714d481091ed209");
}
private void testOutputParameters(final String args, final String md5) {
@@ -363,7 +363,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
" -o %s" +
" -L 1:10,000,000-10,500,000",
1,
- Arrays.asList("39c7a813fd6ee82d3604f2a868b35b2a"));
+ Arrays.asList("8231ae37b52b927db9fc1e5c221b0ba0"));
executeTest(String.format("test indel calling, multiple technologies"), spec);
}
@@ -391,13 +391,13 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
public void testMultiSampleIndels1() {
WalkerTest.WalkerTestSpec spec1 = new WalkerTest.WalkerTestSpec(
baseCommandIndels + " -I " + validationDataLocation + "low_coverage_CEU.chr1.10k-11k.bam -o %s -L 1:10450700-10551000", 1,
- Arrays.asList("5667a699a3a13474f2d1cd2d6b01cd5b"));
+ Arrays.asList("a47810de2f6ef8087f4644064a0814bc"));
List result = executeTest("test MultiSample Pilot1 CEU indels", spec1).getFirst();
WalkerTest.WalkerTestSpec spec2 = new WalkerTest.WalkerTestSpec(
baseCommandIndels + " --genotyping_mode GENOTYPE_GIVEN_ALLELES -alleles " + result.get(0).getAbsolutePath() + " -I " + validationDataLocation +
"low_coverage_CEU.chr1.10k-11k.bam -o %s -L 1:10450700-10551000", 1,
- Arrays.asList("b6c1d5cd28ff584c5f5037afef4e883a"));
+ Arrays.asList("53b8d2b0fa63c5d1019855e8e0db28f0"));
executeTest("test MultiSample Pilot1 CEU indels using GENOTYPE_GIVEN_ALLELES", spec2);
}
@@ -497,18 +497,18 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
public void testReducedBam() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
"-T UnifiedGenotyper -R " + b37KGReference + " --no_cmdline_in_header -I " + privateTestDir + "bamExample.ReducedRead.ADAnnotation.bam -o %s -L 1:67,225,396-67,288,518", 1,
- Arrays.asList("092e42a712afb660ec79ff11c55933e2"));
+ Arrays.asList("8b9a9fc2e7150acbe2dac91b4620f304"));
executeTest("test calling on a ReducedRead BAM", spec);
}
@Test
public void testReducedBamSNPs() {
- testReducedCalling("SNP", "c0de74ab8f4f14eb3a2c5d55c200ac5f");
+ testReducedCalling("SNP", "b5991dddbfb59366614ff8819062649f");
}
@Test
public void testReducedBamINDELs() {
- testReducedCalling("INDEL", "1c9aaf65ffaa12bb766855265a1c3f8e");
+ testReducedCalling("INDEL", "a85c110fcac9574a54c7daccb1e2d5ae");
}
diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java
index 142fa39bf..939b9873c 100644
--- a/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java
+++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java
@@ -50,6 +50,7 @@ import org.broadinstitute.sting.WalkerTest;
import org.testng.annotations.Test;
import java.util.Arrays;
+import java.util.Collections;
public class HaplotypeCallerIntegrationTest extends WalkerTest {
final static String REF = b37KGReference;
@@ -67,7 +68,7 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
@Test
public void testHaplotypeCallerMultiSample() {
- HCTest(CEUTRIO_BAM, "", "35c8425b44429ac7468c2cd26f8f5a42");
+ HCTest(CEUTRIO_BAM, "", "b8f7b741445ce6b6ea491c794ce75c17");
}
@Test
@@ -75,11 +76,33 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
HCTest(NA12878_BAM, "", "a2c63f6e6e51a01019bdbd23125bdb15");
}
- // TODO -- add more tests for GGA mode, especially with input alleles that are complex variants and/or not trimmed
+ @Test(enabled = false)
+ public void testHaplotypeCallerSingleSampleWithDbsnp() {
+ HCTest(NA12878_BAM, "-D " + b37dbSNP132, "");
+ }
+
@Test
public void testHaplotypeCallerMultiSampleGGA() {
HCTest(CEUTRIO_BAM, "--max_alternate_alleles 3 -gt_mode GENOTYPE_GIVEN_ALLELES -out_mode EMIT_ALL_SITES -alleles " + validationDataLocation + "combined.phase1.chr20.raw.indels.sites.vcf",
- "d918d25b22a551cae5d70ea30d7feed1");
+ "c679ae7f04bdfda896b5c046d35e043c");
+ }
+
+ private void HCTestComplexGGA(String bam, String args, String md5) {
+ final String base = String.format("-T HaplotypeCaller -R %s -I %s", REF, bam) + " --no_cmdline_in_header -o %s -minPruning 3 -gt_mode GENOTYPE_GIVEN_ALLELES -out_mode EMIT_ALL_SITES -alleles " + validationDataLocation + "combined.phase1.chr20.raw.indels.sites.vcf";
+ final WalkerTestSpec spec = new WalkerTestSpec(base + " " + args, Arrays.asList(md5));
+ executeTest("testHaplotypeCallerComplexGGA: args=" + args, spec);
+ }
+
+ @Test
+ public void testHaplotypeCallerMultiSampleGGAComplex() {
+ HCTestComplexGGA(NA12878_CHR20_BAM, "-L 20:119673-119823 -L 20:121408-121538",
+ "8730a9ebaeecae913dca2fb5a0d4e946");
+ }
+
+ @Test
+ public void testHaplotypeCallerMultiSampleGGAMultiAllelic() {
+ HCTestComplexGGA(NA12878_CHR20_BAM, "-L 20:133041-133161 -L 20:300207-300337",
+ "d590c8d6d5e58d685401b65a23846893");
}
private void HCTestComplexVariants(String bam, String args, String md5) {
@@ -115,10 +138,15 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
HCTestIndelQualityScores(NA12878_RECALIBRATED_BAM, "", "29f1125df5ab27cc937a144ae08ac735");
}
+ // That problem bam came from a user on the forum and it spotted a problem where the ReadClipper
+ // was modifying the GATKSamRecord and that was screwing up the traversal engine from map call to
+ // map call. So the test is there for consistency but not for correctness. I'm not sure we can trust
+ // any of the calls in that region because it is so messy. The only thing I would maybe be worried about is
+ // that the three calls that are missing happen to all be the left most calls in the region
@Test
public void HCTestProblematicReadsModifiedInActiveRegions() {
final String base = String.format("-T HaplotypeCaller -R %s -I %s", REF, privateTestDir + "haplotype-problem-4.bam") + " --no_cmdline_in_header -o %s -minPruning 3 -L 4:49139026-49139965";
- final WalkerTestSpec spec = new WalkerTestSpec(base, Arrays.asList("2e8e6313228b0387008437feae7f5469"));
+ final WalkerTestSpec spec = new WalkerTestSpec(base, Arrays.asList("31db0a2d9eb07f86e0a89f0d97169072"));
executeTest("HCTestProblematicReadsModifiedInActiveRegions: ", spec);
}
@@ -129,6 +157,14 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
executeTest("HCTestStructuralIndels: ", spec);
}
+ @Test
+ public void HCTestDoesNotFailOnBadRefBase() {
+ // don't care about the output - just want to make sure it doesn't fail
+ final String base = String.format("-T HaplotypeCaller -R %s -I %s", REF, privateTestDir + "NA12878.readsOverBadBase.chr3.bam") + " --no_cmdline_in_header -o /dev/null -L 3:60830000-60840000 --minPruning 3 -stand_call_conf 2 -stand_emit_conf 2";
+ final WalkerTestSpec spec = new WalkerTestSpec(base, Collections.emptyList());
+ executeTest("HCTestDoesNotFailOnBadRefBase: ", spec);
+ }
+
// --------------------------------------------------------------------------------------------------------------
//
// testing reduced reads
@@ -142,4 +178,12 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
Arrays.asList("8a400b0c46f41447fcc35a907e34f384"));
executeTest("HC calling on a ReducedRead BAM", spec);
}
+
+ @Test
+ public void testReducedBamWithReadsNotFullySpanningDeletion() {
+ WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
+ "-T HaplotypeCaller -R " + b37KGReference + " --no_cmdline_in_header -I " + privateTestDir + "reduced.readNotFullySpanningDeletion.bam -o %s -L 1:167871297", 1,
+ Arrays.asList("4e8121dd9dc90478f237bd6ae4d19920"));
+ executeTest("test calling on a ReducedRead BAM where the reads do not fully span a deletion", spec);
+ }
}
diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/variantutils/ConcordanceMetricsUnitTest.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/variantutils/ConcordanceMetricsUnitTest.java
index 28f128dd3..6db44efd5 100644
--- a/protected/java/test/org/broadinstitute/sting/gatk/walkers/variantutils/ConcordanceMetricsUnitTest.java
+++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/variantutils/ConcordanceMetricsUnitTest.java
@@ -111,8 +111,8 @@ public class ConcordanceMetricsUnitTest extends BaseTest {
private Pair getData1() {
- Allele reference_A = Allele.create(BaseUtils.A,true);
- Allele alt_C = Allele.create(BaseUtils.C);
+ Allele reference_A = Allele.create(BaseUtils.Base.A.base,true);
+ Allele alt_C = Allele.create(BaseUtils.Base.C.base);
Genotype sam_1_1_eval = GenotypeBuilder.create("test1_sample1", Arrays.asList(reference_A,reference_A));
Genotype sam_1_2_eval = GenotypeBuilder.create("test1_sample2", Arrays.asList(reference_A,alt_C));
@@ -160,9 +160,9 @@ public class ConcordanceMetricsUnitTest extends BaseTest {
private Pair getData2() {
- Allele reference_A = Allele.create(BaseUtils.A,true);
- Allele alt_C = Allele.create(BaseUtils.C);
- Allele alt_T = Allele.create(BaseUtils.T);
+ Allele reference_A = Allele.create(BaseUtils.Base.A.base,true);
+ Allele alt_C = Allele.create(BaseUtils.Base.C.base);
+ Allele alt_T = Allele.create(BaseUtils.Base.T.base);
Genotype sam_1_1_eval = GenotypeBuilder.create("test1_sample1", Arrays.asList(reference_A,reference_A));
Genotype sam_1_2_eval = GenotypeBuilder.create("test1_sample2", Arrays.asList(reference_A,alt_T));
@@ -213,10 +213,10 @@ public class ConcordanceMetricsUnitTest extends BaseTest {
private Pair getData3() {
- Allele reference_ACT = Allele.create(new byte[]{BaseUtils.A,BaseUtils.C,BaseUtils.T},true);
- Allele alt_AC = Allele.create(new byte[]{BaseUtils.A,BaseUtils.C});
- Allele alt_A = Allele.create(BaseUtils.A);
- Allele alt_ATT = Allele.create(new byte[]{BaseUtils.A,BaseUtils.T,BaseUtils.T});
+ Allele reference_ACT = Allele.create(new byte[]{BaseUtils.Base.A.base,BaseUtils.Base.C.base,BaseUtils.Base.T.base},true);
+ Allele alt_AC = Allele.create(new byte[]{BaseUtils.Base.A.base,BaseUtils.Base.C.base});
+ Allele alt_A = Allele.create(BaseUtils.Base.A.base);
+ Allele alt_ATT = Allele.create(new byte[]{BaseUtils.Base.A.base,BaseUtils.Base.T.base,BaseUtils.Base.T.base});
Genotype sam_1_1_eval = GenotypeBuilder.create("test1_sample1", Arrays.asList(reference_ACT,alt_ATT));
Genotype sam_1_2_eval = GenotypeBuilder.create("test1_sample2", Arrays.asList(alt_A,alt_A));
@@ -267,9 +267,9 @@ public class ConcordanceMetricsUnitTest extends BaseTest {
private Pair getData4() {
- Allele reference_A = Allele.create(BaseUtils.A,true);
- Allele alt_C = Allele.create(BaseUtils.C);
- Allele alt_T = Allele.create(BaseUtils.T);
+ Allele reference_A = Allele.create(BaseUtils.Base.A.base,true);
+ Allele alt_C = Allele.create(BaseUtils.Base.C.base);
+ Allele alt_T = Allele.create(BaseUtils.Base.T.base);
Genotype sam_1_1_eval = GenotypeBuilder.create("test1_sample1", Arrays.asList(reference_A,reference_A));
Genotype sam_1_2_eval = GenotypeBuilder.create("test1_sample2", Arrays.asList(Allele.NO_CALL,Allele.NO_CALL));
@@ -316,9 +316,9 @@ public class ConcordanceMetricsUnitTest extends BaseTest {
private Pair getData5() {
- Allele reference_A = Allele.create(BaseUtils.A,true);
- Allele alt_C = Allele.create(BaseUtils.C);
- Allele alt_T = Allele.create(BaseUtils.T);
+ Allele reference_A = Allele.create(BaseUtils.Base.A.base,true);
+ Allele alt_C = Allele.create(BaseUtils.Base.C.base);
+ Allele alt_T = Allele.create(BaseUtils.Base.T.base);
Genotype sam_1_1_eval = GenotypeBuilder.create("test1_sample1", Arrays.asList(reference_A,reference_A));
Genotype sam_1_2_eval = GenotypeBuilder.create("test1_sample2", new ArrayList(0));
@@ -368,8 +368,8 @@ public class ConcordanceMetricsUnitTest extends BaseTest {
private List> getData6() {
- Allele reference_A = Allele.create(BaseUtils.A,true);
- Allele alt_C = Allele.create(BaseUtils.C);
+ Allele reference_A = Allele.create(BaseUtils.Base.A.base,true);
+ Allele alt_C = Allele.create(BaseUtils.Base.C.base);
// site 1 -
@@ -396,8 +396,8 @@ public class ConcordanceMetricsUnitTest extends BaseTest {
Pair testDataSite1 = new Pair(eval_1_builder.make(),truth_1_builder.make());
- reference_A = Allele.create(BaseUtils.A,true);
- Allele alt_T = Allele.create(BaseUtils.T);
+ reference_A = Allele.create(BaseUtils.Base.A.base,true);
+ Allele alt_T = Allele.create(BaseUtils.Base.T.base);
// site 2 -
// sample 1: no-call/hom-ref
@@ -421,7 +421,7 @@ public class ConcordanceMetricsUnitTest extends BaseTest {
Pair testDataSite2 = new Pair(eval_1_builder.make(),truth_1_builder.make());
- Allele alt_G = Allele.create(BaseUtils.G);
+ Allele alt_G = Allele.create(BaseUtils.Base.G.base);
// site 3 -
// sample 1: alleles do not match
@@ -605,10 +605,10 @@ public class ConcordanceMetricsUnitTest extends BaseTest {
public List> getData7() {
- Allele ref1 = Allele.create(BaseUtils.T,true);
- Allele alt1 = Allele.create(BaseUtils.C);
- Allele alt2 = Allele.create(BaseUtils.G);
- Allele alt3 = Allele.create(BaseUtils.A);
+ Allele ref1 = Allele.create(BaseUtils.Base.T.base,true);
+ Allele alt1 = Allele.create(BaseUtils.Base.C.base);
+ Allele alt2 = Allele.create(BaseUtils.Base.G.base);
+ Allele alt3 = Allele.create(BaseUtils.Base.A.base);
GenomeLoc loc1 = genomeLocParser.createGenomeLoc("chr1",1,1);
VariantContextBuilder site1Eval = new VariantContextBuilder();
diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/variantutils/GenotypeConcordanceIntegrationTest.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/variantutils/GenotypeConcordanceIntegrationTest.java
new file mode 100644
index 000000000..e69d1ee60
--- /dev/null
+++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/variantutils/GenotypeConcordanceIntegrationTest.java
@@ -0,0 +1,118 @@
+/*
+ * Copyright (c) 2010.
+ *
+ * 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.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+ * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+ * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+ * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+ * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
+ * THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+package org.broadinstitute.sting.gatk.walkers.variantutils;
+
+import org.broadinstitute.sting.WalkerTest;
+import org.broadinstitute.sting.utils.exceptions.UserException;
+import org.testng.annotations.Test;
+
+import java.util.Arrays;
+
+public class GenotypeConcordanceIntegrationTest extends WalkerTest {
+
+ protected static final String emptyMd5 = "d41d8cd98f00b204e9800998ecf8427e";
+
+ public static String baseTestString(String eval, String comp) {
+ return "-T GenotypeConcordance -R " + b37KGReference + " --eval " + validationDataLocation + eval + " --comp " + validationDataLocation + comp + " -o %s";
+ }
+
+ @Test
+ public void testIndelConcordance() {
+ WalkerTestSpec spec = new WalkerTestSpec(
+ baseTestString("NA12878.Jan2013.haplotypeCaller.subset.indels.vcf", "NA12878.Jan2013.bestPractices.subset.indels.vcf"),
+ 0,
+ Arrays.asList("0f29a0c6dc44066228c8cb204fd53ec0")
+ );
+
+ executeTest("test indel concordance", spec);
+ }
+
+ @Test
+ public void testNonoverlapingSamples() {
+ WalkerTestSpec spec = new WalkerTestSpec(
+ baseTestString("GenotypeConcordanceNonOverlapTest_Eval.vcf", "GenotypeConcordanceNonOverlapTest_Comp.vcf"),
+ 0,
+ Arrays.asList("fc725022d47b4b5f8a6ef87f0f1ffe89")
+ );
+
+ executeTest("test non-overlapping samples", spec);
+ }
+
+ @Test
+ public void testNonoverlappingSamplesMoltenized() {
+ WalkerTestSpec spec = new WalkerTestSpec(
+ baseTestString("GenotypeConcordanceNonOverlapTest_Eval.vcf", "GenotypeConcordanceNonOverlapTest_Comp.vcf") + " -moltenize",
+ 0,
+ Arrays.asList("370141088362d0ab7054be5249c49c11")
+ );
+
+ executeTest("Test moltenized output",spec);
+ }
+
+ @Test
+ public void testMultipleRecordsPerSite() {
+ WalkerTestSpec spec = new WalkerTestSpec(
+ baseTestString("GenotypeConcordance.multipleRecordsTest1.eval.vcf","GenotypeConcordance.multipleRecordsTest1.comp.vcf"),
+ 0,
+ Arrays.asList("352d59c4ac0cee5eb8ddbc9404b19ce9")
+ );
+
+ executeTest("test multiple records per site",spec);
+ }
+
+ @Test
+ public void testGQFilteringEval() {
+ WalkerTestSpec spec = new WalkerTestSpec(
+ baseTestString("genotypeConcordanceFilterTest.vcf","genotypeConcordanceFilterTest.vcf") + " -gfe 'GQ<30'",
+ 0,
+ Arrays.asList("b7b495ccfa6d50a6be3e095d3f6d3c52")
+ );
+
+ executeTest("Test filtering on the EVAL rod",spec);
+ }
+
+ @Test
+ public void testFloatFilteringComp() {
+ WalkerTestSpec spec = new WalkerTestSpec(
+ baseTestString("genotypeConcordanceFilterTest.vcf","genotypeConcordanceFilterTest.vcf") + " -gfc 'LX<0.50'",
+ 0,
+ Arrays.asList("6406b16cde7960b8943edf594303afd6")
+ );
+
+ executeTest("Test filtering on the COMP rod", spec);
+ }
+
+ @Test
+ public void testCombinedFilters() {
+ WalkerTestSpec spec = new WalkerTestSpec(
+ baseTestString("genotypeConcordanceFilterTest.vcf","genotypeConcordanceFilterTest.vcf") + " -gfc 'LX<0.52' -gfe 'DP<5' -gfe 'GQ<37'",
+ 0,
+ Arrays.asList("26ffd06215b6177acce0ea9f35d73d31")
+ );
+
+ executeTest("Test filtering on both rods",spec);
+ }
+}
diff --git a/public/java/src/org/broadinstitute/sting/commandline/IntervalArgumentCollection.java b/public/java/src/org/broadinstitute/sting/commandline/IntervalArgumentCollection.java
new file mode 100644
index 000000000..3f76ae652
--- /dev/null
+++ b/public/java/src/org/broadinstitute/sting/commandline/IntervalArgumentCollection.java
@@ -0,0 +1,70 @@
+/*
+* Copyright (c) 2012 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.
+*
+* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
+* THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+
+package org.broadinstitute.sting.commandline;
+
+import org.broad.tribble.Feature;
+import org.broadinstitute.sting.utils.interval.IntervalMergingRule;
+import org.broadinstitute.sting.utils.interval.IntervalSetRule;
+
+import java.util.List;
+
+public class IntervalArgumentCollection {
+ /**
+ * Using this option one can instruct the GATK engine to traverse over only part of the genome. This argument can be specified multiple times.
+ * One may use samtools-style intervals either explicitly (e.g. -L chr1 or -L chr1:100-200) or listed in a file (e.g. -L myFile.intervals).
+ * Additionally, one may specify a rod file to traverse over the positions for which there is a record in the file (e.g. -L file.vcf).
+ * To specify the completely unmapped reads in the BAM file (i.e. those without a reference contig) use -L unmapped.
+ */
+ @Input(fullName = "intervals", shortName = "L", doc = "One or more genomic intervals over which to operate. Can be explicitly specified on the command line or in a file (including a rod file)", required = false)
+ public List> intervals = null;
+
+ /**
+ * Using this option one can instruct the GATK engine NOT to traverse over certain parts of the genome. This argument can be specified multiple times.
+ * One may use samtools-style intervals either explicitly (e.g. -XL chr1 or -XL chr1:100-200) or listed in a file (e.g. -XL myFile.intervals).
+ * Additionally, one may specify a rod file to skip over the positions for which there is a record in the file (e.g. -XL file.vcf).
+ */
+ @Input(fullName = "excludeIntervals", shortName = "XL", doc = "One or more genomic intervals to exclude from processing. Can be explicitly specified on the command line or in a file (including a rod file)", required = false)
+ public List> excludeIntervals = null;
+
+ /**
+ * How should the intervals specified by multiple -L or -XL arguments be combined? Using this argument one can, for example, traverse over all of the positions
+ * for which there is a record in a VCF but just in chromosome 20 (-L chr20 -L file.vcf -isr INTERSECTION).
+ */
+ @Argument(fullName = "interval_set_rule", shortName = "isr", doc = "Indicates the set merging approach the interval parser should use to combine the various -L or -XL inputs", required = false)
+ public IntervalSetRule intervalSetRule = IntervalSetRule.UNION;
+
+ /**
+ * Should abutting (but not overlapping) intervals be treated as separate intervals?
+ */
+ @Argument(fullName = "interval_merging", shortName = "im", doc = "Indicates the interval merging rule we should use for abutting intervals", required = false)
+ public IntervalMergingRule intervalMerging = IntervalMergingRule.ALL;
+
+ /**
+ * For example, '-L chr1:100' with a padding value of 20 would turn into '-L chr1:80-120'.
+ */
+ @Argument(fullName = "interval_padding", shortName = "ip", doc = "Indicates how many basepairs of padding to include around each of the intervals specified with the -L/--intervals argument", required = false)
+ public int intervalPadding = 0;
+}
diff --git a/public/java/src/org/broadinstitute/sting/gatk/GenomeAnalysisEngine.java b/public/java/src/org/broadinstitute/sting/gatk/GenomeAnalysisEngine.java
index ba5577730..9b801be7d 100644
--- a/public/java/src/org/broadinstitute/sting/gatk/GenomeAnalysisEngine.java
+++ b/public/java/src/org/broadinstitute/sting/gatk/GenomeAnalysisEngine.java
@@ -55,7 +55,6 @@ import org.broadinstitute.sting.gatk.samples.SampleDBBuilder;
import org.broadinstitute.sting.gatk.walkers.*;
import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.classloader.PluginManager;
-import org.broadinstitute.sting.utils.collections.Pair;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.interval.IntervalUtils;
@@ -361,7 +360,6 @@ public class GenomeAnalysisEngine {
* Returns a list of active, initialized read transformers
*
* @param walker the walker we need to apply read transformers too
- * @return a non-null list of read transformers
*/
public void initializeReadTransformers(final Walker walker) {
final List activeTransformers = new ArrayList();
@@ -672,41 +670,7 @@ public class GenomeAnalysisEngine {
* Setup the intervals to be processed
*/
protected void initializeIntervals() {
- // return if no interval arguments at all
- if ( argCollection.intervals == null && argCollection.excludeIntervals == null )
- return;
-
- // Note that the use of '-L all' is no longer supported.
-
- // if include argument isn't given, create new set of all possible intervals
-
- final Pair includeExcludePair = IntervalUtils.parseIntervalBindingsPair(
- this.referenceDataSource,
- argCollection.intervals,
- argCollection.intervalSetRule, argCollection.intervalMerging, argCollection.intervalPadding,
- argCollection.excludeIntervals);
-
- final GenomeLocSortedSet includeSortedSet = includeExcludePair.getFirst();
- final GenomeLocSortedSet excludeSortedSet = includeExcludePair.getSecond();
-
- // if no exclude arguments, can return parseIntervalArguments directly
- if ( excludeSortedSet == null )
- intervals = includeSortedSet;
-
- // otherwise there are exclude arguments => must merge include and exclude GenomeLocSortedSets
- else {
- intervals = includeSortedSet.subtractRegions(excludeSortedSet);
-
- // logging messages only printed when exclude (-XL) arguments are given
- final long toPruneSize = includeSortedSet.coveredSize();
- final long toExcludeSize = excludeSortedSet.coveredSize();
- final long intervalSize = intervals.coveredSize();
- logger.info(String.format("Initial include intervals span %d loci; exclude intervals span %d loci", toPruneSize, toExcludeSize));
- logger.info(String.format("Excluding %d loci from original intervals (%.2f%% reduction)",
- toPruneSize - intervalSize, (toPruneSize - intervalSize) / (0.01 * toPruneSize)));
- }
-
- logger.info(String.format("Processing %d bp from intervals", intervals.coveredSize()));
+ intervals = IntervalUtils.parseIntervalArguments(this.referenceDataSource, argCollection.intervalArguments);
}
/**
@@ -842,6 +806,8 @@ public class GenomeAnalysisEngine {
if (argCollection.keepProgramRecords)
removeProgramRecords = false;
+ final boolean keepReadsInLIBS = walker instanceof ActiveRegionWalker;
+
return new SAMDataSource(
samReaderIDs,
threadAllocation,
@@ -856,7 +822,8 @@ public class GenomeAnalysisEngine {
readTransformers,
includeReadsWithDeletionAtLoci(),
argCollection.defaultBaseQualities,
- removeProgramRecords);
+ removeProgramRecords,
+ keepReadsInLIBS);
}
/**
diff --git a/public/java/src/org/broadinstitute/sting/gatk/ReadProperties.java b/public/java/src/org/broadinstitute/sting/gatk/ReadProperties.java
index 409b08e5d..1ca0a8a46 100644
--- a/public/java/src/org/broadinstitute/sting/gatk/ReadProperties.java
+++ b/public/java/src/org/broadinstitute/sting/gatk/ReadProperties.java
@@ -61,6 +61,7 @@ public class ReadProperties {
private final ValidationExclusion exclusionList;
private final Collection supplementalFilters;
private final List readTransformers;
+ private final boolean keepUniqueReadListInLIBS;
private final boolean includeReadsWithDeletionAtLoci;
private final boolean useOriginalBaseQualities;
private final byte defaultBaseQualities;
@@ -74,6 +75,10 @@ public class ReadProperties {
return includeReadsWithDeletionAtLoci;
}
+ public boolean keepUniqueReadListInLIBS() {
+ return keepUniqueReadListInLIBS;
+ }
+
/**
* Gets a list of the files acting as sources of reads.
* @return A list of files storing reads data.
@@ -161,6 +166,8 @@ public class ReadProperties {
* will explicitly list reads with deletion over the current reference base; otherwise, only observed
* bases will be seen in the pileups, and the deletions will be skipped silently.
* @param defaultBaseQualities if the reads have incomplete quality scores, set them all to defaultBaseQuality.
+ * @param keepUniqueReadListInLIBS If true, we will tell LocusIteratorByState to track the unique reads it sees
+ * This is really useful for ActiveRegionTraversals
*/
public ReadProperties( Collection samFiles,
SAMFileHeader header,
@@ -172,7 +179,8 @@ public class ReadProperties {
Collection supplementalFilters,
List readTransformers,
boolean includeReadsWithDeletionAtLoci,
- byte defaultBaseQualities) {
+ byte defaultBaseQualities,
+ final boolean keepUniqueReadListInLIBS) {
this.readers = samFiles;
this.header = header;
this.sortOrder = sortOrder;
@@ -184,5 +192,6 @@ public class ReadProperties {
this.includeReadsWithDeletionAtLoci = includeReadsWithDeletionAtLoci;
this.useOriginalBaseQualities = useOriginalBaseQualities;
this.defaultBaseQualities = defaultBaseQualities;
+ this.keepUniqueReadListInLIBS = keepUniqueReadListInLIBS;
}
}
diff --git a/public/java/src/org/broadinstitute/sting/gatk/arguments/GATKArgumentCollection.java b/public/java/src/org/broadinstitute/sting/gatk/arguments/GATKArgumentCollection.java
index ab09064dd..62ca38ad2 100644
--- a/public/java/src/org/broadinstitute/sting/gatk/arguments/GATKArgumentCollection.java
+++ b/public/java/src/org/broadinstitute/sting/gatk/arguments/GATKArgumentCollection.java
@@ -26,11 +26,7 @@
package org.broadinstitute.sting.gatk.arguments;
import net.sf.samtools.SAMFileReader;
-import org.broad.tribble.Feature;
-import org.broadinstitute.sting.commandline.Argument;
-import org.broadinstitute.sting.commandline.Hidden;
-import org.broadinstitute.sting.commandline.Input;
-import org.broadinstitute.sting.commandline.IntervalBinding;
+import org.broadinstitute.sting.commandline.*;
import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
import org.broadinstitute.sting.gatk.downsampling.DownsampleType;
import org.broadinstitute.sting.gatk.downsampling.DownsamplingMethod;
@@ -38,8 +34,6 @@ import org.broadinstitute.sting.gatk.phonehome.GATKRunReport;
import org.broadinstitute.sting.gatk.samples.PedigreeValidationType;
import org.broadinstitute.sting.utils.QualityUtils;
import org.broadinstitute.sting.utils.baq.BAQ;
-import org.broadinstitute.sting.utils.interval.IntervalMergingRule;
-import org.broadinstitute.sting.utils.interval.IntervalSetRule;
import java.io.File;
import java.util.ArrayList;
@@ -100,41 +94,8 @@ public class GATKArgumentCollection {
@Argument(fullName = "read_filter", shortName = "rf", doc = "Specify filtration criteria to apply to each read individually", required = false)
public List readFilters = new ArrayList();
- /**
- * Using this option one can instruct the GATK engine to traverse over only part of the genome. This argument can be specified multiple times.
- * One may use samtools-style intervals either explicitly (e.g. -L chr1 or -L chr1:100-200) or listed in a file (e.g. -L myFile.intervals).
- * Additionally, one may specify a rod file to traverse over the positions for which there is a record in the file (e.g. -L file.vcf).
- * To specify the completely unmapped reads in the BAM file (i.e. those without a reference contig) use -L unmapped.
- */
- @Input(fullName = "intervals", shortName = "L", doc = "One or more genomic intervals over which to operate. Can be explicitly specified on the command line or in a file (including a rod file)", required = false)
- public List> intervals = null;
-
- /**
- * Using this option one can instruct the GATK engine NOT to traverse over certain parts of the genome. This argument can be specified multiple times.
- * One may use samtools-style intervals either explicitly (e.g. -XL chr1 or -XL chr1:100-200) or listed in a file (e.g. -XL myFile.intervals).
- * Additionally, one may specify a rod file to skip over the positions for which there is a record in the file (e.g. -XL file.vcf).
- */
- @Input(fullName = "excludeIntervals", shortName = "XL", doc = "One or more genomic intervals to exclude from processing. Can be explicitly specified on the command line or in a file (including a rod file)", required = false)
- public List> excludeIntervals = null;
-
- /**
- * How should the intervals specified by multiple -L or -XL arguments be combined? Using this argument one can, for example, traverse over all of the positions
- * for which there is a record in a VCF but just in chromosome 20 (-L chr20 -L file.vcf -isr INTERSECTION).
- */
- @Argument(fullName = "interval_set_rule", shortName = "isr", doc = "Indicates the set merging approach the interval parser should use to combine the various -L or -XL inputs", required = false)
- public IntervalSetRule intervalSetRule = IntervalSetRule.UNION;
-
- /**
- * Should abutting (but not overlapping) intervals be treated as separate intervals?
- */
- @Argument(fullName = "interval_merging", shortName = "im", doc = "Indicates the interval merging rule we should use for abutting intervals", required = false)
- public IntervalMergingRule intervalMerging = IntervalMergingRule.ALL;
-
- /**
- * For example, '-L chr1:100' with a padding value of 20 would turn into '-L chr1:80-120'.
- */
- @Argument(fullName = "interval_padding", shortName = "ip", doc = "Indicates how many basepairs of padding to include around each of the intervals specified with the -L/--intervals argument", required = false)
- public int intervalPadding = 0;
+ @ArgumentCollection
+ public IntervalArgumentCollection intervalArguments = new IntervalArgumentCollection();
@Input(fullName = "reference_sequence", shortName = "R", doc = "Reference sequence file", required = false)
public File referenceFile = null;
diff --git a/public/java/src/org/broadinstitute/sting/gatk/datasources/providers/LocusShardDataProvider.java b/public/java/src/org/broadinstitute/sting/gatk/datasources/providers/LocusShardDataProvider.java
index 41fe5a175..45c9af995 100644
--- a/public/java/src/org/broadinstitute/sting/gatk/datasources/providers/LocusShardDataProvider.java
+++ b/public/java/src/org/broadinstitute/sting/gatk/datasources/providers/LocusShardDataProvider.java
@@ -29,7 +29,7 @@ import net.sf.picard.reference.IndexedFastaSequenceFile;
import org.broadinstitute.sting.gatk.ReadProperties;
import org.broadinstitute.sting.gatk.datasources.reads.Shard;
import org.broadinstitute.sting.gatk.datasources.rmd.ReferenceOrderedDataSource;
-import org.broadinstitute.sting.gatk.iterators.LocusIterator;
+import org.broadinstitute.sting.utils.locusiterator.LocusIterator;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.GenomeLocParser;
diff --git a/public/java/src/org/broadinstitute/sting/gatk/datasources/providers/LocusView.java b/public/java/src/org/broadinstitute/sting/gatk/datasources/providers/LocusView.java
index b020a43ba..f77819426 100644
--- a/public/java/src/org/broadinstitute/sting/gatk/datasources/providers/LocusView.java
+++ b/public/java/src/org/broadinstitute/sting/gatk/datasources/providers/LocusView.java
@@ -28,9 +28,10 @@ package org.broadinstitute.sting.gatk.datasources.providers;
import org.broadinstitute.sting.gatk.downsampling.DownsampleType;
import org.broadinstitute.sting.gatk.ReadProperties;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
-import org.broadinstitute.sting.gatk.iterators.LocusIterator;
+import org.broadinstitute.sting.utils.locusiterator.LocusIterator;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.GenomeLocParser;
+import org.broadinstitute.sting.utils.locusiterator.LocusIteratorByState;
import java.util.Arrays;
import java.util.Collection;
@@ -212,4 +213,10 @@ public abstract class LocusView extends LocusIterator implements View {
private boolean isContainedInShard(GenomeLoc location) {
return locus.containsP(location);
}
+
+ // TODO -- remove me
+ @Override
+ public LocusIteratorByState getLIBS() {
+ return loci.getLIBS();
+ }
}
diff --git a/public/java/src/org/broadinstitute/sting/gatk/datasources/providers/ManagingReferenceOrderedView.java b/public/java/src/org/broadinstitute/sting/gatk/datasources/providers/ManagingReferenceOrderedView.java
index 7d3cac33d..09b72f5eb 100644
--- a/public/java/src/org/broadinstitute/sting/gatk/datasources/providers/ManagingReferenceOrderedView.java
+++ b/public/java/src/org/broadinstitute/sting/gatk/datasources/providers/ManagingReferenceOrderedView.java
@@ -77,13 +77,17 @@ public class ManagingReferenceOrderedView implements ReferenceOrderedView {
* @return A tracker containing information about this locus.
*/
public RefMetaDataTracker getReferenceOrderedDataAtLocus( GenomeLoc loc, ReferenceContext referenceContext ) {
- List bindings = states.isEmpty() ? Collections.emptyList() : new ArrayList(states.size());
+ if ( states.isEmpty() )
+ return RefMetaDataTracker.EMPTY_TRACKER;
+ else {
+ List bindings = new ArrayList(states.size());
- for ( ReferenceOrderedDataState state: states )
- // todo -- warning, I removed the reference to the name from states
- bindings.add( state.iterator.seekForward(loc) );
+ for ( ReferenceOrderedDataState state: states )
+ // todo -- warning, I removed the reference to the name from states
+ bindings.add( state.iterator.seekForward(loc) );
- return new RefMetaDataTracker(bindings);
+ return new RefMetaDataTracker(bindings);
+ }
}
/**
diff --git a/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/SAMDataSource.java b/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/SAMDataSource.java
index cb47ffe4c..c9a3b0df0 100644
--- a/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/SAMDataSource.java
+++ b/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/SAMDataSource.java
@@ -158,6 +158,9 @@ public class SAMDataSource {
/**
* Create a new SAM data source given the supplied read metadata.
+ *
+ * For testing purposes
+ *
* @param samFiles list of reads files.
*/
public SAMDataSource(Collection samFiles, ThreadAllocation threadAllocation, Integer numFileHandles, GenomeLocParser genomeLocParser) {
@@ -177,6 +180,8 @@ public class SAMDataSource {
/**
* See complete constructor. Does not enable BAQ by default.
+ *
+ * For testing purposes
*/
public SAMDataSource(
Collection samFiles,
@@ -203,6 +208,7 @@ public class SAMDataSource {
Collections.emptyList(),
includeReadsWithDeletionAtLoci,
(byte) -1,
+ false,
false);
}
@@ -219,6 +225,7 @@ public class SAMDataSource {
* will explicitly list reads with deletion over the current reference base; otherwise, only observed
* bases will be seen in the pileups, and the deletions will be skipped silently.
* @param defaultBaseQualities if the reads have incomplete quality scores, set them all to defaultBaseQuality.
+ * @param keepReadsInLIBS should we keep a unique list of reads in LIBS?
*/
public SAMDataSource(
Collection samFiles,
@@ -234,7 +241,8 @@ public class SAMDataSource {
List readTransformers,
boolean includeReadsWithDeletionAtLoci,
byte defaultBaseQualities,
- boolean removeProgramRecords) {
+ boolean removeProgramRecords,
+ final boolean keepReadsInLIBS) {
this.readMetrics = new ReadMetrics();
this.genomeLocParser = genomeLocParser;
@@ -306,7 +314,8 @@ public class SAMDataSource {
supplementalFilters,
readTransformers,
includeReadsWithDeletionAtLoci,
- defaultBaseQualities);
+ defaultBaseQualities,
+ keepReadsInLIBS);
// cache the read group id (original) -> read group id (merged)
// and read group id (merged) -> read group id (original) mappings.
diff --git a/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/Shard.java b/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/Shard.java
index 2c03363ba..5b4c2afda 100644
--- a/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/Shard.java
+++ b/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/Shard.java
@@ -95,7 +95,10 @@ public abstract class Shard implements HasGenomeLocation {
*/
private final Map fileSpans;
-
+ /**
+ * Lazy-calculated span of all of the genome locs in this shard
+ */
+ private GenomeLoc spanningLocation = null;
/**
* Statistics about which reads in this shards were used and which were filtered away.
@@ -148,27 +151,34 @@ public abstract class Shard implements HasGenomeLocation {
/**
* Returns the span of the genomeLocs comprising this shard
- * @param
- * @return
+ * @return a GenomeLoc that starts as the first position in getGenomeLocs() and stops at the stop of the last
+ * position in getGenomeLocs()
*/
public GenomeLoc getLocation() {
- if ( getGenomeLocs() == null )
- return GenomeLoc.WHOLE_GENOME;
+ if ( spanningLocation == null ) {
+ if ( getGenomeLocs() == null )
+ spanningLocation = GenomeLoc.WHOLE_GENOME;
+ else if ( getGenomeLocs().size() == 0 ) {
+ spanningLocation = getGenomeLocs().get(0);
+ } else {
+ int start = Integer.MAX_VALUE;
+ int stop = Integer.MIN_VALUE;
+ String contig = null;
- int start = Integer.MAX_VALUE;
- int stop = Integer.MIN_VALUE;
- String contig = null;
+ for ( GenomeLoc loc : getGenomeLocs() ) {
+ if ( GenomeLoc.isUnmapped(loc) )
+ // special case the unmapped region marker, just abort out
+ return loc;
+ contig = loc.getContig();
+ if ( loc.getStart() < start ) start = loc.getStart();
+ if ( loc.getStop() > stop ) stop = loc.getStop();
+ }
- for ( GenomeLoc loc : getGenomeLocs() ) {
- if ( GenomeLoc.isUnmapped(loc) )
- // special case the unmapped region marker, just abort out
- return loc;
- contig = loc.getContig();
- if ( loc.getStart() < start ) start = loc.getStart();
- if ( loc.getStop() > stop ) stop = loc.getStop();
+ spanningLocation = parser.createGenomeLoc(contig, start, stop);
+ }
}
- return parser.createGenomeLoc(contig, start, stop);
+ return spanningLocation;
}
diff --git a/public/java/src/org/broadinstitute/sting/gatk/downsampling/PassThroughDownsampler.java b/public/java/src/org/broadinstitute/sting/gatk/downsampling/PassThroughDownsampler.java
index 600834012..b06d5f5b4 100644
--- a/public/java/src/org/broadinstitute/sting/gatk/downsampling/PassThroughDownsampler.java
+++ b/public/java/src/org/broadinstitute/sting/gatk/downsampling/PassThroughDownsampler.java
@@ -27,8 +27,8 @@ package org.broadinstitute.sting.gatk.downsampling;
import net.sf.samtools.SAMRecord;
-import java.util.ArrayList;
import java.util.Collection;
+import java.util.LinkedList;
import java.util.List;
/**
@@ -41,7 +41,7 @@ import java.util.List;
*/
public class PassThroughDownsampler implements ReadsDownsampler {
- private ArrayList selectedReads;
+ private LinkedList selectedReads;
public PassThroughDownsampler() {
clear();
@@ -59,9 +59,13 @@ public class PassThroughDownsampler implements ReadsDownsam
}
public boolean hasFinalizedItems() {
- return selectedReads.size() > 0;
+ return ! selectedReads.isEmpty();
}
+ /**
+ * Note that this list is a linked list and so doesn't support fast random access
+ * @return
+ */
public List consumeFinalizedItems() {
// pass by reference rather than make a copy, for speed
List downsampledItems = selectedReads;
@@ -74,7 +78,7 @@ public class PassThroughDownsampler implements ReadsDownsam
}
public T peekFinalized() {
- return selectedReads.isEmpty() ? null : selectedReads.get(0);
+ return selectedReads.isEmpty() ? null : selectedReads.getFirst();
}
public T peekPending() {
@@ -90,7 +94,7 @@ public class PassThroughDownsampler implements ReadsDownsam
}
public void clear() {
- selectedReads = new ArrayList();
+ selectedReads = new LinkedList();
}
public void reset() {
diff --git a/public/java/src/org/broadinstitute/sting/gatk/downsampling/ReservoirDownsampler.java b/public/java/src/org/broadinstitute/sting/gatk/downsampling/ReservoirDownsampler.java
index 0d7a0dd14..4331fd723 100644
--- a/public/java/src/org/broadinstitute/sting/gatk/downsampling/ReservoirDownsampler.java
+++ b/public/java/src/org/broadinstitute/sting/gatk/downsampling/ReservoirDownsampler.java
@@ -29,9 +29,7 @@ import net.sf.samtools.SAMRecord;
import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
-import java.util.ArrayList;
-import java.util.Collection;
-import java.util.List;
+import java.util.*;
/**
* Reservoir Downsampler: Selects n reads out of a stream whose size is not known in advance, with
@@ -42,10 +40,25 @@ import java.util.List;
* @author David Roazen
*/
public class ReservoirDownsampler implements ReadsDownsampler {
+ private final int targetSampleSize;
- private ArrayList reservoir;
+ /**
+ * if true, this downsampler will be optimized for the case
+ * where most of the time we won't fill up anything like the
+ * targetSampleSize elements. If this is false, we will allocate
+ * internal buffers to targetSampleSize initially, which minimizes
+ * the cost of allocation if we often use targetSampleSize or more
+ * elements.
+ */
+ private final boolean expectFewOverflows;
- private int targetSampleSize;
+ /**
+ * At times this can be a linked list or an array list, depending on how we're accessing the
+ * data and whether or not we're expecting few overflows
+ */
+ private List reservoir;
+
+ private boolean isLinkedList;
private int totalReadsSeen;
@@ -56,17 +69,35 @@ public class ReservoirDownsampler implements ReadsDownsampl
*
* @param targetSampleSize Size of the reservoir used by this downsampler. Number of items retained
* after downsampling will be min(totalReads, targetSampleSize)
+ * @param expectFewOverflows if true, this downsampler will be optimized for the case
+ * where most of the time we won't fill up anything like the
+ * targetSampleSize elements. If this is false, we will allocate
+ * internal buffers to targetSampleSize initially, which minimizes
+ * the cost of allocation if we often use targetSampleSize or more
+ * elements.
*/
- public ReservoirDownsampler ( int targetSampleSize ) {
+ public ReservoirDownsampler ( final int targetSampleSize, final boolean expectFewOverflows) {
if ( targetSampleSize <= 0 ) {
throw new ReviewedStingException("Cannot do reservoir downsampling with a sample size <= 0");
}
this.targetSampleSize = targetSampleSize;
+ this.expectFewOverflows = expectFewOverflows;
clear();
reset();
}
+ /**
+ * Construct a ReservoirDownsampler
+ *
+ * @param targetSampleSize Size of the reservoir used by this downsampler. Number of items retained
+ * after downsampling will be min(totalReads, targetSampleSize)
+ */
+ public ReservoirDownsampler ( int targetSampleSize ) {
+ this(targetSampleSize, false);
+ }
+
+
public void submit ( T newRead ) {
totalReadsSeen++;
@@ -74,7 +105,12 @@ public class ReservoirDownsampler implements ReadsDownsampl
reservoir.add(newRead);
}
else {
- int randomSlot = GenomeAnalysisEngine.getRandomGenerator().nextInt(totalReadsSeen);
+ if ( isLinkedList ) {
+ reservoir = new ArrayList(reservoir);
+ isLinkedList = false;
+ }
+
+ final int randomSlot = GenomeAnalysisEngine.getRandomGenerator().nextInt(totalReadsSeen);
if ( randomSlot < targetSampleSize ) {
reservoir.set(randomSlot, newRead);
}
@@ -93,10 +129,15 @@ public class ReservoirDownsampler implements ReadsDownsampl
}
public List consumeFinalizedItems() {
- // pass by reference rather than make a copy, for speed
- List downsampledItems = reservoir;
- clear();
- return downsampledItems;
+ if ( reservoir.isEmpty() ) {
+ // if there's nothing here, don't both allocating a new list completely
+ return Collections.emptyList();
+ } else {
+ // pass by reference rather than make a copy, for speed
+ List downsampledItems = reservoir;
+ clear();
+ return downsampledItems;
+ }
}
public boolean hasPendingItems() {
@@ -119,9 +160,18 @@ public class ReservoirDownsampler implements ReadsDownsampl
// NO-OP
}
+ /**
+ * Clear the data structures used to hold information
+ */
public void clear() {
- reservoir = new ArrayList(targetSampleSize);
- totalReadsSeen = 0; // an internal stat used by the downsampling process, so not cleared by reset() below
+ // if we aren't expecting many overflows, allocate a linked list not an arraylist
+ reservoir = expectFewOverflows ? new LinkedList() : new ArrayList(targetSampleSize);
+
+ // it's a linked list if we allocate one
+ isLinkedList = expectFewOverflows;
+
+ // an internal stat used by the downsampling process, so not cleared by reset() below
+ totalReadsSeen = 0;
}
public void reset() {
diff --git a/public/java/src/org/broadinstitute/sting/gatk/executive/LinearMicroScheduler.java b/public/java/src/org/broadinstitute/sting/gatk/executive/LinearMicroScheduler.java
index 36d087735..4c0358d40 100644
--- a/public/java/src/org/broadinstitute/sting/gatk/executive/LinearMicroScheduler.java
+++ b/public/java/src/org/broadinstitute/sting/gatk/executive/LinearMicroScheduler.java
@@ -114,7 +114,7 @@ public class LinearMicroScheduler extends MicroScheduler {
}
// Special function call to empty out the work queue. Ugly for now but will be cleaned up when we eventually push this functionality more into the engine
- if( traversalEngine instanceof TraverseActiveRegions ) {
+ if( traversalEngine instanceof TraverseActiveRegions) {
final Object result = ((TraverseActiveRegions) traversalEngine).endTraversal(walker, accumulator.getReduceInit());
accumulator.accumulate(null, result); // Assumes only used with StandardAccumulator
}
diff --git a/public/java/src/org/broadinstitute/sting/gatk/executive/MicroScheduler.java b/public/java/src/org/broadinstitute/sting/gatk/executive/MicroScheduler.java
index a01af80ac..371cce778 100644
--- a/public/java/src/org/broadinstitute/sting/gatk/executive/MicroScheduler.java
+++ b/public/java/src/org/broadinstitute/sting/gatk/executive/MicroScheduler.java
@@ -213,7 +213,7 @@ public abstract class MicroScheduler implements MicroSchedulerMBean {
// Now that we have a progress meter, go through and initialize the traversal engines
for ( final TraversalEngine traversalEngine : allCreatedTraversalEngines )
- traversalEngine.initialize(engine, progressMeter);
+ traversalEngine.initialize(engine, walker, progressMeter);
// JMX does not allow multiple instances with the same ObjectName to be registered with the same platform MXBean.
// To get around this limitation and since we have no job identifier at this point, register a simple counter that
diff --git a/public/java/src/org/broadinstitute/sting/gatk/executive/WindowMaker.java b/public/java/src/org/broadinstitute/sting/gatk/executive/WindowMaker.java
index cbcc4abae..f587442d7 100644
--- a/public/java/src/org/broadinstitute/sting/gatk/executive/WindowMaker.java
+++ b/public/java/src/org/broadinstitute/sting/gatk/executive/WindowMaker.java
@@ -29,13 +29,14 @@ import net.sf.picard.util.PeekableIterator;
import org.broadinstitute.sting.gatk.ReadProperties;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.datasources.reads.Shard;
-import org.broadinstitute.sting.gatk.iterators.LegacyLocusIteratorByState;
-import org.broadinstitute.sting.gatk.iterators.LocusIterator;
-import org.broadinstitute.sting.gatk.iterators.LocusIteratorByState;
+import org.broadinstitute.sting.gatk.iterators.GATKSAMIterator;
import org.broadinstitute.sting.gatk.iterators.StingSAMIterator;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.GenomeLocParser;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
+import org.broadinstitute.sting.utils.locusiterator.LocusIterator;
+import org.broadinstitute.sting.utils.locusiterator.LocusIteratorByState;
+import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
import java.util.Collection;
import java.util.Iterator;
@@ -71,7 +72,7 @@ public class WindowMaker implements Iterable, I
/**
* Hold the read iterator so that it can be closed later.
*/
- private final StingSAMIterator readIterator;
+ private final GATKSAMIterator readIterator;
/**
* The data source for reads. Will probably come directly from the BAM file.
@@ -104,22 +105,23 @@ public class WindowMaker implements Iterable, I
* @param sampleNames The complete set of sample names in the reads in shard
*/
+ private final LocusIteratorByState libs;
+
public WindowMaker(Shard shard, GenomeLocParser genomeLocParser, StingSAMIterator iterator, List intervals, Collection sampleNames) {
this.sourceInfo = shard.getReadProperties();
- this.readIterator = iterator;
+ this.readIterator = new GATKSAMIterator(iterator);
// Use the legacy version of LocusIteratorByState if legacy downsampling was requested:
- this.sourceIterator = sourceInfo.getDownsamplingMethod().useLegacyDownsampler ?
- new PeekableIterator(new LegacyLocusIteratorByState(iterator,sourceInfo,genomeLocParser,sampleNames))
- :
- new PeekableIterator(new LocusIteratorByState(iterator,sourceInfo,genomeLocParser,sampleNames));
-
+ if ( sourceInfo.getDownsamplingMethod().useLegacyDownsampler )
+ throw new IllegalArgumentException("legacy downsampler no longer supported in the window maker");
+ this.libs = new LocusIteratorByState(readIterator,sourceInfo,genomeLocParser,sampleNames);
+ this.sourceIterator = new PeekableIterator(libs);
this.intervalIterator = intervals.size()>0 ? new PeekableIterator(intervals.iterator()) : null;
}
public WindowMaker(Shard shard, GenomeLocParser genomeLocParser, StingSAMIterator iterator, List intervals ) {
- this(shard, genomeLocParser, iterator, intervals, LegacyLocusIteratorByState.sampleListForSAMWithoutReadGroups());
+ this(shard, genomeLocParser, iterator, intervals, LocusIteratorByState.sampleListForSAMWithoutReadGroups());
}
public Iterator iterator() {
@@ -209,5 +211,10 @@ public class WindowMaker implements Iterable, I
throw new ReviewedStingException("BUG: filtering locus does not contain, is not before, and is not past the given alignment context");
}
}
+
+ @Override
+ public LocusIteratorByState getLIBS() {
+ return libs;
+ }
}
}
diff --git a/public/java/src/org/broadinstitute/sting/gatk/iterators/GATKSAMIterator.java b/public/java/src/org/broadinstitute/sting/gatk/iterators/GATKSAMIterator.java
new file mode 100644
index 000000000..30a520e09
--- /dev/null
+++ b/public/java/src/org/broadinstitute/sting/gatk/iterators/GATKSAMIterator.java
@@ -0,0 +1,57 @@
+/*
+ * Copyright (c) 2012 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.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+ * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+ * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+ * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+ * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
+ * THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+package org.broadinstitute.sting.gatk.iterators;
+
+import net.sf.samtools.SAMRecord;
+import net.sf.samtools.util.CloseableIterator;
+import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
+
+import java.util.Iterator;
+
+/**
+ * Temporarily hack to convert SAMRecords to GATKSAMRecords
+ *
+ * User: depristo
+ * Date: 1/11/13
+ * Time: 1:19 PM
+ */
+public class GATKSAMIterator implements CloseableIterator, Iterable {
+ final CloseableIterator it;
+
+ public GATKSAMIterator(final CloseableIterator it) {
+ this.it = it;
+ }
+
+ public GATKSAMIterator(final StingSAMIterator it) {
+ this.it = it;
+ }
+
+ @Override public boolean hasNext() { return it.hasNext(); }
+ @Override public GATKSAMRecord next() { return (GATKSAMRecord)it.next(); }
+ @Override public void remove() { it.remove(); }
+ @Override public void close() { it.close(); }
+ @Override public Iterator iterator() { return this; }
+}
diff --git a/public/java/src/org/broadinstitute/sting/gatk/iterators/LegacyLocusIteratorByState.java b/public/java/src/org/broadinstitute/sting/gatk/iterators/LegacyLocusIteratorByState.java
deleted file mode 100644
index e4d2fcefc..000000000
--- a/public/java/src/org/broadinstitute/sting/gatk/iterators/LegacyLocusIteratorByState.java
+++ /dev/null
@@ -1,962 +0,0 @@
-/*
-* Copyright (c) 2012 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.
-*
-* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
-* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
-* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
-* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
-* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
-* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
-* THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-*/
-
-package org.broadinstitute.sting.gatk.iterators;
-
-import net.sf.picard.util.PeekableIterator;
-import net.sf.samtools.Cigar;
-import net.sf.samtools.CigarElement;
-import net.sf.samtools.CigarOperator;
-import net.sf.samtools.SAMRecord;
-import org.apache.log4j.Logger;
-import org.broadinstitute.sting.gatk.downsampling.DownsampleType;
-import org.broadinstitute.sting.gatk.downsampling.DownsamplingMethod;
-import org.broadinstitute.sting.gatk.ReadProperties;
-import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
-import org.broadinstitute.sting.utils.GenomeLoc;
-import org.broadinstitute.sting.utils.GenomeLocParser;
-import org.broadinstitute.sting.utils.MathUtils;
-import org.broadinstitute.sting.utils.LegacyReservoirDownsampler;
-import org.broadinstitute.sting.utils.exceptions.UserException;
-import org.broadinstitute.sting.utils.pileup.PileupElement;
-import org.broadinstitute.sting.utils.pileup.ReadBackedPileupImpl;
-import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
-import org.broadinstitute.sting.utils.sam.ReadUtils;
-
-import java.util.*;
-
-/**
- * Iterator that traverses a SAM File, accumulating information on a per-locus basis
- */
-public class LegacyLocusIteratorByState extends LocusIterator {
- /**
- * our log, which we want to capture anything from this class
- */
- private static Logger logger = Logger.getLogger(LegacyLocusIteratorByState.class);
-
- // -----------------------------------------------------------------------------------------------------------------
- //
- // member fields
- //
- // -----------------------------------------------------------------------------------------------------------------
-
- /**
- * Used to create new GenomeLocs.
- */
- private final GenomeLocParser genomeLocParser;
- private final ArrayList samples;
- private final ReadStateManager readStates;
-
- static private class SAMRecordState {
- SAMRecord read;
- int readOffset = -1; // how far are we offset from the start of the read bases?
- int genomeOffset = -1; // how far are we offset from the alignment start on the genome?
-
- Cigar cigar = null;
- int cigarOffset = -1;
- CigarElement curElement = null;
- int nCigarElements = 0;
-
- int cigarElementCounter = -1; // how far are we into a single cigarElement
-
- // The logical model for generating extended events is as follows: the "record state" implements the traversal
- // along the reference; thus stepForwardOnGenome() returns on every and only on actual reference bases. This
- // can be a (mis)match or a deletion (in the latter case, we still return on every individual reference base the
- // deletion spans). In the extended events mode, the record state also remembers if there was an insertion, or
- // if the deletion just started *right before* the current reference base the record state is pointing to upon the return from
- // stepForwardOnGenome(). The next call to stepForwardOnGenome() will clear that memory (as we remember only extended
- // events immediately preceding the current reference base).
-
- public SAMRecordState(SAMRecord read) {
- this.read = read;
- cigar = read.getCigar();
- nCigarElements = cigar.numCigarElements();
-
- //System.out.printf("Creating a SAMRecordState: %s%n", this);
- }
-
- public SAMRecord getRead() {
- return read;
- }
-
- /**
- * What is our current offset in the read's bases that aligns us with the reference genome?
- *
- * @return
- */
- public int getReadOffset() {
- return readOffset;
- }
-
- /**
- * What is the current offset w.r.t. the alignment state that aligns us to the readOffset?
- *
- * @return
- */
- public int getGenomeOffset() {
- return genomeOffset;
- }
-
- public int getGenomePosition() {
- return read.getAlignmentStart() + getGenomeOffset();
- }
-
- public GenomeLoc getLocation(GenomeLocParser genomeLocParser) {
- return genomeLocParser.createGenomeLoc(read.getReferenceName(), getGenomePosition());
- }
-
- public CigarOperator getCurrentCigarOperator() {
- return curElement.getOperator();
- }
-
- public String toString() {
- return String.format("%s ro=%d go=%d co=%d cec=%d %s", read.getReadName(), readOffset, genomeOffset, cigarOffset, cigarElementCounter, curElement);
- }
-
- public CigarElement peekForwardOnGenome() {
- return ( cigarElementCounter + 1 > curElement.getLength() && cigarOffset + 1 < nCigarElements ? cigar.getCigarElement(cigarOffset + 1) : curElement );
- }
-
- public CigarElement peekBackwardOnGenome() {
- return ( cigarElementCounter - 1 == 0 && cigarOffset - 1 > 0 ? cigar.getCigarElement(cigarOffset - 1) : curElement );
- }
-
-
- public CigarOperator stepForwardOnGenome() {
- // we enter this method with readOffset = index of the last processed base on the read
- // (-1 if we did not process a single base yet); this can be last matching base, or last base of an insertion
-
-
- if (curElement == null || ++cigarElementCounter > curElement.getLength()) {
- cigarOffset++;
- if (cigarOffset < nCigarElements) {
- curElement = cigar.getCigarElement(cigarOffset);
- cigarElementCounter = 0;
- // next line: guards against cigar elements of length 0; when new cigar element is retrieved,
- // we reenter in order to re-check cigarElementCounter against curElement's length
- return stepForwardOnGenome();
- } else {
- if (curElement != null && curElement.getOperator() == CigarOperator.D)
- throw new UserException.MalformedBAM(read, "read ends with deletion. Cigar: " + read.getCigarString() + ". Although the SAM spec technically permits such reads, this is often indicative of malformed files. If you are sure you want to use this file, re-run your analysis with the extra option: -rf BadCigar");
-
- // Reads that contain indels model the genomeOffset as the following base in the reference. Because
- // we fall into this else block only when indels end the read, increment genomeOffset such that the
- // current offset of this read is the next ref base after the end of the indel. This position will
- // model a point on the reference somewhere after the end of the read.
- genomeOffset++; // extended events need that. Logically, it's legal to advance the genomic offset here:
- // we do step forward on the ref, and by returning null we also indicate that we are past the read end.
-
- return null;
- }
- }
-
- boolean done = false;
- switch (curElement.getOperator()) {
- case H: // ignore hard clips
- case P: // ignore pads
- cigarElementCounter = curElement.getLength();
- break;
- case I: // insertion w.r.t. the reference
- case S: // soft clip
- cigarElementCounter = curElement.getLength();
- readOffset += curElement.getLength();
- break;
- case D: // deletion w.r.t. the reference
- if (readOffset < 0) // we don't want reads starting with deletion, this is a malformed cigar string
- throw new UserException.MalformedBAM(read, "read starts with deletion. Cigar: " + read.getCigarString() + ". Although the SAM spec technically permits such reads, this is often indicative of malformed files. If you are sure you want to use this file, re-run your analysis with the extra option: -rf BadCigar");
- // should be the same as N case
- genomeOffset++;
- done = true;
- break;
- case N: // reference skip (looks and gets processed just like a "deletion", just different logical meaning)
- genomeOffset++;
- done = true;
- break;
- case M:
- case EQ:
- case X:
- readOffset++;
- genomeOffset++;
- done = true;
- break;
- default:
- throw new IllegalStateException("Case statement didn't deal with cigar op: " + curElement.getOperator());
- }
-
- return done ? curElement.getOperator() : stepForwardOnGenome();
- }
- }
-
- //final boolean DEBUG = false;
- //final boolean DEBUG2 = false && DEBUG;
- private ReadProperties readInfo;
- private AlignmentContext nextAlignmentContext;
-
- // -----------------------------------------------------------------------------------------------------------------
- //
- // constructors and other basic operations
- //
- // -----------------------------------------------------------------------------------------------------------------
-
- public LegacyLocusIteratorByState(final Iterator samIterator, ReadProperties readInformation, GenomeLocParser genomeLocParser, Collection