%prog maker maker.gff3 genome.fasta

    Prepare EVM inputs by separating tracks from MAKER.
760066
def maker(args):
    
    from jcvi.formats.base import SetFile, FileShredder

    A, T, P = "ABINITIO_PREDICTION", "TRANSCRIPT", "PROTEIN"
    # Stores default weights and types
    Registry = {\
        "maker": (A, 5),
        "augustus_masked": (A, 1),
        "snap_masked": (A, 1),
        "genemark": (A, 1),
        "est2genome": (T, 5),
        "est_gff": (T, 5),
        "protein2genome": (P, 5),
        "blastx": (P, 1)
    }

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

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

    gffile, fastafile = args

    types = "type.ids"
    if need_update(gffile, types):
        cmd = "cut -f2 -s {0} | sort -u".format(gffile)
        sh(cmd, outfile=types)

    types = SetFile(types)
    reg = defaultdict(list)
    weightsfile = "weights.txt"
    contents = []
    for s in types:
        rs = s.split(":")[0]
        if rs not in Registry:
            continue

        type, weight = Registry[rs]
        reg[type].append(s)
        contents.append("\t".join(str(x) for x in (type, s, weight)))

    contents = "\n".join(sorted(contents))
    write_file(weightsfile, contents)

    evs = [x + ".gff" for x in (A, T, P)]
    FileShredder(evs)

    for type, tracks in reg.items():
        for t in tracks:
            cmd = "grep '\t{0}' {1} | grep -v '_match\t' >> {2}.gff".format(t, gffile, type)
            sh(cmd)

    partition(evs)
    runfile = "run.sh"
    contents = EVMRUN.format(*evs)
    write_file(runfile, contents)
166241
def MakeProto():
  """"""
  # Start running from one directory above the grr directory which is found by
  # this scripts's location as __file__.
  cwd = os.path.dirname(os.path.abspath(__file__))

  # Find all the .proto files.
  protos_to_compile = []
  for (root, _, files) in os.walk(cwd):
    for filename in files:
      full_filename = os.path.join(root, filename)
      if full_filename.endswith(".proto"):
        proto_stat = os.stat(full_filename)

        pb2_path = full_filename.rsplit(".", 1)[0] + "_pb2.py"
        try:
          pb2_stat = os.stat(pb2_path)
          if pb2_stat.st_mtime >= proto_stat.st_mtime:
            continue

        except (OSError, IOError):
          pass

        protos_to_compile.append(full_filename)

  if protos_to_compile:
    # Find the protoc compiler.
    protoc = os.environ.get("PROTOC", "protoc")
    try:
      output = subprocess.check_output([protoc, "--version"])
    except (IOError, OSError):
      raise RuntimeError("Unable to launch %s protoc compiler. Please "
                         "set the PROTOC environment variable.", protoc)

    if b"3.6.1" not in output:
      raise RuntimeError("Incompatible protoc compiler detected. "
                         "We need 3.6.1 not %s" % output)

    for proto in protos_to_compile:
      command = [
          protoc,
          # Write the python files next to the .proto files.
          "--python_out",
          ROOT,
          "--proto_path=%s" % ROOT,
          proto
      ]
      print(
          "Compiling %s with (cwd: %s): %s" % (proto, ROOT, " ".join(command)))
      # The protoc compiler is too dumb to deal with full paths - it expects a
      # relative path from the current working directory.
      subprocess.check_call(command, cwd=ROOT)
193091
def print_generated_sequence(gen, num, *, sep=", ", fmt='', seed=None):
    
    if seed:
       gen.reset(seed)

    elems = [format(next(gen), fmt) for _ in range(num)]
    sep_initial = "\n\n" if '\n' in sep else " "
    print("Generated sequence:{}{}".format(sep_initial, sep.join(elems)))
54099
def _importGenomeObjects(gtfFilePath, chroSet, genome, batchSize, verbose = 0) :
    """"""

    class Store(object) :
        
        def __init__(self, conf) :
            self.conf = conf
            self.chromosomes = {}
            
            self.genes = {}
            self.transcripts = {}
            self.proteins = {}
            self.exons = {}

        def batch_save(self) :
            self.conf.db.beginTransaction()
            
            for c in self.genes.itervalues() :
                c.save()
                conf.removeFromDBRegistery(c)
                
            for c in self.transcripts.itervalues() :
                c.save()
                conf.removeFromDBRegistery(c.exons)
                conf.removeFromDBRegistery(c)
            
            for c in self.proteins.itervalues() :
                c.save()
                conf.removeFromDBRegistery(c)
            
            self.conf.db.endTransaction()
            
            del(self.genes)
            del(self.transcripts)
            del(self.proteins)
            del(self.exons)
            
            self.genes = {}
            self.transcripts = {}
            self.proteins = {}
            self.exons = {}

            gc.collect()

        def save_chros(self) :
            pBar = ProgressBar(nbEpochs = len(self.chromosomes))
            for c in self.chromosomes.itervalues() :
                pBar.update(label = 'Chr %s' % c.number)
                c.save()
            pBar.close()
        
    printf('Importing gene set infos from %s...' % gtfFilePath)
    
    printf('Backuping indexes...')
    indexes = conf.db.getIndexes()
    printf("Droping all your indexes, (don't worry i'll restore them later)...")
    Genome_Raba.flushIndexes()
    Chromosome_Raba.flushIndexes()
    Gene_Raba.flushIndexes()
    Transcript_Raba.flushIndexes()
    Protein_Raba.flushIndexes()
    Exon_Raba.flushIndexes()
    
    printf("Parsing gene set...")
    gtf = GTFFile(gtfFilePath, gziped = True)
    printf('Done. Importation begins!')
    
    store = Store(conf)
    chroNumber = None
    pBar = ProgressBar(nbEpochs = len(gtf))
    for line in gtf :
        chroN = line['seqname']
        pBar.update(label = "Chr %s" % chroN)
        
        if (chroN.upper() in chroSet or chroN.lower() in chroSet):
            strand = line['strand']
            gene_biotype = line['gene_biotype']
            regionType = line['feature']
            frame = line['frame']

            start = int(line['start']) - 1
            end = int(line['end'])
            if start > end :
                start, end = end, start

            chroNumber = chroN.upper()
            if chroNumber not in store.chromosomes :
                store.chromosomes[chroNumber] = Chromosome_Raba()
                store.chromosomes[chroNumber].set(genome = genome, number = chroNumber)
            
            try :
                geneId = line['gene_id']
                geneName =  line['gene_name']
            except KeyError :
                geneId = None
                geneName = None
                if verbose :
                    printf('Warning: no gene_id/name found in line %s' % gtf[i])
            
            if geneId is not None :
                if geneId not in store.genes :
                    if len(store.genes) > batchSize :
                        store.batch_save()
                    
                    if verbose > 0 :
                        printf('\tGene %s, %s...' % (geneId, geneName))
                    store.genes[geneId] = Gene_Raba()
                    store.genes[geneId].set(genome = genome, id = geneId, chromosome = store.chromosomes[chroNumber], name = geneName, strand = strand, biotype = gene_biotype)
                if start < store.genes[geneId].start or store.genes[geneId].start is None :
                        store.genes[geneId].start = start
                if end > store.genes[geneId].end or store.genes[geneId].end is None :
                    store.genes[geneId].end = end
            try :
                transId = line['transcript_id']
                transName = line['transcript_name']
                try :
                    transcript_biotype = line['transcript_biotype']
                except KeyError :
                    transcript_biotype = None
            except KeyError :
                transId = None
                transName = None
                if verbose > 2 :
                    printf('\t\tWarning: no transcript_id, name found in line %s' % gtf[i])
            
            if transId is not None :
                if transId not in store.transcripts :
                    if verbose > 1 :
                        printf('\t\tTranscript %s, %s...' % (transId, transName))
                    store.transcripts[transId] = Transcript_Raba()
                    store.transcripts[transId].set(genome = genome, id = transId, chromosome = store.chromosomes[chroNumber], gene = store.genes.get(geneId, None), name = transName, biotype=transcript_biotype)
                if start < store.transcripts[transId].start or store.transcripts[transId].start is None:
                    store.transcripts[transId].start = start
                if end > store.transcripts[transId].end or store.transcripts[transId].end is None:
                    store.transcripts[transId].end = end
            
                try :
                    protId = line['protein_id']
                except KeyError :
                    protId = None
                    if verbose > 2 :
                        printf('Warning: no protein_id found in line %s' % gtf[i])

                # Store selenocysteine positions in transcript
                if regionType == 'Selenocysteine':
                    store.transcripts[transId].selenocysteine.append(start)
                        
                if protId is not None and protId not in store.proteins :
                    if verbose > 1 :
                        printf('\t\tProtein %s...' % (protId))
                    store.proteins[protId] = Protein_Raba()
                    store.proteins[protId].set(genome = genome, id = protId, chromosome = store.chromosomes[chroNumber], gene = store.genes.get(geneId, None), transcript = store.transcripts.get(transId, None), name = transName)
                    store.transcripts[transId].protein = store.proteins[protId]

                try :
                    exonNumber = int(line['exon_number']) - 1
                    exonKey = (transId, exonNumber)
                except KeyError :
                    exonNumber = None
                    exonKey = None
                    if verbose > 2 :
                        printf('Warning: no exon number or id found in line %s' % gtf[i])
                
                if exonKey is not None :
                    if verbose > 3 :
                        printf('\t\t\texon %s...' % (exonId))
                    
                    if exonKey not in store.exons and regionType == 'exon' :
                        store.exons[exonKey] = Exon_Raba()
                        store.exons[exonKey].set(genome = genome, chromosome = store.chromosomes[chroNumber], gene = store.genes.get(geneId, None), transcript = store.transcripts.get(transId, None), protein = store.proteins.get(protId, None), strand = strand, number = exonNumber, start = start, end = end)
                        store.transcripts[transId].exons.append(store.exons[exonKey])
                    
                    try :
                        store.exons[exonKey].id = line['exon_id']
                    except KeyError :
                        pass
                    
                    if regionType == 'exon' :
                        if start < store.exons[exonKey].start or store.exons[exonKey].start is None:
                            store.exons[exonKey].start = start
                        if end > store.transcripts[transId].end or store.exons[exonKey].end is None:
                            store.exons[exonKey].end = end
                    elif regionType == 'CDS' :
                        store.exons[exonKey].CDS_start = start
                        store.exons[exonKey].CDS_end = end
                        store.exons[exonKey].frame = frame
                    elif regionType == 'stop_codon' :
                        if strand == '+' :
                            if store.exons[exonKey].CDS_end != None :
                                store.exons[exonKey].CDS_end += 3
                                if store.exons[exonKey].end < store.exons[exonKey].CDS_end :
                                    store.exons[exonKey].end = store.exons[exonKey].CDS_end
                                if store.transcripts[transId].end < store.exons[exonKey].CDS_end :
                                    store.transcripts[transId].end = store.exons[exonKey].CDS_end
                                if store.genes[geneId].end < store.exons[exonKey].CDS_end :
                                    store.genes[geneId].end = store.exons[exonKey].CDS_end
                        if strand == '-' :
                            if store.exons[exonKey].CDS_start != None :
                                store.exons[exonKey].CDS_start -= 3
                                if store.exons[exonKey].start > store.exons[exonKey].CDS_start :
                                    store.exons[exonKey].start = store.exons[exonKey].CDS_start
                                if store.transcripts[transId].start > store.exons[exonKey].CDS_start :
                                    store.transcripts[transId].start = store.exons[exonKey].CDS_start
                                if store.genes[geneId].start > store.exons[exonKey].CDS_start :
                                    store.genes[geneId].start = store.exons[exonKey].CDS_start
    pBar.close()
    
    store.batch_save()
    
    conf.db.beginTransaction()
    printf('almost done saving chromosomes...')
    store.save_chros()
    
    printf('saving genome object...')
    genome.save()
    conf.db.endTransaction()
    
    conf.db.beginTransaction()
    printf('restoring core indexes...')
    # Genome.ensureGlobalIndex(('name', 'species'))
    # Chromosome.ensureGlobalIndex('genome')
    # Gene.ensureGlobalIndex('genome')
    # Transcript.ensureGlobalIndex('genome')
    # Protein.ensureGlobalIndex('genome')
    # Exon.ensureGlobalIndex('genome')
    Transcript.ensureGlobalIndex('exons')
    
    printf('commiting changes...')
    conf.db.endTransaction()
    
    conf.db.beginTransaction()
    printf('restoring user indexes')
    pBar = ProgressBar(label = "restoring", nbEpochs = len(indexes))
    for idx in indexes :
        pBar.update()
        conf.db.execute(idx[-1].replace('CREATE INDEX', 'CREATE INDEX IF NOT EXISTS'))
    pBar.close()
    
    printf('commiting changes...')
    conf.db.endTransaction()
    
    return store.chromosomes.values()
1291224
def fixpartials(args):
    
    p = OptionParser(fixpartials.__doc__)
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    gffile, gfasta, partials, = args

    gff = make_index(gffile)
    genome = Fasta(gfasta, index=True)
    partials = LineFile(partials, load=True).lines

    # all_transcripts = [f.id for f in gff.features_of_type("mRNA", \
    #     order_by=("seqid", "start"))]
    seen = set()
    fw = must_open(opts.outfile, "w")
    for gene in gff.features_of_type('gene',  order_by=('seqid', 'start')):
        children = AutoVivification()
        cflag = False
        transcripts = list(gff.children(gene, level=1, order_by=('start')))
        for transcript in transcripts:
            trid, seqid, strand = transcript.id, transcript.seqid, transcript.strand

            for child in gff.children(transcript, order_by=('start')):
                ftype = child.featuretype
                if ftype not in children[trid]: children[trid][ftype] = []
                children[trid][ftype].append(child)

            five_prime, three_prime = True, True
            nstart, nstop = (None, None), (None, None)
            cds_span = [children[trid]['CDS'][0].start, \
                children[trid]['CDS'][-1].stop]
            new_cds_span = [x for x in cds_span]

            start_codon = (cds_span[0], cds_span[0] + 2)
            stop_codon = (cds_span[1] - 2, cds_span[1])
            if strand == '-': start_codon, stop_codon = stop_codon, start_codon

            if trid in partials:
                seen.add(trid)
                start_codon_fasta = _fasta_slice(genome, seqid, \
                    start_codon[0], start_codon[1], strand)
                stop_codon_fasta = _fasta_slice(genome, seqid, \
                    stop_codon[0], stop_codon[1], strand)

                if not is_valid_codon(start_codon_fasta, type='start'):
                    five_prime = False
                    nstart = scan_for_valid_codon(start_codon, strand, \
                        seqid, genome, type='start')

                if not is_valid_codon(stop_codon_fasta, type='stop'):
                    three_prime = False
                    nstop = scan_for_valid_codon(stop_codon, strand, \
                        seqid, genome, type='stop')

                logging.debug("feature={0} ({1})".format(trid, strand) + \
                ", 5'={0}, 3'={1}".format(five_prime, three_prime) + \
                ", {0} <== {1} ==> {2}".format(nstart if strand == '+' else nstop, \
                    cds_span, nstop if strand == '+' else nstart))

            if not five_prime or not three_prime:
                if nstart != (None, None) and (start_codon != nstart):
                    i = 0 if strand == '+' else 1
                    new_cds_span[i] = nstart[i]
                if nstop != (None, None) and (stop_codon != nstop):
                    i = 1 if strand == '+' else 0
                    new_cds_span[i] = nstop[i]
                new_cds_span.sort()

            if set(cds_span) != set(new_cds_span):
                cflag = True
                # if CDS has been extended, appropriately adjust all relevent
                # child feature (CDS, exon, UTR) coordinates
                for ftype in children[trid]:
                    for idx in range(len(children[trid][ftype])):
                        child_span = (children[trid][ftype][idx].start, \
                            children[trid][ftype][idx].stop)
                        if ftype in ('exon', 'CDS'):
                            # if exons/CDSs, adjust start and stop according to
                            # new CDS start and stop, respectively
                            if child_span[0] == cds_span[0]:
                                children[trid][ftype][idx].start = new_cds_span[0]
                            if child_span[1] == cds_span[1]:
                                children[trid][ftype][idx].stop = new_cds_span[1]
                        elif ftype.endswith('UTR'):
                            # if *_prime_UTR, adjust stop according to new CDS start and
                            #                 adjust start according to new CDS stop
                            if child_span[1] == cds_span[0]:
                                children[trid][ftype][idx].stop = new_cds_span[0]
                            if child_span[0] == cds_span[1]:
                                children[trid][ftype][idx].start = new_cds_span[1]

                transcript.start, transcript.stop = \
                    children[trid]['exon'][0].start, children[trid]['exon'][-1].stop

        if cflag:
            _gene_span = range_minmax([(tr.start, tr.stop) for tr in transcripts])
            gene.start, gene.stop = _gene_span[0], _gene_span[1]

        # print gff file
        print(gene, file=fw)
        for transcript in transcripts:
            trid = transcript.id
            print(transcript, file=fw)
            for cftype in children[trid]:
                for child in children[trid][cftype]:
                    print(child, file=fw)

    fw.close()
999395
def prepare(args):
    
    from jcvi.utils.table import write_csv
    from jcvi.formats.base import write_file
    from jcvi.formats.fastq import guessoffset, readlen

    p = OptionParser(prepare.__doc__ + FastqNamings)
    p.add_option("--corr", default=False, action="store_true",
                 help="Extra parameters for corrected data [default: %default]")
    p.add_option("--norun", default=False, action="store_true",
                 help="Don't write `run.sh` script [default: %default]")
    p.add_option("--ploidy", default="2", choices=("1", "2"),
                 help="Ploidy [default: %default]")
    p.set_cpus()
    opts, args = p.parse_args(args)

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

    organism_name = args[0]
    project_name = "".join(x[0] for x in organism_name.split()).upper()
    fnames = sorted(glob("*.fastq*") if len(args) == 1 else args[1:])
    for x in fnames:
        assert op.exists(x), "File `{0}` not found.".format(x)

    groupheader = "group_name library_name file_name".split()
    libheader = "library_name project_name organism_name type paired "\
        "frag_size frag_stddev insert_size insert_stddev read_orientation "\
        "genomic_start genomic_end".split()
    groups_33 = []
    groups_64 = []
    libs = []
    for file_name in fnames:
        offset = guessoffset([file_name])
        group_name = op.basename(file_name).split(".")[0]
        library_name = "-".join(group_name.split("-")[:2])

        # Handle paired files and convert to wildcard
        if ".1." in file_name:
            file_name = file_name.replace(".1.", ".?.")
        elif ".2." in file_name:
            continue

        groupscontents = groups_64 if offset == 64 else groups_33
        groupscontents.append((group_name, library_name, file_name))
        if library_name not in libs:
            libs.append(library_name)

    libcontents = []
    for library_name in libs:
        L = Library(library_name)
        size = L.size
        stddev = L.stddev
        type = L.type
        paired = L.paired
        read_orientation = L.read_orientation

        size = size or ""
        stddev = stddev or ""
        frag_size = size if type == "fragment" else ""
        frag_stddev = stddev if type == "fragment" else ""
        insert_size = size if type != "fragment" else ""
        insert_stddev = stddev if type != "fragment" else ""
        genomic_start, genomic_end = "", ""
        libcontents.append((library_name, project_name, organism_name, type, \
            paired, frag_size, frag_stddev, insert_size, insert_stddev, \
            read_orientation, genomic_start, genomic_end))

    for groups, csvfile in ((groups_33, "in_groups_33.csv"), \
                            (groups_64, "in_groups_64.csv"), \
                            (groups_33 + groups_64, "in_groups.csv")):
        if not groups:
            continue
        write_csv(groupheader, groups, filename=csvfile, tee=True)
        logging.debug("`{0}` created (# of groups = {1}).".\
            format(csvfile, len(groups)))

    write_csv(libheader, libcontents, filename="in_libs.csv", tee=True)
    logging.debug("`in_libs.csv` created (# of libs = {0}).".\
        format(len(libcontents)))

    runfile = "run.sh"

    # ALLPATHS stalls on reads over 250bp <https://www.biostars.org/p/122091/>
    max_rd_len = max(readlen([f]) for f in fnames)
    extra = "CLOSE_UNIPATH_GAPS=False " if max_rd_len > 200 else ""
    if opts.corr:
        extra += "FE_NUM_CYCLES=1 EC_K=28 FE_QUAL_CEIL_RADIUS=0"
        extra += " REMOVE_DODGY_READS_FRAG=False FE_MAX_KMER_FREQ_TO_MARK=1"

    if not opts.norun:
        contents = ALLPATHSRUN.format(opts.ploidy, opts.cpus, extra)
        write_file(runfile, contents)
885499
def demo(args):
    
    p = OptionParser(demo.__doc__)
    opts, args = p.parse_args(args)

    fig = plt.figure(1, (8, 5))
    root = fig.add_axes([0, 0, 1, 1])

    panel_space = .23
    dup_space = .025
    # Draw a gene and two regulatory elements at these arbitrary locations
    locs = [(.5, .9), # ancestral gene
            (.5, .9 - panel_space + dup_space), # identical copies
            (.5, .9 - panel_space - dup_space),
            (.5, .9 - 2 * panel_space + dup_space), # degenerate copies
            (.5, .9 - 2 * panel_space - dup_space),
            (.2, .9 - 3 * panel_space + dup_space), # sub-functionalization
            (.2, .9 - 3 * panel_space - dup_space),
            (.5, .9 - 3 * panel_space + dup_space), # neo-functionalization
            (.5, .9 - 3 * panel_space - dup_space),
            (.8, .9 - 3 * panel_space + dup_space), # non-functionalization
            (.8, .9 - 3 * panel_space - dup_space),
            ]

    default_regulator = "gm"
    regulators = [default_regulator,
            default_regulator, default_regulator,
            "wm", default_regulator,
            "wm", "gw",
            "wb", default_regulator,
            "ww", default_regulator,
            ]

    width = .24
    for i, (xx, yy) in enumerate(locs):
        regulator = regulators[i]
        x1, x2 = xx - .5 * width, xx + .5 * width
        Glyph(root, x1, x2, yy)
        if i == 9:  # upper copy for non-functionalization
            continue

        # coding region
        x1, x2 = xx - .16 * width, xx + .45 * width
        Glyph(root, x1, x2, yy, fc="k")

        # two regulatory elements
        x1, x2 = xx - .4 * width, xx - .28 * width
        for xx, fc in zip((x1, x2), regulator):
            if fc == 'w':
                continue

            DoubleCircle(root, xx, yy, fc=fc)

        rotation = 30
        tip = .02
        if i == 0:
            ya = yy + tip
            root.text(x1, ya, "Flower", rotation=rotation, va="bottom")
            root.text(x2, ya, "Root", rotation=rotation, va="bottom")
        elif i == 7:
            ya = yy + tip
            root.text(x2, ya, "Leaf", rotation=rotation, va="bottom")

    # Draw arrows between panels (center)
    arrow_dist = .08
    ar_xpos = .5
    for ar_ypos in (.3, .53, .76):
        root.annotate(" ", (ar_xpos, ar_ypos),
                (ar_xpos, ar_ypos + arrow_dist),
                arrowprops=arrowprops)

    ar_ypos = .3
    for ar_xpos in (.2, .8):
        root.annotate(" ", (ar_xpos, ar_ypos),
                (.5, ar_ypos + arrow_dist),
                arrowprops=arrowprops)

    # Duplication, Degeneration
    xx = .6
    ys = (.76, .53)
    processes = ("Duplication", "Degeneration")
    for yy, process in zip(ys, processes):
        root.text(xx, yy + .02, process, fontweight="bold")

    # Label of fates
    xs = (.2, .5, .8)
    fates = ("Subfunctionalization", "Neofunctionalization",
            "Nonfunctionalization")
    yy = .05
    for xx, fate in zip(xs, fates):
        RoundLabel(root, xx, yy, fate)

    root.set_xlim(0, 1)
    root.set_ylim(0, 1)
    root.set_axis_off()

    figname = "demo.pdf"
    savefig(figname, dpi=300)
1024169
def setup_jobs(outpath, options, input_files):
    ''' 
    '''

    job_inputs = None
    reverse_mapping = None
    fasta_file_contents = {}

    # Generate FASTA files for PDB inputs
    # fasta_file_contents is a mapping from a file path to a pair (FASTA contents, file type). We remember the file type
    # since we offset residue IDs depending on file type i.e. for FASTA files, we treat each sequence separately and do
    # not renumber the fragments in postprocessing. For PDB files, however, we need to respect the order and length of
    # sequences so that we renumber the fragments appropriately in postprocessing - we assume that if a PDB file is passed in
    # then all chains (protein, RNA, or DNA) will be used in a Rosetta run.
    for input_file in input_files:
        assert(not(fasta_file_contents.get(input_file)))
        if any(fnmatch(input_file, x) for x in pdb_file_wildcards):
            pdb = PDB.from_filepath(input_file, strict=True)
            pdb.pdb_id = os.path.basename(input_file).split('.')[0]            
            if pdb.pdb_id.startswith('pdb') and len(pdb.pdb_id) >= 7:
                # Hack to rename FASTA identifiers for pdb*.ent files which are present in mirrors of the PDB
                pdb.pdb_id = pdb.pdb_id.replace('pdb', '')    
            fasta_file_contents[input_file] = (pdb.create_fasta(prefer_seqres_order = False), 'PDB')
        else:
            fasta_file_contents[input_file] = (read_file(input_file), 'FASTA')

    # Extract sequences from the input FASTA files.
    found_sequences, reverse_mapping, errors = get_sequences(options, fasta_file_contents)
    if found_sequences:
        reformat(found_sequences)
    if errors:
        return None, False, errors

    # Discard sequences that are the wrong chain.
    desired_sequences = {}
    for key, sequence in found_sequences.iteritems():
        pdb_id, chain, file_name = key
        if options.chain is None or chain == options.chain:
            desired_sequences[key] = sequence

    # Create the input FASTA and script files.
    job_inputs, errors = create_inputs(options, outpath, desired_sequences)

    # Create the reverse mapping file
    if reverse_mapping:
        segment_mapping_file = os.path.join(outpath, "segment_map.json")
        colorprinter.message("Creating a reverse mapping file %s." % segment_mapping_file)
        write_file(segment_mapping_file, json.dumps(reverse_mapping))

    # Create the post-processing script file
    post_processing_script = read_file(os.path.join(os.path.split(os.path.realpath(__file__))[0], 'post_processing.py'))
    write_file(os.path.join(outpath, 'post_processing.py'), post_processing_script, 'w')

    # Create the secondary structure filter file
    if options.secondary_structure_file:
        write_file(os.path.join(outpath, 'ss_filter.json'), json.dumps({'secondary_structure_filter' : SecondaryStructureDefinition.from_filepath(options.secondary_structure_file).data}), 'w')

    return job_inputs, reverse_mapping != None, errors
345429
def write_prepareref_fasta_file(outfile, gene_coords, genes_need_upstream, genes_non_upstream, upstream_before=100, upstream_after=100):
    ''''''
    tmp_dict = {}
    fasta_in = os.path.join(data_dir, 'NC_000962.3.fa.gz')
    pyfastaq.tasks.file_to_dict(fasta_in, tmp_dict)
    ref_seq = tmp_dict['NC_000962.3']

    with open(outfile, 'w') as f:
        for gene in genes_non_upstream:
            start = gene_coords[gene]['start']
            end = gene_coords[gene]['end']
            if start < end:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start:end+1])
            else:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[end:start+1])
                gene_fa.revcomp()

            print(gene_fa, file=f)

        for gene in genes_need_upstream:
            start = gene_coords[gene]['start']
            end = gene_coords[gene]['end']
            if start < end:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start - upstream_before:start + upstream_after])
            else:
                gene_fa = pyfastaq.sequences.Fasta(gene, ref_seq[start - upstream_after + 1:start + upstream_before + 1])
                gene_fa.revcomp()

            gene_fa.id += '_upstream'
            print(gene_fa, file=f)
1187584
def infer(args):
    
    from jcvi.apps.grid import WriteJobs
    from jcvi.formats.bed import sort

    p = OptionParser(infer.__doc__)
    p.set_cpus()
    opts, args = p.parse_args(args)

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

    scaffoldsf, genomef = args
    inferbed = "infer-components.bed"
    if need_update((scaffoldsf, genomef), inferbed):
        scaffolds = Fasta(scaffoldsf, lazy=True)
        genome = Fasta(genomef)
        genome = genome.tostring()
        args = [(scaffold_name, scaffold, genome) \
                for scaffold_name, scaffold in scaffolds.iteritems_ordered()]

        pool = WriteJobs(map_one_scaffold, args, inferbed, cpus=opts.cpus)
        pool.run()

    sort([inferbed, "-i"])
    bed = Bed(inferbed)
    inferagpbed = "infer.bed"
    fw = open(inferagpbed, "w")
    seen = []
    for b in bed:
        r = (b.seqid, b.start, b.end)
        if check_seen(r, seen):
            continue
        print("\t".join(str(x) for x in \
                    (b.accn, 0, b.span, b.seqid, b.score, b.strand)), file=fw)
        seen.append(r)
    fw.close()

    frombed([inferagpbed])
655084
def mainRun():
    '''
    
    '''
    #start up the logger
    initLogger()
    
    #attempt to parse the arguments
    p = ap.ArgumentParser(description=util.DESC, formatter_class=ap.RawTextHelpFormatter)
    
    #version data
    p.add_argument('-V', '--version', action='version', version='%(prog)s' + \
                   ' %s in MSBWT %s' % (util.VERSION, util.PKG_VERSION))
    
    #TODO: do we want subparsers groups by type or sorted by name? it's type currently
    
    sp = p.add_subparsers(dest='subparserID')
    p2 = sp.add_parser('cffq', help='create a MSBWT from FASTQ files (pp + cfpp)')
    p2.add_argument('-p', metavar='numProcesses', dest='numProcesses', type=int, default=1, help='number of processes to run (default: 1)')
    p2.add_argument('-u', '--uniform', dest='areUniform', action='store_true', help='the input sequences have uniform length', default=False)
    p2.add_argument('-c', '--compressed', dest='buildCompressed', action='store_true', help='build the RLE BWT (faster, less disk I/O)', default=False)
    p2.add_argument('outBwtDir', type=util.newDirectory, help='the output MSBWT directory')
    p2.add_argument('inputFastqs', nargs='+', type=util.readableFastqFile, help='the input FASTQ files')
    
    p7 = sp.add_parser('pp', help='pre-process FASTQ files before BWT creation')
    p7.add_argument('-u', '--uniform', dest='areUniform', action='store_true', help='the input sequences have uniform length', default=False)
    p7.add_argument('outBwtDir', type=util.newDirectory, help='the output MSBWT directory')
    p7.add_argument('inputFastqs', nargs='+', type=util.readableFastqFile, help='the input FASTQ files')
    
    p3 = sp.add_parser('cfpp', help='create a MSBWT from pre-processed sequences and offsets')
    p3.add_argument('-p', metavar='numProcesses', dest='numProcesses', type=int, default=1, help='number of processes to run (default: 1)')
    p3.add_argument('-u', '--uniform', dest='areUniform', action='store_true', help='the input sequences have uniform length', default=False)
    p3.add_argument('-c', '--compressed', dest='buildCompressed', action='store_true', help='build the RLE BWT (faster, less disk I/O)', default=False)
    p3.add_argument('bwtDir', type=util.existingDirectory, help='the MSBWT directory to process')
    
    p4 = sp.add_parser('merge', help='merge many MSBWTs into a single MSBWT')
    p4.add_argument('-p', metavar='numProcesses', dest='numProcesses', type=int, default=1, help='number of processes to run (default: 1)')
    p4.add_argument('outBwtDir', type=util.newDirectory, help='the output MSBWT directory')
    p4.add_argument('inputBwtDirs', nargs='+', type=util.existingDirectory, help='input BWT directories to merge')
    
    p5 = sp.add_parser('query', help='search for a sequence in an MSBWT, prints sequence and seqID')
    p5.add_argument('inputBwtDir', type=util.existingDirectory, help='the BWT to query')
    p5.add_argument('kmer', type=util.validKmer, help='the input k-mer to search for')
    p5.add_argument('-d', '--dump-seqs', dest='dumpSeqs', action='store_true', help='print all sequences with the given kmer (default=False)', default=False)
    
    p6 = sp.add_parser('massquery', help='search for many sequences in an MSBWT')
    p6.add_argument('inputBwtDir', type=util.existingDirectory, help='the BWT to query')
    p6.add_argument('kmerFile', help='a file with one k-mer per line')
    p6.add_argument('outputFile', help='output file with counts per line')
    p6.add_argument('-r', '--rev-comp', dest='reverseComplement', action='store_true', help='also search for each kmer\'s reverse complement', default=False)
    
    p8 = sp.add_parser('compress', help='compress a MSBWT from byte/base to RLE')
    p8.add_argument('-p', metavar='numProcesses', dest='numProcesses', type=int, default=1, help='number of processes to run (default: 1)')
    p8.add_argument('srcDir', type=util.existingDirectory, help='the source directory for the BWT to compress')
    p8.add_argument('dstDir', type=util.newDirectory, help='the destination directory')
    
    p9 = sp.add_parser('decompress', help='decompress a MSBWT from RLE to byte/base')
    p9.add_argument('-p', metavar='numProcesses', dest='numProcesses', type=int, default=1, help='number of processes to run (default: 1)')
    p9.add_argument('srcDir', type=util.existingDirectory, help='the source directory for the BWT to compress')
    p9.add_argument('dstDir', type=util.newDirectory, help='the destination directory')
    
    p10 = sp.add_parser('convert', help='convert from a raw text input to RLE')
    p10.add_argument('-i', metavar='inputTextFN', dest='inputTextFN', default=None, help='input text filename (default: stdin)')
    p10.add_argument('dstDir', type=util.newDirectory, help='the destination directory')
    
    args = p.parse_args()
    
    if args.subparserID == 'cffq':
        logger.info('Inputs:\t'+str(args.inputFastqs))
        logger.info('Uniform:\t'+str(args.areUniform))
        logger.info('Output:\t'+args.outBwtDir)
        logger.info('Output Compressed:\t'+str(args.buildCompressed))
        logger.info('Processes:\t'+str(args.numProcesses))
        if args.numProcesses > 1:
            logger.warning('Using multi-processing with slow disk accesses can lead to slower build times.')
        print
        if args.areUniform:
            #if they are uniform, use the method developed by Bauer et al., it's likely short Illumina seq
            if args.buildCompressed:
                MultiStringBWT.createMSBWTCompFromFastq(args.inputFastqs, args.outBwtDir, args.numProcesses, args.areUniform, logger)
            else:
                MultiStringBWT.createMSBWTFromFastq(args.inputFastqs, args.outBwtDir, args.numProcesses, args.areUniform, logger)
        else:
            #if they aren't uniform, use the merge method by Holt et al., it's likely longer PacBio seq
            if args.buildCompressed:
                logger.error('No compressed builder for non-uniform datasets, compress after creation.')
            else:
                Multimerge.createMSBWTFromFastq(args.inputFastqs, args.outBwtDir, args.numProcesses, args.areUniform, logger)
        
    elif args.subparserID == 'pp':
        logger.info('Inputs:\t'+str(args.inputFastqs))
        logger.info('Uniform:\t'+str(args.areUniform))
        logger.info('Output:\t'+args.outBwtDir)
        if args.areUniform:
            #preprocess for Bauer et al. method
            MultiStringBWT.preprocessFastqs(args.inputFastqs, args.outBwtDir, args.areUniform, logger)
        else:
            #preprocess for Holt et al. method
            numProcs = 1
            Multimerge.preprocessFastqs(args.inputFastqs, args.outBwtDir, numProcs, args.areUniform, logger)
        
    elif args.subparserID == 'cfpp':
        logger.info('BWT dir:\t'+args.bwtDir)
        logger.info('Uniform:\t'+str(args.areUniform))
        logger.info('Output Compressed:\t'+str(args.buildCompressed))
        logger.info('Processes:\t'+str(args.numProcesses))
        if args.numProcesses > 1:
            logger.warning('Using multi-processing with slow disk accesses can lead to slower build times.')
        print
        seqFN = args.bwtDir+'/seqs.npy'
        offsetFN = args.bwtDir+'/offsets.npy'
        bwtFN = args.bwtDir+'/msbwt.npy'
        
        if args.areUniform:
            #process it using the column-wise Bauer et al. method
            if args.buildCompressed:
                MSBWTCompGenCython.createMsbwtFromSeqs(args.bwtDir, args.numProcesses, logger)
            else:
                MSBWTGenCython.createMsbwtFromSeqs(args.bwtDir, args.numProcesses, logger)
        else:
            #process it using the Holt et al. merge method
            if args.buildCompressed:
                logger.error('No compressed builder for non-uniform datasets, compress after creation.')
            else:
                Multimerge.interleaveLevelMerge(args.bwtDir, args.numProcesses, args.areUniform, logger)
        
    elif args.subparserID == 'compress':
        logger.info('Source Directory:'+args.srcDir)
        logger.info('Dest Directory:'+args.dstDir)
        logger.info('Processes:'+str(args.numProcesses))
        if args.srcDir == args.dstDir:
            raise Exception('Source and destination directories cannot be the same directory.')
        print
        MSBWTGen.compressBWT(args.srcDir+'/msbwt.npy', args.dstDir+'/comp_msbwt.npy', args.numProcesses, logger)
        
    elif args.subparserID == 'decompress':
        logger.info('Source Directory: '+args.srcDir)
        logger.info('Dest Directory: '+args.dstDir)
        logger.info('Processes: '+str(args.numProcesses))
        print
        MSBWTGen.decompressBWT(args.srcDir, args.dstDir, args.numProcesses, logger)
        #TODO: remove if srcdir and dstdir are the same?
        
    elif args.subparserID == 'merge':
        logger.info('Inputs:\t'+str(args.inputBwtDirs))
        logger.info('Output:\t'+args.outBwtDir)
        logger.info('Processes:\t'+str(args.numProcesses))
        if args.numProcesses > 1:
            logger.warning('Multi-processing is not supported at this time, but will be included in a future release.')
            numProcs = 1
            #logger.warning('Using multi-processing with slow disk accesses can lead to slower build times.')
        print
        #MSBWTGen.mergeNewMSBWT(args.outBwtDir, args.inputBwtDirs, args.numProcesses, logger)
        if len(args.inputBwtDirs) > 2:
            #this is a deprecated method, it may still work if you feel daring
            #MSBWTGenCython.mergeMsbwts(args.inputBwtDirs, args.outBwtDir, 1, logger)
            logger.error('Merging more than two MSBWTs at once is not currently supported.')
        else:
            GenericMerge.mergeTwoMSBWTs(args.inputBwtDirs[0], args.inputBwtDirs[1], args.outBwtDir, numProcs, logger)
        
    elif args.subparserID == 'query':
        #this is the easiest thing we can do, don't dump the standard info, just do it
        msbwt = MultiStringBWT.loadBWT(args.inputBwtDir, logger=logger)
        
        #always print how many are found, users can parse it out if they want
        r = msbwt.findIndicesOfStr(args.kmer)
        print r[1]-r[0]
        
        #dump the seqs if request
        if args.dumpSeqs:
            for x in xrange(r[0], r[1]):
                dInd = msbwt.getSequenceDollarID(x)
                print msbwt.recoverString(dInd)[1:]+','+str(dInd)
    
    elif args.subparserID == 'massquery':
        logger.info('Input:\t'+str(args.inputBwtDir))
        logger.info('Queries:\t'+str(args.kmerFile))
        logger.info('Output:\t'+args.outputFile)
        logger.info('Rev-comp:\t'+str(args.reverseComplement))
        print
        msbwt = MultiStringBWT.loadBWT(args.inputBwtDir, logger=logger)
        
        output = open(args.outputFile, 'w+')
        output.write('k-mer,counts')
        if args.reverseComplement:
            output.write(',revCompCounts\n')
        else:
            output.write('\n')
        
        logger.info('Beginning queries...')
        for line in open(args.kmerFile, 'r'):
            kmer = line.strip('\n')
            c = msbwt.countOccurrencesOfSeq(kmer)
            if args.reverseComplement:
                rc = msbwt.countOccurrencesOfSeq(MultiStringBWT.reverseComplement(kmer))
                output.write(kmer+','+str(c)+','+str(rc)+'\n')
            else:
                output.write(kmer+','+str(c)+'\n')
        logger.info('Queries complete.')
    
    elif args.subparserID == 'convert':
        if args.inputTextFN == None:
            logger.info('Input: stdin')
        else:
            logger.info('Input: '+args.inputTextFN)
        logger.info('Output: '+args.dstDir)
        logger.info('Beginning conversion...')
        CompressToRLE.compressInput(args.inputTextFN, args.dstDir)
        logger.info('Finished conversion.')
        
    else:
        print args.subparserID+" is currently not implemented, please wait for a future release."
46355
def fpkm(args):
    
    p = OptionParser(fpkm.__doc__)
    opts, args = p.parse_args(args)

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

    fastafile = args[0]
    bamfiles = args[1:]
    # Create a DUMMY gff file for cuffdiff
    gffile = fastafile.rsplit(".", 1)[0] + ".gff"
    if need_update(fastafile, gffile):
        fw = open(gffile, "w")
        f = Fasta(fastafile, lazy=True)
        for key, size in f.itersizes_ordered():
            print("\t".join(str(x) for x in (key, "dummy", "transcript",\
                1, size, ".", ".", ".", "ID=" + key)), file=fw)
        fw.close()
        logging.debug("Dummy GFF created: {0}".format(gffile))

    cmd = "cuffdiff {0} {1}".format(gffile, " ".join(bamfiles))
    sh(cmd)
240081
def fasta(self):
        
        logging.info('Extracting FASTA sequences matching PointFinder database')
        for sample in self.runmetadata.samples:
            # Ensure that there are sequence data to extract from the GenObject
            if GenObject.isattr(sample[self.analysistype], 'sequences'):
                # Set the name of the FASTA file
                sample[self.analysistype].pointfinderfasta = \
                    os.path.join(sample[self.analysistype].outputdir,
                                 '{seqid}_pointfinder.fasta'.format(seqid=sample.name))
                # Create a list to store all the SeqRecords created
                sequences = list()
                with open(sample[self.analysistype].pointfinderfasta, 'w') as fasta:
                    for gene, sequence in sample[self.analysistype].sequences.items():
                        # Create a SeqRecord using a Seq() of the sequence - both SeqRecord and Seq are from BioPython
                        seq = SeqRecord(seq=Seq(sequence),
                                        id=gene,
                                        name=str(),
                                        description=str())
                        sequences.append(seq)
                    # Write all the SeqRecords to file
                    SeqIO.write(sequences, fasta, 'fasta')
1108010
def vcf_to_df_worker(arg):
    
    canvasvcf, exonbed, i = arg
    logging.debug("Working on job {}: {}".format(i, canvasvcf))
    samplekey = op.basename(canvasvcf).split(".")[0].rsplit('_', 1)[0]
    d = {'SampleKey': samplekey}

    exons = BedTool(exonbed)
    cn = parse_segments(canvasvcf)
    overlaps = exons.intersect(cn, wao=True)
    gcn_store = {}
    for ov in overlaps:
        # Example of ov.fields:
        # [u'chr1', u'11868', u'12227', u'ENSG00000223972.5',
        # u'ENST00000456328.2', u'transcribed_unprocessed_pseudogene',
        # u'DDX11L1', u'.', u'-1', u'-1', u'.', u'0']
        gene_name = "|".join((ov.fields[6], ov.fields[3], ov.fields[5]))
        if gene_name not in gcn_store:
            gcn_store[gene_name] = defaultdict(int)

        cn = ov.fields[-2]
        if cn == ".":
            continue
        cn = int(cn)
        if cn > 10:
            cn = 10
        amt = int(ov.fields[-1])
        gcn_store[gene_name][cn] += amt

    for k, v in sorted(gcn_store.items()):
        v_mean, v_median = counter_mean_and_median(v)
        d[k + ".avgcn"] = v_mean
        d[k + ".medcn"] = v_median
    cleanup()
    return d
174123
function Generator (fn, opts) {
	if (fn instanceof Function) {
		if (typeof opts === 'number') {
			opts = {duration: opts};
		}
		else {
			opts = opts || {};
		}
		opts.generate = fn;
	}
	else {
		opts = fn || {};
	}

	//sort out arguments
	opts = extend({
		//total duration of a stream
		duration: Infinity,

		//time repeat period, in seconds, or 1/frequency
		period: Infinity,

		//inferred from period
		//frequency: 0,

		/**
		 * Generate sample value for a time.
		 * Returns [L, R, ...] or a number for each channel
		 *
		 * @param {number} time current time
		 */
		generate: Math.random
	}, pcm.defaults, opts);

	//align frequency/period
	if (opts.frequency != null) {
		opts.period = 1 / opts.frequency;
	} else {
		opts.frequency = 1 / opts.period;
	}

	let time = 0, count = 0;

	return generate;

	//return sync source/map
	function generate (buffer) {
		if (!buffer) buffer = util.create(opts.samplesPerFrame, opts.channels, opts.sampleRate);

		//get audio buffer channels data in array
		var data = util.data(buffer);

		//enough?
		if (time + buffer.length / opts.sampleRate > opts.duration) return null;

		//generate [channeled] samples
		for (var i = 0; i < buffer.length; i++) {
			var moment = time + i / opts.sampleRate;

			//rotate by period
			if (opts.period !== Infinity) {
				moment %= opts.period;
			}

			var gen = opts.generate(moment);

			//treat null as end
			if (gen === null) {
				return gen;
			}

			//wrap number value
			if (!Array.isArray(gen)) {
				gen = [gen || 0];
			}

			//distribute generated data by channels
			for (var channel = 0; channel < buffer.numberOfChannels; channel++) {
				data[channel][i] = (gen[channel] == null ? gen[0] : gen[channel]);
			}
		}

		//update counters
		count += buffer.length;
		time = count / opts.sampleRate;

		return buffer;
	}
}
1148274
def prepare(args):
    
    valid_enzymes = "ApeKI|ApoI|BamHI|EcoT22I|HinP1I|HpaII|MseI|MspI|" \
                    "NdeI|PasI|PstI|Sau3AI|SbfI|AsiSI-MspI|BssHII-MspI|" \
                    "FseI-MspI|PaeR7I-HhaI|PstI-ApeKI|PstI-EcoT22I|PstI-MspI" \
                    "PstI-TaqI|SalI-MspI|SbfI-MspI".split("|")
    p = OptionParser(prepare.__doc__)
    p.add_option("--enzyme", default="ApeKI", choices=valid_enzymes,
                 help="Restriction enzyme used [default: %default]")
    p.set_home("tassel")
    p.set_aligner(aligner="bwa")
    p.set_cpus()
    opts, args = p.parse_args(args)

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

    barcode, reference = args
    thome = opts.tassel_home
    reference = get_abs_path(reference)
    folders = ("fastq", "tagCounts", "mergedTagCounts", "topm",
               "tbt", "mergedTBT", "hapmap", "hapmap/raw",
               "hapmap/mergedSNPs", "hapmap/filt", "hapmap/bpec")
    for f in folders:
        mkdir(f)

    # Build the pipeline
    runsh = []
    o = "-i fastq -k {0} -e {1} -o tagCounts".format(barcode, opts.enzyme)
    cmd = run_pipeline(thome, "FastqToTagCountPlugin", o)
    runsh.append(cmd)

    o = "-i tagCounts -o mergedTagCounts/myMasterTags.cnt"
    o += " -c 5 -t mergedTagCounts/myMasterTags.cnt.fq"
    cmd = run_pipeline(thome, "MergeMultipleTagCountPlugin", o)
    runsh.append(cmd)
    runsh.append("cd mergedTagCounts")

    cmd = "python -m jcvi.apps.{0} align --cpus {1}".\
                format(opts.aligner, opts.cpus)
    cmd += " {0} myMasterTags.cnt.fq".format(reference)
    runsh.append(cmd)
    runsh.append("cd ..")

    o = "-i mergedTagCounts/*.sam -o topm/myMasterTags.topm"
    cmd = run_pipeline(thome, "SAMConverterPlugin", o)
    runsh.append(cmd)

    o = "-i mergedTBT/myStudy.tbt.byte -y -m topm/myMasterTags.topm"
    o += " -mUpd topm/myMasterTagsWithVariants.topm"
    o += " -o hapmap/raw/myGBSGenos_chr+.hmp.txt"
    o += " -mnF 0.8 -p myPedigreeFile.ped -mnMAF 0.02 -mnMAC 100000"
    o += " -ref {0} -sC 1 -eC 10".format(reference)
    cmd = run_pipeline(thome, "TagsToSNPByAlignmentPlugin", o)
    runsh.append(cmd)

    o = "-hmp hapmap/raw/myGBSGenos_chr+.hmp.txt"
    o += " -o hapmap/mergedSNPs/myGBSGenos_mergedSNPs_chr+.hmp.txt"
    o += " -misMat 0.1 -p myPedigreeFile.ped -callHets -sC 1 -eC 10"
    cmd = run_pipeline(thome, "MergeDuplicateSNPsPlugin", o)
    runsh.append(cmd)

    o = "-hmp hapmap/mergedSNPs/myGBSGenos_mergedSNPs_chr+.hmp.txt"
    o += " -o hapmap/filt/myGBSGenos_mergedSNPsFilt_chr+.hmp.txt"
    o += " -mnTCov 0.01 -mnSCov 0.2 -mnMAF 0.01 -sC 1 -eC 10"
    #o += "-hLD -mnR2 0.2 -mnBonP 0.005"
    cmd = run_pipeline(thome, "GBSHapMapFiltersPlugin", o)
    runsh.append(cmd)

    runfile = "run.sh"
    write_file(runfile, "\n".join(runsh))
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)