Split up GWAS jobs based on portions of a chromosome

    :param args: arguments from parseargs
    :param chrom_file: marker info file
    :return dictionary name=>job_details
544169
def split_chrom_jobs(args, chrom_file):
    
    max_snp_count = args.snps_per_job

    poscol = 3
    if args.map3:
        poscol = 2

    job_list = {}
    cur = None
    last_pos = None
    job_string = ""
    job_name = ""

    for line in sys_call("cut -f 1,%d %s" % (poscol, chrom_file)):
        pos = -1

        values = line.split()
        if len(values) > 0:
            chrom, pos = [int(x) for x in values]

        if cur is None:     # First line observed
            cur = chrom
            job_string = "--chr %d --from-bp %d" % (chrom, pos)
            job_name = "Chr%d_%d-" % (chrom, pos)
            snp_count = 0
        elif cur != cur:    # Changed chromosome
            job_string += " --to-bp %d" % (last_pos)
            job_name += str(last_pos)
            job_list[job_name] = job_string
            cur = chrom
            job_string = "--chr %d --from-bp %d" % (chrom, pos)
            job_name = "Chr%d_%d-" % (chrom, pos)
            snp_count = 0
                            # create new job based on snp count
        elif snp_count < max_snp_count:
            snp_count += 1
        else:
            job_string += " --to-bp %d" % (last_pos)
            job_name += str(last_pos)
            job_list[job_name] = job_string
            job_string = "--chr %d --from-bp" % (chrom, pos)
            job_name = "Chr%d_%d-" % (chrom, pos)
            snp_count = 0

        last_pos = pos
    if job_string != "":
        job_string += " --to-bp %d" % (last_pos)
        job_name += str(last_pos)
        job_list[job_name] = job_string

    return job_list
3608
def _target_chroms_and_header(bam_file, data):
    
    special_remaps = {"chrM": "MT", "MT": "chrM"}
    target_chroms = dict([(x.name, i) for i, x in enumerate(ref.file_contigs(dd.get_ref_file(data)))
                          if chromhacks.is_autosomal_or_sex(x.name)])
    out_chroms = []
    with pysam.Samfile(bam_file, "rb") as bamfile:
        for bami, bam_contig in enumerate([c["SN"] for c in bamfile.header["SQ"]]):
            if bam_contig in target_chroms:
                target_chrom = bam_contig
            elif bam_contig in special_remaps and special_remaps[bam_contig] in target_chroms:
                target_chrom = special_remaps[bam_contig]
            elif bam_contig.startswith("chr") and bam_contig.replace("chr", "") in target_chroms:
                target_chrom = bam_contig.replace("chr", "")
            elif "chr%s" % bam_contig in target_chroms:
                target_chrom = "chr%s" % bam_contig
            else:
                target_chrom = None
            # target_chrom == bam_contig ensures we don't try chr1 -> 1 style remapping
            if target_chrom and target_chrom == bam_contig:
                # Order not required if dealing with SAM file header fixing
                #assert bami == target_chroms[target_chrom], \
                #    ("remove_extracontigs: Non-matching order of standard contig: %s %s (%s vs %s)" %
                #     (bam_file, target_chrom, bami, target_chroms[target_chrom]))
                out_chroms.append(target_chrom)
    assert out_chroms, ("remove_extracontigs: Did not find any chromosomes in reference file: %s %s" %
                        (bam_file, target_chroms))
    return out_chroms
30382
def summary(args):
    
    p = OptionParser(summary.__doc__)
    opts, args = p.parse_args(args)

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

    gff_file, ref = args
    s = Fasta(ref)
    g = make_index(gff_file)
    geneseqs, exonseqs, intronseqs = [], [], []  # Calc % GC
    for f in g.features_of_type("gene"):
        fid = f.id
        fseq = s.sequence({'chr': f.chrom, 'start': f.start, 'stop': f.stop})
        geneseqs.append(fseq)
        exons = set((c.chrom, c.start, c.stop) for c in g.children(fid, 2) \
                    if c.featuretype == "exon")
        exons = list(exons)
        for chrom, start, stop in exons:
            fseq = s.sequence({'chr': chrom, 'start': start, 'stop': stop})
            exonseqs.append(fseq)
        introns = range_interleave(exons)
        for chrom, start, stop in introns:
            fseq = s.sequence({'chr': chrom, 'start': start, 'stop': stop})
            intronseqs.append(fseq)

    r = {}  # Report
    for t, tseqs in zip(("Gene", "Exon", "Intron"), (geneseqs, exonseqs, intronseqs)):
        tsizes = [len(x) for x in tseqs]
        tsummary = SummaryStats(tsizes, dtype="int")
        r[t, "Number"] = tsummary.size
        r[t, "Average size (bp)"] = tsummary.mean
        r[t, "Median size (bp)"] = tsummary.median
        r[t, "Total length (Mb)"] = human_size(tsummary.sum, precision=0, target="Mb")
        r[t, "% of genome"] = percentage(tsummary.sum, s.totalsize, precision=0, mode=-1)
        r[t, "% GC"] = gc(tseqs)

    print(tabulate(r), file=sys.stderr)
866488
def multiprocess_permutation(bed_dict, mut_df, opts):
    
    chroms = sorted(bed_dict.keys())
    multiprocess_flag = opts['processes']>0
    if multiprocess_flag:
        num_processes = opts['processes']
    else:
        num_processes = 1
    num_permutations = opts['num_permutations']
    if not opts['by_sample']:
        obs_result = []
    else:
        uniq_samp = mut_df['Tumor_Sample'].unique()
        obs_result = pd.DataFrame(np.zeros((len(uniq_samp), len(cols))),
                                  index=uniq_samp, columns=cols)

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

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

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

    return result_list, obs_result
45835
def to_standardonly(in_file, ref_file, data):
    
    from bcbio.heterogeneity import chromhacks
    out_file = "%s-stdchrs.vcf.gz" % utils.splitext_plus(in_file)[0]
    if not utils.file_exists(out_file):
        stds = []
        for c in ref.file_contigs(ref_file):
            if chromhacks.is_nonalt(c.name):
                stds.append(c.name)
        if stds:
            with file_transaction(data, out_file) as tx_out_file:
                stds = ",".join(stds)
                in_file = bgzip_and_index(in_file, data["config"])
                cmd = "bcftools view -o {tx_out_file} -O z {in_file} {stds}"
                do.run(cmd.format(**locals()), "Subset to standard chromosomes")
    return bgzip_and_index(out_file, data["config"]) if utils.file_exists(out_file) else in_file
619475
def read_source_alleles(file_name):
    """"""
    alleles = {}
    open_function = open
    if file_name.endswith(".gz"):
        open_function = gzip.open
    with open_function(file_name, 'rb') as input_file:
        for line in input_file:
            row = line.rstrip("\r\n").split("\t")
            marker_name = row[0]
            allele1, allele2 = row[1].split(" ")
            alleles[marker_name] = {allele1, allele2}

    return alleles
670565
def create_dexseq_annotation(gff, count_file):
    
    out_file = count_file + ".ann"
    if file_exists(out_file):
        return out_file
    with file_transaction(out_file) as tx_out:
        with open(tx_out, 'w') as out_handle:
            with open(gff) as in_handle:
                for line in in_handle:
                    cols = line.strip().split("\t")
                    if cols[2] == "exonic_part":
                        exon = [f for f in cols[8].split(";") if f.strip().startswith("exonic_part_number")]
                        gene = [f for f in cols[8].split(";") if f.strip().startswith("gene_id")]
                        exon = exon[0].replace("\"", "").split()[1]
                        gene = gene[0].replace("\"", "").split()[1]
                        length = int(cols[4]) - int(cols[3]) + 1
                        line = "%s:%s\t%s\t%s\t%s\t%s\t%s\t%s\n" % (gene, exon, gene,
                                                                    cols[0], cols[3],
                                                                    cols[4],
                                                                    length, cols[6])
                        out_handle.write(line)
1149397
def concat_variant_files(orig_files, out_file, regions, ref_file, config):
    
    if not utils.file_exists(out_file):
        input_file_list = _get_file_list(orig_files, out_file, regions, ref_file, config)
        try:
            out_file = _run_concat_variant_files_gatk4(input_file_list, out_file, config)
        except subprocess.CalledProcessError as msg:
            if ("We require all VCFs to have complete VCF headers" in str(msg) or
                  "Features added out of order" in str(msg) or
                  "The reference allele cannot be missing" in str(msg)):
                out_file = _run_concat_variant_files_bcftools(input_file_list, out_file, config, naive=True)
            else:
                raise
    if out_file.endswith(".gz"):
        bgzip_and_index(out_file, config)
    return out_file
218612
def _make_ignore_file(work_dir, ref_file, ignore_file, impute_file):
    
    gl_file = os.path.join(work_dir, "gender_loci.txt")
    chroms = set([])
    with open(impute_file) as in_handle:
        for line in in_handle:
            chrom = line.split()[0]
            chroms.add(chrom)
            if not chrom.startswith("chr"):
                chroms.add("chr%s" % chrom)
    with open(ignore_file, "w") as out_handle:
        for contig in ref.file_contigs(ref_file):
            if contig.name not in chroms:
                out_handle.write("%s\n" % contig.name)
    with open(gl_file, "w") as out_handle:
        for contig in ref.file_contigs(ref_file):
            if contig.name in ["Y", "chrY"]:
                # From https://github.com/cancerit/cgpBattenberg/blob/dev/perl/share/gender/GRCh37d5_Y.loci
                positions = [2934912, 4546684, 4549638, 4550107]
                for pos in positions:
                    out_handle.write("%s\t%s\n" % (contig.name, pos))
    return ignore_file, gl_file
1036069
def reporter(self):
        
        logging.info('Preparing reports')
        # Populate self.plusdict in order to reuse parsing code from an assembly-based method
        for sample in self.runmetadata.samples:
            self.plusdict[sample.name] = dict()
            self.matchdict[sample.name] = dict()
            if sample.general.bestassemblyfile != 'NA':
                for gene in sample[self.analysistype].allelenames:
                    self.plusdict[sample.name][gene] = dict()
                    for allele, percentidentity in sample[self.analysistype].results.items():

                        if gene in allele:
                            # Split the allele number from the gene name using the appropriate delimiter
                            if '_' in allele:
                                splitter = '_'
                            elif '-' in allele:
                                splitter = '-'
                            else:
                                splitter = ''
                            self.matchdict[sample.name].update({gene: allele.split(splitter)[-1]})
                            # Create the plusdict dictionary as in the assembly-based (r)MLST method. Allows all the
                            # parsing and sequence typing code to be reused.
                            try:
                                self.plusdict[sample.name][gene][allele.split(splitter)[-1]][percentidentity] \
                                    = sample[self.analysistype].avgdepth[allele]
                            except KeyError:
                                self.plusdict[sample.name][gene][allele.split(splitter)[-1]] = dict()
                                self.plusdict[sample.name][gene][allele.split(splitter)[-1]][percentidentity] \
                                    = sample[self.analysistype].avgdepth[allele]
                    if gene not in self.matchdict[sample.name]:
                        self.matchdict[sample.name].update({gene: 'N'})
        self.profiler()
        self.sequencetyper()
        self.mlstreporter()
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()
550494
def haplotype_caller(align_bams, items, ref_file, assoc_files,
                       region=None, out_file=None):
    
    if out_file is None:
        out_file = "%s-variants.vcf.gz" % utils.splitext_plus(align_bams[0])[0]
    if not utils.file_exists(out_file):
        num_cores = dd.get_num_cores(items[0])
        broad_runner, params = \
            _shared_gatk_call_prep(align_bams, items, ref_file, region, out_file, num_cores)
        gatk_type = broad_runner.gatk_type()
        assert gatk_type in ["restricted", "gatk4"], \
            "Require full version of GATK 2.4+, or GATK4 for haplotype calling"
        with file_transaction(items[0], out_file) as tx_out_file:
            resources = config_utils.get_resources("gatk-spark", items[0]["config"])
            spark_opts = [str(x) for x in resources.get("options", [])]
            if _use_spark(num_cores, gatk_type, items, spark_opts):
                params += ["-T", "HaplotypeCallerSpark"]
                if spark_opts:
                    params += spark_opts
                else:
                    params += ["--spark-master", "local[%s]" % num_cores,
                               "--conf", "spark.local.dir=%s" % os.path.dirname(tx_out_file),
                               "--conf", "spark.driver.host=localhost", "--conf", "spark.network.timeout=800",
                               "--conf", "spark.executor.heartbeatInterval=100"]
            else:
                params += ["-T", "HaplotypeCaller"]
            params += ["--annotation", "ClippingRankSumTest",
                       "--annotation", "DepthPerSampleHC"]
            # Enable hardware based optimizations in GATK 3.1+
            if LooseVersion(broad_runner.gatk_major_version()) >= LooseVersion("3.1"):
                if _supports_avx():
                    # Scale down HMM thread default to avoid overuse of cores
                    # https://github.com/bcbio/bcbio-nextgen/issues/2442
                    if gatk_type == "gatk4":
                        params += ["--native-pair-hmm-threads", "1"]
                    # GATK4 selects the right HMM optimization automatically with FASTEST_AVAILABLE
                    # GATK3 needs to be explicitly set
                    else:
                        params += ["--pair_hmm_implementation", "VECTOR_LOGLESS_CACHING"]
            resources = config_utils.get_resources("gatk-haplotype", items[0]["config"])
            if "options" in resources:
                params += [str(x) for x in resources.get("options", [])]
            # Prepare gVCFs if doing joint calling
            is_joint = False
            if _joint_calling(items) or any("gvcf" in dd.get_tools_on(d) for d in items):
                is_joint = True
                # If joint calling parameters not set in user options
                if not any([x in ["--emit-ref-confidence", "-ERC", "--emitRefConfidence"] for x in params]):
                    if gatk_type == "gatk4":
                        params += ["--emit-ref-confidence", "GVCF"]
                    else:
                        params += ["--emitRefConfidence", "GVCF"]
                        params += ["--variant_index_type", "LINEAR", "--variant_index_parameter", "128000"]
                # Set GQ banding to not be single GQ resolution
                # No recommended default but try to balance resolution and size
                # http://gatkforums.broadinstitute.org/gatk/discussion/7051/recommendation-best-practices-gvcf-gq-bands

                if not any([x in ["-GQB"] for x in params]):
                    for boundary in [10, 20, 30, 40, 60, 80]:
                        params += ["-GQB", str(boundary)]
            # Enable non-diploid calling in GATK 3.3+
            if LooseVersion(broad_runner.gatk_major_version()) >= LooseVersion("3.3"):
                params += ["-ploidy", str(ploidy.get_ploidy(items, region))]
            if gatk_type == "gatk4":
                # GATK4 Spark calling does not support bgzipped output, use plain VCFs
                if is_joint and _use_spark(num_cores, gatk_type, items, spark_opts):
                    tx_out_file = tx_out_file.replace(".vcf.gz", ".vcf")
                params += ["--output", tx_out_file]
            else:
                params += ["-o", tx_out_file]
            broad_runner.new_resources("gatk-haplotype")
            memscale = {"magnitude": 0.9 * num_cores, "direction": "increase"} if num_cores > 1 else None
            try:
                broad_runner.run_gatk(params, os.path.dirname(tx_out_file), memscale=memscale,
                                      parallel_gc=_use_spark(num_cores, gatk_type, items, spark_opts))
            except subprocess.CalledProcessError as msg:
                # Spark failing on regions without any reads, write an empty VCF instead
                # https://github.com/broadinstitute/gatk/issues/4234
                if (_use_spark(num_cores, gatk_type, items, spark_opts) and
                      str(msg).find("java.lang.UnsupportedOperationException: empty collection") >= 0 and
                      str(msg).find("at org.apache.spark.rdd.RDD") >= 0):
                    vcfutils.write_empty_vcf(tx_out_file, samples=[dd.get_sample_name(d) for d in items])
                else:
                    raise
            if tx_out_file.endswith(".vcf"):
                vcfutils.bgzip_and_index(tx_out_file, items[0]["config"])


    # avoid bug in GATK where files can get output as non-compressed
    if out_file.endswith(".gz") and not os.path.exists(out_file + ".tbi"):
        with open(out_file, "r") as in_handle:
            is_plain_text = in_handle.readline().startswith("##fileformat")
        if is_plain_text:
            text_out_file = out_file
            out_file = out_file.replace(".vcf.gz", ".vcf")
            shutil.move(text_out_file, out_file)
    return vcfutils.bgzip_and_index(out_file, items[0]["config"])
341093
def from_command_line():
    
    # Parse options
    parser = argparse.ArgumentParser(
        description='Convert Complete Genomics var files to gVCF format.')
    parser.add_argument(
        '-d', '--refseqdir', metavar='REFSEQDIR', required=True,
        dest='refseqdir',
        help='Directory twobit reference genomes files are stored.')
    parser.add_argument(
        '-i', '--input', metavar='INPUTVARFILE',
        dest='cgivarfile',
        help='Path to Complete Genomics var file to convert. If omitted, data '
        ' also be piped in as standard input.')
    parser.add_argument(
        '-o', '--output', metavar='OUTPUTVCFFILE',
        dest='vcfoutfile',
        help='Path to where to save output VCF file.')
    parser.add_argument(
        '-D', '--download', action='store_true', dest='downloadrefseq',
        help='Download the 2bit file from UCSC to REFSEQDIR, if needed.')
    parser.add_argument(
        '-v', '--var-only', action='store_true', dest='varonly',
        help='Only report variant lines (i.e. VCF, but not gVCF)')
    args = parser.parse_args()

    # Get local twobit file from its directory. Download and store if needed.
    twobit_path, twobit_name = get_reference_genome_file(
        args.refseqdir, build='b37')
    # Handle input
    if sys.stdin.isatty():  # false if data is piped in
        var_input = args.cgivarfile
    else:
        var_input = sys.stdin
    # Handle output
    if args.vcfoutfile:
        convert_to_file(var_input,
                        args.vcfoutfile,
                        twobit_path,
                        twobit_name,
                        args.varonly)
    else:
        for line in convert(
                cgi_input=var_input,
                twobit_ref=twobit_path,
                twobit_name=twobit_name,
                var_only=args.varonly):
            print(line)
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
586621
def splitGenoSlidingWindow(pos,out_file,size=5e4,step=None):
    
    if step is None:    step = 0.5*size
    chroms = SP.unique(pos[:,0])

    RV = []
    wnd_i = 0
    wnd_file = csv.writer(open(out_file,'w'),delimiter='\t')
    nSnps = [] 
    for chrom_i in chroms:
        Ichrom = pos[:,0]==chrom_i
        idx_chrom_start = SP.where(Ichrom)[0][0]
        pos_chr = pos[Ichrom,1]
        start = pos_chr.min()
        pos_chr_max = pos_chr.max()
        while 1:
            if start>pos_chr_max: break
            end = start+size
            Ir = (pos_chr>=start)*(pos_chr<end)
            _nSnps = Ir.sum()
            if _nSnps>0:
                idx_wnd_start = idx_chrom_start+SP.where(Ir)[0][0]
                nSnps.append(_nSnps)
                line = SP.array([wnd_i,chrom_i,start,end,idx_wnd_start,_nSnps],dtype=int)
                wnd_file.writerow(line)
                wnd_i+=1
            start += step
    nSnps = SP.array(nSnps)
    return wnd_i,nSnps
538543
def _indexes_to_secondary_files(gresources, genome_build):
    
    out = {}
    for refname, val in gresources.items():
        if isinstance(val, dict) and "indexes" in val:
            # list of indexes -- aligners
            if len(val.keys()) == 1:
                indexes = sorted(val["indexes"])
                if len(indexes) == 0:
                    if refname not in alignment.allow_noindices():
                        raise ValueError("Did not find indexes for %s: %s" % (refname, val))
                elif len(indexes) == 1:
                    val = {"indexes": indexes[0]}
                else:
                    val = {"indexes": {"base": indexes[0], "indexes": indexes[1:]}}
            # directory plus indexes -- snpEff
            elif "base" in val and os.path.isdir(val["base"]) and len(val["indexes"]) > 0:
                indexes = val["indexes"]
                val = {"base": indexes[0], "indexes": indexes[1:]}
        elif isinstance(val, dict) and genome_build in val:
            val = _indexes_to_secondary_files(val, genome_build)
        out[refname] = val
    return out
254543
def get_coords(data):
    
    for category, vtypes in [("LOH", {"LOSS", "HETEROZYGOSITY"}),
                             ("amplification", {"AMPLIFICATION"})]:
        out = tz.get_in([category, dd.get_genome_build(data)], _COORDS, {})
        priority_file = dd.get_svprioritize(data)
        if priority_file:
            if os.path.basename(priority_file).find("civic") >= 0:
                for chrom, start, end, gene in _civic_regions(priority_file, vtypes, dd.get_disease(data)):
                    out[gene] = (chrom, start, end)
            elif os.path.basename(priority_file).find(".bed") >= 0:
                for line in utils.open_gzipsafe(priority_file):
                    parts = line.strip().split("\t")
                    if len(parts) >= 4:
                        chrom, start, end, gene = parts[:4]
                        out[gene] = (chrom, int(start), int(end))
        yield category, out