Helper function to make a pyensembl.Genome instance.
541445
def infer_genome_from_vcf(genome, vcf_reader, reference_vcf_key):
    
    if genome:
        return infer_genome(genome)
    elif reference_vcf_key not in vcf_reader.metadata:
        raise ValueError("Unable to infer reference genome for %s" % (
            vcf_reader.filename,))
    else:
        reference_path = vcf_reader.metadata[reference_vcf_key]
        return infer_genome(reference_path)
883344
def build(self):
        """"""
        ang_per_res = (2 * numpy.pi) / self.nucleotides_per_turn
        atom_offset_coords = _backbone_properties[self.helix_type]['atoms']
        if self.handedness == 'l':
            handedness = -1
        else:
            handedness = 1
        base_atom_labels = _backbone_properties[self.helix_type]['labels']
        monomers = []
        mol_code_format = _backbone_properties[self.helix_type]['mol_code_format']
        for i, b in enumerate(self.base_sequence):
            nucleotide = Nucleotide(
                mol_code=mol_code_format.format(b), ampal_parent=self)
            atoms_dict = OrderedDict()
            if (i == (len(self.base_sequence) - 1)) and not self.phos_3_prime:
                # Do not include phosphate on last nucleotide
                atom_labels = base_atom_labels[3:] + [_bases[b]['labels'][2]]
                atom_offsets = {k: v for k, v in zip(
                    atom_labels, atom_offset_coords[3:])}
            else:
                atom_labels = base_atom_labels + [_bases[b]['labels'][2]]
                atom_offsets = {k: v for k, v in zip(
                    atom_labels, atom_offset_coords)}
            for atom_label in atom_labels:
                r, zeta, z_shift = atom_offsets[atom_label]
                rot_ang = ((i * ang_per_res) + zeta) * handedness
                z = (self.rise_per_nucleotide * i) + z_shift
                coords = cylindrical_to_cartesian(
                    radius=r, azimuth=rot_ang, z=z, radians=True)
                atom = Atom(
                    coordinates=coords, element=atom_label[0],
                    ampal_parent=nucleotide, res_label=atom_label)
                atoms_dict[atom_label] = atom
            base_ref = _bases[b]['ref_atom']
            rot_adj = _bases[b]['rot_adj']
            base_dict = OrderedDict(
                zip(_bases[b]['labels'], _bases[b]['atoms']))
            translation, angle, axis, point = find_transformations(
                base_dict[base_ref], base_dict["C1'"],
                atoms_dict[base_ref]._vector, atoms_dict["C1'"]._vector)
            q1 = Quaternion.angle_and_axis(angle, axis)
            # Align N9 C1'
            for k, v in base_dict.items():
                base_dict[k] = q1.rotate_vector(v, point) + translation
            # Rotate to align O4'
            axis = numpy.array(base_dict["C1'"]) - base_dict[base_ref]
            angle = dihedral(base_dict["O4'"], base_dict[base_ref],
                             base_dict["C1'"], atoms_dict["O4'"]) - rot_adj
            q2 = Quaternion.angle_and_axis(angle, axis)
            for k, v in list(base_dict.items()):
                if k not in atoms_dict:
                    atom = Atom(q2.rotate_vector(v, base_dict[base_ref]),
                                element=k[0], ampal_parent=nucleotide,
                                res_label=k)
                    atoms_dict[k] = atom
            nucleotide.atoms = atoms_dict
            monomers.append(nucleotide)
        self._monomers = monomers
        self.relabel_monomers()
        self.relabel_atoms()
        return
429906
def importRemoteSNPs(name) :
	""""""
	try :
		dw = listRemoteDatawraps()["Flat"]["SNPs"]
	except AttributeError :
		raise AttributeError("There's no remote genome datawrap by the name of: '%s'" % name)

	finalFile = _DW(name, dw["url"])
	PS.importSNPs(finalFile)
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
529634
def create(genome, config):
        

        # Gather expressed connections.
        connections = [cg.key for cg in itervalues(genome.connections) if cg.enabled]

        layers = feed_forward_layers(config.genome_config.input_keys, config.genome_config.output_keys, connections)
        node_evals = []
        for layer in layers:
            for node in layer:
                inputs = []
                node_expr = [] # currently unused
                for conn_key in connections:
                    inode, onode = conn_key
                    if onode == node:
                        cg = genome.connections[conn_key]
                        inputs.append((inode, cg.weight))
                        node_expr.append("v[{}] * {:.7e}".format(inode, cg.weight))


                ng = genome.nodes[node]
                aggregation_function = config.genome_config.aggregation_function_defs.get(ng.aggregation)
                activation_function = config.genome_config.activation_defs.get(ng.activation)
                node_evals.append((node, activation_function, aggregation_function, ng.bias, ng.response, inputs))

        return FeedForwardNetwork(config.genome_config.input_keys, config.genome_config.output_keys, node_evals)
667708
def sequencetyper(self):
        """"""
        for sample in self.metadata:
            if sample.general.bestassemblyfile != 'NA':
                if type(sample[self.analysistype].allelenames) == list:
                    # Initialise variables
                    header = 0
                    # Iterate through the genomes
                    # for sample in self.metadata:
                    genome = sample.name
                    # Initialise self.bestmatch[genome] with an int that will eventually be replaced by the # of matches
                    self.bestmatch[genome] = defaultdict(int)
                    if sample[self.analysistype].profile != 'NA':
                        # Create the profiledata variable to avoid writing self.profiledata[self.analysistype]
                        profiledata = sample[self.analysistype].profiledata
                        # For each gene in plusdict[genome]
                        for gene in sample[self.analysistype].allelenames:
                            # Clear the appropriate count and lists
                            multiallele = []
                            multipercent = []
                            # Go through the alleles in plusdict
                            for allele in self.plusdict[genome][gene]:
                                percentid = list(self.plusdict[genome][gene][allele].keys())[0]
                                # "N" alleles screw up the allele splitter function
                                if allele != "N":
                                    # Use the alleleSplitter function to get the allele number
                                    # allelenumber, alleleprenumber = allelesplitter(allele)
                                    # Append as appropriate - alleleNumber is treated as an integer for proper sorting
                                    multiallele.append(int(allele))
                                    multipercent.append(percentid)
                                # If the allele is "N"
                                else:
                                    # Append "N" and a percent identity of 0
                                    multiallele.append("N")
                                    multipercent.append(0)
                                if not multiallele:
                                    multiallele.append("N")
                                    multipercent.append(0)
                            # if self.analysistype == 'rmlst':
                            #     # For whatever reason, the rMLST profile scheme treat multiple allele hits as 'N's.
                            #     multiallele = multiallele if len(multiallele) == 1 else ['N']
                            #     if multipercent:
                            #         multipercent = multipercent if len(multiallele) == 1 else [0, 0]
                            #     else:
                            #         multipercent = [0]
                            # Populate self.bestdict with genome, gene, alleles joined with a space (this was made like
                            # this because allele is a list generated by the .iteritems() above
                            self.bestdict[genome][gene][" ".join(str(allele)
                                                                 for allele in sorted(multiallele))] = multipercent[0]
                            # Find the profile with the most alleles in common with the query genome
                            for sequencetype in profiledata:
                                # The number of genes in the analysis
                                header = len(profiledata[sequencetype])
                                # refallele is the allele number of the sequence type
                                refallele = profiledata[sequencetype][gene]
                                # If there are multiple allele matches for a gene in the reference profile e.g. 10 692
                                if len(refallele.split(" ")) > 1:
                                    # Map the split (on a space) alleles as integers - if they are treated as integers,
                                    # the alleles will sort properly
                                    intrefallele = map(int, refallele.split(" "))
                                    # Create a string of the joined, sorted alleles
                                    sortedrefallele = " ".join(str(allele) for allele in sorted(intrefallele))
                                else:
                                    # Use the reference allele as the sortedRefAllele
                                    sortedrefallele = refallele
                                for allele, percentid in self.bestdict[genome][gene].items():
                                    # If the allele in the query genome matches the allele in the reference profile, add
                                    # the result to the bestmatch dictionary. Genes with multiple alleles were sorted
                                    # the same, strings with multiple alleles will match: 10 692 will never be 692 10
                                    if allele == sortedrefallele and float(percentid) == 100.00:
                                        # Increment the number of matches to each profile
                                        self.bestmatch[genome][sequencetype] += 1
                                    # Special handling of BACT000060 and BACT000065 genes. When the reference profile
                                    # has an allele of 'N', and the query allele doesn't, set the allele to 'N', and
                                    # count it as a match
                                    elif gene == 'BACT000060' or gene == 'BACT000065':
                                        if sortedrefallele == 'N' and allele != 'N':
                                            # Increment the number of matches to each profile
                                            self.bestmatch[genome][sequencetype] += 1
                                    elif allele == sortedrefallele and sortedrefallele == 'N':
                                        # Increment the number of matches to each profile
                                        self.bestmatch[genome][sequencetype] += 1
                        # Get the best number of matches
                        # From: https://stackoverflow.com/questions/613183/sort-a-python-dictionary-by-value
                        try:
                            sortedmatches = sorted(self.bestmatch[genome].items(), key=operator.itemgetter(1),
                                                   reverse=True)[0][1]
                        # If there are no matches, set :sortedmatches to zero
                        except IndexError:
                            sortedmatches = 0
                        # Otherwise, the query profile matches the reference profile
                        if int(sortedmatches) == header:
                            # Iterate through best match
                            for sequencetype, matches in self.bestmatch[genome].items():
                                if matches == sortedmatches:
                                    for gene in profiledata[sequencetype]:
                                        # Populate resultProfile with the genome, best match to profile, # of matches
                                        # to the profile, gene, query allele(s), reference allele(s), and % identity
                                        self.resultprofile[genome][sequencetype][sortedmatches][gene][
                                            list(self.bestdict[genome][gene]
                                                .keys())[0]] = str(list(self.bestdict[genome][gene].values())[0])
                                    sample[self.analysistype].sequencetype = sequencetype
                                    sample[self.analysistype].matchestosequencetype = matches
                        # If there are fewer matches than the total number of genes in the typing scheme
                        elif 0 < int(sortedmatches) < header:
                            mismatches = []
                            # Iterate through the sequence types and the number of matches in bestDict for each genome
                            for sequencetype, matches in self.bestmatch[genome].items():
                                # If the number of matches for a profile matches the best number of matches
                                if matches == sortedmatches:
                                    # Iterate through the gene in the analysis
                                    for gene in profiledata[sequencetype]:
                                        # Get the reference allele as above
                                        refallele = profiledata[sequencetype][gene]
                                        # As above get the reference allele split and ordered as necessary
                                        if len(refallele.split(" ")) > 1:
                                            intrefallele = map(int, refallele.split(" "))
                                            sortedrefallele = " ".join(str(allele) for allele in sorted(intrefallele))
                                        else:
                                            sortedrefallele = refallele
                                        # Populate self.mlstseqtype with the genome, best match to profile, # of matches
                                        # to the profile, gene, query allele(s), reference allele(s), and % identity
                                        if self.updateprofile:
                                            self.mlstseqtype[genome][sequencetype][sortedmatches][gene][
                                                str(list(self.bestdict[genome][gene]
                                                    .keys())[0])][sortedrefallele] = str(list(self.bestdict[genome][gene]
                                                                                        .values())[0])
                                        else:
                                            self.resultprofile[genome][sequencetype][sortedmatches][gene][
                                                list(self.bestdict[genome][gene].keys())[0]] \
                                                = str(list(self.bestdict[genome][gene])[0])
                                                #= str(list(self.bestdict[genome][gene].values())[0])
                                            if sortedrefallele != list(self.bestdict[sample.name][gene].keys())[0]:
                                                mismatches.append(
                                                    ({gene: ('{} ({})'.format(list(self.bestdict[sample.name][gene]
                                                                              .keys())[0], sortedrefallele))}))
                                        if not self.updateprofile or self.analysistype == 'mlst':
                                            self.resultprofile[genome][sequencetype][sortedmatches][gene][
                                                list(self.bestdict[genome][gene]
                                                    .keys())[0]] = str(list(self.bestdict[genome][gene].values())[0])
                                            sample[self.analysistype].mismatchestosequencetype = mismatches
                                            sample[self.analysistype].sequencetype = sequencetype
                                            sample[self.analysistype].matchestosequencetype = matches
                            # Add the new profile to the profile file (if the option is enabled)
                            if self.updateprofile and self.analysistype != 'mlst':
                                self.reprofiler(int(header), genome, sample)
                        elif sortedmatches == 0:
                            for gene in sample[self.analysistype].allelenames:
                                # Populate the results profile with negative values for sequence type and sorted matches
                                self.resultprofile[genome]['NA'][sortedmatches][gene]['NA'] = 0
                            # Add the new profile to the profile file (if the option is enabled)
                            if self.updateprofile:
                                self.reprofiler(int(header), genome, sample)
                            sample[self.analysistype].sequencetype = 'NA'
                            sample[self.analysistype].matchestosequencetype = 'NA'
                            sample[self.analysistype].mismatchestosequencetype = 'NA'
                        else:
                            sample[self.analysistype].matchestosequencetype = 'NA'
                            sample[self.analysistype].mismatchestosequencetype = 'NA'
                            sample[self.analysistype].sequencetype = 'NA'
                        dotter()
                else:
                    sample[self.analysistype].matchestosequencetype = 'NA'
                    sample[self.analysistype].mismatchestosequencetype = 'NA'
                    sample[self.analysistype].sequencetype = 'NA'

            else:
                sample[self.analysistype].matchestosequencetype = 'NA'
                sample[self.analysistype].mismatchestosequencetype = 'NA'
                sample[self.analysistype].sequencetype = 'NA'
73938
def build_proteintable(pqdb, headerfields, mergecutoff, isobaric=False,
                       precursor=False, probability=False, fdr=False,
                       pep=False, genecentric=False):
    """"""
    pdmap = create_featuredata_map(pqdb, genecentric=genecentric,
                                   psm_fill_fun=pinfo.add_psms_to_proteindata,
                                   pgene_fill_fun=pinfo.add_protgene_to_protdata,
                                   count_fun=pinfo.count_peps_psms,
                                   get_uniques=True)
    empty_return = lambda x, y, z: {}
    iso_fun = {True: get_isobaric_quant, False: empty_return}[isobaric]
    ms1_fun = {True: get_precursor_quant, False: empty_return}[precursor]
    prob_fun = {True: get_prot_probability,
                False: empty_return}[probability]
    fdr_fun = {True: get_prot_fdr,
               False: empty_return}[fdr]
    pep_fun = {True: get_prot_pep,
               False: empty_return}[pep]
    pdata_fun = {True: get_protein_data_genecentric,
                 False: get_protein_data}[genecentric is not False]
    protein_sql, sqlfieldmap = pqdb.prepare_mergetable_sql(precursor, isobaric,
                                                           probability, fdr,
                                                           pep)
    accession_field = prottabledata.ACCESSIONS[genecentric]
    proteins = pqdb.get_merged_features(protein_sql)
    protein = next(proteins)
    outprotein = {accession_field: protein[sqlfieldmap['p_acc']]}
    check_prot = {k: v for k, v in outprotein.items()}
    if not mergecutoff or protein_pool_fdr_cutoff(protein, sqlfieldmap,
                                                  mergecutoff):
        fill_mergefeature(outprotein, iso_fun, ms1_fun, prob_fun, fdr_fun,
                          pep_fun, pdata_fun, protein, sqlfieldmap,
                          headerfields, pdmap, accession_field)
    for protein in proteins:
        if mergecutoff and not protein_pool_fdr_cutoff(protein, sqlfieldmap,
                                                       mergecutoff):
            continue
        p_acc = protein[sqlfieldmap['p_acc']]
        if p_acc != outprotein[accession_field]:
            # check if protein has been filled, otherwise do not output
            # sometimes proteins have NA in all fields
            if outprotein != check_prot:
                yield outprotein
            outprotein = {accession_field: p_acc}
            check_prot = {k: v for k, v in outprotein.items()}
        fill_mergefeature(outprotein, iso_fun, ms1_fun, prob_fun, fdr_fun,
                          pep_fun, pdata_fun, protein, sqlfieldmap,
                          headerfields, pdmap, accession_field)
    if outprotein != check_prot:
        yield outprotein
673720
def blastnprep(self):
        """"""
        # Populate threads for each gene, genome combination
        for sample in self.metadata:
            if sample.general.bestassemblyfile != 'NA':
                #
                # sample[self.analysistype].alleleresults = GenObject()
                sample[self.analysistype].closealleles = dict()
                sample[self.analysistype].mismatches = dict()
                sample[self.analysistype].alignmentlength = dict()
                sample[self.analysistype].subjectlength = dict()
                sample[self.analysistype].queryid = dict()
                sample[self.analysistype].start = dict()
                sample[self.analysistype].end = dict()
                sample[self.analysistype].queryseq = dict()
                if type(sample[self.analysistype].allelenames) == list:
                    for allele in sample[self.analysistype].combinedalleles:
                        # Add each fasta/allele file combination to the threads
                        self.runblast(sample.general.bestassemblyfile, allele, sample)
458316
def _get_biodata(base_file, args):
    
    with open(base_file) as in_handle:
        config = yaml.safe_load(in_handle)
    config["install_liftover"] = False
    config["genome_indexes"] = args.aligners
    ann_groups = config.pop("annotation_groups", {})
    config["genomes"] = [_setup_genome_annotations(g, args, ann_groups)
                         for g in config["genomes"] if g["dbkey"] in args.genomes]
    return config
670565
def create_dexseq_annotation(gff, count_file):
    
    out_file = count_file + ".ann"
    if file_exists(out_file):
        return out_file
    with file_transaction(out_file) as tx_out:
        with open(tx_out, 'w') as out_handle:
            with open(gff) as in_handle:
                for line in in_handle:
                    cols = line.strip().split("\t")
                    if cols[2] == "exonic_part":
                        exon = [f for f in cols[8].split(";") if f.strip().startswith("exonic_part_number")]
                        gene = [f for f in cols[8].split(";") if f.strip().startswith("gene_id")]
                        exon = exon[0].replace("\"", "").split()[1]
                        gene = gene[0].replace("\"", "").split()[1]
                        length = int(cols[4]) - int(cols[3]) + 1
                        line = "%s:%s\t%s\t%s\t%s\t%s\t%s\t%s\n" % (gene, exon, gene,
                                                                    cols[0], cols[3],
                                                                    cols[4],
                                                                    length, cols[6])
                        out_handle.write(line)
1288342
def build_peptidetable(pqdb, headerfields, isobaric=False,
                       precursor=False, fdr=False, pep=False,
                       genecentric=False):
    """"""
    peptidedatamap = create_featuredata_map(pqdb, genecentric=genecentric, 
            psm_fill_fun=add_psm_to_peptidedata, pgene_fill_fun=add_protgene_to_pepdata,
            is_peptides=True)
    count_psms(peptidedatamap)
    empty_return = lambda x, y, z: {}
    iso_fun = {True: get_isobaric_quant, False: empty_return}[isobaric]
    ms1_fun = {True: get_precursor_quant, False: empty_return}[precursor]
    fdr_fun = {True: get_pep_fdr,
               False: empty_return}[fdr]
    pep_fun = {True: get_peptide_pep,
               False: empty_return}[pep]
    pdata_fun = get_protein_data
    peptide_sql, sqlfieldmap = pqdb.prepare_mergetable_sql(precursor, isobaric,
                                                           probability=False,
                                                           fdr=fdr, pep=pep)
    peptides = pqdb.get_merged_features(peptide_sql)
    peptide = next(peptides)
    outpeptide = {peptabledata.HEADER_PEPTIDE: peptide[sqlfieldmap['p_acc']]}
    check_pep = {k: v for k, v in outpeptide.items()}
    fill_mergefeature(outpeptide, iso_fun, ms1_fun, empty_return, fdr_fun,
                      pep_fun, pdata_fun, peptide, sqlfieldmap,
                      headerfields, peptidedatamap)
    for peptide in peptides:
        p_seq = peptide[sqlfieldmap['p_acc']]
        if p_seq != outpeptide[peptabledata.HEADER_PEPTIDE]:
            if outpeptide != check_pep:
                yield outpeptide
            outpeptide = {peptabledata.HEADER_PEPTIDE: p_seq}
            check_pep = {k: v for k, v in outpeptide.items()}
        fill_mergefeature(outpeptide, iso_fun, ms1_fun, empty_return, fdr_fun,
                          pep_fun, pdata_fun, peptide, sqlfieldmap,
                          headerfields, peptidedatamap)
    yield outpeptide
457146
def build_tree_from_alignment(aln, moltype=DNA, best_tree=False, params=None):
    
    if params is None:
        params = {}

    if moltype == DNA or moltype == RNA:
        params['-nt'] = True
    elif moltype == PROTEIN:
        params['-nt'] = False
    else:
        raise ValueError, \
                "FastTree does not support moltype: %s" % moltype.label

    if best_tree:
        params['-slow'] = True

    #Create mapping between abbreviated IDs and full IDs
    int_map, int_keys = aln.getIntMap()
    #Create SequenceCollection from int_map.
    int_map = SequenceCollection(int_map,MolType=moltype)

    app = FastTree(params=params)

    result = app(int_map.toFasta())
    tree = DndParser(result['Tree'].read(), constructor=PhyloNode)
    #remap tip names
    for tip in tree.tips():
        tip.Name = int_keys[tip.Name]

    return tree
25832
def generate(cls):
        
        return cls(libnacl.randombytes(PrivateKey.SIZE), encoder=encoding.RawEncoder)
1009590
def _sample_template(sample, out_dir):
    """"""
    bam_fn = dd.get_work_bam(sample)
    genome = dd.get_genome_build(sample)
    if genome in supported:
        peaks = sample.get("peaks_files", []).get("main")
        if peaks:
            r_code = ("library(ChIPQC);\n"
                      "sample = ChIPQCsample(\"{bam_fn}\","
                      "\"{peaks}\", "
                      "annotation = \"{genome}\","
                      ");\n"
                      "ChIPQCreport(sample);\n")
            r_code_fn = os.path.join(out_dir, "chipqc.r")
            with open(r_code_fn, 'w') as inh:
                inh.write(r_code.format(**locals()))
            return r_code_fn
1187703
def main():
    
    parser = argparse.ArgumentParser(
        description='monoseq: pretty-printing DNA and protein sequences',
        epilog='If INPUT is in FASTA format, each record is pretty-printed '
        'after printing its name and ANNOTATION (if supplied) is used by '
        'matching chromosome/record name. If INPUT contains a raw sequence, '
        'only the first chromosome in ANNOTATION is used.')
    parser.add_argument(
        'sequence_file', metavar='INPUT', nargs='?', default=sys.stdin,
        type=argparse.FileType('r'), help='file to read sequence(s) from, '
        'can be in FASTA format (default: standard input)')
    parser.add_argument(
        '-b', '--block-length', metavar='LENGTH', dest='block_length',
        type=int, default=10, help='block length in letters (default: 10)')
    parser.add_argument(
        '-l', '--blocks-per-line', metavar='BLOCKS', dest='blocks_per_line',
        type=int, default=6, help='blocks per line (default: 6)')
    parser.add_argument(
        '-a', '--annotation', metavar='POS', dest='annotation', nargs=2,
        action='append', type=int, help='first and last positions of '
        'subsequence to annotate (allowed more than once)')
    parser.add_argument(
        '-e', '--bed', metavar='ANNOTATION', dest='annotation_file',
        type=argparse.FileType('r'), help='file to read annotation from in '
        'BED format')

    args = parser.parse_args()
    pprint(_until_eof(args.sequence_file), annotation=args.annotation,
           annotation_file=args.annotation_file,
           block_length=args.block_length,
           blocks_per_line=args.blocks_per_line)
413031
def genome_coverage(genomes, scaffold_coverage, total_bases):
	
	coverage = {}
	custom = {}
	std = {}
	for genome in genomes:
		for sequence in parse_fasta(genome):
			scaffold = sequence[0].split('>')[1].split()[0]
			coverage, std = sum_coverage(coverage, std, genome, scaffold, sequence, scaffold_coverage)
			custom = calc_custom(custom, genome, scaffold, sequence, scaffold_coverage, total_bases)
	std = calc_std(std)
	custom_std = calc_std(custom)
	custom_av = {}
	for genome in custom:
		custom_av[genome] = []
		for sample in custom[genome]:
			custom_av[genome].append(numpy.mean(sample))
	for genome in coverage:
		print('%s\t%s' % (genome, coverage[genome][0][1]))
	if total_bases is True:
		total_bases = calc_total_mapped_bases(coverage)
	absolute = absolute_abundance(coverage, total_bases)
	for genome in coverage:
		calculated = []
		for calc in coverage[genome]:
			calculated.append(calc[0] / calc[1])
		coverage[genome] = calculated
	relative = relative_abundance(coverage)
	return coverage, std, absolute, relative, custom_av, custom_std
1340134
def _keep_assembled_chrom(bam_file, genome, config):
    """"""
    fai = "%s.fai" % genome
    chrom = []
    with open(fai) as inh:
        for line in inh:
            c = line.split("\t")[0]
            if c.find("_") < 0:
                chrom.append(c)
    chroms = " ".join(chrom)
    out_file = utils.append_stem(bam_file, '_chrom')
    samtools = config_utils.get_program("samtools", config)
    if not utils.file_exists(out_file):
        with file_transaction(out_file) as tx_out:
            cmd = "{samtools} view -b {bam_file} {chroms} > {tx_out}"
            do.run(cmd.format(**locals()), "Remove contigs from %s" % bam_file)
        bam.index(out_file, config)
    return out_file