%prog prepare genomesize *.fastq

    Prepare MERACULOUS configuation file. Genome size should be entered in Mb.
510374
def prepare(args):
    
    p = OptionParser(prepare.__doc__ + FastqNamings)
    p.add_option("-K", default=51, type="int", help="K-mer size")
    p.set_cpus(cpus=32)
    opts, args = p.parse_args(args)

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

    genomesize = float(args[0]) / 1000
    fnames = args[1:]
    for x in fnames:
        assert op.exists(x), "File `{0}` not found.".format(x)

    s = comment_banner("Meraculous params file") + "\n"
    s += comment_banner("Basic parameters") + "\n"
    s += "# Describe the libraries ( one line per library )\n"
    s += "# " + " ".join(header.split()) + "\n"

    libs = get_libs(fnames)
    lib_seqs = []
    rank = 0
    for lib, fs in libs:
        size = lib.size
        if size == 0:
            continue
        rank += 1
        library_name = lib.library_name
        name = library_name.replace("-", "")
        wildcard = "{0}*.1.*,{0}*.2.*".format(library_name)
        rl = max(readlen([x]) for x in fs)
        lib_seq = lib.get_lib_seq(wildcard, name, rl, rank)
        lib_seqs.append(lib_seq)

    s += "\n" + "\n".join(load_csv(None, lib_seqs, sep=" ")) + "\n"
    params = [("genome_size", genomesize),
              ("is_diploid", 0),
              ("mer_size", opts.K),
              ("num_prefix_blocks", 1),
              ("no_read_validation", 0),
              ("local_num_procs", opts.cpus)]
    s += "\n" + "\n".join(load_csv(None, params, sep=" ")) + "\n"

    cfgfile = "meraculous.config"
    write_file(cfgfile, s, tee=True)

    s = "~/export/meraculous/bin/run_meraculous.sh -c {0}"\
                .format(cfgfile)
    runsh = "run.sh"
    write_file(runsh, s)
886330
def fasta(args):
    
    from jcvi.formats.sizes import Sizes

    p = OptionParser(fasta.__doc__)
    p.add_option("--extend", default=1000, type="int",
                 help="Extend seq flanking the gaps [default: %default]")
    opts, args = p.parse_args(args)

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

    mapout, sfasta = args
    Flank = opts.extend
    pf = mapout.split(".")[0]
    mapbed = pf + ".bed"
    bm = BinMap(mapout)
    bm.print_to_bed(mapbed)

    bed = Bed(mapbed, sorted=False)
    markersbed = pf + ".markers.bed"
    fw = open(markersbed, "w")
    sizes = Sizes(sfasta).mapping
    for b in bed:
        accn = b.accn
        scf, pos = accn.split(".")
        pos = int(pos)
        start = max(0, pos - Flank)
        end = min(pos + Flank, sizes[scf])
        print("\t".join(str(x) for x in \
                    (scf, start, end, accn)), file=fw)

    fw.close()

    fastaFromBed(markersbed, sfasta, name=True)
39877
def sequence_prep(self):
        
        # Create a sorted list of all the FASTA files in the sequence path
        strains = sorted(glob(os.path.join(self.fastapath, '*.fa*'.format(self.fastapath))))
        for sample in strains:
            # Create the object
            metadata = MetadataObject()
            # Set the sample name to be the file name of the sequence by removing the path and file extension
            sample_name = os.path.splitext(os.path.basename(sample))[0]
            if sample_name in self.strainset:
                # Extract the OLNID from the dictionary using the SEQID
                samplename = self.straindict[sample_name]
                # samplename = sample_name
                # Set and create the output directory
                outputdir = os.path.join(self.path, samplename)
                make_path(outputdir)
                # Set the name of the JSON file
                json_metadata = os.path.join(outputdir, '{name}.json'.format(name=samplename))
                if not os.path.isfile(json_metadata):
                    # Create the name and output directory attributes
                    metadata.name = samplename
                    metadata.seqid = sample_name
                    metadata.outputdir = outputdir
                    metadata.jsonfile = json_metadata
                    # Set the name of the FASTA file to use in the analyses
                    metadata.bestassemblyfile = os.path.join(metadata.outputdir,
                                                             '{name}.fasta'.format(name=metadata.name))
                    # Symlink the original file to the output directory
                    relative_symlink(sample, outputdir, '{sn}.fasta'.format(sn=metadata.name))
                    # Associate the corresponding FASTQ files with the assembly
                    metadata.fastqfiles = sorted(glob(os.path.join(self.fastqpath,
                                                                   '{name}*.gz'.format(name=metadata.name))))
                    metadata.forward_fastq, metadata.reverse_fastq = metadata.fastqfiles
                    # Write the object to file
                    self.write_json(metadata)
                else:
                    metadata = self.read_json(json_metadata)
                # Add the metadata object to the list of objects
                self.metadata.append(metadata)
1251410
def generate_sizes(name, genome_dir):
    """"""
    fa = os.path.join(genome_dir, name, "{}.fa".format(name))
    sizes = fa + ".sizes"
    g = Fasta(fa)
    with open(sizes, "w") as f:
        for seqname in g.keys():
            f.write("{}\t{}\n".format(seqname, len(g[seqname])))
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)))
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)
732217
def _seqtk_fastq_prep_cl(data, in_file=None, read_num=0):
    
    needs_convert = dd.get_quality_format(data).lower() == "illumina"
    trim_ends = dd.get_trim_ends(data)
    seqtk = config_utils.get_program("seqtk", data["config"])
    if in_file:
        in_file = objectstore.cl_input(in_file)
    else:
        in_file = "/dev/stdin"
    cmd = ""
    if needs_convert:
        cmd += "{seqtk} seq -Q64 -V {in_file}".format(**locals())
    if trim_ends:
        left_trim, right_trim = trim_ends[0:2] if data.get("read_num", read_num) == 0 else trim_ends[2:4]
        if left_trim or right_trim:
            trim_infile = "/dev/stdin" if needs_convert else in_file
            pipe = " | " if needs_convert else ""
            cmd += "{pipe}{seqtk} trimfq -b {left_trim} -e {right_trim} {trim_infile}".format(**locals())
    return cmd
1004431
def _prep_inputs(vrn_info, cnv_info, somatic_info, work_dir, config):
    
    exe = os.path.join(os.path.dirname(sys.executable), "create_phylowgs_inputs.py")
    assert os.path.exists(exe), "Could not find input prep script for PhyloWGS runs."
    ssm_file = os.path.join(work_dir, "ssm_data.txt")
    cnv_file = os.path.join(work_dir, "cnv_data.txt")
    if not utils.file_exists(ssm_file) or not utils.file_exists(cnv_file):
        with file_transaction(somatic_info.tumor_data, ssm_file, cnv_file) as (tx_ssm_file, tx_cnv_file):
            variant_type, input_vcf_file = _prep_vrn_file(vrn_info["vrn_file"], vrn_info["variantcaller"],
                                                          work_dir, somatic_info, cnv_info["ignore"], config)
            input_cnv_file = _prep_cnv_file(cnv_info["subclones"], work_dir, somatic_info)
            cmd = [sys.executable, exe,
                   "--sample-size", str(config["sample_size"]), "--tumor-sample", somatic_info.tumor_name,
                   "--battenberg", input_cnv_file, "--cellularity", _read_contam(cnv_info["contamination"]),
                   "--output-cnvs", tx_cnv_file, "--output-variants", tx_ssm_file,
                   "--variant-type", variant_type, input_vcf_file]
            do.run(cmd, "Prepare PhyloWGS inputs.")
    return ssm_file, cnv_file
888641
def pasteprepare(args):
    
    p = OptionParser(pasteprepare.__doc__)
    p.add_option("--flank", default=5000, type="int",
                 help="Get the seq of size on two ends [default: %default]")
    opts, args = p.parse_args(args)

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

    goodfasta, = args
    flank = opts.flank
    pf = goodfasta.rsplit(".", 1)[0]
    extbed = pf + ".ext.bed"

    sizes = Sizes(goodfasta)
    fw = open(extbed, "w")
    for bac, size in sizes.iter_sizes():
        print("\t".join(str(x) for x in \
                               (bac, 0, min(flank, size), bac + "L")), file=fw)
        print("\t".join(str(x) for x in \
                               (bac, max(size - flank, 0), size, bac + "R")), file=fw)
    fw.close()

    fastaFromBed(extbed, goodfasta, name=True)
1008397
def generate_combined_fasta(self, genome_list, genome_dir):
        '''
        '''

        fasta = []
            
        for genome in genome_list:
            
            full_path = genome_dir + genome
            handle = open(full_path, "rU")
            print 'making combined fasta for', genome
            try:
                seq_record = SeqIO.read(handle, 'genbank')
                org_accession = seq_record.name
                
            except AssertionError,e:
                print str(e), genome
                        
            for feature in seq_record.features:
                if feature.type == 'CDS':
                    try:                                  
                        prot_accession = feature.qualifiers['protein_id'][0]
                        prot_translation = feature.qualifiers['translation'][0]
                        
                        newfast = '>' + org_accession + ',' + prot_accession + \
                         '\n' + prot_translation + '\n'
                        
                        #if newfast not in fasta:
                        fasta.append(newfast)
                            
                    except AttributeError,e:
                        print "organism %s, protein %s did not have \
                        the right attributes" % (org_accession, prot_accession)
                        print str(e)
                    
                    except KeyError,e:
                        print "organism %s, protein %s did not have \
                        the right key" % (org_accession, prot_accession)
                        print str(e)
                          
            handle.close()
            
            print "%s proteins were added" % len(fasta)

            set_fasta = set(fasta)

            print "%s unique proteins were added -- dropping redundant ones" % len(set_fasta)
                                                
        faastring = "".join(set_fasta)

        write_fasta = open('combined_fasta', 'w')
        write_fasta.write(faastring)
        write_fasta.close()

        return set_fasta
1311084
def GC_wide(genome, window=1000):
    

    GC = []
    from Bio import SeqIO
    with open(genome) as handle:
        sequence = "".join([str(record.seq)
                            for record in SeqIO.parse(handle, "fasta")])

    n = len(sequence)
    for i in range(0, n, window):
        portion = sequence[i:min(i + window, n)]
        GC.append(GC_partial(portion))

    return GC
103326
def gc(args):
    
    p = OptionParser(gc.__doc__)
    p.add_option("--binsize", default=500, type="int",
                 help="Bin size to use")
    opts, args = p.parse_args(args)

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

    fastafile, = args
    binsize = opts.binsize
    allbins = []
    for name, seq in parse_fasta(fastafile):
        for i in range(len(seq) / binsize):
            atcnt = gccnt = 0
            for c in seq[i * binsize: (i + 1) * binsize].upper():
                if c in "AT":
                    atcnt += 1
                elif c in "GC":
                    gccnt += 1
            totalcnt = atcnt + gccnt
            if totalcnt == 0:
                continue
            gcpct = gccnt * 100 / totalcnt
            allbins.append(gcpct)

    from jcvi.graphics.base import asciiplot
    from collections import Counter

    title = "Total number of bins={}".format(len(allbins))
    c = Counter(allbins)
    x, y = zip(*sorted(c.items()))
    asciiplot(x, y, title=title)
1228563
def size(args):
    
    from jcvi.utils.cbook import human_size

    p = OptionParser(size.__doc__)
    fp = sys.stdin

    results = []
    for link_name in fp:
        link_name = link_name.strip()
        if not op.islink(link_name):
            continue

        source = get_abs_path(link_name)

        link_name = op.basename(link_name)
        filesize = op.getsize(source)
        results.append((filesize, link_name))

    # sort by descending file size
    for filesize, link_name in sorted(results, reverse=True):
        filesize = human_size(filesize, a_kilobyte_is_1024_bytes=True)
        print("%10s\t%s" % (filesize, link_name), file=sys.stderr)
723772
def prepare_input_data(config):
    

    if not dd.get_disambiguate(config):
        return dd.get_input_sequence_files(config)

    work_bam = dd.get_work_bam(config)
    logger.info("Converting disambiguated reads to fastq...")
    fq_files = convert_bam_to_fastq(
        work_bam, dd.get_work_dir(config), None, None, config
    )
    return fq_files
508546
def _add_umis_with_fastp(read_fq, umi_fq, out_fq, cores):
    
    with utils.open_gzipsafe(umi_fq) as in_handle:
        in_handle.readline()  # name
        umi_size = len(in_handle.readline().strip())
    cmd = ("fastp -Q -A -L -G -w 1 --in1 {read_fq} --in2 {umi_fq} "
           "--umi --umi_prefix UMI --umi_loc read2 --umi_len {umi_size} "
           "--out1 >(bgzip --threads {cores} -c > {out_fq}) --out2 /dev/null "
           "-j /dev/null -h /dev/null")
    do.run(cmd.format(**locals()), "Add UMIs to fastq file with fastp")
518088
def importGenome(name, batchSize = 100) :
	""""""
	path = os.path.join(this_dir, "bootstrap_data", "genomes/" + name)
	PG.importGenome(path, batchSize)
195672
def run_mrtrans(align_fasta, recs, work_dir, outfmt="paml"):
    
    align_file = op.join(work_dir, "prot-align.fasta")
    nuc_file = op.join(work_dir, "nuc.fasta")
    output_file = op.join(work_dir, "nuc-align.mrtrans")

    # make the prot_align file and nucleotide file
    align_h0 = open(align_file + "0", "w")
    align_h0.write(str(align_fasta))
    align_h0.close()
    prot_seqs = {}
    i = 0
    for rec in SeqIO.parse(align_h0.name, "fasta"):
        prot_seqs[i] = rec.seq
        i += 1
    align_h = open(align_file, "w")
    for i, rec in enumerate(recs):
        if len(rec.id) > 30:
            rec.id = rec.id[:28] + "_" + str(i)
            rec.description = ""
        print(">{0}\n{1}".format(rec.id, prot_seqs[i]), file=align_h)
    align_h.close()
    SeqIO.write(recs, file(nuc_file, "w"), "fasta")

    # run the program
    cl = MrTransCommandline(align_file, nuc_file, output_file, outfmt=outfmt)
    r, e = cl.run()
    if e is None:
        print("\tpal2nal:", cl, file=sys.stderr)
        return output_file
    elif e.read().find("could not translate") >= 0:
        print("***pal2nal could not translate", file=sys.stderr)
        return None