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

fixing oddities in duplicates 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@628 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
depristo 2009-05-07 18:55:45 +00:00
parent 4a26f35caa
commit 5a6892900e
3 changed files with 17 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
java/src/org/broadinstitute/sting/gatk/traversals/TraverseDuplicates.java
View File

@ -115,7 +115,12 @@ public class TraverseDuplicates extends TraversalEngine {
final GenomeLoc readLoc = new GenomeLoc(read);
final GenomeLoc readMateLoc = new GenomeLoc(read.getMateReferenceIndex(), read.getMateAlignmentStart(), read.getMateAlignmentStart());
if (DEBUG) logger.debug(String.format("Examining reads at %s vs. %s at %s / %s vs. %s / %s%n", key.getReadName(), read.getReadName(), keyLoc, keyMateLoc, readLoc, readMateLoc));
if ( readLoc.compareTo(keyLoc) == 0 && readMateLoc.compareTo(keyMateLoc) == 0 ) {
// read and key start at the same place, and either the this read and the key
// share a mate location or the read is flagged as a duplicate
if ( readLoc.compareTo(keyLoc) == 0 &&
( readMateLoc.compareTo(keyMateLoc) == 0) ||
read.getDuplicateReadFlag() ) {
// we are at the same position as the dup and have the same mat pos, it's a dup
if (DEBUG) logger.debug(String.format(" => Adding read to dups list: %s%n", read));
dups.add(read);
@ -157,8 +162,11 @@ public class TraverseDuplicates extends TraversalEngine {
List<SAMRecord> uniqueReads = split.getFirst();
List<SAMRecord> duplicateReads = split.getSecond();
logger.debug(String.format("*** TraverseDuplicates.traverse at %s has %d unique and %d duplicate reads",
site, uniqueReads.size(), duplicateReads.size()));
logger.debug(String.format("*** TraverseDuplicates.traverse at %s with %d reads has %d unique and %d duplicate reads",
site, reads.size(), uniqueReads.size(), duplicateReads.size()));
if ( reads.size() != uniqueReads.size() + duplicateReads.size() )
throw new RuntimeException(String.format("Bug occurred spliting reads [N=%d] at loc %s into unique [N=%d] and duplicates [N=%d], sizes don't match",
reads.size(), uniqueReads.size(), duplicateReads.size()));
// Jump forward in the reference to this locus location
LocusContext locus = new LocusContext(site, duplicateReads, Arrays.asList(0));

6
java/src/org/broadinstitute/sting/playground/gatk/walkers/CombineDuplicatesWalker.java
View File

@ -92,7 +92,7 @@ public class CombineDuplicatesWalker extends DuplicateWalker<SAMRecord, SAMFileW
//logger.info(String.format("%s has %d duplicates%n", loc, duplicateReads.size()));
SAMRecord combinedRead = null;
if ( duplicateReads.size() == 1 ) {
if ( duplicateReads.size() == 1 && ! duplicateReads.get(0).getDuplicateReadFlag() ) {
// we are a unique read
combinedRead = duplicateReads.get(0);
} else {
@ -102,6 +102,10 @@ public class CombineDuplicatesWalker extends DuplicateWalker<SAMRecord, SAMFileW
//}
combinedRead = DupUtils.combineDuplicates(duplicateReads, MAX_QUALITY_SCORE);
}
if ( combinedRead.getDuplicateReadFlag() )
throw new RuntimeException(String.format("Combined read %s [of %d] is a duplicate after combination -- this is a bug%n%s",
combinedRead.getReadName(), duplicateReads.size(), combinedRead.format()));
return combinedRead;
}

1
java/src/org/broadinstitute/sting/utils/duplicates/DupUtils.java
View File

@ -114,6 +114,7 @@ public class DupUtils {
// make the combined read by copying the first read and setting the
// bases and quals to new arrays
SAMRecord comb = tmpCopyRead(duplicates.get(0));
comb.setDuplicateReadFlag(false);
int readLen = comb.getReadBases().length;
byte[] bases = new byte[readLen];
byte[] quals = new byte[readLen];