Sequence of number of sites of each type associated with the Poscar.
        Similar to 7th line in vasp 5+ Xdatcar.
168258
def natoms(self):
        
        syms = [site.specie.symbol for site in self.structures[0]]
        return [len(tuple(a[1])) for a in itertools.groupby(syms)]
1197641
def calc_frequencies(genomes, bp_table, min_cov, min_per):
    
    nucs = ['A', 'T', 'G', 'C', 'N']
    if bp_table is not False:
        bp_table = open(bp_table, 'w')
        header = ['#genome', 'contig', 'sample', 'position', \
                    'reference', 'ref. frequency', \
                    'consensus', 'con. frequency', \
                    'A', 'T', 'G', 'C', 'N', '# insertions', '# deletions']
        print('\t'.join(header), file=bp_table)
    for genome, contigs in list(genomes.items()):
        for contig, samples in list(contigs.items()):
            for sample, stats in list(samples.items()):
                for pos, ps in enumerate(stats['bp_stats'], 1):
                    coverage = sum([ps[nuc] for nuc in nucs])
                    for nuc in nucs:
                        # make sure support for base passes thresholds
                        nuc_cov = ps[nuc]
                        if coverage == 0:
                            nuc_per = 0
                        else:
                            nuc_per = (float(nuc_cov)/coverage)*100
                        if nuc_cov < min_cov or nuc_per < min_per:
                            del ps[nuc]
                    ps = find_consensus(ps)
                    genomes[genome][contig][sample][pos] = ps
                    if bp_table is not False:
                        out = [genome, contig, sample, pos]
                        for i in ['ref', 'reference frequency', \
                                  'consensus', 'consensus frequency', \
                                  'A', 'T', 'G', 'C', 'N', \
                                  'In', 'Del']:
                            try:
                                if i == 'consensus':
                                    out.append(ps[i][0])
                                elif i in ['In', 'Del']:
                                    out.append(len(ps[i]))
                                else:
                                    out.append(ps[i])
                            except:
                                out.append('n/a')
                        print('\t'.join([str(i) for i in out]), file=bp_table)
    return genomes
514412
def extract_full(rec, sites, flank, fw):
    
    for s in sites:
        newid = "{0}:{1}".format(rec.name, s)
        left = max(s - flank, 0)
        right = min(s + flank, len(rec))
        frag = rec.seq[left:right].strip("Nn")
        newrec = SeqRecord(frag, id=newid, description="")
        SeqIO.write([newrec], fw, "fasta")
270084
def variants(ctx, variant_id, chromosome, end_chromosome, start, end, variant_type,
             sv_type):
    """"""
    if sv_type:
        variant_type = 'sv'

    adapter = ctx.obj['adapter']

    if (start or end):
        if not (chromosome and start and end):
            LOG.warning("Regions must be specified with chromosome, start and end")
            return

    if variant_id:
        variant = adapter.get_variant({'_id':variant_id})
        if variant:
            click.echo(variant)
        else:
            LOG.info("Variant {0} does not exist in database".format(variant_id))
        return

    if variant_type == 'snv':
        result = adapter.get_variants(
            chromosome=chromosome,
            start=start,
            end=end
        )
    else:
        LOG.info("Search for svs")
        result = adapter.get_sv_variants(
            chromosome=chromosome,
            end_chromosome=end_chromosome,
            sv_type=sv_type,
            pos=start,
            end=end
        )

    i = 0
    for variant in result:
        i += 1
        pp(variant)

    LOG.info("Number of variants found in database: %s", i)
979496
def _parse_coords(coord_lines):
        
        paras = {}
        var_pattern = re.compile(r"^([A-Za-z]+\S*)[\s=,]+([\d\-\.]+)$")
        for l in coord_lines:
            m = var_pattern.match(l.strip())
            if m:
                paras[m.group(1).strip("=")] = float(m.group(2))

        species = []
        coords = []
        # Stores whether a Zmatrix format is detected. Once a zmatrix format
        # is detected, it is assumed for the remaining of the parsing.
        zmode = False
        for l in coord_lines:
            l = l.strip()
            if not l:
                break
            if (not zmode) and GaussianInput._xyz_patt.match(l):
                m = GaussianInput._xyz_patt.match(l)
                species.append(m.group(1))
                toks = re.split(r"[,\s]+", l.strip())
                if len(toks) > 4:
                    coords.append([float(i) for i in toks[2:5]])
                else:
                    coords.append([float(i) for i in toks[1:4]])
            elif GaussianInput._zmat_patt.match(l):
                zmode = True
                toks = re.split(r"[,\s]+", l.strip())
                species.append(toks[0])
                toks.pop(0)
                if len(toks) == 0:
                    coords.append(np.array([0, 0, 0]))
                else:
                    nn = []
                    parameters = []
                    while len(toks) > 1:
                        ind = toks.pop(0)
                        data = toks.pop(0)
                        try:
                            nn.append(int(ind))
                        except ValueError:
                            nn.append(species.index(ind) + 1)
                        try:
                            val = float(data)
                            parameters.append(val)
                        except ValueError:
                            if data.startswith("-"):
                                parameters.append(-paras[data[1:]])
                            else:
                                parameters.append(paras[data])
                    if len(nn) == 1:
                        coords.append(np.array([0, 0, parameters[0]]))
                    elif len(nn) == 2:
                        coords1 = coords[nn[0] - 1]
                        coords2 = coords[nn[1] - 1]
                        bl = parameters[0]
                        angle = parameters[1]
                        axis = [0, 1, 0]
                        op = SymmOp.from_origin_axis_angle(coords1, axis,
                                                           angle, False)
                        coord = op.operate(coords2)
                        vec = coord - coords1
                        coord = vec * bl / np.linalg.norm(vec) + coords1
                        coords.append(coord)
                    elif len(nn) == 3:
                        coords1 = coords[nn[0] - 1]
                        coords2 = coords[nn[1] - 1]
                        coords3 = coords[nn[2] - 1]
                        bl = parameters[0]
                        angle = parameters[1]
                        dih = parameters[2]
                        v1 = coords3 - coords2
                        v2 = coords1 - coords2
                        axis = np.cross(v1, v2)
                        op = SymmOp.from_origin_axis_angle(
                            coords1, axis, angle, False)
                        coord = op.operate(coords2)
                        v1 = coord - coords1
                        v2 = coords1 - coords2
                        v3 = np.cross(v1, v2)
                        adj = get_angle(v3, axis)
                        axis = coords1 - coords2
                        op = SymmOp.from_origin_axis_angle(
                            coords1, axis, dih - adj, False)
                        coord = op.operate(coord)
                        vec = coord - coords1
                        coord = vec * bl / np.linalg.norm(vec) + coords1
                        coords.append(coord)

        def _parse_species(sp_str):
            """
            The species specification can take many forms. E.g.,
            simple integers representing atomic numbers ("8"),
            actual species string ("C") or a labelled species ("C1").
            Sometimes, the species string is also not properly capitalized,
            e.g, ("c1"). This method should take care of these known formats.
            """
            try:
                return int(sp_str)
            except ValueError:
                sp = re.sub(r"\d", "", sp_str)
                return sp.capitalize()

        species = [_parse_species(sp) for sp in species]

        return Molecule(species, coords)
34723
def identify_assembly_points(graph, bgtree, target_multicolor, exclude=None, verbose=False, verbose_destination=None):
    
    if verbose:
        print(">>Identifying assemblies for target multicolor:",
              [e.name for e in target_multicolor.multicolors.elements()], file=verbose_destination)

    guidance = bgtree.consistent_multicolors[:]
    offset = len(Multicolor.split_colors(target_multicolor, guidance=guidance,
                                         account_for_color_multiplicity_in_guidance=False)) - 1

    threshold = 1 if offset == 0 else 2
    assemblies = []  # the overall result
    if exclude is None:
        exclude = []  # a container with single colors of genomes, that are to be considered fully assembled

    p_t_consistent_multicolors_in_target = Multicolor.split_colors(target_multicolor, guidance=guidance,
                                                                   account_for_color_multiplicity_in_guidance=False)
    t_consistent_multicolors_in_target = []
    for tcmc in p_t_consistent_multicolors_in_target:
        t_consistent_multicolors_in_target.append(sorted(color.name for color in tcmc.colors))
    # we work with each connected component separately, as connected components usually preserve fragmentation points
    ################################################################################################
    #
    # its important that we iterate over connected components making each particular one a deepcopy of an
    #   underlying breakpoint graph connected component
    #
    ################################################################################################
    for i, cc in enumerate(graph.connected_components_subgraphs(copy=True)):

        # we filter current connected component of uninteresting / ambiguous edges and retrieve a list of
        #   connected components that are left in the original connected components after filtration
        irregular_subnets = get_irregular_subnets(cc, target_multicolor, exclude)

        if len(irregular_subnets) > 0 and verbose:
            print(">>Processing", str(i) + "th", "connected component", file=verbose_destination)
            print("\tcontains", len(irregular_subnets), "subnet groups", file=verbose_destination)

        # each subnet can be processed separately
        for subnet in irregular_subnets:

            supporting_edge_scores = get_support_edge_scores(graph, subnet, target_multicolor, bgtree)

            # we create a new dummy graph for the purpose of computing maximum weight matching for support edges in it
            new_graph = nx.Graph()
            if verbose:
                print("\tcontains", len(supporting_edge_scores), "possible assembly points", file=verbose_destination)

            # we'll keep track of possible assembly points for future reference
            support_edge_dict = {}
            for (v1, v2), before, after, ex_data in supporting_edge_scores:
                ex_data["tcmc"] = t_consistent_multicolors_in_target
                ##########################################################################################
                #
                # INSERT YOUR CODE ASSEMBLY SCORE THRESHOLD FILTRATION HERE IF NEED BE
                #
                ##########################################################################################
                if before - after - offset < threshold:
                    continue
                ##########################################################################################
                #
                # by default networkx assumes all edges, that have weight >= 0 are good
                #
                ##########################################################################################
                new_graph.add_edge(v1, v2, weight=before - after - offset)
                support_edge_dict[(v1, v2)] = (before, after + offset, ex_data)
                support_edge_dict[(v2, v1)] = (before, after + offset, ex_data)

            maximal_matching = nx.max_weight_matching(new_graph)

            if verbose:
                print("\t", len(maximal_matching) // 2, "assembly points are identified", file=verbose_destination)

            # as networkx provides a maximum matching in a form of adjacency list, every identified edge
            #   (pair of vertices) is present their twice (i.e. matching[v1]=v2 and matching[v2]=v1)
            #   we need to make sure we only add every edge only once
            visited = set()
            for v1, v2 in maximal_matching.items():
                if v1 in visited or v2 in visited:
                    continue
                visited.add(v1)
                visited.add(v2)
                assemblies.append((v1, v2, support_edge_dict[(v1, v2)]))

    # we return the result as a list of assembly points that were identified for the targeted multicolor
    # as a list of tuples (v1, v2, ("before", "after"))
    # where v1 and v2 correspond to assembly point and "before" and "after" are used to compute the assembly score
    return assemblies
905394
def _get_callsites(self, function_address):
        

        all_predecessors = []

        nodes = self.get_all_nodes(function_address)
        for n in nodes:
            predecessors = list(self.get_predecessors(n))
            all_predecessors.extend(predecessors)

        return all_predecessors
1029772
def sequencetyper(self):
        
        for sample in self.metadata.samples:
            if sample.general.bestassemblyfile != 'NA':
                if type(sample[self.analysistype].allelenames) == list:
                    #
                    if sample[self.analysistype].profile != 'NA':
                        # Initialise dictionaries
                        sample[self.analysistype].profilematches = dict()
                        sample[self.analysistype].sequencetypematches = dict()
                        # Create the profiledata variable to avoid writing self.profiledata[self.analysistype]
                        profiledata = sample[self.analysistype].profiledata
                        # For each gene
                        for gene in sorted(sample[self.analysistype].allelenames):
                            try:
                                allelenumber = sample[self.analysistype].allelematches[gene].split('-')[1]
                                # Find the profile with the most alleles in common with the query genome
                                for sequencetype in profiledata:
                                    # refallele is the allele number of the sequence type
                                    refallele = profiledata[sequencetype][gene]
                                    if allelenumber == refallele:
                                        # Add matching alleles
                                        try:
                                            sample[self.analysistype].profilematches[sequencetype] += 1
                                            sample[self.analysistype].sequencetypematches[sequencetype].append(
                                                refallele)
                                        except KeyError:
                                            sample[self.analysistype].profilematches[sequencetype] = 1
                                            sample[self.analysistype].sequencetypematches[sequencetype] = list()
                                            sample[self.analysistype].sequencetypematches[sequencetype].append(
                                                refallele)
                            except KeyError:
                                pass
1302920
def __getSequenceVariants(self, x1, polyStart, polyStop, listSequence) :
		""""""
		if polyStart < len(self.polymorphisms) and polyStart < polyStop: 
			sequence = copy.copy(listSequence)
			ret = []
			
			pk = self.polymorphisms[polyStart]
			posInSequence = pk[0]-x1
			if posInSequence < len(listSequence) : 
				for allele in pk[1] :
					sequence[posInSequence] = allele
					ret.extend(self.__getSequenceVariants(x1, polyStart+1, polyStop, sequence))
			
			return ret
		else :
			return [''.join(listSequence)]
656584
def parse_tstv_by_count(self):
        

        self.vcftools_tstv_by_count = dict()
        for f in self.find_log_files('vcftools/tstv_by_count', filehandles=True):
            d = {}
            for line in f['f'].readlines()[1:]: # don't add the header line (first row)
                key = float(line.split()[0]) # taking the first column (alternative allele count) as key
                val = float(line.split()[3]) # taking Ts/Tv as value
                d[key] = val
            self.vcftools_tstv_by_count[f['s_name']] = d

        # Filter out ignored sample names
        self.vcftools_tstv_by_count = self.ignore_samples(self.vcftools_tstv_by_count)

        if len(self.vcftools_tstv_by_count) == 0:
            return 0

        pconfig = {
            'id': 'vcftools_tstv_by_count',
            'title': 'VCFTools: TsTv by Count',
            'ylab': 'TsTv Ratio',
            'xlab': 'Alternative Allele Count',
            'xmin': 0,
            'ymin': 0,
            'smooth_points': 400, # this limits huge filesizes and prevents browser crashing
            'smooth_points_sumcounts': False
        }

        helptext = '''
        `Transition` is a purine-to-purine or pyrimidine-to-pyrimidine point mutations.
        `Transversion` is a purine-to-pyrimidine or pyrimidine-to-purine point mutation.
        `Alternative allele count` is the number of alternative alleles at the site.
        Note: only bi-allelic SNPs are used (multi-allelic sites and INDELs are skipped.)
        Refer to Vcftools's manual (https://vcftools.github.io/man_latest.html) on `--TsTv-by-count`
        '''

        self.add_section(
            name = 'TsTv by Count',
            anchor = 'vcftools-tstv-by-count',
            description = "Plot of `TSTV-BY-COUNT` - the transition to transversion ratio as a function of alternative allele count from the output of vcftools TsTv-by-count.",
            helptext = helptext,
            plot = linegraph.plot(self.vcftools_tstv_by_count,pconfig)
        )

        return len(self.vcftools_tstv_by_count)
371979
def inferred_var_seqs_plus_flanks(self, ref_seq, flank_length):
        ''''''
        flank_start = max(0, self.POS - flank_length)
        flank_end = min(len(ref_seq) - 1, self.ref_end_pos() + flank_length)
        seqs = [ref_seq[flank_start:self.POS] + self.REF + ref_seq[self.ref_end_pos() + 1: flank_end + 1]]

        for alt in self.ALT:
            seqs.append(ref_seq[flank_start:self.POS] + alt + ref_seq[self.ref_end_pos() + 1: flank_end + 1])

        return flank_start, seqs
61602
def sequence_dir(self):
        """"""
        if self.root_dir:
            return op.join(self.protein_dir, 'sequences')
        else:
            log.debug('Root directory not set')
            return None
210950
def extract_finditer(pos_seq, regex=SimpleNP):
	""""""
	ss = coarse_tag_str(pos_seq)

	def gen():
		for m in re.finditer(regex, ss):
			yield (m.start(), m.end())

	return list(gen())
703495
def cli(ctx, feature_id, organism="", sequence=""):
    
    return ctx.gi.annotations.get_feature_sequence(feature_id, organism=organism, sequence=sequence)
1282511
def _run_collect_allelic_counts(pos_file, pos_name, work_dir, data):
    
    out_dir = utils.safe_makedir(os.path.join(dd.get_work_dir(data), "structural", "counts"))
    out_file = os.path.join(out_dir, "%s-%s-counts.tsv" % (dd.get_sample_name(data), pos_name))
    if not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            params = ["-T", "CollectAllelicCounts", "-L", pos_file, "-I", dd.get_align_bam(data),
                      "-R", dd.get_ref_file(data), "-O", tx_out_file]
            _run_with_memory_scaling(params, tx_out_file, data)
    return out_file
636978
def _cluster(bam_file, ma_file, out_dir, reference, annotation_file=None):
    
    seqcluster = op.join(get_bcbio_bin(), "seqcluster")
    # cl = ["cluster", "-o", out_dir, "-m", ma_file, "-a", bam_file, "-r", reference]
    if annotation_file:
        annotation_file = "-g " + annotation_file
    else:
        annotation_file = ""

    if not file_exists(op.join(out_dir, "counts.tsv")):
        cmd = ("{seqcluster} cluster -o {out_dir} -m {ma_file} -a {bam_file} -r {reference} {annotation_file}")
        do.run(cmd.format(**locals()), "Running seqcluster.")
    counts = op.join(out_dir, "counts.tsv")
    stats = op.join(out_dir, "read_stats.tsv")
    json = op.join(out_dir, "seqcluster.json")
    return {'out_dir': out_dir, 'count_file': counts, 'stat_file': stats, 'json': json}
550494
def haplotype_caller(align_bams, items, ref_file, assoc_files,
                       region=None, out_file=None):
    
    if out_file is None:
        out_file = "%s-variants.vcf.gz" % utils.splitext_plus(align_bams[0])[0]
    if not utils.file_exists(out_file):
        num_cores = dd.get_num_cores(items[0])
        broad_runner, params = \
            _shared_gatk_call_prep(align_bams, items, ref_file, region, out_file, num_cores)
        gatk_type = broad_runner.gatk_type()
        assert gatk_type in ["restricted", "gatk4"], \
            "Require full version of GATK 2.4+, or GATK4 for haplotype calling"
        with file_transaction(items[0], out_file) as tx_out_file:
            resources = config_utils.get_resources("gatk-spark", items[0]["config"])
            spark_opts = [str(x) for x in resources.get("options", [])]
            if _use_spark(num_cores, gatk_type, items, spark_opts):
                params += ["-T", "HaplotypeCallerSpark"]
                if spark_opts:
                    params += spark_opts
                else:
                    params += ["--spark-master", "local[%s]" % num_cores,
                               "--conf", "spark.local.dir=%s" % os.path.dirname(tx_out_file),
                               "--conf", "spark.driver.host=localhost", "--conf", "spark.network.timeout=800",
                               "--conf", "spark.executor.heartbeatInterval=100"]
            else:
                params += ["-T", "HaplotypeCaller"]
            params += ["--annotation", "ClippingRankSumTest",
                       "--annotation", "DepthPerSampleHC"]
            # Enable hardware based optimizations in GATK 3.1+
            if LooseVersion(broad_runner.gatk_major_version()) >= LooseVersion("3.1"):
                if _supports_avx():
                    # Scale down HMM thread default to avoid overuse of cores
                    # https://github.com/bcbio/bcbio-nextgen/issues/2442
                    if gatk_type == "gatk4":
                        params += ["--native-pair-hmm-threads", "1"]
                    # GATK4 selects the right HMM optimization automatically with FASTEST_AVAILABLE
                    # GATK3 needs to be explicitly set
                    else:
                        params += ["--pair_hmm_implementation", "VECTOR_LOGLESS_CACHING"]
            resources = config_utils.get_resources("gatk-haplotype", items[0]["config"])
            if "options" in resources:
                params += [str(x) for x in resources.get("options", [])]
            # Prepare gVCFs if doing joint calling
            is_joint = False
            if _joint_calling(items) or any("gvcf" in dd.get_tools_on(d) for d in items):
                is_joint = True
                # If joint calling parameters not set in user options
                if not any([x in ["--emit-ref-confidence", "-ERC", "--emitRefConfidence"] for x in params]):
                    if gatk_type == "gatk4":
                        params += ["--emit-ref-confidence", "GVCF"]
                    else:
                        params += ["--emitRefConfidence", "GVCF"]
                        params += ["--variant_index_type", "LINEAR", "--variant_index_parameter", "128000"]
                # Set GQ banding to not be single GQ resolution
                # No recommended default but try to balance resolution and size
                # http://gatkforums.broadinstitute.org/gatk/discussion/7051/recommendation-best-practices-gvcf-gq-bands

                if not any([x in ["-GQB"] for x in params]):
                    for boundary in [10, 20, 30, 40, 60, 80]:
                        params += ["-GQB", str(boundary)]
            # Enable non-diploid calling in GATK 3.3+
            if LooseVersion(broad_runner.gatk_major_version()) >= LooseVersion("3.3"):
                params += ["-ploidy", str(ploidy.get_ploidy(items, region))]
            if gatk_type == "gatk4":
                # GATK4 Spark calling does not support bgzipped output, use plain VCFs
                if is_joint and _use_spark(num_cores, gatk_type, items, spark_opts):
                    tx_out_file = tx_out_file.replace(".vcf.gz", ".vcf")
                params += ["--output", tx_out_file]
            else:
                params += ["-o", tx_out_file]
            broad_runner.new_resources("gatk-haplotype")
            memscale = {"magnitude": 0.9 * num_cores, "direction": "increase"} if num_cores > 1 else None
            try:
                broad_runner.run_gatk(params, os.path.dirname(tx_out_file), memscale=memscale,
                                      parallel_gc=_use_spark(num_cores, gatk_type, items, spark_opts))
            except subprocess.CalledProcessError as msg:
                # Spark failing on regions without any reads, write an empty VCF instead
                # https://github.com/broadinstitute/gatk/issues/4234
                if (_use_spark(num_cores, gatk_type, items, spark_opts) and
                      str(msg).find("java.lang.UnsupportedOperationException: empty collection") >= 0 and
                      str(msg).find("at org.apache.spark.rdd.RDD") >= 0):
                    vcfutils.write_empty_vcf(tx_out_file, samples=[dd.get_sample_name(d) for d in items])
                else:
                    raise
            if tx_out_file.endswith(".vcf"):
                vcfutils.bgzip_and_index(tx_out_file, items[0]["config"])


    # avoid bug in GATK where files can get output as non-compressed
    if out_file.endswith(".gz") and not os.path.exists(out_file + ".tbi"):
        with open(out_file, "r") as in_handle:
            is_plain_text = in_handle.readline().startswith("##fileformat")
        if is_plain_text:
            text_out_file = out_file
            out_file = out_file.replace(".vcf.gz", ".vcf")
            shutil.move(text_out_file, out_file)
    return vcfutils.bgzip_and_index(out_file, items[0]["config"])