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

a few changes; checked in to allow debugging. 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@688 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
jmaguire 2009-05-13 15:50:48 +00:00
parent 5f924c46e0
commit 7084ecdeb6
2 changed files with 28 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
java/src/org/broadinstitute/sting/playground/gatk/walkers/PoolCaller.java
View File

@ -19,6 +19,9 @@ import java.util.List;
import java.util.Random;
import java.io.*;
// Draft iterative pooled caller
// j.maguire 4-27-2009
public class PoolCaller extends LocusWalker<AlleleFrequencyEstimate, String>
{
List<SingleSampleGenotyper> callers = null;
@ -28,11 +31,11 @@ public class PoolCaller extends LocusWalker<AlleleFrequencyEstimate, String>
@Argument(required=false, shortName="max_iterations", doc="Maximum number of iterations for EM") public int MAX_ITERATIONS = 10;
@Argument(fullName="lodThreshold", shortName="lod", required=false, doc="lod threshold for outputting individual genotypes") public Double lodThreshold = 2.0;
@Argument(fullName="discovery_output", shortName="discovery_output", required=true, doc="file to write SNP discovery output to") public String DISCOVERY_OUTPUT;
@Argument(fullName="individual_output_prefix", shortName="individual_output_prefix", required=true, doc="prefix to write individual SNP calls to") public String INDIVIDUAL_OUTPUT_PREFIX;
private Random random;
private SAMFileHeader header;
private PrintStream discovery_output_file;
public void initialize()
@ -69,27 +72,26 @@ public class PoolCaller extends LocusWalker<AlleleFrequencyEstimate, String>
System.out.println("SAMPLE: " + sample_name);
SingleSampleGenotyper caller = new SingleSampleGenotyper();
caller.metricsFileName = "/dev/null";
caller.callsFileName = INDIVIDUAL_OUTPUT_PREFIX + "." + sample_name + ".calls";
caller.metricsFileName = INDIVIDUAL_OUTPUT_PREFIX + "." + sample_name + ".metrics";
caller.lodThreshold = lodThreshold;
caller.fourBaseMode = false;
caller.printMetrics = false;
caller.initialize();
//caller.initialize();
callers.add(caller);
}
}
public AlleleFrequencyEstimate map(RefMetaDataTracker tracker, char ref, LocusContext context)
{
if (ref == 'N') { return null; }
// 1. seperate each context.
LocusContext[] contexts = new LocusContext[sample_names.size()];
for (int i = 0; i < sample_names.size(); i++)
{
contexts[i] = filterLocusContext(context, sample_names.get(i), 0);
//System.out.printf("DEPTH %s %d\n", sample_names.get(i), contexts[i].getReads().size());
}
//System.out.printf("DEPTH %s %d\n", "TOTAL", context.getReads().size());
// EM Loop:
AlleleFrequencyEstimate[] calls = null;

42
java/src/org/broadinstitute/sting/playground/gatk/walkers/SingleSampleGenotyper.java
View File

@ -15,6 +15,7 @@ import org.broadinstitute.sting.utils.QualityUtils;
import net.sf.samtools.SAMRecord;
import java.io.*;
import java.util.*;
// Draft single sample genotyper
@ -22,6 +23,7 @@ import java.util.*;
public class SingleSampleGenotyper extends LocusWalker<AlleleFrequencyEstimate, String>
{
@Argument(fullName="calls", shortName="calls", doc="File to write SNP calls to", required=true) public String callsFileName;
@Argument(fullName="metrics", shortName="met", doc="metrics", required=false) public String metricsFileName = "/dev/null";
@Argument(fullName="metInterval", shortName="mi", doc="metInterval", required=false) public Integer metricsInterval = 50000;
@Argument(fullName="printMetrics", shortName="printMetrics", doc="printMetrics", required=false) public Boolean printMetrics = true;
@ -34,38 +36,32 @@ public class SingleSampleGenotyper extends LocusWalker<AlleleFrequencyEstimate,
@Argument(fullName="qHomNonRef", shortName="qHomNonRef", doc="qHomNonRef", required=false) public Double qHomNonRef = 0.97;
public AlleleMetrics metrics;
public PrintStream calls_file;
public String sample_name;
public boolean filter(RefMetaDataTracker tracker, char ref, LocusContext context) { return true; }
public boolean requiresReads() { return true; }
public void initialize() { metrics = new AlleleMetrics(metricsFileName, lodThreshold); sample_name = null; }
public void initialize()
{
try
{
metrics = new AlleleMetrics(metricsFileName, lodThreshold);
sample_name = null;
calls_file = new PrintStream(callsFileName);
}
catch (Exception e)
{
e.printStackTrace();
System.exit(-1);
}
}
public AlleleFrequencyEstimate map(RefMetaDataTracker tracker, char ref, LocusContext context) {
String rodString = getRodString(tracker);
if (ref == 'N') { return null; }
/*
AlleleFrequencyEstimate freq;
if (fourBaseMode) {
if (retest) {
// Compute single quality score genotype likelihoods
freq = getOneProbAlleleFrequency(ref, context, rodString);
if (freq.qhat > 0.1) {
// If we found a putative variant, recompute four-base prob genotype likelihoods
freq = getFourProbAlleleFrequency(ref, context, rodString, tracker);
}
} else {
freq = getFourProbAlleleFrequency(ref, context, rodString, tracker);
}
} else {
// Compute single quality score genotype likelihoods
freq = getOneProbAlleleFrequency(ref, context, rodString);
}
*/
if (context.getReads().size() != 0)
{
SAMRecord read = context.getReads().get(0);
@ -367,7 +363,7 @@ public class SingleSampleGenotyper extends LocusWalker<AlleleFrequencyEstimate,
// No longer hom-ref, so output a ref line.
double lod = confident_ref_interval_LOD_sum / confident_ref_interval_length;
out.format("%s\tCALLER\tREFERENCE\t%d\t%d\t%f\t.\t.\tLENGTH %d\n",
calls_file.format("%s\tCALLER\tREFERENCE\t%d\t%d\t%f\t.\t.\tLENGTH %d\n",
confident_ref_interval_contig,
confident_ref_interval_start,
last_position_considered,
@ -399,7 +395,7 @@ public class SingleSampleGenotyper extends LocusWalker<AlleleFrequencyEstimate,
last_position_considered = current_offset;
if (alleleFreq.lodVsRef >= 5) {
out.print(alleleFreq.asGFFString());
calls_file.print(alleleFreq.asGFFString());
/*
String gtype = genotypeTypeString(alleleFreq.qstar, alleleFreq.N);