>>> encode_genotype("1/1:128,18,0:6:18")  # homozygote B
    'B'
    >>> encode_genotype("0/1:0,0,0:0:3")      # missing data
    '-'
    >>> encode_genotype("0/1:128,0,26:7:22")  # heterozygous A/B
    'X'
777711
def encode_genotype(s, mindepth=3, depth_index=2, nohet=False):
    
    atoms = s.split(":")
    if len(atoms) < 3:
        return g2x[atoms[0]]

    inferred = atoms[0]
    depth = int(atoms[depth_index])
    if depth < mindepth:
        return '-'
    if inferred == '0/0':
        return 'A'
    if inferred == '0/1':
        return '-' if nohet else 'X'
    if inferred == '1/1':
        return 'B'
    return '-'
1151992
def _prep_genome(out_dir, data):
    
    genome_name = utils.splitext_plus(os.path.basename(dd.get_ref_file(data)))[0]
    out_dir = utils.safe_makedir(os.path.join(out_dir, genome_name))
    ref_file = dd.get_ref_file(data)
    utils.symlink_plus(ref_file, os.path.join(out_dir, os.path.basename(ref_file)))
    with open(os.path.join(out_dir, "GenomeSize.xml"), "w") as out_handle:
        out_handle.write('<sequenceSizes genomeName="%s">' % genome_name)
        for c in pysam.AlignmentFile("%s.dict" % utils.splitext_plus(ref_file)[0]).header["SQ"]:
            cur_ploidy = ploidy.get_ploidy([data], region=[c["SN"]])
            out_handle.write('<chromosome fileName="%s" contigName="%s" totalBases="%s" knownBases="%s" '
                             'isCircular="false" ploidy="%s" md5="%s"/>' %
                             (os.path.basename(ref_file), c["SN"], c["LN"], c["LN"], cur_ploidy, c["M5"]))
        out_handle.write('</sequenceSizes>')
    return out_dir
1276336
def _identify_heterogeneity_blocks_seg(in_file, seg_file, params, work_dir, somatic_info):
    
    def _segment_by_cns(target_chrom, freqs, coords):
        with open(seg_file) as in_handle:
            reader = csv.reader(in_handle, dialect="excel-tab")
            next(reader)  # header
            for cur_chrom, start, end in (xs[:3] for xs in reader):
                if cur_chrom == target_chrom:
                    block_freqs = []
                    for i, (freq, coord) in enumerate(zip(freqs, coords)):
                        if coord >= int(start) and coord < int(end):
                            block_freqs.append(freq)
                        elif coord >= int(end):
                            break
                    coords = coords[max(0, i - 1):]
                    freqs = freqs[max(0, i - 1):]
                    if len(block_freqs) > params["hetblock"]["min_alleles"]:
                        yield start, end
    return _identify_heterogeneity_blocks_shared(in_file, _segment_by_cns, params, work_dir, somatic_info)
1071804
def iter_acgt_geno(self):
        
        # Need to iterate over itself, and modify the actual genotypes
        for i, (marker, geno) in enumerate(self.iter_geno()):
            yield marker, self._allele_encoding[i][geno]
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)
91845
def get_variants(data, include_germline=False):
    
    data = utils.deepish_copy(data)
    supported = ["precalled", "vardict", "vardict-java", "vardict-perl",
                 "freebayes", "octopus", "strelka2"]
    # Right now mutect2 and mutect do not provide heterozygous germline calls
    # to be useful https://github.com/bcbio/bcbio-nextgen/issues/2464
    # supported += ["mutect2", "mutect"]
    if include_germline:
        supported.insert(1, "gatk-haplotype")
    out = []
    # CWL based input
    if isinstance(data.get("variants"), dict) and "samples" in data["variants"]:
        cur_vs = []
        # Unpack single sample list of files
        if (isinstance(data["variants"]["samples"], (list, tuple)) and
              len(data["variants"]["samples"]) == 1 and isinstance(data["variants"]["samples"][0], (list, tuple))):
            data["variants"]["samples"] = data["variants"]["samples"][0]
        for fname in data["variants"]["samples"]:
            variantcaller = utils.splitext_plus(os.path.basename(fname))[0]
            variantcaller = variantcaller.replace(dd.get_sample_name(data) + "-", "")
            for batch in dd.get_batches(data):
                variantcaller = variantcaller.replace(batch + "-", "")
            cur_vs.append({"vrn_file": fname, "variantcaller": variantcaller})
        data["variants"] = cur_vs
    for v in data.get("variants", []):
        if v["variantcaller"] in supported and v.get("vrn_file"):
            out.append((supported.index(v["variantcaller"]), v))
    out.sort()
    return [xs[1] for xs in out]
427795
def loss(args):
    
    p = OptionParser(loss.__doc__)
    p.add_option("--bed", default=False, action="store_true",
                 help="Genomic BLAST is in bed format [default: %default]")
    p.add_option("--gdist", default=20, type="int",
                 help="Gene distance [default: %default]")
    p.add_option("--bdist", default=20000, type="int",
                 help="Base pair distance [default: %default]")
    p.set_beds()
    opts, args = p.parse_args(args)

    if len(args) not in (1, 2):
        sys.exit(not p.print_help())

    blocksfile = args[0]
    emptyblast = (len(args) == 1)
    if emptyblast:
        genomicblast = "empty.blast"
        sh("touch {0}".format(genomicblast))
    else:
        genomicblast = args[1]

    gdist, bdist = opts.gdist, opts.bdist
    qbed, sbed, qorder, sorder, is_self = check_beds(blocksfile, p, opts)
    blocks = []
    fp = open(blocksfile)
    genetrack = {}
    proxytrack = {}
    for row in fp:
        a, b = row.split()
        genetrack[a] = b
        blocks.append((a, b))

    data = []
    for key, rows in groupby(blocks, key=lambda x: x[-1]):
        rows = list(rows)
        data.append((key, rows))

    imax = len(data) - 1
    for i, (key, rows) in enumerate(data):
        if i == 0 or i == imax:
            continue
        if key != '.':
            continue

        before, br = data[i - 1]
        after, ar = data[i + 1]
        bi, bx = sorder[before]
        ai, ax = sorder[after]
        dist = abs(bi - ai)
        if bx.seqid != ax.seqid or dist > gdist:
            continue

        start, end = range_minmax(((bx.start, bx.end), (ax.start, ax.end)))
        start, end = max(start - bdist, 1), end + bdist
        proxy = (bx.seqid, start, end)
        for a, b in rows:
            proxytrack[a] = proxy

    tags = {}
    if opts.bed:
        bed = Bed(genomicblast, sorted=False)
        key = lambda x: gene_name(x.accn.rsplit(".", 1)[0])
        for query, bb in groupby(bed, key=key):
            bb = list(bb)
            if query not in proxytrack:
                continue

            proxy = proxytrack[query]
            tag = "NS"
            best_b = bb[0]
            for b in bb:
                hsp = (b.seqid, b.start, b.end)
                if range_overlap(proxy, hsp):
                    tag = "S"
                    best_b = b
                    break

            hsp = (best_b.seqid, best_b.start, best_b.end)
            proxytrack[query] = hsp
            tags[query] = tag

    else:
        blast = Blast(genomicblast)
        for query, bb in blast.iter_hits():
            bb = list(bb)
            query = gene_name(query)
            if query not in proxytrack:
                continue

            proxy = proxytrack[query]
            tag = "NS"
            best_b = bb[0]
            for b in bb:
                hsp = (b.subject, b.sstart, b.sstop)
                if range_overlap(proxy, hsp):
                    tag = "S"
                    best_b = b
                    break

            hsp = (best_b.subject, best_b.sstart, best_b.sstop)
            proxytrack[query] = hsp
            tags[query] = tag

    for b in qbed:
        accn = b.accn
        target_region = genetrack[accn]
        if accn in proxytrack:
            target_region = region_str(proxytrack[accn])
            if accn in tags:
                ptag = "[{0}]".format(tags[accn])
            else:
                ptag = "[NF]"
            target_region = ptag + target_region

        print("\t".join((b.seqid, accn, target_region)))

    if emptyblast:
        sh("rm -f {0}".format(genomicblast))
671759
def _identify_heterogeneity_blocks_shared(in_file, segment_fn, params, work_dir, somatic_info):
    
    out_file = os.path.join(work_dir, "%s-hetblocks.bed" % utils.splitext_plus(os.path.basename(in_file))[0])
    if not utils.file_uptodate(out_file, in_file):
        with file_transaction(somatic_info.tumor_data, out_file) as tx_out_file:
            with open(tx_out_file, "w") as out_handle:
                for chrom, freqs, coords in _freqs_by_chromosome(in_file, params, somatic_info):
                    for start, end in segment_fn(chrom, freqs, coords):
                        out_handle.write("%s\t%s\t%s\n" % (chrom, start, end))
    return out_file
420092
def filter_missing(self):
        """"""

        missing             = None
        locus_count         = 0

        # Filter out individuals according to missingness
        self.genotype_file.seek(0)
        for genotypes in self.genotype_file:
            genotypes = genotypes.split()
            chr, rsid, junk, pos = genotypes[0:4]
            if DataParser.boundary.TestBoundary(chr, pos, rsid):
                locus_count += 1
                allelic_data = numpy.array(genotypes[4:], dtype="S2").reshape(-1, 2)
                if missing is None:
                    missing             = numpy.zeros(allelic_data.shape[0], dtype='int8')
                missing += (numpy.sum(0+(allelic_data==DataParser.missing_representation), axis=1)/2)

        max_missing = DataParser.ind_miss_tol * locus_count
        dropped_individuals = 0+(max_missing<missing)

        self.ind_mask[:,0] = self.ind_mask[:,0]|dropped_individuals
        self.ind_mask[:,1] = self.ind_mask[:,1]|dropped_individuals

        valid_individuals = numpy.sum(self.ind_mask==0)
        max_missing = DataParser.snp_miss_tol * valid_individuals

        self.locus_count = 0
        # We can't merge these two iterations since we need to know which individuals
        # to consider for filtering on MAF
        dropped_snps = []
        self.genotype_file.seek(0)
        for genotypes in self.genotype_file:
            genotypes = genotypes.split()
            chr, rsid, junk, pos = genotypes[0:4]
            chr = int(chr)
            pos = int(pos)
            if DataParser.boundary.TestBoundary(chr, pos, rsid):
                allelic_data = numpy.ma.MaskedArray(numpy.array(genotypes[4:], dtype="S2").reshape(-1, 2), self.ind_mask).compressed()
                missing = numpy.sum(0+(allelic_data==DataParser.missing_representation))
                if missing > max_missing:
                    DataParser.boundary.dropped_snps[int(chr)].add(int(pos))
                    dropped_snps.append(rsid)
                else:
                    self.locus_count += 1
872958
function encode(data, structure, separator) {
	var encoded = data.split('').map(function (val, idx) {
		return BINARIES[structure[idx]];
	}).map(function (val, idx) {
		return val ? val[data[idx]] : '';
	});

	if (separator) {
		var last = data.length - 1;
		encoded = encoded.map(function (val, idx) {
			return idx < last ? val + separator : val;
		});
	}

	return encoded.join('');
}
448978
def n_to_g(self, n_interval):
        """"""

        frs, fre = n_interval.start.base - 1, n_interval.end.base - 1
        start_offset, end_offset = n_interval.start.offset, n_interval.end.offset

        if self.strand == -1:
            fre, frs = self.tgt_len - frs - 1, self.tgt_len - fre - 1
            start_offset, end_offset = -end_offset, -start_offset

        # returns the genomic range start (grs) and end (gre)
        grs, _, grs_cigar = self._map(
            from_pos=self.tgt_pos, to_pos=self.ref_pos, pos=frs, base="start")
        gre, _, gre_cigar = self._map(
            from_pos=self.tgt_pos, to_pos=self.ref_pos, pos=fre, base="end")
        grs, gre = grs + self.gc_offset + 1, gre + self.gc_offset + 1
        gs, ge = grs + start_offset, gre + end_offset

        # The returned interval would be uncertain when locating at alignment gaps
        return hgvs.location.Interval(
            start=hgvs.location.SimplePosition(gs, uncertain=n_interval.start.uncertain),
            end=hgvs.location.SimplePosition(ge, uncertain=n_interval.end.uncertain),
            uncertain=grs_cigar in 'DI' or gre_cigar in 'DI')
459244
def variant_sequences_generator_to_dataframe(variant_sequences_generator):
    
    # TODO: Change VariantSequence.alt to VariantSequence.alt_nucleotides
    # or something else that doesn't clash with a variant's `alt` field
    return dataframe_from_generator(
        VariantSequence,
        variant_sequences_generator,
        rename_dict={"alt": "allele"},
        extra_column_fns={
            "gene": lambda variant, _: ";".join(variant.gene_names),
        })
136930
def parse_sampleinfo(data: dict) -> dict:
    

    genome_build = data['human_genome_build']
    genome_build_str = f"{genome_build['source']}{genome_build['version']}"
    if 'svdb' in data['program']:
        svdb_outpath = (f"{data['program']['svdb']['path']}")
    else:
        svdb_outpath = ''
    outdata = {
        'date': data['analysis_date'],
        'family': data['family'],
        'genome_build': genome_build_str,
        'rank_model_version': data['program']['genmod']['rank_model']['version'],
        'is_finished': True if data['analysisrunstatus'] == 'finished' else False,
        'pedigree_path': data['pedigree_minimal'],
        'peddy': {
            'ped': (data['program']['peddy']['peddy']['path'] if
                    'peddy' in data['program'] else None),
            'ped_check': (data['program']['peddy']['ped_check']['path'] if
                          'peddy' in data['program'] else None),
            'sex_check': (data['program']['peddy']['sex_check']['path'] if
                          'peddy' in data['program'] else None),
        },
        'qcmetrics_path': data['program']['qccollect']['path'],
        'samples': [],
        'snv': {
            'bcf': data['most_complete_bcf']['path'],
            'clinical_vcf': data['vcf_binary_file']['clinical']['path'],
            'gbcf': data['gbcf_file']['path'],
            'research_vcf': data['vcf_binary_file']['research']['path'],
        },
        'svdb_outpath': svdb_outpath,
        'sv': {
            'bcf': data.get('sv_bcf_file', {}).get('path'),
            'clinical_vcf': (data['sv_vcf_binary_file']['clinical']['path'] if
                             'sv_vcf_binary_file' in data else None),
            'merged': svdb_outpath,
            'research_vcf': (data['sv_vcf_binary_file']['research']['path'] if
                             'sv_vcf_binary_file' in data else None),
        },
        'version': data['mip_version'],
    }

    for sample_id, sample_data in data['sample'].items():
        sample = {
            'id': sample_id,
            'bam': sample_data['most_complete_bam']['path'],
            'sambamba': list(sample_data['program']['sambamba_depth'].values())[0]['path'],
            'sex': sample_data['sex'],
            # subsample mt is only for wgs data
            'subsample_mt': (list(sample_data['program']['samtools_subsample_mt'].values())[0]['path'] if
                             'samtools_subsample_mt' in sample_data['program'] else None),
            'vcf2cytosure': list(sample_data['program']['vcf2cytosure'].values())[0]['path'],
        }
        chanjo_sexcheck = list(sample_data['program']['chanjo_sexcheck'].values())[0]
        sample['chanjo_sexcheck'] = chanjo_sexcheck['path']
        outdata['samples'].append(sample)

    return outdata
479355
def group_alleles_by_start_end_Xbp(arr, bp=28):
    
    starts = arr[:,0:bp]
    ends = arr[:,-bp:]
    starts_ends_idxs = defaultdict(list)
    l, seq_len = arr.shape
    for i in range(l):
        start_i = starts[i]
        end_i = ends[i]
        start_i_str = ''.join([str(x) for x in start_i])
        end_i_str = ''.join([str(x) for x in end_i])
        starts_ends_idxs[start_i_str + end_i_str].append(i)
    return starts_ends_idxs
478395
def _compatible_small_variants(data, items):
    
    from bcbio import heterogeneity
    VarFile = collections.namedtuple("VarFile", ["name", "sample", "normal"])
    out = []
    paired = vcfutils.get_paired(items)
    for v in heterogeneity.get_variants(data, include_germline=not paired):
        vrn_file = v["vrn_file"]
        base, ext = utils.splitext_plus(os.path.basename(vrn_file))
        if paired:
            out.append(VarFile(vrn_file, paired.tumor_name, paired.normal_name))
        else:
            out.append(VarFile(vrn_file, dd.get_sample_name(data), None))
    return out
1310248
def _encode_binary(name, value, dummy0, dummy1):
    """"""
    subtype = value.subtype
    if subtype == 2:
        value = _PACK_INT(len(value)) + value
    return b"\x05" + name + _PACK_LENGTH_SUBTYPE(len(value), subtype) + value
197293
def encode(data: Union[str, bytes]) -> str:
        
        return ensure_str(base58.b58encode(ensure_bytes(data)))