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

Minor changes to a qscript and the GQ constants on PrivatePermutations 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@4956 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
chartl 2011-01-07 18:26:21 +00:00
parent 79fcff13ff
commit 3e7802a3e0
2 changed files with 28 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
java/src/org/broadinstitute/sting/oneoffprojects/walkers/varianteval/PrivatePermutations.java
View File

@ -25,8 +25,9 @@ import java.util.Set;
@Analysis(name = "PrivatePermutations", description = "Number of additional mutations from each new sample; random permutations")
public class PrivatePermutations extends VariantEvaluator {
private final int NUM_PERMUTATIONS = 50;
private final double LOW_GQ_PCT = 0.95;
private final double LOW_GQ_THRSH = 30.0;
private final double LOW_GQ_PCT = 0.80;
private final double LOW_GQ_THRSH = 25.0;
private final boolean IGNORE_FILTER_FIELD = true;
private boolean initialized = false;
private long skipped = 0l;
@ -80,18 +81,21 @@ public class PrivatePermutations extends VariantEvaluator {
}
private boolean isGood(VariantContext vc) {
if ( vc == null || vc.isFiltered() || (vc.getHetCount() + vc.getHomVarCount() == 0) ) { // todo -- should be is variant, but need to ensure no alt alleles at ref sites
if ( vc == null || (vc.isFiltered() && ! IGNORE_FILTER_FIELD) || (vc.getHetCount() + vc.getHomVarCount() == 0) ) { // todo -- should be is variant, but need to ensure no alt alleles at ref sites
return false;
} else {
Collection<Genotype> gtypes = vc.getGenotypes().values();
int ngood = 0;
int nbad = 0;
for ( Genotype g : gtypes) {
if ( g.getPhredScaledQual() >= LOW_GQ_THRSH ) {
ngood ++;
} else {
nbad ++;
}
}
return ( (0.0+ngood)/(0.0+gtypes.size()) >= LOW_GQ_PCT );
return ( (0.0+ngood)/(0.0+ngood + nbad) >= LOW_GQ_PCT );
}
}

25
scala/qscript/chartl/expanded_targets.q
View File

@ -12,6 +12,8 @@ class expanded_targets extends QScript {
@Argument(shortName="thisTrigger",doc="The trigger track to use",required=false) var thisTrigger : File = new File("/humgen/gsa-hphome1/chartl/projects/exome/expanded/triggers/joined.omni.hiseq.vcf")
def script = {
// note : bam sorting and indexing handled outside of this script
// GATK hacked by modifying GenomeLocSortedSet not to fuss about getting multiple instances of the same interval
val intervalExpands : List[ExpandIntervals] = (new Range(0,40,1)).toList.map( u => {
new ExpandIntervals(args.projectIntervals,1+5*u,5,new File("./"+args.projectName+"_expanded_%d_%d.interval_list".format(1+5*u,6+5*u)),args.projectRef,"TSV","INTERVALS")
@ -37,14 +39,14 @@ class expanded_targets extends QScript {
rtc.out = swapExt(userDir,u,".bam",".clean.targets.interval_list")
rtc.input_file :+= u.getAbsoluteFile
rtc.intervals :+= cleanIntervals.outList
rtc.memoryLimit = Some(4)
rtc.memoryLimit = Some(6)
rtc
})
val clean : List[IndelRealigner] = realign.map( u => {
var cleaner : IndelRealigner = new IndelRealigner with GATKArgs
cleaner.targetIntervals = u.out
cleaner.input_file = u.input_file
cleaner.memoryLimit = Some(4)
cleaner.memoryLimit = Some(6)
cleaner.out = new File(userDir+"/"+swapExt(u.out,".bam",".expanded.targets.bam").getName)
cleaner.intervals :+= cleanIntervals.outList
cleaner
@ -53,7 +55,7 @@ class expanded_targets extends QScript {
addAll(realign)
addAll(clean)
val callFiles: List[File] = intervalExpands.map(u => makeCalls(u.outList,clean.map(h => h.out)))
val callFiles: List[File] = intervalExpands.map(u => makeCalls(u.outList,clean.map(h => swapExt(h.out,".bam",".sorted.bam"))))
}
@ -81,7 +83,8 @@ class expanded_targets extends QScript {
filter.filterName :+= "HighStrandBias"
filter.out = swapExt(iList,".interval_list",".filtered.vcf")
var callHiseq : UnifiedGenotyper = new UnifiedGenotyper with GATKArgs
callHiseq.input_file = List(new File("/humgen/1kg/analysis/bamsForDataProcessingPapers/NA12878.HiSeq.WGS.bwa.cleaned.recal.bam"))
callHiseq.reference_sequence = new File("/seq/references/Homo_sapiens_assembly19/v1/Homo_sapiens_assembly19.fasta")
callHiseq.input_file = List(new File("/seq/picard_aggregation/EXT1/NA12878/v3/NA12878.bam"))
callHiseq.rodBind :+= new RodBind("trigger","vcf",filter.out)
callHiseq.out = swapExt(iList,".interval_list",".hiSeq.genotypes.vcf")
callHiseq.trig_emit_conf = Some(0.0)
@ -93,12 +96,24 @@ class expanded_targets extends QScript {
eval.rodBind :+= new RodBind("evalInterval","vcf",filter.out)
eval.rodBind :+= new RodBind("compHiSeq","vcf",new File("/humgen/gsa-hpprojects/GATK/data/Comparisons/Unvalidated/NA12878/NA12878.hg19.HiSeq.WGS.cleaned.ug.snpfiltered.indelfiltered.optimized.cut.vcf"))
eval.rodBind :+= new RodBind("compHiSeq_atSites","vcf",callHiseq.out)
eval.rodBind :+= new RodBind("compOMNI","vcf",new File("/humgen/gsa-hpprojects/GATK/data/Comparisons/Validated/Omni2.5_chip/Omni_2.5_764_samples.b37.deduped.annot.vcf"))
eval.rodBind :+= new RodBind("compOMNI","vcf",new File("/humgen/gsa-hpprojects/GATK/data/Comparisons/Validated/Omni2.5_chip/764samples.deduped.b37.annot.vcf"))
eval.out = swapExt(iList,".interval_list",".eval")
eval.reportType = Option(org.broadinstitute.sting.utils.report.VE2ReportFactory.VE2TemplateType.CSV)
eval.memoryLimit = Some(4)
add(eval)
eval.out
}
class B37_to_HG19 extends CommandLineFunction {
@Input(doc="vcf") var vcf : File = _
@Output(doc="out") var outVCF : File = _
def commandLine = "python /humgen/gsa-hphome1/chartl/sting/python/vcf_b36_to_hg18.py %s %s".format(vcf.getAbsolutePath,outVCF.getAbsolutePath)
}
class HG19_to_B37 extends B37_to_HG19 {
override def commandLine = "python /humgen/gsa-hphome1/chartl/sting/python/vcf_b36_to_hg18.py -r %s %s".format(vcf.getAbsolutePath,outVCF.getAbsolutePath)
}
}