zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'master' of ssh://gsa1/humgen/gsa-scr1/gsa-engineering/git/unstable

This commit is contained in:
Mark DePristo 2011-09-22 17:26:21 -04:00
parent dab7232e9a 96c875399c
commit c49cc623de
3 changed files with 38 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
public/java/src/org/broadinstitute/sting/utils/clipreads/ClippingOp.java
View File

@ -432,25 +432,37 @@ public class ClippingOp {
} }
private int calculateAlignmentStartShift(Cigar oldCigar, Cigar newCigar) { private int calculateAlignmentStartShift(Cigar oldCigar, Cigar newCigar) {
int shift = 0; int newShift = 0;
int oldShift = 0;
int deletionShift = 0;
// Rewind to previous start (by counting everything that was already clipped in this read)
for (CigarElement cigarElement : oldCigar.getCigarElements()) {
if (!cigarElement.getOperator().consumesReferenceBases())
shift -= cigarElement.getLength();
else
break;
}
// Advance to new start (by counting everything new that has been clipped )
for (CigarElement cigarElement : newCigar.getCigarElements()) { for (CigarElement cigarElement : newCigar.getCigarElements()) {
if (!cigarElement.getOperator().consumesReferenceBases()) if (cigarElement.getOperator() == CigarOperator.HARD_CLIP || cigarElement.getOperator() == CigarOperator.SOFT_CLIP)
shift += cigarElement.getLength(); newShift += cigarElement.getLength();
else else
break; break;
} }
return shift; for (CigarElement cigarElement : oldCigar.getCigarElements()) {
if (cigarElement.getOperator() == CigarOperator.HARD_CLIP || cigarElement.getOperator() == CigarOperator.SOFT_CLIP )
oldShift += Math.min(cigarElement.getLength(), newShift - oldShift);
else
break;
}
int basesClipped = 0;
for (CigarElement cigarElement : oldCigar.getCigarElements()) {
if (basesClipped > newShift) // are we beyond the clipped region?
break;
else if (cigarElement.getOperator() == CigarOperator.DELETION) // if this is a deletion, we have to adjust the starting shift
deletionShift += cigarElement.getLength();
else
basesClipped += cigarElement.getLength();
}
return newShift - oldShift + deletionShift;
} }
private int calculateHardClippingAlignmentShift(CigarElement cigarElement, int clippedLength) { private int calculateHardClippingAlignmentShift(CigarElement cigarElement, int clippedLength) {
@ -462,8 +474,8 @@ public class ClippingOp {
return -clippedLength; return -clippedLength;
} }
if (cigarElement.getOperator() == CigarOperator.DELETION) // if (cigarElement.getOperator() == CigarOperator.DELETION)
return cigarElement.getLength(); // return cigarElement.getLength();
return 0; return 0;
} }

20
public/java/src/org/broadinstitute/sting/utils/clipreads/ReadClipper.java
View File

@ -71,21 +71,21 @@ public class ReadClipper {
private SAMRecord hardClipByReferenceCoordinates(int refStart, int refStop) { private SAMRecord hardClipByReferenceCoordinates(int refStart, int refStop) {
int start = (refStart < 0) ? 0 : ReadUtils.getReadCoordinateForReferenceCoordinate(read, refStart); int start = (refStart < 0) ? 0 : ReadUtils.getReadCoordinateForReferenceCoordinate(read, refStart);
int stop = (refStop < 0) ? read.getReadLength() - 1: ReadUtils.getReadCoordinateForReferenceCoordinate(read, refStop); int stop = (refStop < 0) ? read.getReadLength() - 1 : ReadUtils.getReadCoordinateForReferenceCoordinate(read, refStop);
if (start < 0 || stop > read.getReadLength() - 1 + numDeletions(read)) if (start < 0 || stop > read.getReadLength() - 1)
throw new ReviewedStingException("Trying to clip before the start or after the end of a read"); throw new ReviewedStingException("Trying to clip before the start or after the end of a read");
// TODO add check in the Hardclip function if ( start > stop ) {
if ( start > stop ) // stop = ReadUtils.getReadCoordinateForReferenceCoordinate(read, ReadUtils.getRefCoordSoftUnclippedEnd(read));
stop = ReadUtils.getReadCoordinateForReferenceCoordinate(read, ReadUtils.getRefCoordSoftUnclippedEnd(read)); throw new ReviewedStingException("START > STOP -- this should never happen -- call Mauricio!");
}
//This tries to fix the bug where the deletion is counted a read base and as a result, the hardCLipper runs into //This tries to fix the bug where the deletion is counted a read base and as a result, the hardCLipper runs into
//an endless loop when hard clipping the cigar string because the read coordinates are not covered by the read //an endless loop when hard clipping the cigar string because the read coordinates are not covered by the read
stop -= numDeletions(read); // stop -= numDeletions(read);
if ( start > stop ) // if ( start > stop )
start -= numDeletions(read); // start -= numDeletions(read);
//System.out.println("Clipping start/stop: " + start + "/" + stop); //System.out.println("Clipping start/stop: " + start + "/" + stop);
@ -145,7 +145,7 @@ public class ReadClipper {
cutLeft = readIndex + cigarElement.getLength() - 1; cutLeft = readIndex + cigarElement.getLength() - 1;
} }
} }
else else if (cigarElement.getOperator() != CigarOperator.HARD_CLIP)
rightTail = true; rightTail = true;
if (cigarElement.getOperator().consumesReadBases()) if (cigarElement.getOperator().consumesReadBases())

3
public/java/test/org/broadinstitute/sting/utils/genotype/vcf/VCFWriterUnitTest.java
View File

@ -105,7 +105,6 @@ public class VCFWriterUnitTest extends BaseTest {
public static VCFHeader createFakeHeader(Set<VCFHeaderLine> metaData, Set<String> additionalColumns) { public static VCFHeader createFakeHeader(Set<VCFHeaderLine> metaData, Set<String> additionalColumns) {
metaData.add(new VCFHeaderLine(VCFHeaderVersion.VCF4_0.getFormatString(), VCFHeaderVersion.VCF4_0.getVersionString())); metaData.add(new VCFHeaderLine(VCFHeaderVersion.VCF4_0.getFormatString(), VCFHeaderVersion.VCF4_0.getVersionString()));
metaData.add(new VCFHeaderLine("two", "2")); metaData.add(new VCFHeaderLine("two", "2"));
additionalColumns.add("FORMAT");
additionalColumns.add("extra1"); additionalColumns.add("extra1");
additionalColumns.add("extra2"); additionalColumns.add("extra2");
return new VCFHeader(metaData, additionalColumns); return new VCFHeader(metaData, additionalColumns);
@ -159,6 +158,6 @@ public class VCFWriterUnitTest extends BaseTest {
Assert.assertTrue(additionalColumns.contains(key)); Assert.assertTrue(additionalColumns.contains(key));
index++; index++;
} }
Assert.assertEquals(index+1, additionalColumns.size() /* for the header field we don't see */); Assert.assertEquals(index, additionalColumns.size());
} }
} }