Ensure frequency of tumor to normal passes a reasonable threshold.

    Avoids calling low frequency tumors also present at low frequency in normals,
    which indicates a contamination or persistent error.
321060
def _check_freqs(parts, indexes):
    
    thresh_ratio = 2.7
    try:  # FreeBayes
        ao_index = parts[8].split(":").index("AO")
        ro_index = parts[8].split(":").index("RO")
    except ValueError:
        ao_index, ro_index = None, None
    try:  # VarDict
        af_index = parts[8].split(":").index("AF")
    except ValueError:
        af_index = None
    if af_index is None and ao_index is None:
        # okay to skip if a gVCF record
        if parts[4].find("<*>") == -1:
            raise NotImplementedError("Unexpected format annotations: %s" % parts[8])
    def _calc_freq(item):
        try:
            if ao_index is not None and ro_index is not None:
                ao = sum([int(x) for x in item.split(":")[ao_index].split(",")])
                ro = int(item.split(":")[ro_index])
                freq = ao / float(ao + ro)
            elif af_index is not None:
                freq = float(item.split(":")[af_index])
            else:
                freq = 0.0
        except (IndexError, ValueError, ZeroDivisionError):
            freq = 0.0
        return freq
    tumor_freq, normal_freq = _calc_freq(parts[indexes["tumor"]]), _calc_freq(parts[indexes["normal"]])
    return normal_freq <= 0.001 or normal_freq <= tumor_freq / thresh_ratio
1124760
def check_frequencies(
    feed: "Feed", *, as_df: bool = False, include_warnings: bool = False
) -> List:
    
    table = "frequencies"
    problems = []

    # Preliminary checks
    if feed.frequencies is None:
        return problems

    f = feed.frequencies.copy()
    problems = check_for_required_columns(problems, table, f)
    if problems:
        return format_problems(problems, as_df=as_df)

    if include_warnings:
        problems = check_for_invalid_columns(problems, table, f)

    # Check trip_id
    problems = check_column_linked_id(
        problems, table, f, "trip_id", feed.trips
    )

    # Check start_time and end_time
    time_cols = ["start_time", "end_time"]
    for col in time_cols:
        problems = check_column(problems, table, f, col, valid_time)

    for col in time_cols:
        f[col] = f[col].map(hp.timestr_to_seconds)

    # Start_time should be earlier than end_time
    cond = f["start_time"] >= f["end_time"]
    problems = check_table(
        problems, table, f, cond, "start_time not earlier than end_time"
    )

    # Headway periods should not overlap
    f = f.sort_values(["trip_id", "start_time"])
    for __, group in f.groupby("trip_id"):
        a = group["start_time"].values
        b = group["end_time"].values
        indices = np.flatnonzero(a[1:] < b[:-1]).tolist()
        if indices:
            problems.append(
                [
                    "error",
                    "Headway periods for the same trip overlap",
                    table,
                    indices,
                ]
            )

    # Check headway_secs
    v = lambda x: x >= 0
    problems = check_column(problems, table, f, "headway_secs", v)

    # Check exact_times
    v = lambda x: x in range(2)
    problems = check_column(
        problems, table, f, "exact_times", v, column_required=False
    )

    return format_problems(problems, as_df=as_df)
16001
def __normalize_progress(self):
        

        progress = self['funding_progress']
        if progress % 10 != 0:
            progress = round(float(progress) / 10)
            progress = int(progress) * 10

            self['funding_progress'] = progress
25564
def threshold(args):
    """"""
    if args.fpr < 0 or args.fpr > 1:
        print("Please specify a FPR between 0 and 1")
        sys.exit(1)

    motifs = read_motifs(args.pwmfile)
    
    s = Scanner()
    s.set_motifs(args.pwmfile)
    s.set_threshold(args.fpr, filename=args.inputfile)

    print("Motif\tScore\tCutoff")
    for motif in motifs:
        min_score = motif.pwm_min_score()
        max_score = motif.pwm_max_score()
        opt_score = s.threshold[motif.id]
        if opt_score is None:
            opt_score = motif.pwm_max_score()
        threshold = (opt_score - min_score) / (max_score - min_score)
        print("{0}\t{1}\t{2}".format(
                motif.id, opt_score, threshold))
125788
def _freqs_by_chromosome(in_file, params, somatic_info):
    
    freqs = []
    coords = []
    cur_chrom = None
    with pysam.VariantFile(in_file) as bcf_in:
        for rec in bcf_in:
            if _is_biallelic_snp(rec) and _passes_plus_germline(rec) and chromhacks.is_autosomal(rec.chrom):
                if cur_chrom is None or rec.chrom != cur_chrom:
                    if cur_chrom and len(freqs) > 0:
                        yield cur_chrom, freqs, coords
                    cur_chrom = rec.chrom
                    freqs = []
                    coords = []
                stats = _tumor_normal_stats(rec, somatic_info)
                if tz.get_in(["tumor", "depth"], stats, 0) > params["min_depth"]:
                    # not a ref only call
                    if len(rec.samples) == 0 or sum(rec.samples[somatic_info.tumor_name].allele_indices) > 0:
                        freqs.append(tz.get_in(["tumor", "freq"], stats))
                        coords.append(rec.start)
        if cur_chrom and len(freqs) > 0:
            yield cur_chrom, freqs, coords
1137402
def ensure_standard_filter_columns(self, low_frequency_cutoff=None):
        

        # Make sure we have a template duration field
        if not hasattr(self.table, 'template_duration'):
            self.table = self.table.add_fields(np.zeros(len(self.table),
                                     dtype=np.float32), 'template_duration')

        # Make sure we have a f_lower field
        if low_frequency_cutoff is not None:
            if not hasattr(self.table, 'f_lower'):
                vec = np.zeros(len(self.table), dtype=np.float32)
                self.table = self.table.add_fields(vec, 'f_lower')
            self.table['f_lower'][:] = low_frequency_cutoff

        self.min_f_lower = min(self.table['f_lower'])
        if self.f_lower is None and self.min_f_lower == 0.:
            raise ValueError('Invalid low-frequency cutoff settings')
799785
def filter_against_normal(self, normal_mutations, maf_min=0.2,
                              maf_count_threshold=20, count_min=1):
        """"""
        assert(normal_mutations.chrom == self.chrom)
        assert(normal_mutations.pos == self.pos)
        assert(normal_mutations.ref == self.ref)

        def passes_normal_criteria(mut):
            return (mut.count >= maf_count_threshold and mut.maf > maf_min) or \
                (mut.count < maf_count_threshold and mut.count > count_min)

        nms = normal_mutations
        muts = MutationsAtSinglePosition(self.chrom, self.pos, self.cov, self.ref)

        for snv in self.snvs:
            if not (snv in nms.snvs and passes_normal_criteria(nms.snvs[snv])):
                muts.add_snv(self.snvs[snv])

        for dlt in self.deletions:
            if not (dlt in nms.deletions and passes_normal_criteria(nms.deletions[dlt])):
                muts.add_deletion(self.deletions[dlt])

        for ins in self.insertions:
            if not (ins in nms.insertions and passes_normal_criteria(nms.insertions[ins])):
                muts.add_insertion(self.insertions[ins])

        return muts
306099
def update(self, forecasts, observations):
        
        for t, threshold in enumerate(self.thresholds[:-1]):
            self.frequencies.loc[t, "Positive_Freq"] += np.count_nonzero((threshold <= forecasts) &
                                                                         (forecasts < self.thresholds[t+1]) &
                                                                         (observations >= self.obs_threshold))
            self.frequencies.loc[t, "Total_Freq"] += np.count_nonzero((threshold <= forecasts) &
                                                                      (forecasts < self.thresholds[t+1]))
962588
private static boolean checkCriticalValue(double score, boolean is_twoTailed, double aLevel) {
        double probability = ContinuousDistributions.gaussCdf(score);
        
        boolean rejectH0=false;
        
        double a=aLevel;
        if(is_twoTailed) { //if to tailed test then split the statistical significance in half
            a=aLevel/2.0;
        }
        if(probability<=a || probability>=(1.0-a)) {
            rejectH0=true;
        }
        
        return rejectH0;
    }
1306518
def validate_candidate(candidate, freq_set, k):
    
    for subcand in combinations(candidate, k-1):
        if subcand not in freq_set:
            return False
    return True
947460
def estimate_percent_contamination(contamination_report_file):
    
    contam_levels = list()
    with open(contamination_report_file) as csvfile:
        reader = csv.DictReader(csvfile)
        for row in reader:
            lowest_count = 99999
            base_counts = row['Bases'].split(';')
            for count in base_counts:
                num_bases = int(count.split(':')[1])
                if num_bases < lowest_count:
                    lowest_count = num_bases
            total_coverage = int(row['Coverage'])
            contam_levels.append(lowest_count*100/total_coverage)
    return '%.2f' % (np.mean(contam_levels)), '%.2f' % np.std(contam_levels)
518077
def _do_prioritize(items):
    
    if not any("tumoronly-prioritization" in dd.get_tools_off(d) for d in items):
        if vcfutils.get_paired_phenotype(items[0]):
            has_tumor = False
            has_normal = False
            for sub_data in items:
                if vcfutils.get_paired_phenotype(sub_data) == "tumor":
                    has_tumor = True
                elif vcfutils.get_paired_phenotype(sub_data) == "normal":
                    has_normal = True
            return has_tumor and not has_normal
96757
function getFrequency(value) {
    if (context.isFake) {
        return 0;
    }
    // get frequency by passing number from 0 to 1
    // Clamp the frequency between the minimum value (40 Hz) and half of the
    // sampling rate.
    const minValue = 40;
    const maxValue = context.sampleRate / 2;
    // Logarithm (base 2) to compute how many octaves fall in the range.
    const numberOfOctaves = Math.log(maxValue / minValue) / Math.LN2;
    // Compute a multiplier from 0 to 1 based on an exponential scale.
    const multiplier = Math.pow(2, numberOfOctaves * (value - 1.0));
    // Get back to the frequency value between min and max.
    return maxValue * multiplier;
}
598589
private void removeSporadic() {
		//System.out.println("REMOVE SPORADIC CALLED");
		// 1. For each grid g in grid_list
		//    a. If g is sporadic
		//       i. If currTime - tg > gap, delete g from grid_list
		//       ii. Else if (S1 && S2), mark as sporadic
		//       iii. Else, mark as normal
		//    b. Else
		//       i. If (S1 && S2), mark as sporadic
		
		// For each grid g in grid_list
		Iterator<Map.Entry<DensityGrid, CharacteristicVector>> glIter = this.grid_list.entrySet().iterator();
		HashMap<DensityGrid, CharacteristicVector> newGL = new HashMap<DensityGrid, CharacteristicVector>();
		ArrayList<DensityGrid> remGL = new ArrayList<DensityGrid>();
				
		while(glIter.hasNext())
		{
			Map.Entry<DensityGrid, CharacteristicVector> grid = glIter.next();
			DensityGrid dg = grid.getKey();
			CharacteristicVector cv = grid.getValue();
			
			// If g is sporadic
			if (cv.isSporadic())
			{
				// If currTime - tg > gap, delete g from grid_list
				if ((this.getCurrTime() - cv.getUpdateTime()) >= gap)
				{
					int dgClass = cv.getLabel();
					
					if (dgClass != -1)
						this.cluster_list.get(dgClass).removeGrid(dg);
					
					remGL.add(dg);
					//System.out.println("Removed "+dg.toString()+" from cluster "+dgClass);
				}
				// Else if (S1 && S2), mark as sporadic - Else mark as normal
				else
				{
					cv.setSporadic(checkIfSporadic(cv));
					//System.out.println("within gap" + dg.toString() + " sporadicity assessed "+cv.isSporadic());
					newGL.put(dg, cv);
				}
				
			}
			// Else if (S1 && S2), mark as sporadic
			else
			{
				cv.setSporadic(checkIfSporadic(cv));
				//System.out.println(dg.toString() + " sporadicity assessed "+cv.isSporadic());
				newGL.put(dg, cv);
			}
		}
		
		this.grid_list.putAll(newGL);
		
		//System.out.println(" - Removed "+remGL.size()+" grids from grid_list.");
		Iterator<DensityGrid> remIter = remGL.iterator();
		
		while(remIter.hasNext())
		{
			DensityGrid sporadicDG = remIter.next();
			//System.out.println("Removing sporadic grid "+sporadicDG.toString()+" at time "+this.getCurrTime()+".");
			this.deleted_grids.put(sporadicDG, new Integer(this.getCurrTime()));
			this.grid_list.remove(sporadicDG);
		}
		
	}
157088
def filter_repeated(self):
        '''
        
        '''
        defaults = {
            'occurrences': 1,
            'interval': 30,
            'refresh': 1800
        }

        # Override defaults with anything defined in the settings
        if isinstance(self.settings['sensu_plugin'], dict):
            defaults.update(self.settings['sensu_plugin'])

        occurrences = int(self.event['check'].get(
            'occurrences', defaults['occurrences']))
        interval = int(self.event['check'].get(
            'interval', defaults['interval']))
        refresh = int(self.event['check'].get(
            'refresh', defaults['refresh']))

        if self.event['occurrences'] < occurrences:
            self.bail('not enough occurrences')

        if (self.event['occurrences'] > occurrences and
                self.event['action'] == 'create'):
            return

        number = int(refresh / interval)
        if (number == 0 or
                (self.event['occurrences'] - occurrences) % number == 0):
            return

        self.bail('only handling every ' + str(number) + ' occurrences')
428434
def run(align_bams, items, ref_file, assoc_files, region=None, out_file=None):
    
    paired = vcfutils.get_paired_bams(align_bams, items)
    assert paired and not paired.normal_bam, ("smCounter2 supports tumor-only variant calling: %s" %
                                              (",".join([dd.get_sample_name(d) for d in items])))
    vrs = bedutils.population_variant_regions(items)
    target = shared.subset_variant_regions(vrs, region,
                                            out_file, items=items, do_merge=True)
    out_file = out_file.replace(".vcf.gz", ".vcf")
    out_prefix = utils.splitext_plus(os.path.basename(out_file))[0]
    if not utils.file_exists(out_file) and not utils.file_exists(out_file + ".gz"):
        with file_transaction(paired.tumor_data, out_file) as tx_out_file:
            cmd = ["smCounter2", "--runPath", os.path.dirname(tx_out_file),
                   "--outPrefix", out_prefix,
                   "--bedTarget", target, "--refGenome", ref_file,
                   "--bamFile", paired.tumor_bam, "--bamType", "consensus",
                   "--nCPU", dd.get_num_cores(paired.tumor_data)]
            do.run(cmd, "smcounter2 variant calling")
            for fname in glob.glob(os.path.join(os.path.dirname(tx_out_file), "*.smCounter*")):
                shutil.move(fname, os.path.join(os.path.dirname(out_file), os.path.basename(fname)))
            utils.symlink_plus(os.path.join(os.path.dirname(out_file),
                                            "%s.smCounter.cut.vcf" % out_prefix),
                               out_file)
    return vcfutils.bgzip_and_index(out_file, paired.tumor_data["config"], remove_orig=False,
                                    prep_cmd="sed 's#FORMAT\t%s#FORMAT\t%s#' | %s" %
                                    (out_prefix, dd.get_sample_name(paired.tumor_data),
                                     vcfutils.add_contig_to_header_cl(dd.get_ref_file(paired.tumor_data), out_file)))
904069
function checkThreshold(item, ctx) {
		if (item.count >= item.opts.minRequestCount) {
			const rate = item.failures / item.count;
			if (rate >= item.opts.threshold)
				trip(item, ctx);
		}
	}