r"""
    Parameters
    ----------
    bim : pandas.DataFrame
        Variant annotation
    geno_range : tuple
        (chrom, pos_start, pos_end)

    Returns
    -------
    Isnp : bool array
        Variant filter
495833
def is_in(bim, geno_range):
    
    """
    Ichrom = bim.chrom == geno_range[0]
    Ipos1 = bim.pos >= geno_range[1]
    Ipos2 = bim.pos < geno_range[2]
    return Ichrom & Ipos1 & Ipos2
1098679
def _get_variant_regions(items):
    
    return list(filter(lambda x: x is not None,
                       [tz.get_in(("config", "algorithm", "variant_regions"), data)
                        for data in items
                        if tz.get_in(["config", "algorithm", "coverage_interval"], data) != "genome"]))
95450
def postprocess_variants(items):
    
    vrn_key = "vrn_file"
    if not isinstance(items, dict):
        items = [utils.to_single_data(x) for x in items]
        if "vrn_file_joint" in items[0]:
            vrn_key = "vrn_file_joint"
    data, items = _get_batch_representative(items, vrn_key)
    items = cwlutils.unpack_tarballs(items, data)
    data = cwlutils.unpack_tarballs(data, data)
    cur_name = "%s, %s" % (dd.get_sample_name(data), get_variantcaller(data, require_bam=False))
    logger.info("Finalizing variant calls: %s" % cur_name)
    orig_vrn_file = data.get(vrn_key)
    data = _symlink_to_workdir(data, [vrn_key])
    data = _symlink_to_workdir(data, ["config", "algorithm", "variant_regions"])
    if data.get(vrn_key):
        logger.info("Calculating variation effects for %s" % cur_name)
        ann_vrn_file, vrn_stats = effects.add_to_vcf(data[vrn_key], data)
        if ann_vrn_file:
            data[vrn_key] = ann_vrn_file
        if vrn_stats:
            data["vrn_stats"] = vrn_stats
        orig_items = _get_orig_items(items)
        logger.info("Annotate VCF file: %s" % cur_name)
        data[vrn_key] = annotation.finalize_vcf(data[vrn_key], get_variantcaller(data, require_bam=False),
                                                orig_items)
        if cwlutils.is_cwl_run(data):
            logger.info("Annotate with population level variation data")
            ann_file = population.run_vcfanno(data[vrn_key], data)
            if ann_file:
                data[vrn_key] = ann_file
        logger.info("Filtering for %s" % cur_name)
        data[vrn_key] = variant_filtration(data[vrn_key], dd.get_ref_file(data),
                                           tz.get_in(("genome_resources", "variation"), data, {}),
                                           data, orig_items)
        logger.info("Prioritization for %s" % cur_name)
        prio_vrn_file = prioritize.handle_vcf_calls(data[vrn_key], data, orig_items)
        if prio_vrn_file != data[vrn_key]:
            data[vrn_key] = prio_vrn_file
            logger.info("Germline extraction for %s" % cur_name)
            data = germline.extract(data, orig_items)

        if dd.get_align_bam(data):
            data = damage.run_filter(data[vrn_key], dd.get_align_bam(data), dd.get_ref_file(data),
                                     data, orig_items)
    if orig_vrn_file and os.path.samefile(data[vrn_key], orig_vrn_file):
        data[vrn_key] = orig_vrn_file
    return [[data]]
1280576
def subset_variant_regions(variant_regions, region, out_file, items=None, do_merge=True, data=None):
    
    if region is None:
        return variant_regions
    elif variant_regions is None:
        return region
    elif not isinstance(region, (list, tuple)) and region.find(":") > 0:
        raise ValueError("Partial chromosome regions not supported")
    else:
        merge_text = "-unmerged" if not do_merge else ""
        subset_file = "{0}".format(utils.splitext_plus(out_file)[0])
        subset_file += "%s-regions.bed" % (merge_text)
        if not os.path.exists(subset_file):
            data = items[0] if items else data
            with file_transaction(data, subset_file) as tx_subset_file:
                if isinstance(region, (list, tuple)):
                    _subset_bed_by_region(variant_regions, tx_subset_file, to_multiregion(region),
                                          dd.get_ref_file(data), do_merge=do_merge)
                else:
                    _rewrite_bed_with_chrom(variant_regions, tx_subset_file, region)
        if os.path.getsize(subset_file) == 0:
            return region
        else:
            return subset_file
640550
def filter_genes_cv_deprecated(X, Ecutoff, cvFilter):
    
    if issparse(X):
        raise ValueError('Not defined for sparse input. See `filter_genes_dispersion`.')
    mean_filter = np.mean(X, axis=0) > Ecutoff
    var_filter = np.std(X, axis=0) / (np.mean(X, axis=0) + .0001) > cvFilter
    gene_subset = np.nonzero(np.all([mean_filter, var_filter], axis=0))[0]
    return gene_subset
225673
def has_variant_regions(items, base_file, chrom=None):
    
    if chrom:
        all_vrs = _get_variant_regions(items)
        if len(all_vrs) > 0:
            test = shared.subset_variant_regions(tz.first(all_vrs), chrom, base_file, items)
            if test == chrom:
                return False
    return True
953136
def compare_variant_type_plot(data):
    
    keys = OrderedDict()
    keys['snps'] = {'name': 'SNPs', 'color': '#7cb5ec'}
    keys['indels'] = {'name': 'InDels', 'color': '#90ed7d'}
    keys['multiallelic_snps'] = {'name': 'multi-allelic SNP', 'color': 'orange'}
    keys['complex_indels'] = {'name': 'Complex InDels', 'color': '#8085e9'}

    total_variants = dict()
    known_variants = dict()
    novel_variants = dict()
    for s_name, values in data.items():
        total_variants[s_name] = {
            'snps': values['TOTAL_SNPS'],
            'indels': values['TOTAL_INDELS'],
            'multiallelic_snps': values['TOTAL_MULTIALLELIC_SNPS'],
            'complex_indels': values['TOTAL_COMPLEX_INDELS'],
        }

        known_variants[s_name] = {
            'snps': values['NUM_IN_DB_SNP'],
            'indels': int(values['TOTAL_INDELS'])-int(values['NOVEL_INDELS']),
            'multiallelic_snps': values['NUM_IN_DB_SNP_MULTIALLELIC'],
            'complex_indels': values['NUM_IN_DB_SNP_COMPLEX_INDELS'],
        }

        novel_variants[s_name] = {
            'snps': values['NOVEL_SNPS'],
            'indels': values['NOVEL_INDELS'],
            'multiallelic_snps': int(values['TOTAL_MULTIALLELIC_SNPS'])-int(values['NUM_IN_DB_SNP_MULTIALLELIC']),
            'complex_indels': int(values['TOTAL_COMPLEX_INDELS'])-int(values['NUM_IN_DB_SNP_COMPLEX_INDELS']),
        }

    plot_conf = {
        'id': 'picard_variantCallingMetrics_variant_type',
        'title': 'Picard: Variants Called',
        'ylab': 'Counts of Variants',
        'hide_zero_cats': False,
        'data_labels': [
            {'name': 'Total'},
            {'name': 'Known'},
            {'name': 'Novel'}
        ],
    }
    return bargraph.plot(data=[total_variants, known_variants, novel_variants],
                         cats=[keys, keys, keys],
                         pconfig=plot_conf)
432595
def validate_mutations(vcffile, bams, reffa, chr_reffa, samples, output_format, outfile,
                       mpileup_parameters=mpileup_parser.MPILEUP_DEFAULT_PARAMS):
    """"""
    output_header = "sample chrom pos ref cov A C G T * - + " \
        "val_ref val_alt val_al_type val_al_count val_maf "\
        "most_common_indel most_common_indel_count most_common_indel_maf most_common_indel_type most_common_al " \
        "most_common_al_count most_common_al_maf most_common_count most_common_maf".split()

    # for backwards compatibility
    # if bam or samples is a string, convert to list instead
    if isinstance(samples, six.string_types):
        samples = [samples]
    if isinstance(bams, six.string_types):
        bams = [bams]

    if output_format == 'vcf':
        vcf_reader = vcf.Reader(open(vcffile))
        vcf_reader.samples = samples
        vcf_reader.formats['GT'] = vcf.parser._Format(id='GT', num=1, type='String', desc="Genotype")
        vcf_reader.formats['AD'] = vcf.parser._Format(id='AD', num='R', type='Integer', desc="Allelic depth")
        vcf_reader.formats['DP'] = vcf.parser._Format(id='DP', num=1, type='Integer', desc="Depth")
        vcf_writer = vcf.Writer(outfile, vcf_reader)
    else:
        vcf_reader = open(vcffile)


    if output_format == "sufam":
        outfile.write("\t".join(output_header))
        outfile.write("\n")
    for record in vcf_reader:
        if output_format != 'vcf':
            line = record
            if line.startswith("#CHROM"):
                header = line[1:].rstrip('\n').split("\t")
                # create spoof pyvcf record if vcf_reader is not used
                _Record = namedtuple('Record', header)
            if line.startswith("#"):
                continue
            if len(header) == 0:
                raise(Exception("No header found in vcf file #CHROM not found"))
            # zip all column values, except alt (needs to be list in pyvcf)
            record_args = dict(zip(header, line.rstrip('\n').split("\t")))
            record_args['ALT'] = [record_args['ALT']]
            record = _Record(**record_args)

        # determine type of mutation
        record_type = "snv"
        if len(record.ALT) > 1:
            warnings.warn("Multiple ALT in one record is not implemented - using first")
        if len(record.REF) > len(record.ALT[0]):
            record_type = "deletion"
        elif len(record.ALT[0]) > len(record.REF):
            record_type = "insertion"

        # no coverage results
        no_cov = pd.Series({
            "chrom": str(record.CHROM), "pos": str(record.POS),
            "ref": str(record.REF),
            "cov": 0, "A": 0, "C": 0, "G": 0, "T": 0,
            "val_ref": str(record.REF), "val_alt": str(record.ALT[0]),
            "val_al_type": record_type, "val_al_count": 0, "val_maf": 0})

        # collect mpileup baseparser results per bam
        bps = []
        for i, bam in enumerate(bams):
            sample = samples[i]
            no_cov['sample'] = sample
            bp_lines = mpileup_parser.run_and_parse(bam, str(record.CHROM), str(record.POS), str(record.POS), reffa,
                                                    chr_reffa, mpileup_parameters)
            bpdf = get_baseparser_extended_df(sample, bp_lines, str(record.REF), str(record.ALT[0]))
            if bpdf is None:
                bp = no_cov
            else:
                bp = bpdf.ix[0, :]
            bps += [bp]

        # output call
        if output_format == "vcf":
            _write_bp_vcf(outfile, bps, vcf_writer, record)
        else:
            # only one bam file supported for outputs other than vcf
            _write_bp(outfile, bps[0], output_header, output_format)
28154
def get_variant_slice(self, package_name, range_):
        
        variant_list = self.variant_lists.get(package_name)

        if variant_list is None:
            variant_list = _PackageVariantList(package_name, self.solver)
            self.variant_lists[package_name] = variant_list

        entries = variant_list.get_intersection(range_)
        if not entries:
            return None

        slice_ = _PackageVariantSlice(package_name,
                                      entries=entries,
                                      solver=self.solver)
        return slice_
599150
def get_type(data):
    
    if data["analysis"].lower().startswith("var") or dd.get_variantcaller(data):
        return tz.get_in(("config", "algorithm", "effects"), data, "snpeff")
299884
def _load_gtf_as_dataframe(self, usecols=None, features=None):
        
        logger.info("Reading GTF from %s", self.gtf_path)
        df = read_gtf(
            self.gtf_path,
            column_converters={
                "seqname": normalize_chromosome,
                "strand": normalize_strand,
            },
            infer_biotype_column=True,
            usecols=usecols,
            features=features)

        column_names = set(df.keys())
        expect_gene_feature = features is None or "gene" in features
        expect_transcript_feature = features is None or "transcript" in features
        observed_features = set(df["feature"])

        # older Ensembl releases don't have "gene" or "transcript"
        # features, so fill in those rows if they're missing
        if expect_gene_feature and "gene" not in observed_features:
            # if we have to reconstruct gene feature rows then
            # fill in values for 'gene_name' and 'gene_biotype'
            # but only if they're actually present in the GTF
            logger.info("Creating missing gene features...")
            df = create_missing_features(
                dataframe=df,
                unique_keys={"gene": "gene_id"},
                extra_columns={
                    "gene": {
                        "gene_name",
                        "gene_biotype"
                    }.intersection(column_names),
                },
                missing_value="")
            logger.info("Done.")

        if expect_transcript_feature and "transcript" not in observed_features:
            logger.info("Creating missing transcript features...")
            df = create_missing_features(
                dataframe=df,
                unique_keys={"transcript": "transcript_id"},
                extra_columns={
                    "transcript": {
                        "gene_id",
                        "gene_name",
                        "gene_biotype",
                        "transcript_name",
                        "transcript_biotype",
                        "protein_id",
                    }.intersection(column_names)
                },
                missing_value="")
            logger.info("Done.")

        return df
92572
def construction_error(tot_items, block_shape, axes, e=None):
    
    passed = tuple(map(int, [tot_items] + list(block_shape)))
    # Correcting the user facing error message during dataframe construction
    if len(passed) <= 2:
        passed = passed[::-1]

    implied = tuple(len(ax) for ax in axes)
    # Correcting the user facing error message during dataframe construction
    if len(implied) <= 2:
        implied = implied[::-1]

    if passed == implied and e is not None:
        raise e
    if block_shape[0] == 0:
        raise ValueError("Empty data passed with indices specified.")
    raise ValueError("Shape of passed values is {0}, indices imply {1}".format(
        passed, implied))
1113020
def genotype_to_df(g, samples, as_string=False):
    
    name = g.variant.name if g.variant.name else "genotypes"
    df = pd.DataFrame(g.genotypes, index=samples, columns=[name])

    if as_string:
        df["alleles"] = None

        hard_calls = df[name].round()
        df.loc[hard_calls == 0, "alleles"] = "{0}/{0}".format(g.reference)
        df.loc[hard_calls == 1, "alleles"] = "{0}/{1}".format(g.reference,
                                                              g.coded)
        df.loc[hard_calls == 2, "alleles"] = "{0}/{0}".format(g.coded)

        df = df[["alleles"]]
        df.columns = [name]

    return df
1345613
def gff(args):
    
    from jcvi.utils.range import range_minmax
    from jcvi.formats.gff import valid_gff_parent_child, valid_gff_type

    p = OptionParser(gff.__doc__)
    p.add_option("--source", default=None, help="Specify GFF source." +
                " By default, it picks algorithm used to generate btab file." +
                " [default: %default]")
    p.add_option("--type", default="protein_match", choices=valid_gff_type,
                help="GFF feature type [default: %default]")

    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(not p.print_help())

    btabfile, = args
    btabdict = {}
    btab = Btab(btabfile, aat_dialect=True)
    osource = opts.source or "aat"
    otype = opts.type
    octype = valid_gff_parent_child[otype]
    for b in btab:
        nargs = b.nargs
        id = b.query + "-" + otype + "{0:05d}".format(b.chainNum)
        key = b.key
        if key not in btabdict:
            btabdict[key] = { 'id': id,
                              'method': b.method,
                              'query': b.query,
                              'subject': b.subject,
                              'strand': b.qStrand,
                              'sDesc': b.sDesc,
                              'coords': [],
                              'children': []
                            }

        btabdict[key]['coords'].append((b.qStart, b.qStop))
        btabdict[key]['children'].append(b.gffline(source=osource, type=octype, id=id))

    for v in btabdict.itervalues():
        b = BtabLine("\t".join(str(x) for x in [0] * nargs), aat_dialect=True)
        id = v['id']
        b.query = v['query']
        b.method = v['method']
        b.subject = v['subject']
        b.qStrand = v['strand']
        b.sDesc = v['sDesc']
        b.qStart, b.qStop = range_minmax(v['coords'])
        print(b.gffline(source=osource, type=otype, primary_tag="ID", id=id))
        print("\n".join(v['children']))
842874
def _check_variantcaller(item):
    
    allowed = set(list(genotype.get_variantcallers().keys()) + [None, False])
    vcs = item["algorithm"].get("variantcaller")
    if not isinstance(vcs, dict):
        vcs = {"variantcaller": vcs}
    for vc_set in vcs.values():
        if not isinstance(vc_set, (tuple, list)):
            vc_set = [vc_set]
        problem = [x for x in vc_set if x not in allowed]
        if len(problem) > 0:
            raise ValueError("Unexpected algorithm 'variantcaller' parameter: %s\n"
                             "Supported options: %s\n" % (problem, sorted(list(allowed))))
    # Ensure germline somatic calling only specified with tumor/normal samples
    if "germline" in vcs or "somatic" in vcs:
        paired = vcfutils.get_paired_phenotype(item)
        if not paired:
            raise ValueError("%s: somatic/germline calling in 'variantcaller' "
                             "but tumor/normal metadata phenotype not specified" % dd.get_sample_name(item))
763211
def true_range(arg, high_col='high', low_col='low', close_col='close', skipna=0):
    
    _ensure_col(arg, high_col=high_col, low_col=low_col, close_col=close_col)
    yclose = arg[close_col].shift(1)
    low, high = arg[low_col], arg[high_col]
    mx = pd.DataFrame({'a': high, 'b': yclose}).max(axis=1, skipna=skipna)
    mn = pd.DataFrame({'a': low, 'b': yclose}).min(axis=1, skipna=skipna)
    result = mx - mn
    return pd.Series(result, index=arg.index, name='true_range')
1305445
def annotate_variant(variant, var_obj=None):
    
    if var_obj:
    
        variant.INFO['Obs'] = var_obj['observations']
        if var_obj.get('homozygote'):
            variant.INFO['Hom'] = var_obj['homozygote']
        if var_obj.get('hemizygote'):
            variant.INFO['Hem'] = var_obj['hemizygote']
    
    return variant