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

additional output tweaking 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1028 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
kcibul 2009-06-17 15:37:38 +00:00
parent 7d281296a7
commit 673205ed5f
1 changed files with 37 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

66
java/src/org/broadinstitute/sting/playground/gatk/walkers/SomaticMutationWalker.java
View File

@ -75,12 +75,13 @@ public class SomaticMutationWalker extends LocusWalker<Integer, Integer> {
public int MIN_MUTANT_SUM = 100; public int MIN_MUTANT_SUM = 100;
@Argument(fullName = "mode", required = false, doc="Mode of operation (detect, full)")
public String mode = "full";
// @Argument(fullName = "output_failures", required = false, doc="produce output for failed sites") // @Argument(fullName = "output_failures", required = false, doc="produce output for failed sites")
public boolean OUTPUT_FAILURES = true; public boolean OUTPUT_FAILURES = true;
public boolean bedOutput = true;
public void initialize() { public void initialize() {
} }
@ -283,9 +284,9 @@ public class SomaticMutationWalker extends LocusWalker<Integer, Integer> {
// than one alternate allele (hints at being an alignment artifact) // than one alternate allele (hints at being an alignment artifact)
ReferenceSequence refSeq; ReferenceSequence refSeq;
// TODO: don't hardcode. Make this the max read length in the pile // TODO: don't hardcode. Make this the max read length in the pile
long refStart = context.getLocation().getStart() - 100; long refStart = context.getLocation().getStart() - 150;
//tumorReadPile.offsets.get(0); //tumorReadPile.offsets.get(0);
long refStop = context.getLocation().getStart() + 100; long refStop = context.getLocation().getStart() + 150;
try { try {
IndexedFastaSequenceFile seqFile = new IndexedFastaSequenceFile(getToolkit().getArguments().referenceFile); IndexedFastaSequenceFile seqFile = new IndexedFastaSequenceFile(getToolkit().getArguments().referenceFile);
refSeq = seqFile.getSubsequenceAt(context.getContig(),refStart, refStop); refSeq = seqFile.getSubsequenceAt(context.getContig(),refStart, refStop);
@ -304,13 +305,21 @@ public class SomaticMutationWalker extends LocusWalker<Integer, Integer> {
// at least 100 or the LOD score for mutant:ref + mutant:mutant vs ref:ref must // at least 100 or the LOD score for mutant:ref + mutant:mutant vs ref:ref must
// be at least 6.3; // be at least 6.3;
double tumorLod = t2.getAltVsRef(altAllele); double tumorLod = t2.getAltVsRef(altAllele);
if (t2.qualitySums.get(altAllele) < MIN_MUTANT_SUM && tumorLod < TUMOR_LOD_THRESHOLD) { if (mode.equals("full") && t2.qualitySums.get(altAllele) < MIN_MUTANT_SUM && tumorLod < TUMOR_LOD_THRESHOLD) {
if (OUTPUT_FAILURES) { if (OUTPUT_FAILURES) {
System.out.println("FAILED " + context.getContig() + ":" + context.getPosition() + " due to MAX MM QSCORE TEST." +
String msg = "FAILED due to MAX MM QSCORE TEST." +
" LOD was " + tumorReadPile.getAltVsRef(altAllele) + " LOD was " + tumorReadPile.getAltVsRef(altAllele) +
" LOD is now " + t2.getAltVsRef(altAllele) + " LOD is now " + t2.getAltVsRef(altAllele) +
" QSUM was " + tumorReadPile.qualitySums.get(altAllele) + " QSUM was " + tumorReadPile.qualitySums.get(altAllele) +
" QSUM is now " + t2.qualitySums.get(altAllele)); " QSUM is now " + t2.qualitySums.get(altAllele);
System.out.println(
context.getContig() + "\t" +
context.getPosition() + "\t" +
context.getPosition() + "\t"
+ msg.replace(' ','_')
);
} }
continue; continue;
} }
@ -320,9 +329,17 @@ public class SomaticMutationWalker extends LocusWalker<Integer, Integer> {
boolean shouldDisalign = boolean shouldDisalign =
disaligner(context.getPosition(), tumorReadPile, StringUtil.bytesToString(refSeq.getBases()), refStart); disaligner(context.getPosition(), tumorReadPile, StringUtil.bytesToString(refSeq.getBases()), refStart);
if (shouldDisalign) { if (mode.equals("full") && shouldDisalign) {
if (OUTPUT_FAILURES) { if (OUTPUT_FAILURES) {
System.out.println("FAILED " + context.getContig() + ":" + context.getPosition() + " due to DISALIGNMENT TEST."); String msg = "FAILED due to DISALIGNMENT TEST.";
System.out.println(
context.getContig() + "\t" +
context.getPosition() + "\t" +
context.getPosition() + "\t"
+ msg.replace(' ','_')
);
} }
continue; continue;
} }
@ -341,33 +358,24 @@ public class SomaticMutationWalker extends LocusWalker<Integer, Integer> {
// if we're still here... we've got a somatic mutation! Output the results // if we're still here... we've got a somatic mutation! Output the results
// and stop looking for mutants! // and stop looking for mutants!
if (bedOutput) { String msg =
out.println( (failedMidpointCheck?"__FAILED-MPCHECK":"") +
context.getContig() + "\t" +
context.getPosition() + "\t" +
context.getPosition() + "\t" +
"TScore:" + tumorLod + "TScore:" + tumorLod +
"__TRefSum: " + tumorReadPile.qualitySums.get(ref) + "__TRefSum: " + tumorReadPile.qualitySums.get(ref) +
"__TAltSum: " + tumorReadPile.qualitySums.get(altAllele) + "__TAltSum: " + tumorReadPile.qualitySums.get(altAllele) +
"__NScore:" + normalLod + "__NScore:" + normalLod +
"__NRefSum: " + normalReadPile.qualitySums.get(ref) + "__NRefSum: " + normalReadPile.qualitySums.get(ref) +
"__NAltSum: " + normalReadPile.qualitySums.get(altAllele) + "__NAltSum: " + normalReadPile.qualitySums.get(altAllele) +
"__MIDP: " + midp.get(altAllele) + "__MIDP: " + midp.get(altAllele);
(failedMidpointCheck?"__FAILED-MPCHECK":"")
); System.out.println(
context.getContig() + "\t" +
context.getPosition() + "\t" +
context.getPosition() + "\t"
+ msg.replace(' ','_')
);
} else {
out.println(context.getLocation() + " " + upRef + " " + altAllele +
" TScore:" + tumorLod +
" TRefSum: " + tumorReadPile.qualitySums.get(ref) +
" TAltSum: " + tumorReadPile.qualitySums.get(altAllele) +
" NScore:" + normalLod +
" NRefSum: " + normalReadPile.qualitySums.get(ref) +
" NAltSum: " + normalReadPile.qualitySums.get(altAllele) + " " +
tumorReadPile.getLocusBases().toString() + " " +
normalReadPile.getLocusBases().toString()
);
}
return 1; return 1;