Reconciles the given gene tree with the species tree and
    report the number of needed duplications and losses
777964
def calculateDupsAndLossesByReconcilingTrees(speciesTree, geneTree, processID):
      
    checkGeneTreeMatchesSpeciesTree(speciesTree, geneTree, processID)   
    def fn(tree, m):  
        if tree.internal:
            nodes = fn(tree.left, m)
            nodes = nodes.union(fn(tree.right, m))
            m[tree.traversalID.mid] = nodes
        else:
            m[tree.traversalID.mid] = set((processID(tree.iD),))
        return m[tree.traversalID.mid]
    a = {}
    fn(speciesTree, a)
    b = {}
    fn(geneTree, b)
    def fn2(nodes, speciesTree):
        assert nodes.issubset(a[speciesTree.traversalID.mid])
        if speciesTree.internal:
            if nodes.issubset(a[speciesTree.left.traversalID.mid]):
                return fn2(nodes, speciesTree.left)
            if nodes.issubset(a[speciesTree.right.traversalID.mid]):
                return fn2(nodes, speciesTree.right)
        return speciesTree.traversalID.mid
    for iD in b.keys():
        nodes = b[iD]
        b[iD] = fn2(nodes, speciesTree)
    dups = []
    def fn3(geneTree):
        if geneTree.internal:
            i = b[geneTree.traversalID.mid]
            if b[geneTree.left.traversalID.mid] == i or b[geneTree.right.traversalID.mid] == i:
                dups.append(geneTree.traversalID.mid)
            fn3(geneTree.left)
            fn3(geneTree.right)
    fn3(geneTree)
    lossMap = {}
    def fn4(speciesTree):
        nodes = [(speciesTree.traversalID.mid, -1)]
        lossMap[(speciesTree.traversalID.mid, speciesTree.traversalID.mid)] = 0
        if speciesTree.internal:
            for node, losses in fn4(speciesTree.left) + fn4(speciesTree.right):
                lossMap[(speciesTree.traversalID.mid, node)] = losses+1
                nodes.append((node, losses+1))
        return nodes
    for node, losses in fn4(speciesTree):
        lossMap[(sys.maxint, node)] = losses+1
    losses = [0]
    def fn5(geneTree, ancestor):
        if geneTree.internal:
            i = b[geneTree.traversalID.mid]
            if geneTree.traversalID.mid in dups:
                losses[0] += lossMap[(ancestor, b[geneTree.left.traversalID.mid])]
                losses[0] += lossMap[(ancestor, b[geneTree.right.traversalID.mid])]
            else:
                losses[0] += lossMap[(i, b[geneTree.left.traversalID.mid])]
                losses[0] += lossMap[(i, b[geneTree.right.traversalID.mid])]
            fn5(geneTree.left, i)
            fn5(geneTree.right, i)
    ancestorHolder = [None]
    def fn6(speciesTree, ancestor, node):
        if speciesTree.traversalID.mid == node:
            ancestorHolder[0] = ancestor
        if speciesTree.internal:
            fn6(speciesTree.left, speciesTree.traversalID.mid, node)
            fn6(speciesTree.right, speciesTree.traversalID.mid, node)
    ancestor = fn6(speciesTree, sys.maxint, b[geneTree.traversalID.mid])
    assert ancestorHolder[0] is not None
    fn5(geneTree, ancestorHolder[0])
    return len(dups), losses[0]
908881
def merge(args):
    
    from jcvi.formats.base import DictFile

    p = OptionParser(merge.__doc__)
    opts, args = p.parse_args(args)

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

    quartets, registry, lost = args
    qq = DictFile(registry, keypos=1, valuepos=3)
    lost = DictFile(lost, keypos=1, valuepos=0, delimiter='|')
    qq.update(lost)
    fp = open(quartets)
    cases = {
        "AN,CN": 4,
        "BO,AN,CN": 8,
        "BO,CN": 2,
        "BR,AN": 1,
        "BR,AN,CN": 6,
        "BR,BO": 3,
        "BR,BO,AN": 5,
        "BR,BO,AN,CN": 9,
        "BR,BO,CN": 7,
    }
    ip = {
        "syntenic_model": "Syntenic_model_excluded_by_OMG",
        "complete": "Predictable",
        "partial": "Truncated",
        "pseudogene": "Pseudogene",
        "random": "Match_random",
        "real_ns": "Transposed",
        "gmap_fail": "GMAP_fail",
        "AN LOST": "AN_LOST",
        "CN LOST": "CN_LOST",
        "BR LOST": "BR_LOST",
        "BO LOST": "BO_LOST",
        "outside": "Outside_synteny_blocks",
        "[NF]": "Not_found",
    }
    for row in fp:
        atoms = row.strip().split("\t")
        genes = atoms[:4]
        tag = atoms[4]
        a, b, c, d = [qq.get(x, ".").rsplit("-", 1)[-1] for x in genes]
        qqs = [c, d, a, b]
        for i, q in enumerate(qqs):
            if atoms[i] != '.':
                qqs[i] = "syntenic_model"
        # Make comment
        comment = "Case{0}".format(cases[tag])
        dots = sum([1 for x in genes if x == '.'])
        if dots == 1:
            idx = genes.index(".")
            status = qqs[idx]
            status = ip[status]
            comment += "-" + status
        print(row.strip() + "\t" + "\t".join(qqs + [comment]))
1019750
def apply(self, read_tuple_name, read_tuple_id=None, synchronize_widths=True):
        
        parts = read_tuple_name.split("__")
        parts[0] = self._fill_right(parts[0], "-", self.prefix_width)
        if read_tuple_id is not None:
            parts[1] = "{:x}".format(read_tuple_id)
        parts[1] = self._fill_left(parts[1], "0", self.read_tuple_id_width)

        if synchronize_widths:
            new_segments = []
            segments = parts[2][1:-1].split("),(")
            for segment in segments:
                values = segment.split(",")
                values[0] = values[0].zfill(self.genome_id_width)
                values[1] = values[1].zfill(self.chr_id_width)
                values[3] = values[3].zfill(self.coor_width)
                values[4] = values[4].zfill(self.coor_width)
                new_segments.append("(" + ",".join(values) + ")")
            parts[2] = ",".join(new_segments)

        return "__".join(parts)
1149759
def _iterative_bridged_contig_pair_merge(self, outprefix):
        ''''''
        if self.reads is None:
            if self.verbose:
                print('Skipping iterative contig merging because no reads given (see --reads option)')
            return self.original_fasta, None, None

        log_file = outprefix + '.iterations.log'
        log_fh = pyfastaq.utils.open_file_write(log_file)
        genome_fasta = self.original_fasta
        nucmer_coords = outprefix + '.iter.1.coords'
        reads_to_map = self.reads
        act_prefix = None
        made_a_join = True
        iteration = 1

        while made_a_join:
            this_log_prefix = '[' + self.log_prefix + ' iterative_merge ' + str(iteration) + ']'
            print(this_log_prefix, '\tUsing nucmer matches from ', nucmer_coords, sep='', file=log_fh)
            self._run_nucmer(genome_fasta, self.reassembly.contigs_fasta, nucmer_coords)
            act_prefix = outprefix + '.iter.' + str(iteration)
            print(this_log_prefix, '\tYou can view the nucmer matches with ACT using: ./', act_prefix, '.start_act.sh', sep='', file=log_fh)
            self._write_act_files(genome_fasta, self.reassembly.contigs_fasta, nucmer_coords, act_prefix)
            nucmer_hits_by_ref = self._load_nucmer_hits(nucmer_coords)
            made_a_join = self._merge_all_bridged_contigs(nucmer_hits_by_ref, self.original_contigs, self.reassembly_contigs, log_fh, this_log_prefix)
            iteration += 1

            if made_a_join:
                print(this_log_prefix, '\tMade at least one merge. Remapping reads and reassembling',sep='', file=log_fh)
                nucmer_coords = outprefix + '.iter.' + str(iteration) + '.coords'
                genome_fasta = outprefix + '.iter.' + str(iteration) + '.merged.fasta'
                self._contigs_dict_to_file(self.original_contigs, genome_fasta)
                bam = outprefix + '.iter.' + str(iteration) + '.bam'

                circlator.mapping.bwa_mem(
                  genome_fasta,
                  reads_to_map,
                  bam,
                  threads=self.threads,
                  verbose=self.verbose,
                )

                reads_prefix = outprefix + '.iter.' + str(iteration) + '.reads'
                reads_to_map =  reads_prefix + ('.fasta' if self.spades_only_assembler else '.fastq')
                bam_filter = circlator.bamfilter.BamFilter(bam, reads_prefix, fastq_out=not self.spades_only_assembler, split_all_reads=self.split_all_reads)
                bam_filter.run()
                assembler_dir = outprefix + '.iter.' + str(iteration) + '.assembly'
                a = circlator.assemble.Assembler(
                    reads_prefix + ('.fasta' if self.spades_only_assembler else '.fastq'),
                    assembler_dir,
                    threads=self.threads,
                    careful=self.spades_careful,
                    only_assembler=self.spades_only_assembler,
                    verbose=self.verbose,
                    spades_kmers=self.spades_kmers,
                    spades_use_first_success=self.spades_use_first_success,
                    assembler=self.assembler,
                    genomeSize=self.length_cutoff,
                    data_type=self.data_type
                )
                a.run()
                self.reassembly = circlator.assembly.Assembly(assembler_dir, assembler=self.assembler)
                self.reassembly_contigs = self.reassembly.get_contigs()
            elif iteration <= 2:
                print(this_log_prefix, '\tNo contig merges were made',sep='', file=log_fh)

        pyfastaq.utils.close(log_fh)
        return genome_fasta, nucmer_coords, act_prefix + '.start_act.sh'
37998
def merge(cls, trees):
        
        first = trees[0]
        for tree in trees:
            first.update(tree)
        return first
193058
def dump_tree(node, name=None, initial_indent='', indentation='    ',
              maxline=120, maxmerged=80,
              # Runtime optimization
              iter_node=iter_node, special=ast.AST,
              list=list, isinstance=isinstance, type=type, len=len):
    
    def dump(node, name=None, indent=''):
        level = indent + indentation
        name = name and name + '=' or ''
        values = list(iter_node(node))
        if isinstance(node, list):
            prefix, suffix = '%s[' % name, ']'
        elif values:
            prefix, suffix = '%s%s(' % (name, type(node).__name__), ')'
        elif isinstance(node, special):
            prefix, suffix = name + type(node).__name__, ''
        else:
            return '%s%s' % (name, repr(node))
        node = [dump(a, b, level) for a, b in values if b != 'ctx']
        oneline = '%s%s%s' % (prefix, ', '.join(node), suffix)
        if len(oneline) + len(indent) < maxline:
            return '%s' % oneline
        if node and len(prefix) + len(node[0]) < maxmerged:
            prefix = '%s%s,' % (prefix, node.pop(0))
        node = (',\n%s' % level).join(node).lstrip()
        return '%s\n%s%s%s' % (prefix, level, node, suffix)
    return dump(node, name, initial_indent)
481935
def _importSNPs_dbSNPSNP(setName, species, genomeSource, snpsFile) :
	""
	snpData = VCFFile(snpsFile, gziped = True, stream = True)
	dbSNPSNP.dropIndex(('start', 'chromosomeNumber', 'setName'))
	conf.db.beginTransaction()
	pBar = ProgressBar()
	pLabel = ''
	for snpEntry in snpData :
		pBar.update(label = 'Chr %s, %s...' %  (snpEntry['#CHROM'], snpEntry['ID']))
		
		snp = dbSNPSNP()
		for f in snp.getFields() :
			try :
				setattr(snp, f, snpEntry[f])
			except KeyError :
				pass
		snp.chromosomeNumber = snpEntry['#CHROM']
		snp.species = species
		snp.setName = setName
		snp.start = snpEntry['POS']-1
		snp.alt = snpEntry['ALT']
		snp.ref = snpEntry['REF']
		snp.end = snp.start+len(snp.alt)
		snp.save()
	
	pBar.close()
	
	snpMaster = SNPMaster()
	snpMaster.set(setName = setName, SNPType = 'dbSNPSNP', species = species)
	snpMaster.save()
	
	printf('saving...')
	conf.db.endTransaction()
	printf('creating indexes...')
	dbSNPSNP.ensureGlobalIndex(('start', 'chromosomeNumber', 'setName'))
	printf('importation of SNP set %s for species %s done.' %(setName, species))

	return True
266042
def configure_crossover(self, genome1, genome2, config):
        
        assert isinstance(genome1.fitness, (int, float))
        assert isinstance(genome2.fitness, (int, float))
        if genome1.fitness > genome2.fitness:
            parent1, parent2 = genome1, genome2
        else:
            parent1, parent2 = genome2, genome1

        # Inherit connection genes
        for key, cg1 in iteritems(parent1.connections):
            cg2 = parent2.connections.get(key)
            if cg2 is None:
                # Excess or disjoint gene: copy from the fittest parent.
                self.connections[key] = cg1.copy()
            else:
                # Homologous gene: combine genes from both parents.
                self.connections[key] = cg1.crossover(cg2)

        # Inherit node genes
        parent1_set = parent1.nodes
        parent2_set = parent2.nodes

        for key, ng1 in iteritems(parent1_set):
            ng2 = parent2_set.get(key)
            assert key not in self.nodes
            if ng2 is None:
                # Extra gene: copy from the fittest parent
                self.nodes[key] = ng1.copy()
            else:
                # Homologous gene: combine genes from both parents.
                self.nodes[key] = ng1.crossover(ng2)
923295
def fix_vcf_line(parts, ref_base):
    
    swap = {"1/1": "0/0", "0/1": "0/1", "0/0": "1/1", "./.": "./."}
    complements = {"G": "C", "A": "T", "C": "G", "T": "A", "N": "N"}
    varinfo, genotypes = fix_line_problems(parts)
    ref, var = varinfo[3:5]
    # non-reference regions or non-informative, can't do anything
    if ref_base in [None, "N"] or set(genotypes) == set(["./."]):
        varinfo = None
    # matching reference, all good
    elif ref_base == ref:
        assert ref_base == ref, (ref_base, parts)
    # swapped reference and alternate regions
    elif ref_base == var or ref in ["N", "0"]:
        varinfo[3] = var
        varinfo[4] = ref
        genotypes = [swap[x] for x in genotypes]
    # reference is on alternate strand
    elif ref_base != ref and complements.get(ref) == ref_base:
        varinfo[3] = complements[ref]
        varinfo[4] = ",".join([complements[v] for v in var.split(",")])
    # unspecified alternative base
    elif ref_base != ref and var in ["N", "0"]:
        varinfo[3] = ref_base
        varinfo[4] = ref
        genotypes = [swap[x] for x in genotypes]
    # swapped and on alternate strand
    elif ref_base != ref and complements.get(var) == ref_base:
        varinfo[3] = complements[var]
        varinfo[4] = ",".join([complements[v] for v in ref.split(",")])
        genotypes = [swap[x] for x in genotypes]
    else:
        print "Did not associate ref {0} with line: {1}".format(
            ref_base, varinfo)
    if varinfo is not None:
        return varinfo + genotypes
357778
def parse(self,tolerance=0,downsample=None,evidence=2,use_gene_names=False):
      """"""
      g = Graph()
      nodes = [Node(x) for x in self._transcripts]
      for n in nodes: g.add_node(n)
      for i in range(0,len(nodes)):
         for j in range(0,len(nodes)):
            if i == j: continue
            jov = nodes[i].payload.junction_overlap(nodes[j].payload,tolerance)
            sub = jov.is_subset()
            if not sub: continue
            if sub == 1:
               g.add_edge(Edge(nodes[i],nodes[j]))
               g.add_edge(Edge(nodes[j],nodes[i]))
            if sub == 2:
               g.add_edge(Edge(nodes[i],nodes[j]))
      g.merge_cycles()
      roots = g.roots
      groups = []
      for r in roots:
         g2 = g.get_root_graph(r)
         c = CompatibleGraph(g2,tolerance,downsample,evidence,use_gene_names=use_gene_names)
         groups.append(c)
      return groups
660866
public void genStats(List<? extends JCTree> trees, Env<GenContext> env) {
        for (List<? extends JCTree> l = trees; l.nonEmpty(); l = l.tail)
            genStat(l.head, env, CRT_STATEMENT);
    }
669003
def print_genome_matrix(hits, fastas, id2desc, file_name):
    
    out = open(file_name, 'w')
    fastas = sorted(fastas)
    print('## percent identity between genomes', file=out)
    print('# - \t %s' % ('\t'.join(fastas)), file=out)
    for fasta in fastas:
        line = [fasta]
        for other in fastas:
            if other == fasta:
                average = '-'
            else:
                average = numpy.average([hits[fasta][other][i][3] for i in hits[fasta][other]])
            line.append(str(average))
        print('\t'.join(line), file=out)
    print('', file=out)
    print('## percent of orfs that are orthologous between genomes', file=out)
    print('# - \t %s' % ('\t'.join(fastas)), file=out)
    for fasta in fastas:
        line = [fasta]
        for other in fastas:
            if other == fasta:
                percent = '-'
            else:
                orthologs = float(len(hits[fasta][other]))
                orfs = float(len([i for i in id2desc if id2desc[i][0] == fasta]))
                percent = float(orthologs / orfs) * 100
            line.append(str(percent))
        print('\t'.join(line), file=out)
442648
def _importSNPs_AgnosticSNP(setName, species, genomeSource, snpsFile) :
	""
	printf('importing SNP set %s for species %s...' % (setName, species))

	snpData = CSVFile()
	snpData.parse(snpsFile, separator = "\t")

	AgnosticSNP.dropIndex(('start', 'chromosomeNumber', 'setName'))
	conf.db.beginTransaction()
	
	pBar = ProgressBar(len(snpData))
	pLabel = ''
	currChrNumber = None
	for snpEntry in snpData :
		tmpChr = snpEntry['chromosomeNumber']
		if tmpChr != currChrNumber :
			currChrNumber = tmpChr
			pLabel = 'Chr %s...' % currChrNumber

		snp = AgnosticSNP()
		snp.species = species
		snp.setName = setName
		for f in snp.getFields() :
			try :
				setattr(snp, f, snpEntry[f])
			except KeyError :
				if f != 'species' and f != 'setName' :
					printf("Warning filetype as no key %s", f)
		snp.quality = float(snp.quality)
		snp.start = int(snp.start)
		snp.end = int(snp.end)
		snp.save()
		pBar.update(label = pLabel)

	pBar.close()
	
	snpMaster = SNPMaster()
	snpMaster.set(setName = setName, SNPType = 'AgnosticSNP', species = species)
	snpMaster.save()

	printf('saving...')
	conf.db.endTransaction()
	printf('creating indexes...')
	AgnosticSNP.ensureGlobalIndex(('start', 'chromosomeNumber', 'setName'))
	printf('importation of SNP set %s for species %s done.' %(setName, species))
	
	return True
25434
def merge_chromosome_dfs(df_tuple):
    # type: (Tuple[pd.DataFrame, pd.DataFrame]) -> pd.DataFrame
    """"""

    plus_df, minus_df = df_tuple
    index_cols = "Chromosome Bin".split()
    count_column = plus_df.columns[0]

    if plus_df.empty:
        return return_other(minus_df, count_column, index_cols)
    if minus_df.empty:
        return return_other(plus_df, count_column, index_cols)

    # sum duplicate bins
    # TODO: why are there duplicate bins here in the first place?
    plus_df = plus_df.groupby(index_cols).sum()
    minus_df = minus_df.groupby(index_cols).sum()

    # first sum the two bins from each strand
    df = pd.concat([plus_df, minus_df], axis=1).fillna(0).sum(axis=1)
    df = df.reset_index().sort_values(by="Bin")

    df.columns = ["Chromosome", "Bin", count_column]

    df = df.sort_values(["Chromosome", "Bin"])
    df[["Bin", count_column]] = df[["Bin", count_column]].astype(int32)
    df = df[[count_column, "Chromosome", "Bin"]]
    return df.reset_index(drop=True)
947517
def p_term_mult(self, tree):
        ''' '''
        tree.value = convertstr(tree[0].value) * convertstr(tree[2].value)
957299
function(root,new$,old,caret) {
	var k = new$.length;
	var i = k;
	var last = new$[k - 1];
	
	
	if (k == old.length && new$[0] === old[0]) {
		// running through to compare
		while (i--){
			if (new$[i] !== old[i]) { break; };
		};
	};
	
	if (i == -1) {
		return last && last._slot_ || last || caret;
	} else {
		return reconcileCollectionChanges(root,new$,old,caret);
	};
}
27299
def balanced_binary_tree(n_leaves):
    
    def _balanced_subtree(leaves):
        if len(leaves) == 1:
            return leaves[0]
        elif len(leaves) == 2:
            return (leaves[0], leaves[1])
        else:
            split = len(leaves) // 2
            return (_balanced_subtree(leaves[:split]),
                    _balanced_subtree(leaves[split:]))

    return _balanced_subtree(np.arange(n_leaves))