Get a list of chromosomes to target and new updated ref_file header.

    Could potentially handle remapping from chr1 -> 1 but currently disabled due
    to speed issues.
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
174729
def ref_string_matches_dict_of_ref_sequences(self, ref_sequences):
        ''''''
        return self.CHROM in ref_sequences and self.ref_string_matches_ref_sequence(ref_sequences[self.CHROM])
1024956
def createBamHeader(self, baseHeader):
        
        header = dict(baseHeader)
        newSequences = []
        for index, referenceInfo in enumerate(header['SQ']):
            if index < self.numChromosomes:
                referenceName = referenceInfo['SN']
                # The sequence dictionary in the BAM file has to match up
                # with the sequence ids in the data, so we must be sure
                # that these still match up.
                assert referenceName == self.chromosomes[index]
                newReferenceInfo = {
                    'AS': self.referenceSetName,
                    'SN': referenceName,
                    'LN': 0,  # FIXME
                    'UR': 'http://example.com',
                    'M5': 'dbb6e8ece0b5de29da56601613007c2a',  # FIXME
                    'SP': 'Human'
                }
                newSequences.append(newReferenceInfo)
        header['SQ'] = newSequences
        return header
956567
def get_symbol_id_map(self):
        

        symbol_id_map = {}
        f = '/'.join((self.rawdir, self.files['genes']['file']))
        with open(f, 'r', encoding="utf8") as csvfile:
            filereader = csv.reader(csvfile, delimiter='\t', quotechar='\"')
            for row in filereader:
                (hgnc_id, symbol, name, locus_group, locus_type, status,
                 location, location_sortable, alias_symbol, alias_name,
                 prev_symbol, prev_name, gene_family, gene_family_id,
                 date_approved_reserved, date_symbol_changed,
                 date_name_changed, date_modified, entrez_id, ensembl_gene_id,
                 vega_id, ucsc_id, ena, refseq_accession, ccds_id, uniprot_ids,
                 pubmed_id, mgd_id, rgd_id, lsdb, cosmic, omim_id, mirbase,
                 homeodb, snornabase, bioparadigms_slc, orphanet,
                 pseudogene_org, horde_id, merops, imgt, iuphar,
                 kznf_gene_catalog, mamit_trnadb, cd, lncrnadb, enzyme_id,
                 intermediate_filament_db) = row

                symbol_id_map[symbol.strip()] = hgnc_id

        return symbol_id_map
125788
def _freqs_by_chromosome(in_file, params, somatic_info):
    
    freqs = []
    coords = []
    cur_chrom = None
    with pysam.VariantFile(in_file) as bcf_in:
        for rec in bcf_in:
            if _is_biallelic_snp(rec) and _passes_plus_germline(rec) and chromhacks.is_autosomal(rec.chrom):
                if cur_chrom is None or rec.chrom != cur_chrom:
                    if cur_chrom and len(freqs) > 0:
                        yield cur_chrom, freqs, coords
                    cur_chrom = rec.chrom
                    freqs = []
                    coords = []
                stats = _tumor_normal_stats(rec, somatic_info)
                if tz.get_in(["tumor", "depth"], stats, 0) > params["min_depth"]:
                    # not a ref only call
                    if len(rec.samples) == 0 or sum(rec.samples[somatic_info.tumor_name].allele_indices) > 0:
                        freqs.append(tz.get_in(["tumor", "freq"], stats))
                        coords.append(rec.start)
        if cur_chrom and len(freqs) > 0:
            yield cur_chrom, freqs, coords
275199
def _refresh_file_mapping(self):
        '''
        
        '''
        # map of suffix to description for imp
        if self.opts.get('cython_enable', True) is True:
            try:
                global pyximport
                pyximport = __import__('pyximport')  # pylint: disable=import-error
                pyximport.install()
                # add to suffix_map so file_mapping will pick it up
                self.suffix_map['.pyx'] = tuple()
            except ImportError:
                log.info('Cython is enabled in the options but not present '
                    'in the system path. Skipping Cython modules.')
        # Allow for zipimport of modules
        if self.opts.get('enable_zip_modules', True) is True:
            self.suffix_map['.zip'] = tuple()
        # allow for module dirs
        if USE_IMPORTLIB:
            self.suffix_map[''] = ('', '', MODULE_KIND_PKG_DIRECTORY)
        else:
            self.suffix_map[''] = ('', '', imp.PKG_DIRECTORY)

        # create mapping of filename (without suffix) to (path, suffix)
        # The files are added in order of priority, so order *must* be retained.
        self.file_mapping = salt.utils.odict.OrderedDict()

        opt_match = []

        def _replace_pre_ext(obj):
            '''
            Hack so we can get the optimization level that we replaced (if
            any) out of the re.sub call below. We use a list here because
            it is a persistent data structure that we will be able to
            access after re.sub is called.
            '''
            opt_match.append(obj)
            return ''

        for mod_dir in self.module_dirs:
            try:
                # Make sure we have a sorted listdir in order to have
                # expectable override results
                files = sorted(
                    x for x in os.listdir(mod_dir) if x != '__pycache__'
                )
            except OSError:
                continue  # Next mod_dir
            if six.PY3:
                try:
                    pycache_files = [
                        os.path.join('__pycache__', x) for x in
                        sorted(os.listdir(os.path.join(mod_dir, '__pycache__')))
                    ]
                except OSError:
                    pass
                else:
                    files.extend(pycache_files)

            for filename in files:
                try:
                    dirname, basename = os.path.split(filename)
                    if basename.startswith('_'):
                        # skip private modules
                        # log messages omitted for obviousness
                        continue  # Next filename
                    f_noext, ext = os.path.splitext(basename)
                    if six.PY3:
                        f_noext = PY3_PRE_EXT.sub(_replace_pre_ext, f_noext)
                        try:
                            opt_level = int(
                                opt_match.pop().group(1).rsplit('-', 1)[-1]
                            )
                        except (AttributeError, IndexError, ValueError):
                            # No regex match or no optimization level matched
                            opt_level = 0
                        try:
                            opt_index = self.opts['optimization_order'].index(opt_level)
                        except KeyError:
                            log.trace(
                                'Disallowed optimization level %d for module '
                                'name \'%s\', skipping. Add %d to the '
                                '\'optimization_order\' config option if you '
                                'do not want to ignore this optimization '
                                'level.', opt_level, f_noext, opt_level
                            )
                            continue
                    else:
                        # Optimization level not reflected in filename on PY2
                        opt_index = 0

                    # make sure it is a suffix we support
                    if ext not in self.suffix_map:
                        continue  # Next filename
                    if f_noext in self.disabled:
                        log.trace(
                            'Skipping %s, it is disabled by configuration',
                            filename
                        )
                        continue  # Next filename
                    fpath = os.path.join(mod_dir, filename)
                    # if its a directory, lets allow us to load that
                    if ext == '':
                        # is there something __init__?
                        subfiles = os.listdir(fpath)
                        for suffix in self.suffix_order:
                            if '' == suffix:
                                continue  # Next suffix (__init__ must have a suffix)
                            init_file = '__init__{0}'.format(suffix)
                            if init_file in subfiles:
                                break
                        else:
                            continue  # Next filename

                    try:
                        curr_ext = self.file_mapping[f_noext][1]
                        curr_opt_index = self.file_mapping[f_noext][2]
                    except KeyError:
                        pass
                    else:
                        if '' in (curr_ext, ext) and curr_ext != ext:
                            log.error(
                                'Module/package collision: \'%s\' and \'%s\'',
                                fpath,
                                self.file_mapping[f_noext][0]
                            )

                        if six.PY3 and ext == '.pyc' and curr_ext == '.pyc':
                            # Check the optimization level
                            if opt_index >= curr_opt_index:
                                # Module name match, but a higher-priority
                                # optimization level was already matched, skipping.
                                continue
                        elif not curr_ext or self.suffix_order.index(ext) >= self.suffix_order.index(curr_ext):
                            # Match found but a higher-priorty match already
                            # exists, so skip this.
                            continue

                    if six.PY3 and not dirname and ext == '.pyc':
                        # On Python 3, we should only load .pyc files from the
                        # __pycache__ subdirectory (i.e. when dirname is not an
                        # empty string).
                        continue

                    # Made it this far - add it
                    self.file_mapping[f_noext] = (fpath, ext, opt_index)

                except OSError:
                    continue
        for smod in self.static_modules:
            f_noext = smod.split('.')[-1]
            self.file_mapping[f_noext] = (smod, '.o', 0)
277055
def apply_all_link_refs(
    bytecode: bytes, link_refs: List[Dict[str, Any]], attr_dict: Dict[str, str]
) -> bytes:
    
    if link_refs is None:
        return bytecode
    link_fns = (
        apply_link_ref(offset, ref["length"], attr_dict[ref["name"]])
        for ref in link_refs
        for offset in ref["offsets"]
    )
    linked_bytecode = pipe(bytecode, *link_fns)
    return linked_bytecode
243056
def save_snps(self, filename=None):
        
        comment = (
            "# Source(s): {}\n"
            "# Assembly: {}\n"
            "# SNPs: {}\n"
            "# Chromosomes: {}\n".format(
                self.source, self.assembly, self.snp_count, self.chromosomes_summary
            )
        )

        if filename is None:
            filename = self.get_var_name() + "_lineage_" + self.assembly + ".csv"

        return lineage.save_df_as_csv(
            self._snps,
            self._output_dir,
            filename,
            comment=comment,
            header=["chromosome", "position", "genotype"],
        )
1188030
def multiprocess_permutation(bed_dict, mut_df, opts, indel_df=None):
    
    chroms = sorted(bed_dict.keys(), key=lambda x: len(bed_dict[x]), reverse=True)
    multiprocess_flag = opts['processes']>0
    if multiprocess_flag:
        num_processes = opts['processes']
    else:
        num_processes = 1
    #file_handle = open(opts['output'], 'w')
    file_handle = opts['handle']
    mywriter = csv.writer(file_handle, delimiter='\t', lineterminator='\n')
    if opts['maf'] and opts['num_iterations']:
        header = ['Gene', 'strand', 'Chromosome', 'Start_Position',
                  'End_Position', 'Reference_Allele', 'Tumor_Allele',
                  'Context', 'DNA_Change', 'Protein_Change', 'Variant_Classification']
    elif opts['maf']:
        header = ['Gene', 'strand', 'Chromosome', 'Start_Position',
                  'End_Position', 'Reference_Allele', 'Tumor_Allele',
                  'DNA_Change', 'Protein_Change', 'Variant_Classification',
                  'Tumor_Sample', 'Tumor_Type']
    else:
        header = ['Gene', 'ID', 'gene length', 'non-silent snv', 'silent snv', 'nonsense', 'lost stop',
                  'splice site', 'lost start', 'missense', 'recurrent missense',
                  'normalized missense position entropy',]
        # add column header for scores, is user provided one
        if opts['score_dir']:
            header += ['Total Missense MGAEntropy', 'Total Missense VEST Score']
        # add indel columns
        header += ['frameshift indel', 'inframe indel', 'normalized mutation entropy']
    mywriter.writerow(header)
    num_iterations = opts['num_iterations']

    # simulate indel counts
    if opts['summary'] and num_iterations:
        fs_cts, inframe_cts, gene_names = indel.simulate_indel_counts(indel_df,
                                                                      bed_dict,
                                                                      num_iterations,
                                                                      opts['seed'])
        name2ix = {gene_names[z]: z for z in range(len(gene_names))}
    # just count observed indels
    elif opts['summary']:
        # get gene names
        gene_names = [mybed.gene_name
                      for chrom in bed_dict
                      for mybed in bed_dict[chrom]]
        name2ix = {gene_names[z]: z for z in range(len(gene_names))}

        # initiate count vectors
        inframe_cts = np.zeros((1, len(gene_names)))
        fs_cts = np.zeros((1, len(gene_names)))

        # populate observed counts
        indel_cts_dict = indel_df['Gene'].value_counts().to_dict()
        fs_cts_dict = indel_df[indel.is_frameshift_annotation(indel_df)]['Gene'].value_counts().to_dict()
        for mygene in indel_cts_dict:
            if mygene in name2ix:
                # gene should be found in BED file annotation
                ix = name2ix[mygene]
                fs_cts[0, ix] = 0 if mygene not in fs_cts_dict else fs_cts_dict[mygene]
                inframe_cts[0, ix] = indel_cts_dict[mygene] - fs_cts[0, ix]

    # simulate snvs
    obs_result = []
    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:
                # iterate through each chromosome result
                for chrom_result in process_results:
                    # add columns for indels
                    if opts['summary']:
                        tmp_chrom_result = []
                        for gname, grp in it.groupby(chrom_result, lambda x: x[0]):
                            for l, row in enumerate(grp):
                                gene_ix = name2ix[gname]
                                fs_count = fs_cts[l, gene_ix]
                                inframe_count = inframe_cts[l, gene_ix]
                                missense_pos_ct = list(row.pop(-1).values())  # missense codon counts
                                silent_pos_ct = [1 for l in range(row[4])]
                                inactivating_ct = sum(row[5:9]) + fs_count
                                tmp_count_list = missense_pos_ct + silent_pos_ct + [inactivating_ct, inframe_count]
                                norm_ent = math.normalized_mutation_entropy(tmp_count_list)
                                tmp_chrom_result.append(row+[fs_count, inframe_count, norm_ent])
                        chrom_result = tmp_chrom_result

                    # write output to file
                    mywriter.writerows(chrom_result)
            except KeyboardInterrupt:
                pool.close()
                pool.join()
                logger.info('Exited by user. ctrl-c')
                sys.exit(0)
            pool.close()
            pool.join()
        else:
            # perform simulation
            info = (bed_dict[chroms[i]], mut_df, opts)
            chrom_results = singleprocess_permutation(info)

            # add indel columns
            if opts['summary']:
                tmp_chrom_result = []
                for gname, grp in it.groupby(chrom_results, lambda x: x[0]):
                    for l, row in enumerate(grp):
                        gene_ix = name2ix[gname]
                        fs_count = fs_cts[l, gene_ix]
                        inframe_count = inframe_cts[l, gene_ix]
                        missense_pos_ct = list(row.pop(-1).values())  # missense codon counts
                        silent_pos_ct = [1 for l in range(row[4])]
                        inactivating_ct = sum(row[5:9]) + fs_count
                        tmp_count_list = missense_pos_ct + silent_pos_ct + [inactivating_ct, inframe_count]
                        norm_ent = math.normalized_mutation_entropy(tmp_count_list)
                        tmp_chrom_result.append(row+[fs_count, inframe_count, norm_ent])
                chrom_results = tmp_chrom_result

            # write to file
            mywriter.writerows(chrom_results)
848454
def get_reference_genome_file(refseqdir, build):
    
    if not os.path.exists(refseqdir) or not os.path.isdir(refseqdir):
        raise ValueError("No directory at {}".format(refseqdir))
    twobit_name = ''
    if build in ['b37', 'build 37', 'build37', '37', 'hg19']:
        twobit_name = 'hg19.2bit'
        build = 'build37'
    if not twobit_name:
        raise ValueError('Genome bulid "{}" not supported.'.format(build))
    twobit_path = os.path.join(refseqdir, twobit_name)
    if not os.path.exists(twobit_path):
        twobitdownload.save_genome('hg19', destdir=refseqdir)
    return twobit_path, twobit_name
1000796
def getPysamVariants(self, referenceName, startPosition, endPosition):
        
        if referenceName in self._chromFileMap:
            varFileName = self._chromFileMap[referenceName]
            referenceName, startPosition, endPosition = \
                self.sanitizeVariantFileFetch(
                    referenceName, startPosition, endPosition)
            cursor = self.getFileHandle(varFileName).fetch(
                referenceName, startPosition, endPosition)
            for record in cursor:
                yield record
130419
def remap_index_fn(ref_file):
    
    return os.path.join(os.path.dirname(os.path.dirname(ref_file)), "star")
398390
def process_headers(data, filename):
    """"""
    headers = []
    if 'toc' in data:
        for element in PyQuery(data['toc'])('a'):
            headers.append(recurse_while_none(element))
        if None in headers:
            log.info('Unable to index file headers for: %s', filename)
    return headers
857926
def _add_to_ref(self, rec_curr, line):
        """"""
        # Examples of record lines containing ':' include:
        #   id: GO:0000002
        #   name: mitochondrial genome maintenance
        #   namespace: biological_process
        #   def: "The maintenance of ...
        #   is_a: GO:0007005 ! mitochondrion organization
        if line[:4] == "id: ":
            assert not rec_curr.item_id
            item_id = line[4:]
            rec_curr.item_id = item_id
            rec_curr.id = item_id
        elif line[:8] == "alt_id: ":
            rec_curr.alt_ids.add(line[8:])
        elif line[:6] == "name: ":
            assert not rec_curr.name
            rec_curr.name = line[6:]
        elif line[:11] == "namespace: ":
            assert not rec_curr.namespace
            rec_curr.namespace = line[11:]
        elif line[:6] == "is_a: ":
            rec_curr._parents.add(line[6:].split()[0])
        elif line[:13] == "is_obsolete: " and line[13:] == "true":
            rec_curr.is_obsolete = True
        elif self.optobj and ':' in line:
            self.optobj.update_rec(rec_curr, line)
1203474
def get_chromosomes(snps):
    

    if isinstance(snps, pd.DataFrame):
        return list(pd.unique(snps["chrom"]))
    else:
        return []
799136
def sample_chromosomes(job, genome_fai_file):
    
    work_dir = os.getcwd()
    genome_fai = untargz(job.fileStore.readGlobalFile(genome_fai_file), work_dir)
    return chromosomes_from_fai(genome_fai)
674332
def _read_reference_information(self):
    """"""
    ref_lengths = {}
    ref_names = []
    for n in range(self._n_ref):
      l_name = struct.unpack('<i',self._fh.read(4))[0]
      name = self._fh.read(l_name).rstrip('\0')
      l_ref = struct.unpack('<i',self._fh.read(4))[0]
      ref_lengths[name] = l_ref
      ref_names.append(name)
    return ref_lengths, ref_names