Permalink
Browse files

Merge pull request #199 from broadinstitute/md_clipped_reduced_reads

Bugfix for ReadClipper with ReducedReads
  • Loading branch information...
2 parents 5dd73ba + 0387ea8 commit c5701a9adef72fef08b5ab57f0d66205a9afe4de @eitanbanks eitanbanks committed Apr 29, 2013
@@ -188,7 +188,7 @@ public void HCTestDoesNotFailOnBadRefBase() {
public void HCTestReducedBam() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
"-T HaplotypeCaller --disableDithering -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("3c87eb93ffe3a0166aca753050b981e1"));
+ Arrays.asList("0df626cd0d76aca8a05a545d0b36bf23"));
executeTest("HC calling on a ReducedRead BAM", spec);
}
@@ -378,7 +378,13 @@ private GATKSAMRecord hardClip(GATKSAMRecord read, int start, int stop) {
hardClippedRead.setBaseQualities(newBaseInsertionQuals, EventType.BASE_INSERTION);
hardClippedRead.setBaseQualities(newBaseDeletionQuals, EventType.BASE_DELETION);
}
-
+
+ if (read.isReducedRead()) {
+ final byte[] reducedCounts = new byte[newLength];
+ System.arraycopy(read.getReducedReadCounts(), copyStart, reducedCounts, 0, newLength);
+ hardClippedRead.setReducedReadCounts(reducedCounts);
+ }
+
return hardClippedRead;
}
@@ -394,6 +394,22 @@ public void setReducedReadCountsTag(final byte[] counts) {
}
/**
+ * Set the reduced read counts tag for this record to counts
+ *
+ * Note that this function does not set the REDUCED_READ_CONSENSUS_TAG value, it's purely for manipulating
+ * the underlying reduced reads count
+ *
+ * TODO -- this function needs to be fixed when the RR spec is set to 2.0
+ *
+ * @param counts the count array
+ */
+ public void setReducedReadCounts(final byte[] counts) {
+ if ( counts.length != getReadBases().length ) throw new IllegalArgumentException("Reduced counts length " + counts.length + " != bases length " + getReadBases().length);
+ retrievedReduceReadCounts = true;
+ reducedReadCounts = counts;
+ }
+
+ /**
* The number of bases corresponding the i'th base of the reduced read.
*
* @param i the read based coordinate inside the read
@@ -678,6 +694,7 @@ public static GATKSAMRecord emptyRead(GATKSAMRecord read) {
emptyRead.setCigarString("");
emptyRead.setReadBases(new byte[0]);
emptyRead.setBaseQualities(new byte[0]);
+ if ( read.isReducedRead() ) emptyRead.setReducedReadCounts(new byte[0]);
SAMReadGroupRecord samRG = read.getReadGroup();
emptyRead.clearAttributes();
@@ -35,6 +35,7 @@
import org.testng.annotations.BeforeClass;
import org.testng.annotations.Test;
+import java.util.Arrays;
import java.util.HashMap;
import java.util.List;
@@ -348,4 +349,20 @@ public boolean assertHardClippingSoftClips(CigarCounter clipped) {
}
+ @Test(enabled = true)
+ public void testHardClipReducedRead() {
+ GATKSAMRecord read = ReadClipperTestUtils.makeReadFromCigar("10M");
+ final byte[] counts = new byte[read.getReadLength()];
+ for ( int i = 0; i < counts.length; i++ ) counts[i] = (byte)i;
+ read.setReducedReadCounts(counts);
+ int alnStart = read.getAlignmentStart();
+ int alnEnd = read.getAlignmentEnd();
+ int readLength = read.getReadLength();
+ for (int i = 0; i < readLength / 2; i++) {
+ GATKSAMRecord clippedRead = ReadClipper.hardClipBothEndsByReferenceCoordinates(read, alnStart + i, alnEnd - i);
+ final byte[] expectedReducedCounts = Arrays.copyOfRange(counts, i + 1, readLength - i - 1);
+ Assert.assertEquals(clippedRead.getReducedReadCounts(), expectedReducedCounts);
+ }
+ }
+
}

0 comments on commit c5701a9

Please sign in to comment.