Ensure likelihoods for tumor and normal pass thresholds.

    Skipped if no FreeBayes GL annotations available.
435075
def _check_lods(parts, tumor_thresh, normal_thresh, indexes):
    
    try:
        gl_index = parts[8].split(":").index("GL")
    except ValueError:
        return True
    try:
        tumor_gls = [float(x) for x in parts[indexes["tumor"]].strip().split(":")[gl_index].split(",") if x != "."]
        if tumor_gls:
            tumor_lod = max(tumor_gls[i] - tumor_gls[0] for i in range(1, len(tumor_gls)))
        else:
            tumor_lod = -1.0
    # No GL information, no tumor call (so fail it)
    except IndexError:
        tumor_lod = -1.0
    try:
        normal_gls = [float(x) for x in parts[indexes["normal"]].strip().split(":")[gl_index].split(",") if x != "."]
        if normal_gls:
            normal_lod = min(normal_gls[0] - normal_gls[i] for i in range(1, len(normal_gls)))
        else:
            normal_lod = normal_thresh
    # No GL inofmration, no normal call (so pass it)
    except IndexError:
        normal_lod = normal_thresh
    return normal_lod >= normal_thresh and tumor_lod >= tumor_thresh
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
339880
def calculate(bam_file, data, sv_bed):
    
    params = {"min": dd.get_coverage_depth_min(data)}
    variant_regions = dd.get_variant_regions_merged(data)
    if not variant_regions:
        variant_regions = _create_genome_regions(data)
    # Back compatible with previous pre-mosdepth callable files
    callable_file = os.path.join(utils.safe_makedir(os.path.join(dd.get_work_dir(data), "align",
                                                                 dd.get_sample_name(data))),
                                 "%s-coverage.callable.bed" % (dd.get_sample_name(data)))
    if not utils.file_uptodate(callable_file, bam_file):
        vr_quantize = ("0:1:%s:" % (params["min"]), ["NO_COVERAGE", "LOW_COVERAGE", "CALLABLE"])
        to_calculate = [("variant_regions", variant_regions,
                         vr_quantize, None, "coverage_perbase" in dd.get_tools_on(data)),
                        ("sv_regions", bedutils.clean_file(sv_bed, data, prefix="svregions-"),
                         None, None, False),
                        ("coverage", bedutils.clean_file(dd.get_coverage(data), data, prefix="cov-"),
                         None, DEPTH_THRESHOLDS, False)]
        depth_files = {}
        for target_name, region_bed, quantize, thresholds, per_base in to_calculate:
            if region_bed:
                cur_depth = {}
                depth_info = run_mosdepth(data, target_name, region_bed, quantize=quantize, thresholds=thresholds,
                                          per_base=per_base)
                for attr in ("dist", "regions", "thresholds", "per_base"):
                    val = getattr(depth_info, attr, None)
                    if val:
                        cur_depth[attr] = val
                depth_files[target_name] = cur_depth
                if target_name == "variant_regions":
                    callable_file = depth_info.quantize
    else:
        depth_files = {}
    final_callable = _subset_to_variant_regions(callable_file, variant_regions, data)
    return final_callable, depth_files
1125277
function assertThresholds(thresholds) {
  if (!elv(thresholds)) return null;

  if (typeof thresholds !== 'object') {
    throw new TypeError(msg.argThresholds);
  }

  const avg = elv.coalesce(thresholds.avg, false);
  const max = elv.coalesce(thresholds.max, false);
  const min = elv.coalesce(thresholds.min, false);
  const { target } = thresholds;

  if (avg !== false && typeof avg !== 'number') {
    throw new TypeError(msg.argAvg);
  }

  if (max !== false && typeof max !== 'number') {
    throw new TypeError(msg.argMax);
  }

  if (min !== false && typeof min !== 'number') {
    throw new TypeError(msg.argMin);
  }

  if (!Number.isInteger(target) || target < 0) {
    throw new TypeError(msg.argTarget);
  }

  return {
    avg,
    max,
    min,
    target,
  };
}
421769
def _SID_call_prep(align_bams, items, ref_file, assoc_files, region=None, out_file=None):
    
    base_config = items[0]["config"]
    for x in align_bams:
        bam.index(x, base_config)

    params = ["-R", ref_file, "-T", "SomaticIndelDetector", "-U", "ALLOW_N_CIGAR_READS"]
    # Limit per base read start count to between 200-10000, i.e. from any base
    # can no more 10000 new reads begin.
    # Further, limit maxNumberOfReads accordingly, otherwise SID discards
    # windows for high coverage panels.
    paired = vcfutils.get_paired_bams(align_bams, items)
    params += ["--read_filter", "NotPrimaryAlignment"]
    params += ["-I:tumor", paired.tumor_bam]
    min_af = float(get_in(paired.tumor_config, ("algorithm", "min_allele_fraction"), 10)) / 100.0
    if paired.normal_bam is not None:
        params += ["-I:normal", paired.normal_bam]
        # notice there must be at least 4 reads of coverage in normal
        params += ["--filter_expressions", "T_COV<6||N_COV<4||T_INDEL_F<%s||T_INDEL_CF<0.7" % min_af]
    else:
        params += ["--unpaired"]
        params += ["--filter_expressions", "COV<6||INDEL_F<%s||INDEL_CF<0.7" % min_af]
    if region:
        params += ["-L", bamprep.region_to_gatk(region), "--interval_set_rule",
                   "INTERSECTION"]
    return params
849701
def _variant_checkpoints(samples):
    
    checkpoints = {}
    checkpoints["vc"] = any([dd.get_variantcaller(d) or d.get("vrn_file") for d in samples])
    checkpoints["sv"] = any([dd.get_svcaller(d) for d in samples])
    checkpoints["jointvc"] = any([(dd.get_jointcaller(d) or "gvcf" in dd.get_tools_on(d))
                                  for d in samples])
    checkpoints["hla"] = any([dd.get_hlacaller(d) for d in samples])
    checkpoints["align"] = any([(dd.get_aligner(d) or dd.get_bam_clean(d)) for d in samples])
    checkpoints["align_split"] = not all([(dd.get_align_split_size(d) is False or
                                           not dd.get_aligner(d))
                                          for d in samples])
    checkpoints["archive"] = any([dd.get_archive(d) for d in samples])
    checkpoints["umi"] = any([dd.get_umi_consensus(d) for d in samples])
    checkpoints["ensemble"] = any([dd.get_ensemble(d) for d in samples])
    checkpoints["cancer"] = any(dd.get_phenotype(d) in ["tumor"] for d in samples)
    return checkpoints
44185
def validate(self):
        
        if self.stoichiometry_matrix.cols != self.propensities.rows:
            raise ValueError('There must be a column in stoichiometry matrix '
                             'for each row in propensities matrix. '
                             'S ({0.rows}x{0.cols}): {0!r} , '
                             'propensities ({1.rows}x{1.cols}): {1!r}'.format(self.stoichiometry_matrix,
                                                                              self.propensities))

        if self.stoichiometry_matrix.rows != len(self.species):
            raise ValueError('There must be a row in stoichiometry matrix for each variable. '
                             'S ({0.rows}x{0.cols}): {0!r}, variables: {1!r}'.format(self.stoichiometry_matrix,
                                                                                     self.species))

        seen_free_symbols = set()
        parameters = set(self.parameters)
        species = set(self.species)

        # Check if there are any parameters in both lists
        intersection = parameters & species
        if intersection:
            raise ValueError("Some symbols are in both parameters and species lists")

        both = parameters | species
        for row in self.propensities:
            free_symbols = row.free_symbols
            # Do not check the seen symbols twice
            free_symbols = free_symbols - seen_free_symbols
            for symbol in free_symbols:
                if symbol not in both:
                    raise ValueError('Propensity {0!r} '
                                     'contains a free symbol {1!r} '
                                     'that is not in listed in parameters or species lists '
                                     'Parameters: {2!r}; '
                                     'Species: {3!r}'.format(row, symbol,
                                                             self.parameters,
                                                             self.species))

            seen_free_symbols.update(free_symbols)
1223963
def heldout_log_likelihood(self, test_mask=None):
        
        if test_mask is None:
            # If a test mask is not supplied, use the negation of this object's mask
            if self.mask is None:
                return 0
            else:
                test_mask = ~self.mask

        xs = np.hstack((self.gaussian_states, self.inputs))
        if self.single_emission:
            return self.emission_distns[0].\
                log_likelihood((xs, self.data), mask=test_mask).sum()
        else:
            hll = 0
            z = self.stateseq
            for idx, ed in enumerate(self.emission_distns):
                hll += ed.log_likelihood((xs[z == idx], self.data[z == idx]),
                                         mask=test_mask[z == idx]).sum()
949493
def get_global_threshold(threshold_method, image, mask = None, **kwargs):
    """"""
    if mask is not None and not np.any(mask):
        return 1
    
    if threshold_method == TM_OTSU:
        fn = get_otsu_threshold
    elif threshold_method == TM_MOG:
        fn = get_mog_threshold
    elif threshold_method == TM_BACKGROUND:
        fn = get_background_threshold
    elif threshold_method == TM_ROBUST_BACKGROUND:
        fn = get_robust_background_threshold
    elif threshold_method == TM_RIDLER_CALVARD:
        fn = get_ridler_calvard_threshold
    elif threshold_method == TM_KAPUR:
        fn = get_kapur_threshold
    elif threshold_method == TM_MCT:
        fn = get_maximum_correlation_threshold
    else:
        raise NotImplementedError("%s algorithm not implemented"%(threshold_method))
    kwargs = dict([(k, v) for k, v in kwargs.items()
                   if k in fn.args])
    return fn(image, mask, **kwargs)
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))
29896
def verify_consistency(predictions, scores, proba, opt_class):
    
    if predictions.shape != scores.shape:
        raise ValueError("predictions and scores arrays have inconsistent " +
                         "dimensions.")

    n_class = scores.shape[1] if scores.ndim > 1 else 1

    # If proba not given, default to False for all classifiers
    if proba is None:
        proba = np.repeat(False, n_class)

    # If opt_class is not given, default to True for all classifiers
    if opt_class is None:
        opt_class = np.repeat(True, n_class)

    # Convert to numpy arrays if necessary
    proba = np.array(proba, dtype=bool, ndmin=1)
    opt_class = np.array(opt_class, dtype=bool, ndmin=1)

    if np.sum(opt_class) < 1:
        raise ValueError("opt_class should contain at least one True value.")

    if predictions.shape[1] != len(proba):
        raise ValueError("mismatch in shape of proba and predictions.")
    if predictions.shape[1] != len(opt_class):
        raise ValueError("mismatch in shape of opt_class and predictions.")

    for m in range(n_class):
        if (np.any(np.logical_or(scores[:,m] < 0, scores[:,m] > 1)) and proba[m]):
            warnings.warn("scores fall outside the [0,1] interval for " +
                          "classifier {}. Setting proba[m]=False.".format(m))
            proba[m] = False

    return proba, opt_class
913220
def calc_thresholds(rbh, file_name, thresholds = [False, False, False, False], stdevs = 2):
    
    calc_threshold = thresholds[-1]
    norm_threshold = {}
    for pair in itertools.permutations([i for i in rbh], 2):
        if pair[0] not in norm_threshold:
            norm_threshold[pair[0]] = {}
        norm_threshold[pair[0]][pair[1]] = {}
    out = open(file_name, 'w')
    print('#### summary of rbh comparisons\n', file=out)
    comparisons = []
    for genome in rbh:
        for compare in rbh[genome]:
            pair = ''.join(sorted([genome, compare]))
            if pair in comparisons:
                continue
            comparisons.append(pair)
            scores = {'percent identity': [], 'e-value': [], 'bit score': [], 'normalized bit score': [], 'alignment length fraction': []}
            print('### blast between %s and %s\n' % (genome, compare), file=out)
            for id in rbh[genome][compare]:
                pident, length_fraction, e, bit, norm_bit = rbh[genome][compare][id][3:]
                scores['percent identity'].append(pident)
                scores['alignment length fraction'].append(length_fraction)
                scores['e-value'].append(e)
                scores['bit score'].append(bit)
                scores['normalized bit score'].append(norm_bit)
            if calc_threshold is True:
                norms = scores['normalized bit score']
                average = numpy.average(norms) 
                std = numpy.std(norms)
                normal_thresh = average - (std * stdevs)
                print('## average normalized bit score: %s' % average, file=out)
                print('## standard deviation of normalized bit scores: %s' % std, file=out)
                print('## normalized bit score threshold set to: %s\n' % (normal_thresh), file=out)
                norm_threshold[genome][compare], norm_threshold[compare][genome] = normal_thresh, normal_thresh
            for score in scores:
                print('## %s' % (score), file=out)
                if len(scores[score]) > 0:
                    print('## average: %s' % numpy.average(scores[score]), file=out)
#                    hist = histogram(scores[score], [])
#                    for line in hist:
#                        print >> out, line
                print('', file=out)
    out.close()
    if calc_threshold is True:
        return thresholds[0:-1] + [norm_threshold]
    else:
        return thresholds
1016037
def max_normal_germline_depth(in_file, params, somatic_info):
    
    bcf_in = pysam.VariantFile(in_file)
    depths = []
    for rec in bcf_in:
        stats = _is_possible_loh(rec, bcf_in, params, somatic_info)
        if tz.get_in(["normal", "depth"], stats):
            depths.append(tz.get_in(["normal", "depth"], stats))
    if depths:
        return np.median(depths) * NORMAL_FILTER_PARAMS["max_depth_percent"]
981012
def main(argv=None):
    

    if argv is None:
        argv = sys.argv

    # setup command line parser
    parser = U.OptionParser(version="%prog version: $Id$",
                            usage=usage,
                            description=globals()["__doc__"])

    group = U.OptionGroup(parser, "whitelist-specific options")

    group.add_option("--plot-prefix",
                     dest="plot_prefix", type="string",
                     help=("Prefix for plots to visualise the automated "
                           "detection of the number of 'true' cell barcodes"))
    group.add_option("--subset-reads",
                     dest="subset_reads", type="int",
                     help=("Use the first N reads to automatically identify "
                           "the true cell barcodes. If N is greater than the "
                           "number of reads, all reads will be used. "
                           "Default is 100,000,000"))
    group.add_option("--error-correct-threshold",
                     dest="error_correct_threshold",
                     type="int",
                     help=("Hamming distance for correction of barcodes to "
                           "whitelist barcodes. This value will also be used "
                           "for error detection above the knee if required "
                           "(--ed-above-threshold)"))
    group.add_option("--method",
                     dest="method",
                     choices=["reads", "umis"],
                     help=("Use reads or unique umi counts per cell"))
    group.add_option("--knee-method",
                     dest="knee_method",
                     choices=["distance", "density"],
                     help=("Use distance or density methods for detection of knee"))
    group.add_option("--expect-cells",
                     dest="expect_cells",
                     type="int",
                     help=("Prior expectation on the upper limit on the "
                           "number of cells sequenced"))
    group.add_option("--allow-threshold-error",
                     dest="allow_threshold_error", action="store_true",
                     help=("Don't select a threshold. Will still "
                           "output the plots if requested (--plot-prefix)"))
    group.add_option("--set-cell-number",
                     dest="cell_number",
                     type="int",
                     help=("Specify the number of cell barcodes to accept"))

    parser.add_option("--ed-above-threshold",
                      dest="ed_above_threshold", type="choice",
                      choices=["discard", "correct"],
                      help=("Detect CBs above the threshold which may be "
                            "sequence errors from another CB and either "
                            "'discard' or 'correct'. Default=discard"))
    parser.add_option_group(group)

    parser.set_defaults(method="reads",
                        knee_method="distance",
                        extract_method="string",
                        filter_cell_barcodes=False,
                        whitelist_tsv=None,
                        blacklist_tsv=None,
                        error_correct_threshold=1,
                        pattern=None,
                        pattern2=None,
                        read2_in=None,
                        plot_prefix=None,
                        subset_reads=100000000,
                        expect_cells=False,
                        allow_threshold_error=False,
                        cell_number=False,
                        ed_above_threshold=None)

    # add common options (-h/--help, ...) and parse command line

    (options, args) = U.Start(parser, argv=argv,
                              add_extract_options=True,
                              add_group_dedup_options=False,
                              add_umi_grouping_options=False,
                              add_sam_options=False)

    if options.expect_cells:
        if options.knee_method == "distance":
            U.error("Cannot use --expect-cells with 'distance' knee "
                    "method. Switch to --knee-method=density if you want to "
                    "provide an expectation for the number of "
                    "cells. Alternatively, if you know the number of cell "
                    "barcodes, use --cell-number")
        if options.cell_number:
            U.error("Cannot supply both --expect-cells and "
                    "--cell-number options")

    extract_cell, extract_umi = U.validateExtractOptions(options)

    if not extract_cell:
        if options.extract_method == "string":
            U.error("barcode pattern(s) do not include any cell bases "
                    "(marked with 'Cs') %s, %s" % (
                        options.pattern, options.pattern2))
        elif options.extract_method == "regex":
            U.error("barcode regex(es) do not include any cell groups "
                    "(starting with 'cell_') %s, %s" (
                        options.pattern, options.pattern2))

    read1s = umi_methods.fastqIterate(options.stdin)

    # set up read extractor
    ReadExtractor = extract_methods.ExtractFilterAndUpdate(
        method=options.extract_method,
        pattern=options.pattern,
        pattern2=options.pattern2,
        prime3=options.prime3,
        extract_cell=extract_cell)

    cell_barcode_counts = collections.Counter()

    n_reads = 0
    n_cell_barcodes = 0

    # if using the umis method, need to keep a set of umis observed
    if options.method == "umis":
        cell_barcode_umis = collections.defaultdict(set)

    # variables for progress monitor
    displayMax = 100000
    U.info("Starting barcode extraction")

    if not options.read2_in:
        for read1 in read1s:

            # Update display in every 100kth iteration
            if n_reads % displayMax == 0:
                U.info("Parsed {} reads".format(n_reads))

            n_reads += 1
            barcode_values = ReadExtractor.getBarcodes(read1)
            if barcode_values is None:
                continue
            else:
                cell, umi, _, _, _, _, _ = barcode_values
                if options.method == "umis":
                    cell_barcode_umis[cell].add(umi)
                else:
                    cell_barcode_counts[cell] += 1
                n_cell_barcodes += 1

            if options.subset_reads:
                if n_cell_barcodes > options.subset_reads:
                    break
    else:
        read2s = umi_methods.fastqIterate(U.openFile(options.read2_in))
        for read1, read2 in izip(read1s, read2s):

            # Update display in every 100kth iteration
            if n_reads % displayMax == 0:
                U.info("Parsed {} reads".format(n_reads))

            n_reads += 1

            barcode_values = ReadExtractor.getBarcodes(read1, read2)
            if barcode_values is None:
                continue
            else:
                cell, umi, _, _, _, _, _ = barcode_values
                if options.method == "umis":
                    cell_barcode_umis[cell].add(umi)
                else:
                    cell_barcode_counts[cell] += 1
                n_cell_barcodes += 1

            if options.subset_reads:
                if n_reads > options.subset_reads:
                    break

    U.info("Starting - whitelist determination")

    if options.method == "umis":
        for cell in cell_barcode_umis:
            cell_barcode_counts[cell] = len(cell_barcode_umis[cell])

    if options.cell_number and options.cell_number > len(cell_barcode_counts):
        raise ValueError(
            "--set-cell-barcode option specifies more cell barcodes than the "
            "number of observed cell barcodes. This may be because "
            "--subset-reads was set to a value too low to capture reads from "
            "all cells. %s cell barcodes observed from %s parsed reads. "
            "Expected>= %s cell barcodes" % (
                len(cell_barcode_counts),
                options.subset_reads,
                options.cell_number))

    cell_whitelist, true_to_false_map = whitelist_methods.getCellWhitelist(
        cell_barcode_counts,
        options.knee_method,
        options.expect_cells,
        options.cell_number,
        options.error_correct_threshold,
        options.plot_prefix)

    if cell_whitelist:
        U.info("Top %s cell barcodes passed the selected threshold" %
               len(cell_whitelist))

    if options.ed_above_threshold:
        cell_whitelist, true_to_false_map = whitelist_methods.errorDetectAboveThreshold(
            cell_barcode_counts,
            cell_whitelist,
            true_to_false_map,
            errors=options.error_correct_threshold,
            resolution_method=options.ed_above_threshold)

    if cell_whitelist:
        U.info("Writing out whitelist")
        total_correct_barcodes = 0
        total_corrected_barcodes = 0
        for barcode in sorted(list(cell_whitelist)):

            total_correct_barcodes += cell_barcode_counts[barcode]

            if true_to_false_map:
                corrected_barcodes = ",".join(
                    sorted(true_to_false_map[barcode]))

                correct_barcode_counts = [cell_barcode_counts[x] for x in
                                          sorted(true_to_false_map[barcode])]
                total_corrected_barcodes += sum(correct_barcode_counts)

                corrected_barcode_counts = ",".join(
                    map(str, correct_barcode_counts))
            else:
                corrected_barcodes, corrected_barcode_counts = "", ""

            options.stdout.write("%s\t%s\t%s\t%s\n" % (
                barcode, corrected_barcodes, cell_barcode_counts[barcode],
                corrected_barcode_counts))
    else:
        msg = ("No local minima was accepted. Recommend checking the plot "
               "output and counts per local minima (requires `--plot-prefix`"
               "option) and then re-running with manually selected threshold "
               "(`--set-cell-number` option)")

        if options.allow_threshold_error:
            U.info(msg)
        else:
            U.error(msg)

    U.info("Parsed %i reads" % n_reads)
    U.info("%i reads matched the barcode pattern" % n_cell_barcodes)
    U.info("Found %i unique cell barcodes" % len(cell_barcode_counts))

    if cell_whitelist:
        U.info("Found %i total reads matching the selected cell barcodes" %
               total_correct_barcodes)
        U.info("Found %i total reads which can be error corrected to the "
               "selected cell barcodes" % total_corrected_barcodes)

    U.Stop()
1308634
def _run_freebayes_caller(align_bams, items, ref_file, assoc_files,
                          region=None, out_file=None, somatic=None):
    
    config = items[0]["config"]
    if out_file is None:
        out_file = "%s-variants.vcf.gz" % os.path.splitext(align_bams[0])[0]
    if not utils.file_exists(out_file):
        if not utils.file_exists(out_file):
            with file_transaction(items[0], out_file) as tx_out_file:
                freebayes = config_utils.get_program("freebayes", config)
                input_bams = " ".join("-b %s" % x for x in align_bams)
                opts, no_target_regions = _freebayes_options_from_config(items, config, out_file, region)
                if no_target_regions:
                    vcfutils.write_empty_vcf(tx_out_file, config, samples=[dd.get_sample_name(d) for d in items])
                else:
                    opts = " ".join(opts)
                    # Recommended options from 1000 genomes low-complexity evaluation
                    # https://groups.google.com/d/msg/freebayes/GvxIzjcpbas/1G6e3ArxQ4cJ
                    opts += " --min-repeat-entropy 1"
                    # Remove partial observations, which cause a preference for heterozygote calls
                    # https://github.com/ekg/freebayes/issues/234#issuecomment-205331765
                    opts += " --no-partial-observations"
                    if somatic:
                        opts = _add_somatic_opts(opts, somatic)
                    compress_cmd = "| bgzip -c" if out_file.endswith("gz") else ""
                    # For multi-sample outputs, ensure consistent order
                    samples = ("-s" + ",".join([dd.get_sample_name(d) for d in items])) if len(items) > 1 else ""
                    fix_ambig = vcfutils.fix_ambiguous_cl()
                    py_cl = config_utils.get_program("py", config)
                    cmd = ("{freebayes} -f {ref_file} {opts} {input_bams} "
                           """| bcftools filter -i 'ALT="<*>" || QUAL > 5' """
                           "| {fix_ambig} | bcftools view {samples} -a - | "
                           "{py_cl} -x 'bcbio.variation.freebayes.remove_missingalt(x)' | "
                           "vcfallelicprimitives -t DECOMPOSED --keep-geno | vcffixup - | vcfstreamsort | "
                           "vt normalize -n -r {ref_file} -q - | vcfuniqalleles | vt uniq - 2> /dev/null "
                           "{compress_cmd} > {tx_out_file}")
                    do.run(cmd.format(**locals()), "Genotyping with FreeBayes", {})
    return out_file
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
148954
def auto_qc(dset,inside_perc=60,atlas=None,p=0.001):
    ''''''
    with nl.notify('Running quality check on %s:' % dset):
        if not os.path.exists(dset):
            nl.notify('Error: cannot find the file!',level=nl.level.error)
            return False
        
        info = nl.dset_info(dset)
        if not info:
            nl.notify('Error: could not read the dataset!',level=nl.level.error)
        
        if any(['stat' in x for x in info.subbricks]):
            with nl.notify('Statistical results detected...'):
                inside = inside_brain(dset,atlas=atlas,p=p)
                nl.notify('%.1f significant voxels inside brain')
                if inside<inside_perc:
                    nl.notify('Warning: below quality threshold!',level=nl.level.warning)
#                    return False
                nl.notify('Looks ok')
                return inside
        
        if len(info.subbricks)>1:
            with nl.notify('Time-series detected...'):
                return_val = True
                (cost,overlap) = atlas_overlap(dset)
                if cost>0.15 or overlap<80:
                    nl.notify('Warning: does not appear to conform to brain dimensions',level=nl.level.warning)
                    return_val = False
                if len(info.subbricks)>5:
                    (oc,perc_outliers) = outcount(dset)
                    if perc_outliers>0.1:
                        nl.notify('Warning: large amount of outlier time points',level=nl.level.warning)
                        return_val = False
            if return_val:
                nl.notify('Looks ok')
                return min(100*(1-cost),overlap,100*perc_outliers)
            return False
        
        with nl.notify('Single brain image detected...'):
            (cost,overlap) = atlas_overlap(dset)
            # Be more lenient if it's not an EPI, cuz who knows what else is in this image
            if cost>0.45 or overlap<70:
                nl.notify('Warning: does not appear to conform to brain dimensions',level=nl.level.warning)
                return False
            nl.notify('Looks ok')
            return min(100*(1-cost),overlap)