What are all the protein IDs
        (optionally restrict to a given chromosome and/or strand)
519499
def protein_ids(self, contig=None, strand=None):
        
        protein_ids = self._all_feature_values(
            column="protein_id",
            feature="CDS",
            contig=contig,
            strand=strand,
            distinct=True)
        # drop None values
        return [protein_id for protein_id in protein_ids if protein_id]
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
1251120
def reporter(self):
        
        # Create a set of all the gene names without alleles or accessions e.g. sul1_18_AY260546 becomes sul1
        genedict = dict()
        # Load the notes file to a dictionary
        notefile = os.path.join(self.targetpath, 'notes.txt')
        with open(notefile, 'r') as notes:
            for line in notes:
                # Ignore comment lines - they will break the parsing
                if line.startswith('#'):
                    continue
                # Split the line on colons e.g. stx1Aa:  Shiga toxin 1, subunit A, variant a: has three variables after
                # the split: gene(stx1Aa), description(Shiga toxin 1, subunit A, variant a), and _(\n)
                try:
                    gene, description, _ = line.split(':')
                # There are exceptions to the parsing. Some lines only have one :, while others have three. Allow for
                # these possibilities.
                except ValueError:
                    try:
                        gene, description = line.split(':')
                    except ValueError:
                        gene, description, _, _ = line.split(':')
                # Set up the description dictionary
                genedict[gene] = description.replace(', ', '_').strip()
        # Find unique gene names with the highest percent identity
        for sample in self.runmetadata.samples:
            try:
                if sample[self.analysistype].results:
                    # Initialise a dictionary to store the unique genes, and their percent identities
                    sample[self.analysistype].uniquegenes = dict()
                    for name, identity in sample[self.analysistype].results.items():
                        # Split the name of the gene from the string e.g. stx1:11:Z36899:11 yields stx1
                        if ':' in name:
                            sample[self.analysistype].delimiter = ':'
                        else:
                            sample[self.analysistype].delimiter = '_'
                        genename = name.split(sample[self.analysistype].delimiter)[0]
                        # Set the best observed percent identity for each unique gene
                        try:
                            # Pull the previous best identity from the dictionary
                            bestidentity = sample[self.analysistype].uniquegenes[genename]
                            # If the current identity is better than the old identity, save it
                            if float(identity) > float(bestidentity):
                                sample[self.analysistype].uniquegenes[genename] = float(identity)
                        # Initialise the dictionary if necessary
                        except KeyError:
                            sample[self.analysistype].uniquegenes[genename] = float(identity)
            except AttributeError:
                raise
        # Create the path in which the reports are stored
        make_path(self.reportpath)
        # Initialise strings to store the results
        data = 'Strain,Gene,Subtype/Allele,Description,Accession,PercentIdentity,FoldCoverage\n'
        with open(os.path.join(self.reportpath, self.analysistype + '.csv'), 'w') as report:
            for sample in self.runmetadata.samples:
                try:
                    if sample[self.analysistype].results:
                        # If there are many results for a sample, don't write the sample name in each line of the report
                        for name, identity in sorted(sample[self.analysistype].results.items()):
                            # Check to see which delimiter is used to separate the gene name, allele, accession, and
                            # subtype information in the header
                            if len(name.split(sample[self.analysistype].delimiter)) == 4:
                                # Split the name on the delimiter: stx2A:63:AF500190:d; gene: stx2A, allele: 63,
                                # accession: AF500190, subtype: d
                                genename, allele, accession, subtype = name.split(sample[self.analysistype].delimiter)
                            elif len(name.split(sample[self.analysistype].delimiter)) == 3:
                                # Treat samples without a subtype e.g. icaC:intercellular adhesion protein C: differently.
                                # Extract the allele as the 'subtype', and the gene name, and accession as above
                                genename, subtype, accession = name.split(sample[self.analysistype].delimiter)
                            else:
                                genename = name
                                subtype = ''
                                accession = ''
                            # Retrieve the best identity for each gene
                            percentid = sample[self.analysistype].uniquegenes[genename]
                            # If the percent identity of the current gene matches the best percent identity, add it to
                            # the report - there can be multiple occurrences of genes e.g.
                            # sul1,1,AY224185,100.00,840 and sul1,2,CP002151,100.00,927 are both included because they
                            # have the same 100% percent identity
                            if float(identity) == percentid:
                                # Treat the initial vs subsequent results for each sample slightly differently - instead
                                # of including the sample name, use an empty cell instead
                                try:
                                    description = genedict[genename]
                                except KeyError:
                                    description = 'na'
                                # Populate the results
                                data += '{samplename},{gene},{subtype},{description},{accession},{identity},{depth}\n'\
                                    .format(samplename=sample.name,
                                            gene=genename,
                                            subtype=subtype,
                                            description=description,
                                            accession=accession,
                                            identity=identity,
                                            depth=sample[self.analysistype].avgdepth[name])
                    else:
                        data += sample.name + '\n'
                except (KeyError, AttributeError):
                    data += sample.name + '\n'
            # Write the strings to the file
            report.write(data)
350304
def _make_variant_locus_id(gene_id, allele_id):
        

        varloci = '-'.join((gene_id, allele_id))
        varloci = '_:' + re.sub(r'(ZFIN)?:', '', varloci)

        return varloci
687469
def _validate_id_types(self):
        ''''''

        for sequences in [self.uniparc_sequences, self.fasta_sequences, self.seqres_sequences, self.rosetta_sequences]:
            for chain_id, sequence in sequences.iteritems():
                sequence_id_types = set(map(type, sequence.ids()))
                if sequence_id_types:
                    assert(len(sequence_id_types) == 1)
                    assert(sequence_id_types.pop() == types.IntType)

        for chain_id, sequence in self.atom_sequences.iteritems():
            sequence_id_types = set(map(type, sequence.ids()))
            assert(len(sequence_id_types) == 1)
            sequence_id_type = sequence_id_types.pop()
            assert(sequence_id_type == types.StringType or sequence_id_type == types.UnicodeType)
800909
def process_mgi_relationship_transgene_genes(self, limit=None):
        
        if self.test_mode:
            graph = self.testgraph
        else:
            graph = self.graph
        LOG.info("getting transgene genes")
        raw = '/'.join((self.rawdir, 'mgi_relationship_transgene_genes'))
        geno = Genotype(graph)
        col = [
            'rel_key', 'allele_key', 'allele_id', 'allele_label', 'category_key',
            'category_name', 'property_key', 'property_name', 'gene_num'
        ]
        with open(raw, 'r', encoding="utf8") as csvfile:
            filereader = csv.reader(csvfile, delimiter='\t', quotechar='\"')
            header = next(filereader)
            if header != col:
                LOG.error('expected columns:  %s\n\tBut got:\n%s', col, header)
            for row in filereader:
                # rel_key,
                allele_key = int(row[col.index('allele_key')])
                allele_id = row[col.index('allele_id')]
                # allele_label,
                # category_key,
                # category_name,
                # property_key,
                # property_name,
                gene_num = int(row[col.index('gene_num')])

                if self.test_mode and allele_key not in self.test_keys.get('allele')\
                        and gene_num not in self.test_ids:
                    continue

                gene_id = 'NCBIGene:' + str(gene_num)

                # geno.addParts(gene_id, allele_id, self.globaltt['has_variant_part'])
                seqalt_id = self.idhash['seqalt'].get(allele_key)
                if seqalt_id is None:
                    seqalt_id = allele_id
                geno.addSequenceDerivesFrom(seqalt_id, gene_id)

                if not self.test_mode and limit is not None and \
                        filereader.line_num > limit:
                    break

        return
356469
def signature(self):
        
        sig_elems = [self.seqid, self.source, self.type, \
                    self.start, self.end, self.strand, \
                    self.phase]
        if re.search("exon|CDS|UTR", self.type):
            parent = self.get_attr("Parent")
            if parent:
                (locus, iso) = atg_name(parent, retval="locus,iso", \
                        trimpad0=False)
                if locus:
                    sig_elems.append(locus)
        else:
            sig_elems.extend([self.accn])

        return ",".join(str(elem) for elem in sig_elems)
392289
private int getGroupTypIndicator(PolymerType polymerType) {
		if(PolymerType.PROTEIN_ONLY.contains(polymerType)){
			return 1;
		}
		if(PolymerType.POLYNUCLEOTIDE_ONLY.contains(polymerType)){
			return 2;
		}
		return 0;
	}
1248065
def filter(self, chromosome, **kwargs) :
		""""""
		def appendAllele(alleles, sources, snp) :
			pos = snp.start
			if snp.alt[0] == '-' :
				pass
				# print warn % ('DELETION', snpSet, snp.start, snp.chromosomeNumber)
			elif snp.ref[0] == '-' :
				pass
				# print warn % ('INSERTION', snpSet, snp.start, snp.chromosomeNumber)
			else :
				sources[snpSet] = snp
				alleles.append(snp.alt) #if not an indel append the polymorphism
			
			refAllele = chromosome.refSequence[pos]
			alleles.append(refAllele)
			sources['ref'] = refAllele
			return alleles, sources
				
		warn = 'Warning: the default snp filter ignores indels. IGNORED %s of SNP set: %s at pos: %s of chromosome: %s'
		sources = {}
		alleles = []
		for snpSet, data in kwargs.iteritems() :
			if type(data) is list :
				for snp in data :
					alleles, sources = appendAllele(alleles, sources, snp)
			else :
				allels, sources = appendAllele(alleles, sources, data)

		#appends the refence allele to the lot

		#optional we keep a record of the polymorphisms that were used during the process
		return SequenceSNP(alleles, sources = sources)
1275414
def get_ncbi_id_from_symbol(gene_symbol):
        
        monarch_url = 'https://solr.monarchinitiative.org/solr/search/select'
        params = DipperUtil._get_solr_weight_settings()
        params["q"] = "{0} \"{0}\"".format(gene_symbol)
        params["fq"] = ["taxon:\"NCBITaxon:9606\"", "category:\"gene\""]
        gene_id = None
        try:
            monarch_request = requests.get(monarch_url, params=params)
            response = monarch_request.json()
            count = response['response']['numFound']
            if count > 0:
                gene_id = response['response']['docs'][0]['id']
        except requests.ConnectionError:
            print("error fetching {0}".format(monarch_url))

        return gene_id
58448
def get_chromosomes(self, sv=False):
        
        if sv:
            res = self.db.structural_variant.distinct('chrom')
        else:
            res = self.db.variant.distinct('chrom')
            
        return res
1057932
def print_table(seqs, id2name, name):
    
    itable = open('%s.itable' % (name.rsplit('.', 1)[0]), 'w')
    print('\t'.join(['#sequence', 'gene', 'model', 'insertion', 'gene position', 'model position', 'length', 'orf?', 'intron?', 'orf?intron?', 'insertion', 'orf', 'intron']), file=itable)
    for seq, info in list(seqs.items()):
        gene, model, insertions = info
        name = id2name[seq]
        for i, ins in enumerate(insertions, 1):
            gene_pos, model_pos, length, iseq, \
                    orfs, introns, orfs_b, introns_b, orf_annotations = ins
            # append annotation to orf header
            for orf in orfs:
                parts = orf[0].split()
                annotation = orf_annotations[parts[0].split('>')[1]]
                orf[0] = '%s %s %s' % (parts[0], annotation, ' '.join(parts[1:]))
            # get orf position
            gene_pos = '-'.join([str(j) for j in gene_pos])
            # check if orf, intron is present
            if orfs_b is True or introns_b is True:
                orfs_introns_b = True
            else:
                orfs_introns_b = False
            out = [name, gene, model, i, gene_pos, model_pos, length, orfs_b, introns_b, orfs_introns_b]
            out.append('|'.join(iseq))
            out.append('|'.join(['|'.join(orf) for orf in orfs]))
            out.append('|'.join(['|'.join(intron) for intron in introns]))
            print('\t'.join([str(i) for i in out]), file=itable)
    itable.close()
785655
def gene_names(self, contig=None, strand=None):
        
        return self._all_feature_values(
            column="gene_name",
            feature="gene",
            contig=contig,
            strand=strand)
1317921
def genes_with_structures(self):
        """"""
        return DictList(x for x in self.genes if x.protein.num_structures > 0)
976365
def get_transcript_ids_for_ensembl_gene_ids(self, gene_ids, hgnc_symbols):
        
        
        chroms = {"1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", \
             "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", \
              "X", "Y"}
        
        headers = {"content-type": "application/json"}
        
        transcript_ids = []
        for gene_id in gene_ids:
            self.attempt = 0
            ext = "/overlap/id/{}?feature=transcript".format(gene_id)
            r = self.ensembl_request(ext, headers)
            
            for item in json.loads(r):
                # ignore non-coding transcripts
                if item["biotype"] not in ["protein_coding", "polymorphic_pseudogene"]:
                    continue
                
                # ignore transcripts not on the standard chromosomes
                # (non-default chroms fail to map the known de novo variants
                # to the gene location
                if item["Parent"] != gene_id or item["seq_region_name"] not in \
                        chroms or \
                        all([symbol not in item["external_name"] for symbol in hgnc_symbols]):
                    continue
                transcript_ids.append(item["id"])
        
        return transcript_ids
1203474
def get_chromosomes(snps):
    

    if isinstance(snps, pd.DataFrame):
        return list(pd.unique(snps["chrom"]))
    else:
        return []
750108
def gene_features(self):
        
        nm, strand = self.gene_name, self.strand
        feats = [(self.chrom, self.start, self.end, nm, strand, 'gene')]
        for feat in ('introns', 'exons', 'utr5', 'utr3', 'cdss'):
            fname = feat[:-1] if feat[-1] == 's' else feat
            res = getattr(self, feat)
            if res is None or all(r is None for r in res): continue
            if not isinstance(res, list): res = [res]
            feats.extend((self.chrom, s, e, nm, strand, fname) for s, e in res)

        tss = self.tss(down=1)
        if tss is not None:
            feats.append((self.chrom, tss[0], tss[1], nm, strand, 'tss'))
            prom = self.promoter()
            feats.append((self.chrom, prom[0], prom[1], nm, strand, 'promoter'))

        return sorted(feats, key=itemgetter(1))