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

Merge pull request #200 from broadinstitute/gda_ug_rr_bug_48742591 

Fix for indel calling with UG in presence of reduced reads: When a read ...
This commit is contained in:
Eric Banks 2013-04-29 17:54:43 -07:00
parent c5701a9ade 20d3137928
commit a3a2ec5a1c
4 changed files with 15 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
protected/java/src/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModel.java
View File

@ -245,8 +245,13 @@ public class PairHMMIndelErrorModel {
}
}
else {
final int refWindowStart = ref.getWindow().getStart();
final int refWindowStop = ref.getWindow().getStop();
// extra padding on candidate haplotypes to make sure reads are always strictly contained
// in them - a value of 1 will in theory do but we use a slightly higher one just for safety sake, mostly
// in case bases at edge of reads have lower quality.
final int trailingBases = 3;
final int extraOffset = Math.abs(eventLength);
final int refWindowStart = ref.getWindow().getStart()+(trailingBases+extraOffset);
final int refWindowStop = ref.getWindow().getStop()-(trailingBases+extraOffset);
if (DEBUG) {
System.out.format("Read Name:%s, aln start:%d aln stop:%d orig cigar:%s\n",p.getRead().getReadName(), p.getRead().getAlignmentStart(), p.getRead().getAlignmentEnd(), p.getRead().getCigarString());
@ -255,10 +260,10 @@ public class PairHMMIndelErrorModel {
GATKSAMRecord read = ReadClipper.hardClipAdaptorSequence(p.getRead());
if (!read.isEmpty() && (read.getSoftEnd() > refWindowStop && read.getSoftStart() < refWindowStop))
read = ReadClipper.hardClipByReferenceCoordinatesRightTail(read, ref.getWindow().getStop());
read = ReadClipper.hardClipByReferenceCoordinatesRightTail(read, refWindowStop);
if (!read.isEmpty() && (read.getSoftStart() < refWindowStart && read.getSoftEnd() > refWindowStart))
read = ReadClipper.hardClipByReferenceCoordinatesLeftTail (read, ref.getWindow().getStart());
read = ReadClipper.hardClipByReferenceCoordinatesLeftTail (read, refWindowStart);
if (read.isEmpty())
continue;
@ -270,7 +275,6 @@ public class PairHMMIndelErrorModel {
continue;
// get bases of candidate haplotypes that overlap with reads
final int trailingBases = 3;
final long readStart = read.getSoftStart();
final long readEnd = read.getSoftEnd();
@ -286,7 +290,6 @@ public class PairHMMIndelErrorModel {
final int numEndSoftClippedBases = softClips ? read.getSoftEnd()- read.getAlignmentEnd() : 0 ;
final byte [] unclippedReadBases = read.getReadBases();
final byte [] unclippedReadQuals = read.getBaseQualities();
final int extraOffset = Math.abs(eventLength);
/**
* Compute genomic locations that candidate haplotypes will span.
@ -313,6 +316,7 @@ public class PairHMMIndelErrorModel {
startLocationInRefForHaplotypes = ref.getWindow().getStop(); // read starts after haplotype: read will have to be clipped completely;
}
// candidate haplotype cannot go beyond reference context
if (stopLocationInRefForHaplotypes > ref.getWindow().getStop()) {
stopLocationInRefForHaplotypes = ref.getWindow().getStop(); // check also if end of read will go beyond reference context
}

4
protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIndelCallingIntegrationTest.java
View File

@ -101,7 +101,7 @@ public class UnifiedGenotyperIndelCallingIntegrationTest extends WalkerTest {
" -o %s" +
" -L 1:10,000,000-10,500,000",
1,
Arrays.asList("54e13f696f56eb742bf449ad11d0dc5f"));
Arrays.asList("8d9b8f8a1479322961c840e461b6dba8"));
executeTest(String.format("test indel calling, multiple technologies"), spec);
}
@ -136,7 +136,7 @@ public class UnifiedGenotyperIndelCallingIntegrationTest extends WalkerTest {
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 " + result.get(0).getAbsolutePath(), 1,
Arrays.asList("e87f5c76661527ef7aa44e528fe19573"));
Arrays.asList("3d4d66cc253eac55f16e5b0a36f17d8d"));
executeTest("test MultiSample Pilot1 CEU indels using GENOTYPE_GIVEN_ALLELES", spec2);
}

4
protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java
View File

@ -240,7 +240,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
" -o %s" +
" -L 1:10,000,000-10,100,000",
1,
Arrays.asList("d995b76adc3766889f7c2c88da14055c"));
Arrays.asList("31be725b2a7c15e9769391ad940c0587"));
executeTest(String.format("test multiple technologies"), spec);
}
@ -259,7 +259,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
" -L 1:10,000,000-10,100,000" +
" -baq CALCULATE_AS_NECESSARY",
1,
Arrays.asList("9669e1643d22d5ad047b3941aeefd6db"));
Arrays.asList("dcc5cec42730567982def16da4a7f286"));
executeTest(String.format("test calling with BAQ"), spec);
}

2
protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperNormalCallingIntegrationTest.java
View File

@ -64,7 +64,7 @@ public class UnifiedGenotyperNormalCallingIntegrationTest extends WalkerTest{
public void testMultiSamplePilot1() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
baseCommand + " -I " + validationDataLocation + "low_coverage_CEU.chr1.10k-11k.bam -o %s -L 1:10,022,000-10,025,000", 1,
Arrays.asList("e3efd1917192ea743ac1e9958aa0a98f"));
Arrays.asList("a6c224235c21b4af816b1512eb0624df"));
executeTest("test MultiSample Pilot1", spec);
}