Releasable version of the Pileup walker

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@786 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
depristo 2009-05-21 22:25:37 +00:00
parent dc17a5661d
commit 65995887fc
1 changed files with 40 additions and 18 deletions

View File

@ -6,43 +6,65 @@ import org.broadinstitute.sting.gatk.refdata.rodDbSNP;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker; import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.cmdLine.Argument; import org.broadinstitute.sting.utils.cmdLine.Argument;
import org.broadinstitute.sting.utils.ReadBackedPileup; import org.broadinstitute.sting.utils.ReadBackedPileup;
import org.broadinstitute.sting.utils.Utils;
import net.sf.samtools.SAMRecord;
import java.util.List;
/** /**
* Created by IntelliJ IDEA. * samtools pileup [-f in.ref.fasta] [-t in.ref_list] [-l in.site_list] [-iscg] [-T theta] [-N nHap] [-r pairDiffRate] <in.alignment>
* User: mdepristo *
* Date: Feb 22, 2009 * Print the alignment in the pileup format. In the pileup format, each line represents a genomic position,
* Time: 3:22:14 PM * consisting of chromosome name, coordinate, reference base, read bases, read qualities and alignment mapping
* To change this template use File | Settings | File Templates. * qualities. Information on match, mismatch, indel, strand, mapping quality and start and end of a read are all
* encoded at the read base column. At this column, a dot stands for a match to the reference base on the forward strand,
* a comma for a match on the reverse strand, ÔACGTNÕ for a mismatch on the forward strand and ÔacgtnÕ for a mismatch on the
* reverse strand.
*
* A pattern Ô\+[0-9]+[ACGTNacgtn]+Õ indicates there is an insertion between this reference position and the next
* reference position. The length of the insertion is given by the integer in the pattern, followed by the inserted sequence.
* Similarly, a pattern Ô-[0-9]+[ACGTNacgtn]+Õ represents a deletion from the reference.
* Also at the read base column, a symbol Ô^Õ marks the start of a read segment which is a contiguous subsequence on the read
* separated by ÔN/S/HÕ CIGAR operations. The ASCII of the character following Ô^Õ minus 33 gives the mapping quality.
* A symbol Ô$Õ marks the end of a read segment.
*/ */
public class PileupWalker extends LocusWalker<Integer, Integer> implements TreeReducible<Integer> { public class PileupWalker extends LocusWalker<Integer, Integer> implements TreeReducible<Integer> {
@Argument(fullName="verbose",doc="verbose",required=false) @Argument(fullName="alwaysShowSecondBase",doc="If true, prints dummy bases for the second bases in the BAM file where they are missing",required=false)
public boolean VERBOSE = false; public boolean alwaysShowSecondBase = false;
@Argument(fullName="showSecondBaseQuals",doc="If true, prints out second base qualities in the pileup",required=false)
public boolean showSecondBaseQuals = false;
@Argument(fullName="extended",shortName="ext",doc="extended",required=false) @Argument(fullName="extended",shortName="ext",doc="extended",required=false)
public boolean EXTENDED = false; public boolean EXTENDED = false;
public boolean FLAG_UNCOVERED_BASES = true; // todo: how do I make this a command line argument? public boolean FLAG_UNCOVERED_BASES = true; // todo: how do I make this a command line argument?
public void initialize() { public void initialize() {
} }
// Do we actually want to operate on the context?
public boolean filter(RefMetaDataTracker tracker, char ref, LocusContext context) {
return true; // We are keeping all the reads
}
public Integer map(RefMetaDataTracker tracker, char ref, LocusContext context) { public Integer map(RefMetaDataTracker tracker, char ref, LocusContext context) {
ReadBackedPileup pileup = new ReadBackedPileup(ref, context); ReadBackedPileup pileup = new ReadBackedPileup(ref, context);
String bases = pileup.getBases(); String bases = pileup.getBases();
if ( bases.equals("") && FLAG_UNCOVERED_BASES ) { if ( bases.equals("") && FLAG_UNCOVERED_BASES ) {
bases = "*** UNCOVERED SITE ***"; bases = "***UNCOVERED_SITE***";
} }
String extras = ""; StringBuilder extras = new StringBuilder();
if ( VERBOSE ) {
extras += " BQ=" + pileup.getQualsAsInts(); String secondBasePileup = pileup.getSecondaryBasePileup();
extras += " MQ=" + pileup.getMappingQualsAsInts(); if ( secondBasePileup == null && alwaysShowSecondBase ) {
secondBasePileup = Utils.dupString('N', bases.length());
}
if ( secondBasePileup != null ) extras.append(" ").append(secondBasePileup);
if ( showSecondBaseQuals ) {
String secondQualPileup = pileup.getSecondaryQualPileup();
if ( secondQualPileup == null )
secondQualPileup = Utils.dupString((char)(33), bases.length());
extras.append(" ").append(secondQualPileup);
} }
String rodString = ""; String rodString = "";