%prog parallel genome.fasta N

    Partition the genome into parts and run separately. This is useful if MAKER
    is to be run on the grid.
177282
def parallel(args):
    
    from jcvi.formats.base import split

    p = OptionParser(parallel.__doc__)
    p.set_home("maker")
    p.set_tmpdir(tmpdir="tmp")
    p.set_grid_opts(array=True)
    opts, args = p.parse_args(args)

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

    genome, NN = args
    threaded = opts.threaded or 1
    tmpdir = opts.tmpdir

    mkdir(tmpdir)
    tmpdir = get_abs_path(tmpdir)

    N = int(NN)
    assert 1 <= N < 1000, "Required: 1 < N < 1000!"

    outdir = "outdir"
    fs = split([genome, outdir, NN])

    c = CTLFile("maker_opts.ctl")
    c.update_abs_path()
    if threaded > 1:
        c.update_tag("cpus", threaded)

    cwd = os.getcwd()
    dirs = []
    for name in fs.names:
        fn = get_abs_path(name)
        bn = op.basename(name)
        dirs.append(bn)
        c.update_tag("genome", fn)
        mkdir(bn)
        sh("cp *.ctl {0}".format(bn))

        os.chdir(bn)
        c.write_file("maker_opts.ctl")
        os.chdir(cwd)

    jobs = "jobs"
    fw = open(jobs, "w")
    print("\n".join(dirs), file=fw)
    fw.close()

    # Submit to grid
    ncmds = len(dirs)
    runfile = "array.sh"
    cmd = op.join(opts.maker_home, "bin/maker")
    if tmpdir:
        cmd += " -TMP {0}".format(tmpdir)

    engine = get_grid_engine()
    contents = arraysh.format(jobs, cmd) if engine == "SGE" \
                else arraysh_ua.format(N, threaded, jobs, cmd)
    write_file(runfile, contents)

    if engine == "PBS":
        return

    # qsub script
    outfile = "maker.\$TASK_ID.out"
    p = GridProcess(runfile, outfile=outfile, errfile=outfile,
                    arr=ncmds, grid_opts=opts)
    qsubfile = "qsub.sh"
    qsub = p.build()
    write_file(qsubfile, qsub)
866488
def multiprocess_permutation(bed_dict, mut_df, opts):
    
    chroms = sorted(bed_dict.keys())
    multiprocess_flag = opts['processes']>0
    if multiprocess_flag:
        num_processes = opts['processes']
    else:
        num_processes = 1
    num_permutations = opts['num_permutations']
    if not opts['by_sample']:
        obs_result = []
    else:
        uniq_samp = mut_df['Tumor_Sample'].unique()
        obs_result = pd.DataFrame(np.zeros((len(uniq_samp), len(cols))),
                                  index=uniq_samp, columns=cols)

    # initialize list containing output
    if not opts['score_dir']:
        result_list = [[0, 0, 0, 0, 0, 0, 0] for k in range(num_permutations)]
    else:
        result_list = [[0, 0, 0, 0, 0, 0, 0, 0, 0] for k in range(num_permutations)]

    # iterate over each chromosome
    for i in range(0, len(chroms), num_processes):
        if multiprocess_flag:
            pool = Pool(processes=num_processes)
            tmp_num_proc = len(chroms) - i if i + num_processes > len(chroms) else num_processes
            info_repeat = ((bed_dict[chroms[tmp_ix]], mut_df, opts)
                            for tmp_ix in range(i, i+tmp_num_proc))
            process_results = pool.imap(singleprocess_permutation, info_repeat)
            process_results.next = utils.keyboard_exit_wrapper(process_results.next)
            try:
                for chrom_result, obs_mutations in process_results:
                    for j in range(num_permutations):
                        result_list[j][0] += chrom_result[j][0]
                        result_list[j][1] += chrom_result[j][1]
                        result_list[j][2] += chrom_result[j][2]
                        result_list[j][3] += chrom_result[j][3]
                        result_list[j][4] += chrom_result[j][4]
                        result_list[j][5] += chrom_result[j][5]
                        result_list[j][6] += chrom_result[j][6]
                        if opts['score_dir']:
                            result_list[j][7] += chrom_result[j][7]
                            result_list[j][8] += chrom_result[j][8]

                    if not opts['by_sample']:
                        obs_result.append(obs_mutations)
                    else:
                        obs_result = obs_result + obs_mutations
            except KeyboardInterrupt:
                pool.close()
                pool.join()
                logger.info('Exited by user. ctrl-c')
                sys.exit(0)
            pool.close()
            pool.join()
        else:
            info = (bed_dict[chroms[i]], mut_df, opts)
            chrom_result, obs_mutations = singleprocess_permutation(info)
            for j in range(num_permutations):
                result_list[j][0] += chrom_result[j][0]
                result_list[j][1] += chrom_result[j][1]
                result_list[j][2] += chrom_result[j][2]
                result_list[j][3] += chrom_result[j][3]
                result_list[j][4] += chrom_result[j][4]
                result_list[j][5] += chrom_result[j][5]
                result_list[j][6] += chrom_result[j][6]
                if opts['score_dir']:
                    result_list[j][7] += chrom_result[j][7]
                    result_list[j][8] += chrom_result[j][8]
            if not opts['by_sample']:
                obs_result.append(obs_mutations)
            else:
                obs_result = obs_result + obs_mutations

    return result_list, obs_result
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}')
689777
def main(argv=None):
    

    if argv is None:
        argv = sys.argv

    # setup command line parser
    parser = U.OptionParser(version="%prog version: $Id$",
                            usage=usage,
                            description=globals()["__doc__"])

    group = U.OptionGroup(parser, "count-specific options")

    parser.add_option("--wide-format-cell-counts", dest="wide_format_cell_counts",
                      action="store_true",
                      default=False,
                      help=("output the cell counts in a wide format "
                            "(rows=genes, columns=cells)"))

    parser.add_option_group(group)

    # add common options (-h/--help, ...) and parse command line
    (options, args) = U.Start(parser, argv=argv, add_group_dedup_options=False)

    options.per_gene = True  # hardcodes counting to per-gene only

    U.validateSamOptions(options, group=False)

    if options.random_seed:
        np.random.seed(options.random_seed)

    if options.stdin != sys.stdin:
        in_name = options.stdin.name
        options.stdin.close()
    else:
        raise ValueError("Input on standard in not currently supported")

    if options.in_sam:
        in_mode = "r"
    else:
        in_mode = "rb"

    infile = pysam.Samfile(in_name, in_mode)

    # write out to tempfile and then sort to stdout
    tmpfilename = U.getTempFilename(dir=options.tmpdir)
    tmpfile = U.openFile(tmpfilename, mode="w")

    nInput, nOutput, input_reads = 0, 0, 0

    gene_tag = options.gene_tag
    metacontig2contig = None

    if options.chrom:
        inreads = infile.fetch(reference=options.chrom)
    else:
        if options.gene_transcript_map:
            metacontig2contig = sam_methods.getMetaContig2contig(
                infile, options.gene_transcript_map)
            metatag = "MC"
            inreads = sam_methods.metafetcher(infile, metacontig2contig, metatag)
            gene_tag = metatag
        else:
            inreads = infile.fetch()

    bundle_iterator = sam_methods.get_bundles(
        options,
        only_count_reads=True,
        metacontig_contig=metacontig2contig)

    # set up UMIClusterer functor with methods specific to
    # specified options.method

    processor = network.UMIClusterer(options.method)

    for bundle, key, status in bundle_iterator(inreads):
        if status == "single_read":
            continue

        gene, cell = key

        umis = bundle.keys()
        counts = {umi: bundle[umi]["count"] for umi in umis}

        nInput += sum(counts.values())

        while nInput >= input_reads + 1000000:
            input_reads += 1000000
            U.info("Parsed %i input reads" % input_reads)

        # group the umis
        groups = processor(
            counts,
            threshold=options.threshold)

        gene_count = len(groups)

        if options.per_cell:
            tmpfile.write("%s\n" % "\t".join((gene, cell.decode(), str(gene_count))))
        else:
            tmpfile.write("%s\n" % "\t".join((gene, str(gene_count))))

        nOutput += gene_count

    tmpfile.close()

    if options.per_cell:

        gene_counts_dict = {}

        with U.openFile(tmpfilename, mode="r") as inf:
            genes = set()
            cells = set()
            for line in inf:
                gene, cell, gene_count = line.strip().split("\t")
                genes.add(gene)
                cells.add(cell)

                if gene not in gene_counts_dict:
                    gene_counts_dict[gene] = {}

                gene_counts_dict[gene][cell] = gene_count

        if options.wide_format_cell_counts:  # write out in wide format

            options.stdout.write(
                "%s\t%s\n" % ("gene", "\t".join(sorted(cells))))

            for gene in sorted(genes):
                counts = []
                for cell in sorted(cells):
                    if cell in gene_counts_dict[gene]:
                        counts.append(gene_counts_dict[gene][cell])
                    else:
                        counts.append(0)
                options.stdout.write(
                    "%s\t%s\n" % (gene, "\t".join(map(str, counts))))

        else:  # write out in long format
            options.stdout.write("%s\t%s\t%s\n" % ("gene", "cell", "count"))
            for gene in sorted(genes):
                for cell in sorted(list(gene_counts_dict[gene].keys())):
                    options.stdout.write("%s\t%s\t%s\n" % (
                        gene, cell, gene_counts_dict[gene][cell]))
    else:
        options.stdout.write("%s\t%s\n" % ("gene", "count"))

        with U.openFile(tmpfilename, mode="r") as inf:
            for line in inf:
                options.stdout.write(line)

    os.unlink(tmpfilename)

    # output reads events and benchmark information.
    for event in bundle_iterator.read_events.most_common():
        U.info("%s: %s" % (event[0], event[1]))

    U.info("Number of (post deduplication) reads counted: %i" % nOutput)

    U.Stop()
1024609
def mask(args):
    
    p = OptionParser(mask.__doc__)
    p.add_option("--hard", dest="hard", default=False, action="store_true",
            help="Hard mask the low-complexity bases [default: %default]")
    opts, args = p.parse_args(args)

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

    genomefile, = args

    # entire pipeline
    countsfile = genomefile + ".counts"
    wm_mk_counts(infile=genomefile, outfile=countsfile)

    maskedfastafile = "%s.masked%s" % op.splitext(genomefile)
    wm_mk_masks(infile=countsfile, outfile=maskedfastafile, genomefile=genomefile)

    if opts.hard:
        hardmask(maskedfastafile)
1363805
def main():
    """"""
    parser = argparse.ArgumentParser()
    parser.add_argument("-f", "--file",
                        required=True,
                        help="input file",
                        type=str)

    parser.add_argument("-l", "--locus",
                        required=True,
                        help="Locus",
                        type=str)

    parser.add_argument("-k", "--kir",
                        help="Option for running with KIR",
                        action='store_true')

    parser.add_argument("-s", "--server",
                        help="Option for running with a server",
                        action='store_true')

    parser.add_argument("-v", "--verbose",
                        help="Option for running in verbose",
                        action='store_true')

    args = parser.parse_args()
    fastafile = args.file
    locus = args.locus

    verbose = False
    if args.verbose:
        verbose = True

    verbose = False
    if args.verbose:
        verbose = True

    kir = False
    if args.kir:
        kir = True

    serv = False
    if args.server:
        serv = True

    if verbose:
        logging.basicConfig(format='%(asctime)s - %(name)-35s - %(levelname)-5s - %(message)s',
                            datefmt='%m/%d/%Y %I:%M:%S %p',
                            level=logging.INFO)

    server = None
    if serv:
        server = BioSeqDatabase.open_database(driver="pymysql", user="root",
                                              passwd="", host="localhost",
                                              db="bioseqdb")

    seqann = BioSeqAnn(verbose=True, kir=kir)
    for seq in SeqIO.parse(fastafile, "fasta"):
        ann = seqann.annotate(seq, locus=locus)
        print('{:*^20} {:^20} {:*^20}'.format("", str(seq.description), ""))
        l = 0
        for f in ann.annotation:
            if isinstance(ann.annotation[f], DBSeq):
                print(f, ann.method, str(ann.annotation[f]), sep="\t")
                l += len(ann.annotation[f])
            else:
                print(f, ann.method, str(ann.annotation[f].seq), sep="\t")
                l += len(ann.annotation[f].seq)
        print("")

    if serv:
        server.close()
307452
def validate(args):
    
    import pyfasta

    p = OptionParser(validate.__doc__)
    p.add_option("--prefix", help="Add prefix to seqid")
    opts, args = p.parse_args(args)

    vcffile, fastafile = args
    pf = opts.prefix
    genome = pyfasta.Fasta(fastafile, record_class=pyfasta.MemoryRecord)
    fp = must_open(vcffile)
    match_ref = match_alt = total = 0
    for row in fp:
        if row[0] == '#':
            continue
        seqid, pos, id, ref, alt = row.split()[:5]
        total += 1
        if pf:
            seqid = pf + seqid
        pos = int(pos)
        if seqid not in genome:
            continue
        true_ref = genome[seqid][pos - 1]
        if total % 100000 == 0:
            print(total, "sites parsed", file=sys.stderr)
        if ref == true_ref:
            match_ref += 1
        elif alt == true_ref:
            match_alt += 1

    logging.debug("Match REF: {}".format(percentage(match_ref, total)))
    logging.debug("Match ALT: {}".format(percentage(match_alt, total)))
1170163
def _fasta_iter(fasta):
    
    # Adapted from Brent Pedersen.
    # http://www.biostars.org/p/710/#1412
    groups = (group for _, group in
              itertools.groupby(fasta, lambda line: line.startswith('>')))
    for group in groups:
        header = next(group)[1:].strip()
        sequence = ''.join(line.strip() for line in next(groups))
        yield header, sequence
1118843
def evaluate(self, genomes, config):
        
        if self.mode != MODE_PRIMARY:
            raise ModeError("Not in primary mode!")
        tasks = [(genome_id, genome, config) for genome_id, genome in genomes]
        id2genome = {genome_id: genome for genome_id, genome in genomes}
        tasks = chunked(tasks, self.secondary_chunksize)
        n_tasks = len(tasks)
        for task in tasks:
            self.inqueue.put(task)
        tresults = []
        while len(tresults) < n_tasks:
            try:
                sr = self.outqueue.get(block=True, timeout=0.2)
            except (queue.Empty, managers.RemoteError):
                continue
            tresults.append(sr)
        results = []
        for sr in tresults:
            results += sr
        for genome_id, fitness in results:
            genome = id2genome[genome_id]
            genome.fitness = fitness
310936
def snap(args):
    
    p = OptionParser(snap.__doc__)
    p.set_home("maker")
    opts, args = p.parse_args(args)

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

    species, gffile, fastafile = args
    mhome = opts.maker_home
    snapdir = "snap"
    mkdir(snapdir)

    cwd = os.getcwd()
    os.chdir(snapdir)

    newgffile = "training.gff3"
    logging.debug("Construct GFF file combined with sequence ...")
    sh("cat ../{0} > {1}".format(gffile, newgffile))
    sh('echo "##FASTA" >> {0}'.format(newgffile))
    sh("cat ../{0} >> {1}".format(fastafile, newgffile))

    logging.debug("Make models ...")
    sh("{0}/src/bin/maker2zff training.gff3".format(mhome))
    sh("{0}/exe/snap/fathom -categorize 1000 genome.ann genome.dna".format(mhome))
    sh("{0}/exe/snap/fathom -export 1000 -plus uni.ann uni.dna".format(mhome))
    sh("{0}/exe/snap/forge export.ann export.dna".format(mhome))
    sh("{0}/exe/snap/hmm-assembler.pl {1} . > {1}.hmm".format(mhome, species))

    os.chdir(cwd)
    logging.debug("SNAP matrix written to `{0}/{1}.hmm`".format(snapdir, species))
92186
def parallel_fn(f):
    

    def simple_parallel(func, sequence, **args):
        """ f takes an element of sequence as input and the keyword args in **args"""
        if 'processes' in args:
            processes = args.get('processes')
            del args['processes']
        else:
            processes = multiprocessing.cpu_count()

        pool = multiprocessing.Pool(processes)  # depends on available cores

        result = pool.map_async(universal_worker, pool_args(func, sequence, args))
        pool.close()
        pool.join()
        cleaned = [x for x in result.get() if x is not None]  # getting results
        cleaned = asarray(cleaned)
        return cleaned
    return partial(simple_parallel, f)
62004
def _secondary_loop(self, reconnect=False):
        """"""
        if self.num_workers > 1:
            pool = multiprocessing.Pool(self.num_workers)
        else:
            pool = None
        should_reconnect = True
        while should_reconnect:
            i = 0
            running = True
            try:
                self._reset_em()
            except (socket.error, EOFError, IOError, OSError, socket.gaierror, TypeError):
                continue
            while running:
                i += 1
                if i % 5 == 0:
                    # for better performance, only check every 5 cycles
                    try:
                        state = self.em.secondary_state
                    except (socket.error, EOFError, IOError, OSError, socket.gaierror, TypeError):
                        if not reconnect:
                            raise
                        else:
                            break
                    if state == _STATE_FORCED_SHUTDOWN:
                        running = False
                        should_reconnect = False
                    elif state == _STATE_SHUTDOWN:
                        running = False
                    if not running:
                        continue
                try:
                    tasks = self.inqueue.get(block=True, timeout=0.2)
                except queue.Empty:
                    continue
                except (socket.error, EOFError, IOError, OSError, socket.gaierror, TypeError):
                    break
                except (managers.RemoteError, multiprocessing.ProcessError) as e:
                    if ('Empty' in repr(e)) or ('TimeoutError' in repr(e)):
                        continue
                    if (('EOFError' in repr(e)) or ('PipeError' in repr(e)) or
                        ('AuthenticationError' in repr(e))): # Second for Python 3.X, Third for 3.6+
                        break
                    raise
                if pool is None:
                    res = []
                    for genome_id, genome, config in tasks:
                        fitness = self.eval_function(genome, config)
                        res.append((genome_id, fitness))
                else:
                    genome_ids = []
                    jobs = []
                    for genome_id, genome, config in tasks:
                        genome_ids.append(genome_id)
                        jobs.append(
                            pool.apply_async(
                                self.eval_function, (genome, config)
                                )
                            )
                    results = [
                        job.get(timeout=self.worker_timeout) for job in jobs
                        ]
                    res = zip(genome_ids, results)
                try:
                    self.outqueue.put(res)
                except (socket.error, EOFError, IOError, OSError, socket.gaierror, TypeError):
                    break
                except (managers.RemoteError, multiprocessing.ProcessError) as e:
                    if ('Empty' in repr(e)) or ('TimeoutError' in repr(e)):
                        continue
                    if (('EOFError' in repr(e)) or ('PipeError' in repr(e)) or
                        ('AuthenticationError' in repr(e))): # Second for Python 3.X, Third for 3.6+
                        break
                    raise

            if not reconnect:
                should_reconnect = False
                break
        if pool is not None:
            pool.terminate()
503242
def runner(queue, parallel):
    
    def run(fn_name, items):
        logger.info("clusterk: %s" % fn_name)
        assert "wrapper" in parallel, "Clusterk requires bcbio-nextgen-vm wrapper"
        fn = getattr(__import__("{base}.clusterktasks".format(base=parallel["module"]),
                                fromlist=["clusterktasks"]),
                     parallel["wrapper"])
        wrap_parallel = {k: v for k, v in parallel.items() if k in set(["fresources", "pack"])}
        out = []
        for data in [fn(fn_name, queue, parallel.get("wrapper_args"), wrap_parallel, x) for x in items]:
            if data:
                out.extend(data)
        return out
    return run
1151992
def _prep_genome(out_dir, data):
    
    genome_name = utils.splitext_plus(os.path.basename(dd.get_ref_file(data)))[0]
    out_dir = utils.safe_makedir(os.path.join(out_dir, genome_name))
    ref_file = dd.get_ref_file(data)
    utils.symlink_plus(ref_file, os.path.join(out_dir, os.path.basename(ref_file)))
    with open(os.path.join(out_dir, "GenomeSize.xml"), "w") as out_handle:
        out_handle.write('<sequenceSizes genomeName="%s">' % genome_name)
        for c in pysam.AlignmentFile("%s.dict" % utils.splitext_plus(ref_file)[0]).header["SQ"]:
            cur_ploidy = ploidy.get_ploidy([data], region=[c["SN"]])
            out_handle.write('<chromosome fileName="%s" contigName="%s" totalBases="%s" knownBases="%s" '
                             'isCircular="false" ploidy="%s" md5="%s"/>' %
                             (os.path.basename(ref_file), c["SN"], c["LN"], c["LN"], cur_ploidy, c["M5"]))
        out_handle.write('</sequenceSizes>')
    return out_dir
74784
def _pprint_fasta(fasta, annotations=None, annotation_file=None,
                  block_length=10, blocks_per_line=6):
    
    annotations = annotations or []

    # Annotations by chromosome.
    as_by_chrom = collections.defaultdict(lambda: [a for a in annotations] or [])
    if annotation_file:
        for chrom, chrom_iter in _bed_iter(annotation_file):
            as_by_chrom[chrom].append(list(chrom_iter))

    for header, sequence in _fasta_iter(fasta):
        print(header)
        print(pprint_sequence(sequence,
                              annotations=as_by_chrom[header.split()[0]],
                              block_length=block_length,
                              blocks_per_line=blocks_per_line,
                              format=AnsiFormat))
95639
def effective_genome_size(fasta, read_length, nb_cores, tmpdir="/tmp"):
    # type: (str, int, int, str) -> None
    """"""

    idx = Fasta(fasta)

    genome_length = sum([len(c) for c in idx])

    logging.info("Temporary directory: " + tmpdir)
    logging.info("File analyzed: " + fasta)
    logging.info("Genome length: " + str(genome_length))
    print("File analyzed: ", fasta)
    print("Genome length: ", genome_length)

    chromosomes = ", ".join([c.name for c in idx])

    if "_" in chromosomes:
        print("Warning. The following chromosomes are part of your genome:\n",
              chromosomes.replace(">", "") + "\n",
              file=sys.stderr)
        print(
            "You probably want to remove all chromosomes in your fasta containing '_' for the effective genome size computation to be accurate.",
            file=sys.stderr)

    # if tmpdir is None:
    #     try:
    #         tmpdir = os.environ['TMPDIR']
    #     except KeyError:
    #         tmpdir = '/tmp'
    output_file = os.path.join(tmpdir, '{1}.jf'.format(read_length,
                                                       basename(fasta)))
    atexit.register(
        lambda: call("rm {output_file}".format(output_file=output_file), shell=True))

    call(
        "jellyfish count -t {nb_cores} -m {read_length} -s {genome_length} -L 1 -U 1 --out-counter-len 1 --counter-len 1 {fasta} -o {output_file}".format(
            **vars()),
        shell=True)

    stats = check_output("jellyfish stats {output_file}".format(
        output_file=output_file),
                         shell=True)

    unique_kmers = int(stats.split()[1])

    effective_genome_size = unique_kmers / genome_length

    logging.info("Number unique {read_length}-mers: ".format(
        read_length=read_length) + str(unique_kmers))
    logging.info("Effective genome size: " + str(effective_genome_size))
    print("Number unique {read_length}-mers: ".format(read_length=read_length),
          unique_kmers)
    print("Effective genome size: ", effective_genome_size)
    assert effective_genome_size < 1, "Something wrong happened, effective genome size over 1!"
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)