Minor refactoring of the cleaner allows me to add a trivial walker that left aligns the indels present in reads.

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@4581 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
ebanks 2010-10-27 03:39:10 +00:00
parent 04e38929f0
commit ec442086ec
2 changed files with 109 additions and 18 deletions

View File

@ -1071,7 +1071,7 @@ public class IndelRealigner extends ReadWalker<Integer, Integer> {
return reduces;
}
private static Cigar unclipCigar(Cigar cigar) {
protected static Cigar unclipCigar(Cigar cigar) {
ArrayList<CigarElement> elements = new ArrayList<CigarElement>(cigar.numCigarElements());
for ( CigarElement ce : cigar.getCigarElements() ) {
if ( !isClipOperator(ce.getOperator()) )
@ -1084,6 +1084,26 @@ public class IndelRealigner extends ReadWalker<Integer, Integer> {
return op == CigarOperator.S || op == CigarOperator.H || op == CigarOperator.P;
}
protected static Cigar reclipCigar(Cigar cigar, SAMRecord read) {
ArrayList<CigarElement> elements = new ArrayList<CigarElement>();
int i = 0;
int n = read.getCigar().numCigarElements();
while ( i < n && isClipOperator(read.getCigar().getCigarElement(i).getOperator()) )
elements.add(read.getCigar().getCigarElement(i++));
elements.addAll(cigar.getCigarElements());
i++;
while ( i < n && !isClipOperator(read.getCigar().getCigarElement(i).getOperator()) )
i++;
while ( i < n && isClipOperator(read.getCigar().getCigarElement(i).getOperator()) )
elements.add(read.getCigar().getCigarElement(i++));
return new Cigar(elements);
}
private class AlignedRead {
private final SAMRecord read;
private byte[] readBases = null;
@ -1191,23 +1211,7 @@ public class IndelRealigner extends ReadWalker<Integer, Integer> {
// pull out the bases that aren't clipped out
private Cigar reclipCigar(Cigar cigar) {
ArrayList<CigarElement> elements = new ArrayList<CigarElement>();
int i = 0;
int n = read.getCigar().numCigarElements();
while ( i < n && isClipOperator(read.getCigar().getCigarElement(i).getOperator()) )
elements.add(read.getCigar().getCigarElement(i++));
elements.addAll(cigar.getCigarElements());
i++;
while ( i < n && !isClipOperator(read.getCigar().getCigarElement(i).getOperator()) )
i++;
while ( i < n && isClipOperator(read.getCigar().getCigarElement(i).getOperator()) )
elements.add(read.getCigar().getCigarElement(i++));
return new Cigar(elements);
return IndelRealigner.reclipCigar(cigar, read);
}
// tentatively sets the new start, but it needs to be confirmed later

View File

@ -0,0 +1,87 @@
/*
* Copyright (c) 2010.
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
* THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.gatk.walkers.indels;
import net.sf.samtools.*;
import org.broadinstitute.sting.commandline.*;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.io.StingSAMFileWriter;
import org.broadinstitute.sting.gatk.refdata.ReadMetaDataTracker;
import org.broadinstitute.sting.gatk.walkers.ReadWalker;
import org.broadinstitute.sting.utils.sam.AlignmentUtils;
/**
* Left aligns indels in reads.
*/
public class LeftAlignIndels extends ReadWalker<Integer, Integer> {
@Output(required=false, doc="Output bam")
protected StingSAMFileWriter writer = null;
@Argument(fullName="maxReadsInRam", shortName="maxInRam", doc="max reads allowed to be kept in memory at a time by the SAMFileWriter. "+
"If too low, the tool may run out of system file descriptors needed to perform sorting; if too high, the tool may run out of memory.", required=false)
protected int MAX_RECORDS_IN_RAM = 500000;
public void initialize() {
// set up the output writer
if ( writer != null )
writer.setMaxRecordsInRam(MAX_RECORDS_IN_RAM);
}
private void emit(final SAMRecord read) {
if ( writer != null )
writer.addAlignment(read);
}
public Integer map(ReferenceContext ref, SAMRecord read, ReadMetaDataTracker metaDataTracker) {
// we can not deal with screwy records
if ( read.getCigar().numCigarElements() == 0 ) {
emit(read);
return 0;
}
// move existing indels (for 1 indel reads only) to leftmost position within identical sequence
int numBlocks = AlignmentUtils.getNumAlignmentBlocks(read);
if ( numBlocks == 2 ) {
Cigar newCigar = AlignmentUtils.leftAlignIndel(IndelRealigner.unclipCigar(read.getCigar()), ref.getBases(), read.getReadBases(), 0, 0);
newCigar = IndelRealigner.reclipCigar(newCigar, read);
read.setCigar(newCigar);
}
emit(read);
return 1;
}
public Integer reduceInit() {
return 0;
}
public Integer reduce(Integer value, Integer sum) {
return sum + value;
}
public void onTraversalDone(Integer result) {}
}