%prog stutter a.vcf.gz

    Extract info from lobSTR vcf file. Generates a file that has the following
    fields:

    CHR, POS, MOTIF, RL, ALLREADS, Q
970696
def stutter(args):
    
    p = OptionParser(stutter.__doc__)
    opts, args = p.parse_args(args)

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

    vcf, = args
    pf = op.basename(vcf).split(".")[0]
    execid, sampleid = pf.split("_")

    C = "vcftools --remove-filtered-all --min-meanDP 10"
    C += " --gzvcf {} --out {}".format(vcf, pf)
    C += " --indv {}".format(sampleid)

    info = pf + ".INFO"
    if need_update(vcf, info):
        cmd = C + " --get-INFO MOTIF --get-INFO RL"
        sh(cmd)

    allreads = pf + ".ALLREADS.FORMAT"
    if need_update(vcf, allreads):
        cmd = C + " --extract-FORMAT-info ALLREADS"
        sh(cmd)

    q = pf + ".Q.FORMAT"
    if need_update(vcf, q):
        cmd = C + " --extract-FORMAT-info Q"
        sh(cmd)

    outfile = pf + ".STUTTER"
    if need_update((info, allreads, q), outfile):
        cmd = "cut -f1,2,5,6 {}".format(info)
        cmd += r" | sed -e 's/\t/_/g'"
        cmd += " | paste - {} {}".format(allreads, q)
        cmd += " | cut -f1,4,7"
        sh(cmd, outfile=outfile)
858553
def main(args=None):
    """"""

    if args is None:
        parser = get_argument_parser()
        args = parser.parse_args()

    input_file = args.annotation_file
    output_file = args.output_file
    species = args.species
    chrom_pat = args.chromosome_pattern
    field_name = args.field_name

    log_file = args.log_file
    quiet = args.quiet
    verbose = args.verbose

    # configure root logger
    log_stream = sys.stdout
    if output_file == '-':
        # if we print output to stdout, redirect log messages to stderr
        log_stream = sys.stderr

    logger = misc.get_logger(log_stream = log_stream, log_file = log_file,
            quiet = quiet, verbose = verbose)

    if chrom_pat is None:
        chrom_pat = re.compile(ensembl.SPECIES_CHROMPAT[species])
    else:
        chrom_pat = re.compile(chrom_pat)

    logger.info('Regular expression used for filtering chromosome names: "%s"',
            chrom_pat.pattern)

    chromosomes = set()
    excluded_chromosomes = set()
    i = 0
    exons = 0
    logger.info('Parsing data...')
    if input_file == '-':
        input_file = None
    with misc.smart_open_read(input_file, mode = 'rb', try_gzip = True) as fh, \
            misc.smart_open_write(output_file) as ofh:
        #if i >= 500000: break
        reader = csv.reader(fh, dialect = 'excel-tab')
        writer = csv.writer(ofh, dialect = 'excel-tab', lineterminator = os.linesep,
                quoting = csv.QUOTE_NONE , quotechar = '|')
        for l in reader:
            i += 1
            #if i % int(1e5) == 0:
            #   print '\r%d...' %(i), ; sys.stdout.flush() # report progress
            if len(l) > 1 and l[2] == field_name:
                attr = parse_attributes(l[8])
                type_ = attr['gene_biotype']
                if type_ in ['protein_coding','polymorphic_pseudogene']:

                    # test whether chromosome is valid
                    chrom = l[0]
                    m = chrom_pat.match(chrom)
                    if m is None:
                        excluded_chromosomes.add(chrom)
                        continue

                    chromosomes.add(chrom)
                    writer.writerow(l)
                    exons += 1

    logger.info('Done! (Parsed %d lines.)', i)

    logger.info('')
    logger.info('Gene chromosomes (%d):', len(chromosomes))
    logger.info('\t' + ', '.join(sorted(chromosomes)))
    logger.info('')
    logger.info('Excluded chromosomes (%d):', len(excluded_chromosomes))
    logger.info('\t' + ', '.join(sorted(excluded_chromosomes)))
    logger.info('')
    logger.info('Total no. of exons: %d' %(exons))

    return 0
191326
def cli(variant_file, vep, split):
    
    from datetime import datetime
    from pprint import pprint as pp
    if variant_file == '-':
        my_parser = VCFParser(fsock=sys.stdin, split_variants=split)
    else:
        my_parser = VCFParser(infile = variant_file, split_variants=split)
    start = datetime.now()
    nr_of_variants = 0
    for line in my_parser.metadata.print_header():
        print(line)
    for variant in my_parser:
        pp(variant)
        nr_of_variants += 1
    print('Number of variants: %s' % nr_of_variants)
1245523
def to_vcf(in_file, caller, header_fn, vcf_fn, data, sep="\t"):
    
    out_file = "%s.vcf" % utils.splitext_plus(in_file)[0]
    if not utils.file_exists(out_file + ".gz") and not utils.file_exists(out_file):
        with file_transaction(data, out_file) as tx_out_file:
            with open(in_file) as in_handle:
                with open(tx_out_file, "w") as out_handle:
                    out_handle.write(_vcf_header.format(caller=caller))
                    out_handle.write("\t".join(["#CHROM", "POS", "ID", "REF", "ALT", "QUAL",
                                                "FILTER", "INFO", "FORMAT", dd.get_sample_name(data)]) + "\n")
                    header, in_handle = header_fn(in_handle)
                    for line in in_handle:
                        out = vcf_fn(dict(zip(header, line.strip().split(sep))))
                        if out:
                            out_handle.write("\t".join(out) + "\n")
    out_file = vcfutils.bgzip_and_index(out_file, data["config"])
    effects_vcf, _ = effects.add_to_vcf(out_file, data, "snpeff")
    return effects_vcf or out_file
15529
def split_vcf(in_file, ref_file, config, out_dir=None):
    
    if out_dir is None:
        out_dir = os.path.join(os.path.dirname(in_file), "split")
    out_files = []
    with open(ref.fasta_idx(ref_file, config)) as in_handle:
        for line in in_handle:
            chrom, size = line.split()[:2]
            out_file = os.path.join(out_dir,
                                    os.path.basename(replace_suffix(append_stem(in_file, "-%s" % chrom), ".vcf")))
            subset_vcf(in_file, (chrom, 0, size), out_file, config)
            out_files.append(out_file)
    return out_files
172966
def run(call_file, ref_file, vrn_files, data):
    
    algs = [data["config"]["algorithm"]] * len(data.get("vrn_files", [1]))
    if includes_missingalt(data):
        logger.info("Removing variants with missing alts from %s." % call_file)
        call_file = gatk_remove_missingalt(call_file, data)

    if "gatkcnn" in dd.get_tools_on(data):
        return _cnn_filter(call_file, vrn_files, data)
    elif config_utils.use_vqsr(algs, call_file):
        if vcfutils.is_gvcf_file(call_file):
            raise ValueError("Cannot force gVCF output with joint calling using tools_on: [gvcf] and use VQSR. "
                             "Try using cutoff-based soft filtering with tools_off: [vqsr]")
        snp_file, indel_file = vcfutils.split_snps_indels(call_file, ref_file, data["config"])
        snp_filter_file = _variant_filtration(snp_file, ref_file, vrn_files, data, "SNP",
                                              vfilter.gatk_snp_cutoff)
        indel_filter_file = _variant_filtration(indel_file, ref_file, vrn_files, data, "INDEL",
                                                vfilter.gatk_indel_cutoff)
        orig_files = [snp_filter_file, indel_filter_file]
        out_file = "%scombined.vcf.gz" % os.path.commonprefix(orig_files)
        combined_file = vcfutils.combine_variant_files(orig_files, out_file, ref_file, data["config"])
        return combined_file
    else:
        snp_filter = vfilter.gatk_snp_cutoff(call_file, data)
        indel_filter = vfilter.gatk_indel_cutoff(snp_filter, data)
        return indel_filter
313304
def main(argv, stdout=sys.stdout, stderr=sys.stderr):
    """"""
    # Wrap sys stdout for python 2, so print can understand unicode.
    linesep = os.linesep
    if sys.version_info[0] < 3:
        if stdout == sys.stdout:
            stdout = codecs.getwriter("utf-8")(stdout)
        if stderr == sys.stderr:
            stderr = codecs.getwriter("utf-8")(stderr)
        linesep = unicode(os.linesep)  # pylint: disable=undefined-variable

    arguments = docopt.docopt(
        __doc__, argv=argv[1:], version='git-lint v%s' % __VERSION__)

    json_output = arguments['--json']

    vcs, repository_root = get_vcs_root()

    if vcs is None:
        stderr.write('fatal: Not a git repository' + linesep)
        return 128

    commit = None
    if arguments['--last-commit']:
        commit = vcs.last_commit()

    if arguments['FILENAME']:
        invalid_filenames = find_invalid_filenames(arguments['FILENAME'],
                                                   repository_root)
        if invalid_filenames:
            invalid_filenames.append(('', ''))
            stderr.write(
                linesep.join(invalid[1] for invalid in invalid_filenames))
            return 2

        changed_files = vcs.modified_files(
            repository_root,
            tracked_only=arguments['--tracked'],
            commit=commit)
        modified_files = {}
        for filename in arguments['FILENAME']:
            normalized_filename = os.path.abspath(filename)
            modified_files[normalized_filename] = changed_files.get(
                normalized_filename)
    else:
        modified_files = vcs.modified_files(
            repository_root,
            tracked_only=arguments['--tracked'],
            commit=commit)

    linter_not_found = False
    files_with_problems = 0
    gitlint_config = get_config(repository_root)
    json_result = {}

    with futures.ThreadPoolExecutor(max_workers=multiprocessing.cpu_count())\
            as executor:
        processfile = functools.partial(process_file, vcs, commit,
                                        arguments['--force'], gitlint_config)
        for filename, result in executor.map(
                processfile, [(filename, modified_files[filename])
                              for filename in sorted(modified_files.keys())]):

            rel_filename = os.path.relpath(filename)

            if not json_output:
                stdout.write('Linting file: %s%s' % (termcolor.colored(
                    rel_filename, attrs=('bold', )), linesep))

            output_lines = []
            if result.get('error'):
                output_lines.extend('%s: %s' % (ERROR, reason)
                                    for reason in result.get('error'))
                linter_not_found = True
            if result.get('skipped'):
                output_lines.extend('%s: %s' % (SKIPPED, reason)
                                    for reason in result.get('skipped'))
            if not result.get('comments', []):
                if not output_lines:
                    output_lines.append(OK)
            else:
                files_with_problems += 1
                for data in result['comments']:
                    formatted_message = format_comment(data)
                    output_lines.append(formatted_message)
                    data['formatted_message'] = formatted_message

            if json_output:
                json_result[filename] = result
            else:
                output = linesep.join(output_lines)
                stdout.write(output)
                stdout.write(linesep + linesep)

    if json_output:
        # Hack to convert to unicode, Python3 returns unicode, wheres Python2
        # returns str.
        stdout.write(
            json.dumps(json_result,
                       ensure_ascii=False).encode('utf-8').decode('utf-8'))

    if files_with_problems > 0:
        return 1
    if linter_not_found:
        return 4
    return 0
1111560
def mitocompile(args):
    
    from jcvi.formats.vcf import VcfLine
    from six.moves.urllib.parse import parse_qsl

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

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

    vcfs = args
    print("\t".join("vcf samplekey depth seqid pos alt svlen pe sr".split()))
    for i, vcf in enumerate(vcfs):
        if (i + 1) % 100 == 0:
            logging.debug("Process `{}` [{}]".
                          format(vcf, percentage(i + 1, len(vcfs))))
        depthfile = vcf.replace(".sv.vcf.gz", ".depth")
        fp = must_open(depthfile)
        chrm, depth = fp.next().split()
        depth = int(float(depth))
        samplekey = op.basename(vcf).split("_")[0]

        fp = must_open(vcf)
        for row in fp:
            if row[0] == '#':
                continue
            v = VcfLine(row)
            info = dict(parse_qsl(v.info))
            print("\t".join(str(x) for x in (vcf,
                                             samplekey, depth,
                                             v.seqid, v.pos, v.alt,
                                             info.get("SVLEN"), info["PE"], info["SR"])))
1191582
def _extract_germline(in_file, data):
    
    out_file = "%s-germline.vcf" % utils.splitext_plus(in_file)[0]
    if not utils.file_uptodate(out_file, in_file) and not utils.file_uptodate(out_file + ".gz", in_file):
        with file_transaction(data, out_file) as tx_out_file:
            reader = cyvcf2.VCF(str(in_file))
            reader.add_filter_to_header({'ID': 'Somatic', 'Description': 'Variant called as Somatic'})
            #with contextlib.closing(cyvcf2.Writer(tx_out_file, reader)) as writer:
            with open(tx_out_file, "w") as out_handle:
                out_handle.write(reader.raw_header)
                for rec in reader:
                    rec = _update_germline_filters(rec)
                    out_handle.write(str(rec))
                    #writer.write_record(rec)
    return out_file
843632
def _print_header(data):
    
    print("##fileformat=VCFv4.2", file=STDOUT, end="")
    print("##source=seqbuster2.3", file=STDOUT, end="")
    print("##reference=mirbase", file=STDOUT, end="")
    for pos in data:
        print("##contig=<ID=%s>" % pos["chrom"], file=STDOUT, end="")
    print('##INFO=<ID=ID,Number=1,Type=String,Description="miRNA name">', file=STDOUT, end="")
    print('##FORMAT=<ID=GT,Number=1,Type=Integer,Description="Genotype">', file=STDOUT, end="")
    print('##FORMAT=<ID=NR,Number=A,Type=Integer,Description="Total reads supporting the variant">', file=STDOUT, end="")
    print('##FORMAT=<ID=NS,Number=A,Type=Float,Description="Total number of different sequences supporting the variant">', file=STDOUT, end="")
    print("#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tSAMP001", file=STDOUT, end="")
670706
def vcf(args):
    
    from jcvi.apps.grid import Jobs

    valid_callers = ("mpileup", "freebayes")
    p = OptionParser(vcf.__doc__)
    p.set_outfile(outfile="out.vcf.gz")
    p.add_option("--nosort", default=False, action="store_true",
                 help="Do not sort the BAM files")
    p.add_option("--caller", default="mpileup", choices=valid_callers,
                 help="Use variant caller [default: %default]")
    opts, args = p.parse_args(args)

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

    fastafile = args[0]
    bamfiles = args[1:]
    caller = opts.caller

    unsorted = [x for x in bamfiles if ".sorted." not in x]
    if opts.nosort:
        bamfiles = unsorted
    else:
        jargs = [[[x, "--unique"]] for x in unsorted]
        jobs = Jobs(index, args=jargs)
        jobs.run()
        bamfiles = [x.replace(".sorted.bam", ".bam") for x in bamfiles]
        bamfiles = [x.replace(".bam", ".sorted.bam") for x in bamfiles]

    if caller == "mpileup":
        cmd = "samtools mpileup -E -uf"
        cmd += " {0} {1}".format(fastafile, " ".join(bamfiles))
        cmd += " | bcftools call -vmO v"
    elif caller == "freebayes":
        cmd = "freebayes -f"
        cmd += " {0} {1}".format(fastafile, " ".join(bamfiles))
    sh(cmd, outfile=opts.outfile)
984738
def main():
    """"""

    parser = argparse.ArgumentParser(description='Extract parameters from an ULog file')
    parser.add_argument('filename', metavar='file.ulg', help='ULog input file')

    parser.add_argument('-d', '--delimiter', dest='delimiter', action='store',
                        help='Use delimiter in CSV (default is \',\')', default=',')

    parser.add_argument('-i', '--initial', dest='initial', action='store_true',
                        help='Only extract initial parameters', default=False)

    parser.add_argument('-o', '--octave', dest='octave', action='store_true',
                        help='Use Octave format', default=False)

    parser.add_argument('-t', '--timestamps', dest='timestamps', action='store_true',
                        help='Extract changed parameters with timestamps', default=False)

    parser.add_argument('output_filename', metavar='params.txt',
                        type=argparse.FileType('w'), nargs='?',
                        help='Output filename (default=stdout)', default=sys.stdout)

    args = parser.parse_args()
    ulog_file_name = args.filename

    message_filter = []
    if not args.initial: message_filter = None

    ulog = ULog(ulog_file_name, message_filter)

    param_keys = sorted(ulog.initial_parameters.keys())
    delimiter = args.delimiter
    output_file = args.output_filename

    if not args.octave:
        for param_key in param_keys:
            output_file.write(param_key)
            if args.initial:
                output_file.write(delimiter)
                output_file.write(str(ulog.initial_parameters[param_key]))
                output_file.write('\n')
            elif args.timestamps:
                output_file.write(delimiter)
                output_file.write(str(ulog.initial_parameters[param_key]))
                for t, name, value in ulog.changed_parameters:
                    if name == param_key:
                        output_file.write(delimiter)
                        output_file.write(str(value))

                output_file.write('\n')
                output_file.write("timestamp")
                output_file.write(delimiter)
                output_file.write('0')
                for t, name, value in ulog.changed_parameters:
                    if name == param_key:
                        output_file.write(delimiter)
                        output_file.write(str(t))

                output_file.write('\n')
            else:
                for t, name, value in ulog.changed_parameters:
                    if name == param_key:
                        output_file.write(delimiter)
                        output_file.write(str(value))

                output_file.write('\n')

    else:

        for param_key in param_keys:
            output_file.write('# name ')
            output_file.write(param_key)
            values = [ulog.initial_parameters[param_key]]

            if not args.initial:
                for t, name, value in ulog.changed_parameters:
                    if name == param_key:
                        values += [value]

            if len(values) > 1:
                output_file.write('\n# type: matrix\n')
                output_file.write('# rows: 1\n')
                output_file.write('# columns: ')
                output_file.write(str(len(values)) + '\n')
                for value in values:
                    output_file.write(str(value) + ' ')

            else:
                output_file.write('\n# type: scalar\n')
                output_file.write(str(values[0]))

            output_file.write('\n')
865809
def pasa(args):
    
    from jcvi.formats.base import SetFile
    from jcvi.formats.gff import Gff

    p = OptionParser(pasa.__doc__)
    p.set_home("pasa")
    opts, args = p.parse_args(args)

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

    fastafile, gffile = args
    transcodergff = fastafile + ".transdecoder.gff3"
    transcodergenomegff = fastafile + ".transdecoder.genome.gff3"
    if need_update((fastafile, gffile), (transcodergff, transcodergenomegff)):
        cmd = "{0}/scripts/pasa_asmbls_to_training_set.dbi".format(opts.pasa_home)
        cmd += " --pasa_transcripts_fasta {0} --pasa_transcripts_gff3 {1}".\
                format(fastafile, gffile)
        sh(cmd)

    completeids = fastafile.rsplit(".", 1)[0] + ".complete.ids"
    if need_update(transcodergff, completeids):
        cmd = "grep complete {0} | cut -f1 | sort -u".format(transcodergff)
        sh(cmd, outfile=completeids)

    complete = SetFile(completeids)
    seen = set()
    completegff = transcodergenomegff.rsplit(".", 1)[0] + ".complete.gff3"
    fw = open(completegff, "w")
    gff = Gff(transcodergenomegff)
    for g in gff:
        a = g.attributes
        if "Parent" in a:
            id = a["Parent"][0]
        else:
            id = a["ID"][0]
        asmbl_id = id.split("|")[0]
        if asmbl_id not in complete:
            continue
        print(g, file=fw)
        if g.type == "gene":
            seen.add(id)

    fw.close()
    logging.debug("A total of {0} complete models extracted to `{1}`.".\
                    format(len(seen), completegff))
1327489
def info(args):
    
    p = OptionParser(info.__doc__)
    p.add_option("--gaps", default=False, action="store_true",
                 help="Count number of gaps [default: %default]")
    p.set_table()
    p.set_outfile()
    opts, args = p.parse_args(args)

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

    fastafiles = args
    data = []
    for f in fastafiles:
        s = SequenceInfo(f, gapstats=opts.gaps)
        data.append(s.data)
    write_csv(s.header, data, sep=opts.sep,
              filename=opts.outfile, align=opts.align)
801315
def _get_parser(description):
    """"""
    parser = argparse.ArgumentParser(description=description)
    parser.add_argument('key', help="Camellia key.")
    parser.add_argument('input_file', nargs='*',
                        help="File(s) to read as input data. If none are "
                        "provided, assume STDIN.")
    parser.add_argument('-o', '--output_file',
                        help="Output file. If not provided, assume STDOUT.")
    parser.add_argument('-l', '--keylen', type=int, default=128,
                        help="Length of 'key' in bits, must be in one of %s "
                        "(default 128)." % camcrypt.ACCEPTABLE_KEY_LENGTHS)
    parser.add_argument('-H', '--hexkey', action='store_true',
                        help="Treat 'key' as a hex string rather than binary.")

    return parser
419160
def run_pyvcf(args):
    """"""
    # open VCF reader
    reader = vcf.Reader(filename=args.input_vcf)
    # optionally, open VCF writer
    writer = None
    # read through input VCF file, optionally also writing out
    start = time.clock()
    num = 0
    for num, r in enumerate(reader):
        if num % 10000 == 0:
            print(num, "".join(map(str, [r.CHROM, ":", r.POS])), sep="\t", file=sys.stderr)
        if writer:
            writer.write_record(r)
        if args.max_records and num >= args.max_records:
            break
    end = time.clock()
    print("Read {} records in {} seconds".format(num, (end - start)), file=sys.stderr)
679855
def _seg_to_vcf(vals):
    
    call_to_cn = {"+": 3, "-": 1}
    call_to_type = {"+": "DUP", "-": "DEL"}
    if vals["CALL"] not in ["0"]:
        info = ["FOLD_CHANGE_LOG=%s" % vals["MEAN_LOG2_COPY_RATIO"],
                "PROBES=%s" % vals["NUM_POINTS_COPY_RATIO"],
                "SVTYPE=%s" % call_to_type[vals["CALL"]],
                "SVLEN=%s" % (int(vals["END"]) - int(vals["START"])),
                "END=%s" % vals["END"],
                "CN=%s" % call_to_cn[vals["CALL"]]]
        return [vals["CONTIG"], vals["START"], ".", "N", "<%s>" % call_to_type[vals["CALL"]], ".",
                ".", ";".join(info), "GT", "0/1"]