Run skesa to assemble genomes
4794
def skesa_assemble(self):
        
        with progressbar(self.metadata) as bar:
            for sample in bar:
                # Initialise the assembly command
                sample.commands.assemble = str()
                try:
                    if sample.general.trimmedcorrectedfastqfiles:
                        # If the sample is a pure isolate, assemble it. Otherwise, run the pre-metagenome pipeline
                        try:
                            status = sample.run.Description
                        except AttributeError:
                            status = 'unknown'
                        if status == 'metagenome':
                            self.merge(sample)
                        else:
                            # Set the output directory
                            sample.general.assembly_output = os.path.join(sample.general.outputdirectory,
                                                                          'assembly_output')
                            make_path(sample.general.assembly_output)
                            sample.general.assemblyfile = os.path.join(sample.general.assembly_output,
                                                                       '{name}_unfiltered.fasta'
                                                                       .format(name=sample.name))
                            sample.general.bestassemblyfile = os.path.join(sample.general.assembly_output,
                                                                           '{name}.fasta'
                                                                           .format(name=sample.name))
                            fastqfiles = sample.general.trimmedcorrectedfastqfiles

                            # Set the the forward fastq files
                            sample.general.assemblyfastq = fastqfiles
                            forward = fastqfiles[0]
                            gz = True if '.gz' in forward else False
                            # If there are two fastq files
                            if len(fastqfiles) == 2:
                                # Set the reverse fastq name https://github.com/ncbi/SKESA/issues/7
                                sample.commands.assemble = 'skesa --fastq {fastqfiles} --cores {threads} ' \
                                                           '--use_paired_ends --vector_percent 1 ' \
                                                           '--contigs_out {contigs}'\
                                    .format(fastqfiles=','.join(fastqfiles),
                                            threads=self.cpus,
                                            contigs=sample.general.assemblyfile)
                            # Same as above, but use single read settings for the assembler
                            else:
                                sample.commands.assemble = 'skesa --fastq {fastqfiles} --cores {threads} ' \
                                                           '--vector_percent 1 --contigs_out {contigs}'\
                                    .format(fastqfiles=','.join(fastqfiles),
                                            threads=self.cpus,
                                            contigs=sample.general.assemblyfile)
                    # If there are no fastq files, populate the metadata appropriately
                    else:
                        sample.general.assembly_output = 'NA'
                        sample.general.assemblyfastq = 'NA'
                        sample.general.bestassemblyfile = 'NA'
                except AttributeError:
                    sample.general.assembly_output = 'NA'
                    sample.general.assemblyfastq = 'NA'
                    sample.general.trimmedcorrectedfastqfiles = 'NA'
                    sample.general.bestassemblyfile = 'NA'
                if sample.commands.assemble and not os.path.isfile(sample.general.assemblyfile):
                    # Run the assembly
                    out, err = run_subprocess(sample.commands.assemble)
                    write_to_logfile(sample.commands.assemble,
                                     sample.commands.assemble,
                                     self.logfile,
                                     sample.general.logout,
                                     sample.general.logerr,
                                     None,
                                     None)
                    write_to_logfile(out,
                                     err,
                                     self.logfile,
                                     sample.general.logout,
                                     sample.general.logerr,
                                     None,
                                     None)
332298
def annotatethreads(self):
        
        # Move the files to subfolders and create objects
        self.runmetadata = createobject.ObjectCreation(self)
        # Fix headers
        self.headers()
        printtime('Performing prokka analyses', self.start)
        # Create and start threads
        for i in range(self.cpus):
            # Send the threads to the appropriate destination function
            threads = Thread(target=self.annotate, args=())
            # Set the daemon to true - something to do with thread management
            threads.setDaemon(True)
            # Start the threading
            threads.start()
        for sample in self.metadata.samples:
            # Create the prokka attribute in the metadata object
            setattr(sample, 'prokka', GenObject())
            sample.prokka.outputdir = os.path.join(sample.general.outputdirectory, 'prokka')
            if not os.path.isdir(sample.prokka.outputdir):
                os.makedirs(sample.prokka.outputdir)
            # TODO Incorporate MASH/rMLST/user inputted genus, species results in the system call
            # Create the system call
            # prokka 2014-SEQ-0275.fasta --force --genus Escherichia --species coli --usegenus --addgenes
            # --prefix 2014-SEQ-0275 --locustag EC0275 --outputdir /path/to/sequences/2014-SEQ-0275/prokka
            sample.prokka.command = 'prokka {} ' \
                                    '--force ' \
                                    '--genus {} ' \
                                    '--species {} ' \
                                    '--usegenus ' \
                                    '--addgenes ' \
                                    '--prefix {} ' \
                                    '--locustag {} ' \
                                    '--outdir {}' \
                .format(sample.general.fixedheaders,
                        self.genus, self.species, sample.name, sample.name, sample.prokka.outputdir)
            self.queue.put(sample)
        self.queue.join()
43344
def get_resources(genome, ref_file, data):
    
    base_dir = os.path.normpath(os.path.dirname(ref_file))
    resource_file = os.path.join(base_dir, "%s-resources.yaml" % genome.replace("-test", ""))
    if not os.path.exists(resource_file):
        raise IOError("Did not find resource file for %s: %s\n"
                      "To update bcbio_nextgen.py with genome resources for standard builds, run:\n"
                      "bcbio_nextgen.py upgrade -u skip"
                      % (genome, resource_file))
    with open(resource_file) as in_handle:
        resources = yaml.safe_load(in_handle)

    def resource_file_path(x):
        if isinstance(x, six.string_types) and os.path.exists(os.path.join(base_dir, x)):
            return os.path.normpath(os.path.join(base_dir, x))
        return x
    cleaned = utils.dictapply(resources, resource_file_path)
    return ensure_annotations(cleaned, data)
574041
def index(ref_file, out_dir, data):
    
    (ref_dir, local_file) = os.path.split(ref_file)
    gtf_file = dd.get_gtf_file(data)
    if not utils.file_exists(gtf_file):
        raise ValueError("%s not found, could not create a star index." % (gtf_file))
    if not utils.file_exists(out_dir):
        with tx_tmpdir(data, os.path.dirname(out_dir)) as tx_out_dir:
            num_cores = dd.get_cores(data)
            cmd = ("STAR --genomeDir {tx_out_dir} --genomeFastaFiles {ref_file} "
                   "--runThreadN {num_cores} "
                   "--runMode genomeGenerate --sjdbOverhang 99 --sjdbGTFfile {gtf_file}")
            do.run(cmd.format(**locals()), "Index STAR")
            if os.path.exists(out_dir):
                shutil.rmtree(out_dir)
            shutil.move(tx_out_dir, out_dir)
    return out_dir
92729
def run(_, data, out_dir):
    
    # logger.info("Number of aligned reads < than 0.60 in %s: %s" % (dd.get_sample_name(data), ratio))
    logger.info("Running kraken to determine contaminant: %s" % dd.get_sample_name(data))
    # ratio = bam.get_aligned_reads(bam_file, data)
    out = out_stats = None
    db = tz.get_in(["config", "algorithm", "kraken"], data)
    if db and isinstance(db, (list, tuple)):
        db = db[0]
    kraken_cmd = config_utils.get_program("kraken", data["config"])
    if db == "minikraken":
        db = os.path.join(install._get_data_dir(), "genomes", "kraken", "minikraken")

    if not os.path.exists(db):
        logger.info("kraken: no database found %s, skipping" % db)
        return {"kraken_report": "null"}

    if not os.path.exists(os.path.join(out_dir, "kraken_out")):
        work_dir = os.path.dirname(out_dir)
        utils.safe_makedir(work_dir)
        num_cores = data["config"]["algorithm"].get("num_cores", 1)
        fn_file = data["files_orig"][0] if dd.get_save_diskspace(data) else data["files"][0]
        if fn_file.endswith("bam"):
            logger.info("kraken: need fastq files as input")
            return {"kraken_report": "null"}
        with tx_tmpdir(data) as tx_tmp_dir:
            with utils.chdir(tx_tmp_dir):
                out = os.path.join(tx_tmp_dir, "kraken_out")
                out_stats = os.path.join(tx_tmp_dir, "kraken_stats")
                cat = "zcat" if fn_file.endswith(".gz") else "cat"
                cl = ("{cat} {fn_file} | {kraken_cmd} --db {db} --quick "
                      "--preload --min-hits 2 "
                      "--threads {num_cores} "
                      "--output {out} --fastq-input /dev/stdin  2> {out_stats}").format(**locals())
                do.run(cl, "kraken: %s" % dd.get_sample_name(data))
                if os.path.exists(out_dir):
                    shutil.rmtree(out_dir)
                shutil.move(tx_tmp_dir, out_dir)
    metrics = _parse_kraken_output(out_dir, db, data)
    return metrics
1159976
def liftover_to_genome(pass_pos, gtf):
    """"""

    fixed_pos = []
    for pos in pass_pos:
        if pos["chrom"] not in gtf:
            continue
        db_pos = gtf[pos["chrom"]][0]
        mut = _parse_mut(pos["sv"])
        print([db_pos, pos])
        if db_pos[3] == "+":
            pos['pre_pos'] = db_pos[1] + pos["pre_pos"] + 1
        else:
            pos['pre_pos'] = db_pos[2] - (pos["pre_pos"] - 1)
        pos['chrom'] = db_pos[0]
        pos['nt'] = list(mut[0])
        fixed_pos.append(pos)

    _print_header(fixed_pos)
    for pos in fixed_pos:
        print_vcf(pos)
421813
def _importSNPs_CasavaSNP(setName, species, genomeSource, snpsFile) :
	""
	printf('importing SNP set %s for species %s...' % (setName, species))

	snpData = SNPsTxtFile(snpsFile)
	
	CasavaSNP.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 = CasavaSNP()
		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.start -= 1
		snp.end -= 1
		snp.save()
		pBar.update(label = pLabel)

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

	printf('saving...')
	conf.db.endTransaction()
	printf('creating indexes...')
	CasavaSNP.ensureGlobalIndex(('start', 'chromosomeNumber', 'setName'))
	printf('importation of SNP set %s for species %s done.' %(setName, species))
	
	return True
586078
def main(input_bed, output_file, output_features=False, genome=None,
         only_canonical=False, short=False, extended=False, high_confidence=False,
         ambiguities_method=False, coding_only=False, collapse_exons=False, work_dir=False, is_debug=False):
    
    logger.init(is_debug_=is_debug)

    if not genome:
        raise click.BadParameter('Error: please, specify genome build name with -g (e.g. `-g hg19`)', param='genome')

    if short:
        if extended:        raise click.BadParameter('--short and --extended can\'t be set both', param='extended')
        if output_features: raise click.BadParameter('--short and --output-features can\'t be set both', param='output_features')
    elif output_features or extended:
        extended = True
        short    = False

    if not verify_file(input_bed):
        click.BadParameter(f'Usage: {__file__} Input_BED_file -g hg19 -o Annotated_BED_file [options]', param='input_bed')
    input_bed = verify_file(input_bed, is_critical=True, description=f'Input BED file for {__file__}')

    if work_dir:
        work_dir = join(adjust_path(work_dir), os.path.splitext(basename(input_bed))[0])
        safe_mkdir(work_dir)
        info(f'Created work directory {work_dir}')
    else:
        work_dir = mkdtemp('bed_annotate')
        debug('Created temporary work directory {work_dir}')

    input_bed = clean_bed(input_bed, work_dir)
    input_bed = verify_bed(input_bed, is_critical=True, description=f'Input BED file for {__file__} after cleaning')

    output_file = adjust_path(output_file)

    output_file = annotate(
        input_bed, output_file, work_dir, genome=genome,
        only_canonical=only_canonical, short=short, extended=extended,
        high_confidence=high_confidence, collapse_exons=collapse_exons,
        output_features=output_features,
        ambiguities_method=ambiguities_method, coding_only=coding_only,
        is_debug=is_debug)

    if not work_dir:
        debug(f'Removing work directory {work_dir}')
        shutil.rmtree(work_dir)

    info(f'Done, saved to {output_file}')
1032106
def main():
    
    # Define Parser object and add to toil
    url_prefix = 'https://s3-us-west-2.amazonaws.com/cgl-pipeline-inputs/'
    parser = argparse.ArgumentParser(description=main.__doc__, formatter_class=argparse.RawTextHelpFormatter)
    parser.add_argument('--config', required=True,
                        help='Path to configuration file for samples, one per line. UUID,URL_to_bamfile. '
                             'The URL may be a standard "http://", a "file://<abs_path>", or "s3://<bucket>/<key>"')
    parser.add_argument('--gtf', help='URL to annotation GTF file',
                        default=url_prefix + 'rnaseq_cgl/gencode.v23.annotation.gtf')
    parser.add_argument('--gtf-pickle', help='Pickled GTF file',
                        default=url_prefix + 'spladder/gencode.v23.annotation.gtf.pickle')
    parser.add_argument('--gtf-m53', help='M53 preprocessing annotation table',
                        default=url_prefix + 'spladder/gencode.v23.annotation.gtf.m53')
    parser.add_argument('--positions', help='URL to SNP positions over genes file (TSV)',
                        default=url_prefix + 'spladder/positions_fixed.tsv')
    parser.add_argument('--genome', help='URL to Genome fasta',
                        default=url_prefix + 'rnaseq_cgl/hg38_no_alt.fa')
    parser.add_argument('--genome-index', help='Index file (fai) of genome',
                        default=url_prefix + 'spladder/hg38_no_alt.fa.fai')
    parser.add_argument('--ssec', default=None, help='Path to master key used for downloading encrypted files.')
    parser.add_argument('--output-s3-dir', default=None, help='S3 Directory of the form: s3://bucket/directory')
    parser.add_argument('--output-dir', default=None, help='full path where final results will be output')
    parser.add_argument('--sudo', action='store_true', default=False,
                        help='Set flag if sudo is required to run Docker.')
    parser.add_argument('--star-index', help='URL to download STAR Index built from HG38/gencodev23 annotation.',
                        default=url_prefix + 'rnaseq_cgl/starIndex_hg38_no_alt.tar.gz')
    parser.add_argument('--fwd-3pr-adapter', help="Sequence for the FWD 3' Read Adapter.", default='AGATCGGAAGAG')
    parser.add_argument('--rev-3pr-adapter', help="Sequence for the REV 3' Read Adapter.", default='AGATCGGAAGAG')
    Job.Runner.addToilOptions(parser)
    args = parser.parse_args()
    # Sanity Checks
    if args.config:
        assert os.path.isfile(args.config), 'Config not found at: {}'.format(args.config)
    if args.ssec:
        assert os.path.isfile(args.ssec), 'Encryption key not found at: {}'.format(args.config)
    if args.output_s3_dir:
        assert args.output_s3_dir.startswith('s3://'), 'Wrong format for output s3 directory'
    # Program checks
    for program in ['curl', 'docker']:
        assert which(program), 'Program "{}" must be installed on every node.'.format(program)

    Job.Runner.startToil(Job.wrapJobFn(parse_input_samples, args), args)
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
26632
def run
      if @options[:list_assemblers]
        puts @assemblerman.list_assemblers
        return
      elsif @options[:install_assemblers]
        @assemblerman.install_assemblers(@options[:install_assemblers])
        return
      end

      if (@options[:left].nil? || @options[:right].nil?)
        logger.error "Reads must be provided with --left and --right"
        logger.error "Try --help for command-line help"
        exit(1)
      end

      unless (@options[:timelimit].nil?)
        logger.info "Time limit set to #{@options[:timelimit]}"
      end

      subsample_input

      res = @assemblerman.run_all_assemblers @options
      write_metadata res
      merge_assemblies res

    end
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
195645
def run_assembly(stmts, folder, pmcid, background_assertions=None):
    ''''''
    # Folder for index card output (scored submission)
    indexcard_prefix = folder + '/index_cards/' + pmcid
    # Folder for other outputs (for analysis, debugging)
    otherout_prefix = folder + '/other_outputs/' + pmcid

    # Do grounding mapping here
    # Load the TRIPS-specific grounding map and add to the default
    # (REACH-oriented) grounding map:
    trips_gm = load_grounding_map('trips_grounding_map.csv')
    default_grounding_map.update(trips_gm)
    gm = GroundingMapper(default_grounding_map)

    mapped_agent_stmts = gm.map_agents(stmts)
    renamed_agent_stmts = gm.rename_agents(mapped_agent_stmts)

    # Filter for grounding
    grounded_stmts = []
    for st in renamed_agent_stmts:
        if all([is_protein_or_chemical(a) for a in st.agent_list()]):
            grounded_stmts.append(st)

    # Instantiate the Preassembler
    pa = Preassembler(hierarchies)
    pa.add_statements(grounded_stmts)
    print('== %s ====================' % pmcid)
    print('%d statements collected in total.' % len(pa.stmts))

    # Combine duplicates
    unique_stmts = pa.combine_duplicates()
    print('%d statements after combining duplicates.' % len(unique_stmts))

    # Run BeliefEngine on unique statements
    epe = BeliefEngine()
    epe.set_prior_probs(pa.unique_stmts)

    # Build statement hierarchy
    related_stmts = pa.combine_related()
    # Run BeliefEngine on hierarchy
    epe.set_hierarchy_probs(related_stmts)
    print('%d statements after combining related.' % len(related_stmts))

    # Instantiate the mechanism linker
    # Link statements
    linked_stmts = MechLinker.infer_active_forms(related_stmts)
    linked_stmts += MechLinker.infer_modifications(related_stmts)
    linked_stmts += MechLinker.infer_activations(related_stmts)
    # Run BeliefEngine on linked statements
    epe.set_linked_probs(linked_stmts)
    # Print linked statements for debugging purposes
    print('Linked\n=====')
    for ls in linked_stmts:
        print(ls.inferred_stmt.belief, ls.inferred_stmt)
    print('=============')

    # Combine all statements including linked ones
    all_statements = related_stmts + [ls.inferred_stmt for ls in linked_stmts]

    # Instantiate a new preassembler
    pa = Preassembler(hierarchies, all_statements)
    # Build hierarchy again
    pa.combine_duplicates()
    # Choose the top-level statements
    related_stmts = pa.combine_related()

    # Remove top-level statements that came only from the prior
    if background_assertions is not None:
        nonbg_stmts = [stmt for stmt in related_stmts
                       if stmt not in background_assertions]
    else:
        nonbg_stmts = related_stmts

    # Dump top-level statements in a pickle
    with open(otherout_prefix + '.pkl', 'wb') as fh:
        pickle.dump(nonbg_stmts, fh)

    # Flatten evidence for statements
    flattened_evidence_stmts = flatten_evidence(nonbg_stmts)

    # Start a card counter
    card_counter = 1
    # We don't limit the number of cards reported in this round
    card_lim = float('inf')
    top_stmts = []
    ###############################################
    # The belief cutoff for statements
    belief_cutoff = 0.3
    ###############################################
    # Sort by amount of evidence
    for st in sorted(flattened_evidence_stmts,
                     key=lambda x: x.belief, reverse=True):
        if st.belief >= belief_cutoff:
            print(st.belief, st)
        if st.belief < belief_cutoff:
            print('SKIP', st.belief, st)

        # If it's background knowledge, we skip the statement
        if is_background_knowledge(st):
            print('This statement is background knowledge - skipping.')
            continue

        # Assemble IndexCards
        ia = IndexCardAssembler([st], pmc_override=pmcid)
        ia.make_model()
        # If the index card was actually made 
        # (not all statements can be assembled into index cards to
        # this is often not the case)
        if ia.cards:
            # Save the index card json
            ia.save_model(indexcard_prefix + '-%d.json' % card_counter)
            card_counter += 1
            top_stmts.append(st)
            if card_counter > card_lim:
                break

    # Print the English-assembled model for debugging purposes
    ea = EnglishAssembler(top_stmts)
    print('=======================')
    print(ea.make_model().encode('utf-8'))
    print('=======================')

    # Print the statement graph
    graph = render_stmt_graph(nonbg_stmts)
    graph.draw(otherout_prefix + '_graph.pdf', prog='dot')
    # Print statement diagnostics
    print_stmts(pa.stmts, otherout_prefix + '_statements.tsv')
    print_stmts(related_stmts, otherout_prefix + '_related_statements.tsv')
1318665
def save_refs(sam, fastas, genomes, s2b):
    
    if s2b is False:
        s2b = {}
    for fasta in fastas:
        # if no reference sequence supplied, get length from sam file
        if fasta is False:
            for line in sam:
                if line.startswith('@PG'):
                    break
                if line.startswith('@SQ') is False:
                    continue
                line = line.strip().split()
                contig, length = line[1].split(':', 1)[1], int(line[2].split(':', 1)[1])
                if contig not in s2b:
                    genome = 'n/a'
                    s2b[contig] = genome
                else:
                    genome = s2b[contig]
                if genome not in genomes:
                    genomes[genome] = {}
                if contig not in genomes[genome]:
                    genomes[genome][contig] = {}
                bp_stats = []
                for i in range(0, length):
                    bp_stats.append({'A':0, 'T':0, 'G':0, 'C':0, 'N':0, \
                                        'In':[], 'Del':[], 'ref':'N'})
                genomes[genome][contig][sam.name] = {'bp_stats':bp_stats}
                return genomes, s2b
        # save reference sequences, if available
        genome = fasta.rsplit('.', 1)[0]
        if genome not in genomes:
            genomes[genome] = {}
        bp_stats = []
        for seq in parse_fasta(fasta):
            contig = seq[0].split('>')[1].split()[0]
            s2b[contig] = genome
            if contig not in genomes[genome]:
                genomes[genome][contig] = {}
            bp_stats = []
            for base in seq[1]:
                bp_stats.append({'A':0, 'T':0, 'G':0, 'C':0, 'N':0, \
                                    'In':[], 'Del':[], 'ref':base.upper()})
            genomes[genome][contig][sam.name] = {'bp_stats':bp_stats}
        return genomes, s2b
587865
def generate_bio_assembly(data_api, struct_inflator):
    """"""
    bioassembly_count = 0
    for bioassembly in data_api.bio_assembly:
        bioassembly_count += 1
        for transform in bioassembly["transformList"]:
            struct_inflator.set_bio_assembly_trans(bioassembly_count,
                                                   transform["chainIndexList"],
                                                   transform["matrix"])
585055
def run(name, chip_bam, input_bam, genome_build, out_dir, method, resources, data):
    
    # output file name need to have the caller name
    config = dd.get_config(data)
    out_file = os.path.join(out_dir, name + "_peaks_macs2.xls")
    macs2_file = os.path.join(out_dir, name + "_peaks.xls")
    if utils.file_exists(out_file):
        _compres_bdg_files(out_dir)
        return _get_output_files(out_dir)
    macs2 = config_utils.get_program("macs2", config)
    options = " ".join(resources.get("macs2", {}).get("options", ""))
    genome_size = bam.fasta.total_sequence_length(dd.get_ref_file(data))
    genome_size = "" if options.find("-g") > -1 else "-g %s" % genome_size
    paired = "-f BAMPE" if bam.is_paired(chip_bam) else ""
    with utils.chdir(out_dir):
        cmd = _macs2_cmd(method)
        try:
            do.run(cmd.format(**locals()), "macs2 for %s" % name)
            utils.move_safe(macs2_file, out_file)
        except subprocess.CalledProcessError:
            raise RuntimeWarning("macs2 terminated with an error.\n"
                                 "Please, check the message and report "
                                 "error if it is related to bcbio.\n"
                                 "You can add specific options for the sample "
                                 "setting resources as explained in docs: "
                                 "https://bcbio-nextgen.readthedocs.org/en/latest/contents/configuration.html#sample-specific-resources")
    _compres_bdg_files(out_dir)
    return _get_output_files(out_dir)
636235
def parse_cmdline():
    """"""
    parser = ArgumentParser(prog="genbank_get_genomes_by_taxon.py")
    parser.add_argument(
        "-o",
        "--outdir",
        dest="outdirname",
        required=True,
        action="store",
        default=None,
        help="Output directory (required)")
    parser.add_argument(
        "-t",
        "--taxon",
        dest="taxon",
        action="store",
        default=None,
        help="NCBI taxonomy ID")
    parser.add_argument(
        "-v",
        "--verbose",
        dest="verbose",
        action="store_true",
        default=False,
        help="Give verbose output")
    parser.add_argument(
        "-f",
        "--force",
        dest="force",
        action="store_true",
        default=False,
        help="Force file overwriting")
    parser.add_argument(
        "--noclobber",
        dest="noclobber",
        action="store_true",
        default=False,
        help="Don't nuke existing files")
    parser.add_argument(
        "-l",
        "--logfile",
        dest="logfile",
        action="store",
        default=None,
        help="Logfile location")
    parser.add_argument(
        "--format",
        dest="format",
        action="store",
        default="fasta",
        help="Output file format [gbk|fasta]")
    parser.add_argument(
        "--email",
        dest="email",
        required=True,
        action="store",
        default=None,
        help="Email associated with NCBI queries (required)")
    parser.add_argument(
        "--retries",
        dest="retries",
        action="store",
        default=20,
        type=int,
        help="Number of Entrez retry attempts per request.")
    parser.add_argument(
        "--batchsize",
        dest="batchsize",
        action="store",
        default=10000,
        type=int,
        help="Entrez record return batch size")
    parser.add_argument(
        "--timeout",
        dest="timeout",
        action="store",
        default=10,
        type=int,
        help="Timeout for URL connection (s)")
    return parser.parse_args()