Retrieve gene regions and names filtered by variant_types and diseases.
542201
def _civic_regions(civic_file, variant_types=None, diseases=None, drugs=None):
    
    if isinstance(diseases, six.string_types):
        diseases = [diseases]
    with utils.open_gzipsafe(civic_file) as in_handle:
        reader = csv.reader(in_handle, delimiter="\t")
        for chrom, start, end, info_str in reader:
            info = edn_loads(info_str)
            if not variant_types or _matches(info["support"]["variants"], variant_types):
                if not diseases or _matches(info["support"]["diseases"], diseases):
                    if not drugs or _matches(info["support"]["drugs"], drugs):
                        yield (chrom, int(start), int(end), list(info["name"])[0])
1280576
def subset_variant_regions(variant_regions, region, out_file, items=None, do_merge=True, data=None):
    
    if region is None:
        return variant_regions
    elif variant_regions is None:
        return region
    elif not isinstance(region, (list, tuple)) and region.find(":") > 0:
        raise ValueError("Partial chromosome regions not supported")
    else:
        merge_text = "-unmerged" if not do_merge else ""
        subset_file = "{0}".format(utils.splitext_plus(out_file)[0])
        subset_file += "%s-regions.bed" % (merge_text)
        if not os.path.exists(subset_file):
            data = items[0] if items else data
            with file_transaction(data, subset_file) as tx_subset_file:
                if isinstance(region, (list, tuple)):
                    _subset_bed_by_region(variant_regions, tx_subset_file, to_multiregion(region),
                                          dd.get_ref_file(data), do_merge=do_merge)
                else:
                    _rewrite_bed_with_chrom(variant_regions, tx_subset_file, region)
        if os.path.getsize(subset_file) == 0:
            return region
        else:
            return subset_file
704191
def _variants_fields(fields, exclude_fields, info_ids):
    """"""
    if fields is None:
        # no fields specified by user
        # by default extract all standard and INFO fields
        fields = config.STANDARD_VARIANT_FIELDS + info_ids
    else:
        # fields have been specified
        for f in fields:
            # check for non-standard fields not declared in INFO header
            if f not in config.STANDARD_VARIANT_FIELDS and f not in info_ids:
                # support extracting INFO even if not declared in header,
                # but warn...
                print('WARNING: no INFO definition found for field %s' % f,
                      file=sys.stderr)
    # process any exclusions
    if exclude_fields is not None:
        fields = [f for f in fields if f not in exclude_fields]
    return tuple(f for f in fields)
989011
def gene_by_id(self, gene_id):
        
        if gene_id not in self._genes:
            field_names = [
                "seqname",
                "start",
                "end",
                "strand",
            ]
            optional_field_names = [
                "gene_name",
                "gene_biotype",
            ]
            # Do not look for gene_name and gene_biotype if they are
            # not in the database.
            field_names.extend([
                name for name in optional_field_names
                if self.db.column_exists("gene", name)
            ])
            result = self.db.query_one(
                field_names,
                filter_column="gene_id",
                filter_value=gene_id,
                feature="gene")
            if not result:
                raise ValueError("Gene not found: %s" % (gene_id,))

            gene_name, gene_biotype = None, None
            assert len(result) >= 4 and len(result) <= 6, \
                "Result is not the expected length: %d" % len(result)
            contig, start, end, strand = result[:4]
            if len(result) == 5:
                if "gene_name" in field_names:
                    gene_name = result[4]
                else:
                    gene_biotype = result[4]
            elif len(result) == 6:
                gene_name, gene_biotype = result[4:]

            self._genes[gene_id] = Gene(
                gene_id=gene_id,
                gene_name=gene_name,
                contig=contig,
                start=start,
                end=end,
                strand=strand,
                biotype=gene_biotype,
                genome=self)

        return self._genes[gene_id]
91845
def get_variants(data, include_germline=False):
    
    data = utils.deepish_copy(data)
    supported = ["precalled", "vardict", "vardict-java", "vardict-perl",
                 "freebayes", "octopus", "strelka2"]
    # Right now mutect2 and mutect do not provide heterozygous germline calls
    # to be useful https://github.com/bcbio/bcbio-nextgen/issues/2464
    # supported += ["mutect2", "mutect"]
    if include_germline:
        supported.insert(1, "gatk-haplotype")
    out = []
    # CWL based input
    if isinstance(data.get("variants"), dict) and "samples" in data["variants"]:
        cur_vs = []
        # Unpack single sample list of files
        if (isinstance(data["variants"]["samples"], (list, tuple)) and
              len(data["variants"]["samples"]) == 1 and isinstance(data["variants"]["samples"][0], (list, tuple))):
            data["variants"]["samples"] = data["variants"]["samples"][0]
        for fname in data["variants"]["samples"]:
            variantcaller = utils.splitext_plus(os.path.basename(fname))[0]
            variantcaller = variantcaller.replace(dd.get_sample_name(data) + "-", "")
            for batch in dd.get_batches(data):
                variantcaller = variantcaller.replace(batch + "-", "")
            cur_vs.append({"vrn_file": fname, "variantcaller": variantcaller})
        data["variants"] = cur_vs
    for v in data.get("variants", []):
        if v["variantcaller"] in supported and v.get("vrn_file"):
            out.append((supported.index(v["variantcaller"]), v))
    out.sort()
    return [xs[1] for xs in out]
426084
def missing_uniprot_mapping(self):
        """"""
        uniprot_missing = []
        for g in self.genes:
            ups = g.protein.filter_sequences(UniProtProp)

            no_sequence_file_available = True
            for u in ups:
                if u.sequence_file or u.metadata_file:
                    no_sequence_file_available = False
                    break
            if no_sequence_file_available:
                uniprot_missing.append(g.id)

        return list(set(uniprot_missing))
249872
def from_gene_ids_and_names(cls, gene_names: Dict[str, str]):
        """"""
        genes = [ExpGene(id_, name=name) for id_, name in gene_names.items()]
        return cls.from_genes(genes)
475557
def get_gene_info(ensembl_ids=None, hgnc_symbols=None):
    
    uniq_ensembl_ids = set(ensembl_id for ensembl_id in (ensembl_ids or []))
    uniq_hgnc_symbols = set(hgnc_symbol for hgnc_symbol in (hgnc_symbols or []))
    genes = []
    gene_data = []
    
    if uniq_ensembl_ids:
        for ensembl_id in uniq_ensembl_ids:
            for res in query_gene(ensembl_id=ensembl_id):
                gene_data.append(res) 

    elif uniq_hgnc_symbols:
        for hgnc_symbol in uniq_hgnc_symbols:
            query_res = query_gene(hgnc_symbol=hgnc_symbol)
            if query_res:
                for res in query_res:
                    gene_data.append(res)
            else:
                # If no result we add just the symbol
                gene_data.append({
                    'hgnc_symbol': hgnc_symbol,
                    'hgnc_id': None,
                    'ensembl_id': None,
                    'description': None,
                    'chrom': 'unknown',
                    'start': 0,
                    'stop': 0,
                    'hi_score': None,
                    'constraint_score': None,
                })
    for gene in gene_data:
        genes.append(Gene(
            symbol=gene ['hgnc_symbol'],
            hgnc_id=gene['hgnc_id'],
            ensembl_id=gene['ensembl_id'],
            description=gene['description'],
            chrom=gene['chrom'],
            start=gene['start'],
            stop=gene['stop'],
            location=get_cytoband_coord(gene['chrom'], gene['start']),
            hi_score=gene['hi_score'],
            constraint_score=gene['constraint_score'],
            omim_number=get_omim_number(gene['hgnc_symbol'])
            ))

    return genes
393312
def from_genes(cls, genes: List[ExpGene]):
        """"""
        data = [g.to_dict() for g in genes]
        index = [d.pop('ensembl_id') for d in data]
        table = cls(data, index=index)
        return table
315445
def get_variants_in_region(self, chrom, start, end):
        """"""
        if not self.has_index:
            raise NotImplementedError("Not implemented when IMPUTE2 file is "
                                      "not indexed (see genipe)")

        if not self._index_has_location:
            raise NotImplementedError("Not implemented when index doesn't "
                                      "have location information.")

        # Getting the required variants
        required = self._impute2_index.loc[
            (self._impute2_index.chrom == CHROM_STR_TO_INT[chrom]) &
            (start <= self._impute2_index.pos) &
            (self._impute2_index.pos <= end)
        ]

        for name, variant_info in required.iterrows():
            for genotypes in self.get_variant_by_name(name, variant_info):
                self._fix_genotypes_object(genotypes, variant_info)
                yield genotypes
1271327
def _process_diseasegene(self, limit):
        
        if self.test_mode:
            graph = self.testgraph
        else:
            graph = self.graph
        line_counter = 0

        model = Model(graph)

        myfile = '/'.join((self.rawdir, self.files['disease-gene']['file']))

        for event, elem in ET.iterparse(myfile):
            if elem.tag == 'Disorder':
                # get the element name and id, ignore element name
                # id = elem.get('id') # some internal identifier
                disorder_num = elem.find('OrphaNumber').text
                disorder_id = 'ORPHA:' + str(disorder_num)

                if self.test_mode and disorder_id not in self.all_test_ids['disease']:
                    continue
                disorder_label = elem.find('Name').text

                # assuming that these are in the ontology (...any particular one?)
                model.addClassToGraph(disorder_id, disorder_label)
                assoc_list = elem.find('DisorderGeneAssociationList')
                expected_genes = assoc_list.get('count')
                LOG.info(
                    'Expecting %s genes associated with disorder %s.',
                    expected_genes, disorder_id)
                processed_genes = 0
                for assoc in assoc_list.findall('DisorderGeneAssociation'):
                    processed_genes += 1
                    gene = assoc.find('Gene')

                    # get gene's curie  HGNC or Ensembl ...

                    lclid = gene.find('OrphaNumber').text
                    gene_curie = 'ORPHA:' + lclid
                    gene_set = {'ORPHA': lclid}
                    for gene_ref in gene.findall(
                            './ExternalReferenceList/ExternalReference'):
                        gene_set[gene_ref.find('Source').text] = \
                            gene_ref.find('Reference').text

                    # set priority (clique leader if available) but default to OPRHA
                    for pfx in ('HGNC', 'Ensembl', 'SwissProt'):
                        if pfx in gene_set:
                            if pfx in self.localtt:
                                pfx = self.localtt[pfx]
                            gene_curie = pfx + ':' + gene_set[pfx]
                            gene_set.pop(pfx)
                            model.addClassToGraph(gene_curie, None)
                            break

                    # TEC have reservations w.r.t aggerator links being gene classes
                    for prefix in gene_set:
                        lclid = gene_set[prefix]
                        if prefix in self.localtt:
                            prefix = self.localtt[prefix]

                        dbxref = prefix + ':' + lclid

                        if gene_curie != dbxref:
                            model.addClassToGraph(dbxref, None)
                            model.addEquivalentClass(gene_curie, dbxref)

                    # TEC. would prefer this not happen here. let HGNC handle it
                    # except there are some w/o explicit external links ...

                    gene_symbol = gene.find('Symbol').text

                    syn_list = gene.find('./SynonymList')
                    if int(syn_list.get('count')) > 0:
                        for syn in syn_list.findall('./Synonym'):
                            model.addSynonym(gene_curie, syn.text)

                    dg_label = assoc.find('./DisorderGeneAssociationType/Name').text

                    # use dg association status to issue an evidence code
                    # FIXME I think that these codes are sub-optimal
                    eco_id = self.resolve(
                        assoc.find('DisorderGeneAssociationStatus/Name').text)

                    rel_id = self.resolve(dg_label)
                    
                    g2p_assoc = G2PAssoc(self.graph, self.name, gene_curie, disorder_id, rel_id)
                    g2p_assoc.add_evidence(eco_id)
                    g2p_assoc.add_association_to_graph()

                elem.clear()  # empty the element
                if int(expected_genes) != processed_genes:
                    LOG.warning(
                        '% expected %s associated genes but we processed %i',
                        disorder_id, expected_genes, processed_genes)

            if self.test_mode and limit is not None and line_counter > limit:
                return

        return
1152135
def genes_with_a_representative_sequence(self):
        """"""
        return DictList(x for x in self.genes if x.protein.representative_sequence)
10599
def population_variant_regions(items, merged=False):
    
    def _get_variant_regions(data):
        out = dd.get_variant_regions(data) or dd.get_sample_callable(data)
        # Only need to merge for variant region inputs, not callable BED regions which don't overlap
        if merged and dd.get_variant_regions(data):
            merged_out = dd.get_variant_regions_merged(data)
            if merged_out:
                out = merged_out
            else:
                out = merge_overlaps(out, data)
        return out
    import pybedtools
    if len(items) == 1:
        return _get_variant_regions(items[0])
    else:
        paired = vcfutils.get_paired(items)
        if paired:
            return _get_variant_regions(paired.tumor_data)
        else:
            vrs = []
            for data in items:
                vr_bed = _get_variant_regions(data)
                if vr_bed:
                    vrs.append((pybedtools.BedTool(vr_bed).total_coverage(), vr_bed))
            vrs.sort(reverse=True)
            if vrs:
                return vrs[0][1]
92200
def _get_regions(self):
        """"""
        if self._specs_in[_REGIONS_STR] == 'all':
            return [_get_all_objs_of_type(
                Region, getattr(self._obj_lib, 'regions', self._obj_lib)
            )]
        else:
            return [set(self._specs_in[_REGIONS_STR])]
1091344
def get_gene_names(self):
        """"""
        # Collect all gene names in network
        gene_names = []
        for node in self._nodes:
            members = node['data'].get('members')
            if members:
                gene_names += list(members.keys())
            else:
                if node['data']['name'].startswith('Group'):
                    continue
                gene_names.append(node['data']['name'])
        self._gene_names = gene_names
763855
def get_disease_mappings(self, att_ind_start):
        
        all_disease_ids = self.get_all_unique_diseases()
        disease_enum = enumerate(all_disease_ids, start=att_ind_start)
        disease_mappings = {}
        for num, dis in disease_enum:
            disease_mappings[dis] = num
        return disease_mappings
860718
def _get_genes(self, variant):
        
        ensembl_ids = []
        hgnc_symbols = []
        
        for transcript in variant.transcripts:
            if transcript.ensembl_id:
                ensembl_ids.append(transcript.ensembl_id)
            if transcript.hgnc_symbol:
                hgnc_symbols.append(transcript.hgnc_symbol)
        
        genes = get_gene_info(
                        ensembl_ids=ensembl_ids, 
                        hgnc_symbols=hgnc_symbols
                        )
        return genes