Quick fix for Danny Lieber: flesh out the additional functionality required

to align to a reference other than what's specified in the header.


git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@5133 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
hanna 2011-01-31 05:28:37 +00:00
parent b5d1aab8dc
commit 5e7a5cf924
3 changed files with 19 additions and 30 deletions

View File

@ -195,6 +195,12 @@ public class Alignment {
if(newSAMHeader != null) if(newSAMHeader != null)
read.setHeader(newSAMHeader); read.setHeader(newSAMHeader);
// If we're realigning a previously aligned record, strip out the placement of the alignment.
read.setReferenceName(SAMRecord.NO_ALIGNMENT_REFERENCE_NAME);
read.setAlignmentStart(SAMRecord.NO_ALIGNMENT_START);
read.setMateReferenceName(SAMRecord.NO_ALIGNMENT_REFERENCE_NAME);
read.setMateAlignmentStart(SAMRecord.NO_ALIGNMENT_START);
if(alignment != null) { if(alignment != null) {
read.setReadUmappedFlag(false); read.setReadUmappedFlag(false);
read.setReferenceIndex(alignment.getContigIndex()); read.setReferenceIndex(alignment.getContigIndex());

View File

@ -25,6 +25,7 @@
package org.broadinstitute.sting.alignment; package org.broadinstitute.sting.alignment;
import org.broadinstitute.sting.gatk.io.StingSAMFileWriter;
import org.broadinstitute.sting.gatk.refdata.ReadMetaDataTracker; import org.broadinstitute.sting.gatk.refdata.ReadMetaDataTracker;
import org.broadinstitute.sting.commandline.Argument; import org.broadinstitute.sting.commandline.Argument;
import org.broadinstitute.sting.commandline.Output; import org.broadinstitute.sting.commandline.Output;
@ -49,28 +50,19 @@ import java.io.PrintStream;
*/ */
@WalkerName("Align") @WalkerName("Align")
public class AlignmentWalker extends ReadWalker<Integer,Integer> { public class AlignmentWalker extends ReadWalker<Integer,Integer> {
@Argument(fullName="BWTPrefix",shortName="BWT",doc="Index files generated by bwa index -d bwtsw",required=false) @Argument(fullName="target_reference",shortName="target_ref",doc="The reference to which reads in the source file should be aligned. Alongside this reference should sit index files " +
private String prefix = null; "generated by bwa index -d bwtsw. If unspecified, will default " +
"to the reference specified via the -R argument.",required=false)
@Argument(fullName = "outputBam", shortName = "ob", doc = "Write output to this BAM filename instead of STDOUT", required = false) private File targetReferenceFile = null;
private String outputBamFile = null;
@Output @Output
private PrintStream out = null; private StingSAMFileWriter out = null;
@Argument(fullName = "bam_compression", shortName = "compress", doc = "Compression level to use for writing BAM files", required = false)
private Integer bamCompression = 5;
/** /**
* The actual aligner. * The actual aligner.
*/ */
private BWACAligner aligner = null; private BWACAligner aligner = null;
/**
* Target for reads to output.
*/
private SAMFileWriter outputBam = null;
/** /**
* New header to use, if desired. * New header to use, if desired.
*/ */
@ -81,23 +73,20 @@ public class AlignmentWalker extends ReadWalker<Integer,Integer> {
*/ */
@Override @Override
public void initialize() { public void initialize() {
if(prefix == null) if(targetReferenceFile == null)
prefix = getToolkit().getArguments().referenceFile.getAbsolutePath(); targetReferenceFile = getToolkit().getArguments().referenceFile;
BWTFiles bwtFiles = new BWTFiles(prefix); BWTFiles bwtFiles = new BWTFiles(targetReferenceFile.getAbsolutePath());
BWAConfiguration configuration = new BWAConfiguration(); BWAConfiguration configuration = new BWAConfiguration();
aligner = new BWACAligner(bwtFiles,configuration); aligner = new BWACAligner(bwtFiles,configuration);
// Take the header of the SAM file, tweak it by adding in the reference dictionary and specifying that the target file is unsorted. // Take the header of the SAM file, tweak it by adding in the reference dictionary and specifying that the target file is unsorted.
header = getToolkit().getSAMFileHeader().clone(); header = getToolkit().getSAMFileHeader().clone();
SAMSequenceDictionary referenceDictionary = SAMSequenceDictionary referenceDictionary =
ReferenceSequenceFileFactory.getReferenceSequenceFile(getToolkit().getArguments().referenceFile).getSequenceDictionary(); ReferenceSequenceFileFactory.getReferenceSequenceFile(targetReferenceFile).getSequenceDictionary();
header.setSequenceDictionary(referenceDictionary); header.setSequenceDictionary(referenceDictionary);
header.setSortOrder(SAMFileHeader.SortOrder.unsorted); header.setSortOrder(SAMFileHeader.SortOrder.unsorted);
if (outputBamFile != null) { out.writeHeader(header);
SAMFileWriterFactory factory = new SAMFileWriterFactory();
outputBam = factory.makeBAMWriter(header, false, new File(outputBamFile), bamCompression);
}
} }
/** /**
@ -109,11 +98,7 @@ public class AlignmentWalker extends ReadWalker<Integer,Integer> {
@Override @Override
public Integer map(ReferenceContext ref, SAMRecord read, ReadMetaDataTracker metaDataTracker) { public Integer map(ReferenceContext ref, SAMRecord read, ReadMetaDataTracker metaDataTracker) {
SAMRecord alignedRead = aligner.align(read,header); SAMRecord alignedRead = aligner.align(read,header);
if (outputBam != null) { out.addAlignment(alignedRead);
outputBam.addAlignment(alignedRead);
} else {
out.println(alignedRead.format());
}
return 1; return 1;
} }
@ -142,8 +127,6 @@ public class AlignmentWalker extends ReadWalker<Integer,Integer> {
@Override @Override
public void onTraversalDone(Integer result) { public void onTraversalDone(Integer result) {
aligner.close(); aligner.close();
if (outputBam != null)
outputBam.close();
super.onTraversalDone(result); super.onTraversalDone(result);
} }

View File

@ -19,7 +19,7 @@ public class AlignerIntegrationTest extends WalkerTest {
"-R /humgen/gsa-scr1/GATK_Data/bwa/human_b36_both.fasta" + "-R /humgen/gsa-scr1/GATK_Data/bwa/human_b36_both.fasta" +
" -T Align" + " -T Align" +
" -I " + validationDataLocation + "NA12878_Pilot1_20.trimmed.unmapped.bam" + " -I " + validationDataLocation + "NA12878_Pilot1_20.trimmed.unmapped.bam" +
" -ob %s", " -o %s",
1, // just one output file 1, // just one output file
Arrays.asList(md5)); Arrays.asList(md5));
executeTest("testBasicAlignment", spec); executeTest("testBasicAlignment", spec);