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

Added in a check for what would be an empty allele after trimming.

This commit is contained in:
Eric Banks 2014-01-14 12:56:54 -05:00
parent 201ad398ac
commit 9f1ab0087a
2 changed files with 37 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -205,10 +205,21 @@ public class PairHMMIndelErrorModel {
}
}
private LinkedHashMap<Allele, Haplotype> trimHaplotypes(final LinkedHashMap<Allele, Haplotype> haplotypeMap,
long startLocationInRefForHaplotypes,
long stopLocationInRefForHaplotypes,
final ReferenceContext ref){
/**
* Trims the haplotypes in the given map to the provided start/stop.
*
* @param haplotypeMap the input map
* @param startLocationInRefForHaplotypes the start location of the trim
* @param stopLocationInRefForHaplotypes the stop location of the trim
* @param ref the reference context (used for debugging only, so can be null)
* @return a non-null mapping corresponding to the trimmed version of the original;
* some elements may be lost if trimming cannot be performed on them (e.g. they fall outside of the region to keep)
*/
protected static Map<Allele, Haplotype> trimHaplotypes(final Map<Allele, Haplotype> haplotypeMap,
long startLocationInRefForHaplotypes,
long stopLocationInRefForHaplotypes,
final ReferenceContext ref) {
if ( haplotypeMap == null ) throw new IllegalArgumentException("The input allele to haplotype map cannot be null");
final LinkedHashMap<Allele, Haplotype> trimmedHaplotypeMap = new LinkedHashMap<>();
for (final Allele a: haplotypeMap.keySet()) {
@ -225,10 +236,13 @@ public class PairHMMIndelErrorModel {
final long indStart = startLocationInRefForHaplotypes - haplotype.getStartPosition();
final long indStop = stopLocationInRefForHaplotypes - haplotype.getStartPosition();
if ( indStart >= indStop )
continue;
if (DEBUG)
System.out.format("indStart: %d indStop: %d WinStart:%d WinStop:%d start: %d stop: %d\n",
indStart, indStop, ref.getWindow().getStart(), ref.getWindow().getStop(), startLocationInRefForHaplotypes, stopLocationInRefForHaplotypes);
// commented out here because we need to make this method static for unit testing
//if (DEBUG)
// System.out.format("indStart: %d indStop: %d WinStart:%d WinStop:%d start: %d stop: %d\n",
// indStart, indStop, ref.getWindow().getStart(), ref.getWindow().getStop(), startLocationInRefForHaplotypes, stopLocationInRefForHaplotypes);
// get the trimmed haplotype-bases array and create a new haplotype based on it. Pack this into the new map
final byte[] trimmedHaplotypeBases = Arrays.copyOfRange(haplotype.getBases(), (int)indStart, (int)indStop);

17
protected/java/test/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModelUnitTest.java
View File

@ -50,9 +50,12 @@ package org.broadinstitute.sting.gatk.walkers.indels;
import net.sf.picard.reference.IndexedFastaSequenceFile;
import net.sf.samtools.SAMFileHeader;
import org.broadinstitute.sting.BaseTest;
import org.broadinstitute.sting.utils.UnvalidatingGenomeLoc;
import org.broadinstitute.sting.utils.fasta.CachingIndexedFastaSequenceFile;
import org.broadinstitute.sting.utils.haplotype.Haplotype;
import org.broadinstitute.sting.utils.sam.ArtificialSAMUtils;
import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
import org.broadinstitute.variant.variantcontext.Allele;
import org.testng.Assert;
import org.testng.annotations.BeforeClass;
import org.testng.annotations.DataProvider;
@ -138,4 +141,16 @@ public class PairHMMIndelErrorModelUnitTest extends BaseTest {
Assert.assertEquals(PairHMMIndelErrorModel.mustClipDownstream(read, 13), true);
Assert.assertEquals(PairHMMIndelErrorModel.mustClipDownstream(read, 14), false);
}
}
@Test
public void trimHaplotypesToNullAlleleTest() {
// we need a case where start and stop > haplotype coordinates
final int start = 100, stop = 100;
final Haplotype h = new Haplotype(new byte[]{(byte)'A'}, new UnvalidatingGenomeLoc("1", 0, 10, 10));
final Map<Allele, Haplotype> input = new HashMap<Allele, Haplotype>(1);
input.put(Allele.create("A"), h);
final Map<Allele, Haplotype> output = PairHMMIndelErrorModel.trimHaplotypes(input, start, stop, null);
Assert.assertTrue(output.isEmpty());
}
}