What are all the gene IDs
        (optionally restrict to a given chromosome/contig and/or strand)
499037
def gene_ids(self, contig=None, strand=None):
        
        return self._all_feature_values(
            column="gene_id",
            feature="gene",
            contig=contig,
            strand=strand)
458110
def get_genomes(cfg):
    
    if cfg['host']=='homo_sapiens':
        contigs=['1',
                 '2',
                 '3',
                 '4',
                 '5',
                 '6',
                 '7',
                 '8',
                 '9',
                 '10',
                 '11',
                 '12',
                 '13',
                 '14',
                 '15',
                 '16',
                 '17',
                 '18',
                 '19',
                 '20',
                 '21',
                 '22',
                 'X','Y']
    else:
        runbashcmd(f"pyensembl install --reference-name {cfg['genomeassembly']} --release {cfg['genomerelease']} --species {cfg['host']}")

        import pyensembl
        ensembl = pyensembl.EnsemblRelease(species=pyensembl.species.Species.register(
                                            latin_name=cfg['host'],
                                            synonyms=[cfg['host']],
                                            reference_assemblies={
                                                cfg['genomeassembly']: (cfg['genomerelease'], cfg['genomerelease']),
                                            }),release=cfg['genomerelease'])
        contig_mito=['MTDNA','MITO','MT']
        contigs=[c for c in ensembl.contigs() if ((not '.' in c) and (len(c)<5) and (c not in contig_mito))]    
    if len(contigs)==0:
        contigs=[c for c in ensembl.contigs()]
        # logging.error([c for c in ensembl.contigs()])
        # logging.error('no contigs identified by pyensembl; aborting')
        # sys.exit(0)
    logging.info(f"{len(contigs)} contigs/chromosomes in the genome")
    logging.info(contigs)
    # raw genome next
    if 'human' in cfg['host'].lower():
        cfg['host']='homo_sapiens'
    if 'yeast' in cfg['host'].lower():
        cfg['host']='saccharomyces_cerevisiae'
    host_="_".join(s for s in cfg['host'].split('_')).capitalize()
    if 'GRCh37' in cfg['genomeassembly']:    
        ensembl_fastad=f"pub/grch37/update/fasta/{cfg['host']}/dna/"
    else:
        ensembl_fastad=f"pub/release-{cfg['genomerelease']}/fasta/{cfg['host']}/dna/"        
    genome_fastad=f"{cfg['genomed']}/{ensembl_fastad}"
    cfg['genomep']=f'{genome_fastad}/genome.fa'.format()
    if not exists(cfg['genomep']):
        logging.error(f"not found: {cfg['genomep']}")
        ifdlref = input(f"Genome file are not there at {genome_fastad}.\n Download?[Y/n]: ")
        if ifdlref=='Y':
        # #FIXME download contigs and cat and get index, sizes
            for contig in contigs:
                fn=f"{cfg['host'].capitalize()}.{cfg['genomeassembly']}.dna_sm.chromosome.{contig}.fa.gz"
                fp=f'{ensembl_fastad}/{fn}'
                if not exists(f"{cfg['genomed']}{fp.replace('.gz','')}"):
                    if not exists(f"{cfg['genomed']}{fp}"):
                        logging.info(f'downloading {fp}')
                        cmd=f"wget -q -x -nH ftp://ftp.ensembl.org/{fp} -P {cfg['genomed']}"
                        try:
                            runbashcmd(cmd,test=cfg['test'])
                        except:
                            fn=f"{cfg['host'].capitalize()}.{cfg['genomeassembly']}.dna_sm.toplevel.fa.gz"
                            fp='{}/{}'.format(ensembl_fastad,fn)                        
                            if not exists(fp):
                                cmd=f"wget -q -x -nH ftp://ftp.ensembl.org/{fp} -P {cfg['genomed']}"
                                # print(cmd)
                                runbashcmd(cmd,test=cfg['test'])
                                break
                #                 break
            # make the fa ready
            if not exists(cfg['genomep']):
                cmd=f"gunzip {genome_fastad}*.fa.gz;cat {genome_fastad}/*.fa > {genome_fastad}/genome.fa;"
                runbashcmd(cmd,test=cfg['test'])
        else:
            logging.error('abort')
            sys.exit(1)
    if not exists(cfg['genomep']+'.bwt'):
        cmd=f"{cfg['bwa']} index {cfg['genomep']}"
        runbashcmd(cmd,test=cfg['test'])
    else:        
        logging.info('bwa index is present')
    if not exists(cfg['genomep']+'.fai'):
        cmd=f"{cfg['samtools']} faidx {cfg['genomep']}"
        runbashcmd(cmd,test=cfg['test'])
    else:
        logging.info('samtools index is present')
    if not exists(cfg['genomep']+'.sizes'):
        cmd=f"cut -f1,2 {cfg['genomep']}.fai > {cfg['genomep']}.sizes"            
        runbashcmd(cmd,test=cfg['test'])
    else:
        logging.info('sizes of contigs are present')

    if 'GRCh37' in cfg['genomeassembly']:    
#     ftp://ftp.ensembl.org/pub/grch37/update/gff3/homo_sapiens/Homo_sapiens.GRCh37.87.gff3.gz
        ensembl_gff3d=f"pub/grch37/update/gff3/{cfg['host']}/"    
    else:
        ensembl_gff3d=f"pub/release-{cfg['genomerelease']}/gff3/{cfg['host']}/"    
    
    genome_gff3d=f"{cfg['genomed']}/{ensembl_gff3d}"
    cfg['genomegffp']=f'{genome_gff3d}/genome.gff3'
    if not exists(cfg['genomegffp']):
        logging.error('not found: {}'.format(cfg['genomegffp']))
        ifdlref = input("Download genome annotations at {}?[Y/n]: ".format(genome_gff3d))
        if ifdlref=='Y':
        # #FIXME download contigs and cat and get index, sizes
            fn=f"{cfg['host'].capitalize()}.{cfg['genomeassembly']}.{cfg['genomerelease']}.gff3.gz"
            fp=f"{ensembl_gff3d}/{fn}"
            if not exists(fp):
                cmd=f"wget -x -nH ftp://ftp.ensembl.org/{fp} -P {cfg['genomed']}"
                runbashcmd(cmd,test=cfg['test'])
            # move to genome.gff3
                cmd=f"cp {genome_gff3d}/{fn} {cfg['genomegffp']}"
                runbashcmd(cmd,test=cfg['test'])

        else:
            logging.error('abort')
            sys.exit(1)
    logging.info('genomes are installed!')
    return cfg
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)
683988
def gffatributes2ids(s):
    
    Name,gene_id,transcript_id,protein_id,exon_id=np.nan,np.nan,np.nan,np.nan,np.nan
    if '=' in s:
        d=dict([i.split('=') for i in s.split(';')])
        if 'Parent' in d:
            d[d['Parent'].split(':')[0]+'_id']=d['Parent'].split(':')[1]
        Name,gene_id,transcript_id,protein_id,exon_id=np.nan,np.nan,np.nan,np.nan,np.nan
        if 'Name' in d:    
            Name=d['Name']
        if 'gene_id' in d:    
            gene_id=d['gene_id']
        if 'transcript_id' in d:    
            transcript_id=d['transcript_id']
        if 'protein_id' in d:    
            protein_id=d['protein_id']
        if 'exon_id' in d:    
            exon_id=d['exon_id']
    return Name,gene_id,transcript_id,protein_id,exon_id
1136850
def _try_parse_gene_association(self, reaction_id, s):
        
        s = s.strip()
        if s == '':
            return None

        try:
            return boolean.Expression(s)
        except boolean.ParseError as e:
            msg = 'Failed to parse gene association for {}: {}'.format(
                reaction_id, text_type(e))
            if e.indicator is not None:
                msg += '\n{}\n{}'.format(s, e.indicator)
            logger.warning(msg)

        return s
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)
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)
706195
def seq_info(names, id2names, insertions, sequences):
    
    seqs = {} # seqs[id] = [gene, model, [[i-gene_pos, i-model_pos, i-length, iseq, [orfs], [introns]], ...]]
    for name in names:
        id = id2names[name]
        gene = name.split('fromHMM::', 1)[0].rsplit(' ', 1)[1]
        model = name.split('fromHMM::', 1)[1].split('=', 1)[1].split()[0]
        i_gene_pos = insertions[id] # coordinates of each insertion wrt gene
        i_model_pos = name.split('fromHMM::', 1)[1].split('model-pos(ins-len)=')[1].split()[0].split(';') # model overlap
        i_info = []
        for i, ins in enumerate(i_gene_pos):
            model_pos = i_model_pos[i].split('-')[1].split('(')[0]
            length = i_model_pos[i].split('(')[1].split(')')[0]
            iheader = '>%s_%s insertion::seq=%s type=insertion strand=n/a gene-pos=%s-%s model-pos=%s'\
                    % (id, (i + 1), (i + 1), ins[0], ins[1], model_pos)
            iseq = sequences[id][1][ins[0]:(ins[1] + 1)]
            iseq = [iheader, iseq]
            info = [ins, model_pos, length, iseq, [], []]
            i_info.append(info)
        seqs[id] = [gene, model, i_info]
    return seqs
930151
def idsKEGG(organism):
    
    ORG=urlopen("http://rest.kegg.jp/list/"+organism).read()
    ORG=ORG.split("\n")
    final=[]
    for k in ORG:
        final.append(k.split("\t"))
    df=pd.DataFrame(final[0:len(final)-1])[[0,1]]
    df.columns=['KEGGid','description']
    field = pd.DataFrame(df['description'].str.split(';',1).tolist())[0]
    field = pd.DataFrame(field)
    df = pd.concat([df[['KEGGid']],field],axis=1)
    df.columns=['KEGGid','gene_name']
    df=df[['gene_name','KEGGid']]
    return df
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))
223849
def from_gene_ids(cls, gene_ids: List[str]):
        """"""
        genes = [ExpGene(id_) for id_ in gene_ids]
        return cls.from_genes(genes)
316366
def get_assc_pruned(assc_geneid2gos, min_genecnt=None, max_genecnt=None, prt=sys.stdout):
    """"""
    # DEFN WAS: get_assc_pruned(assc_geneid2gos, max_genecnt=None, prt=sys.stdout):
    #      ADDED min_genecnt argument and functionality
    if max_genecnt is None and min_genecnt is None:
        return assc_geneid2gos, set()
    go2genes_orig = get_b2aset(assc_geneid2gos)
    # go2genes_prun = {go:gs for go, gs in go2genes_orig.items() if len(gs) <= max_genecnt}
    go2genes_prun = {}
    for goid, genes in go2genes_orig.items():
        num_genes = len(genes)
        if (min_genecnt is None or num_genes >= min_genecnt) and \
           (max_genecnt is None or num_genes <= max_genecnt):
            go2genes_prun[goid] = genes
    num_was = len(go2genes_orig)
    num_now = len(go2genes_prun)
    gos_rm = set(go2genes_orig.keys()).difference(set(go2genes_prun.keys()))
    assert num_was-num_now == len(gos_rm)
    if prt is not None:
        if min_genecnt is None:
            min_genecnt = 1
        if max_genecnt is None:
            max_genecnt = "Max"
        prt.write("{N:4} GO IDs pruned. Kept {NOW} GOs assc w/({m} to {M} genes)\n".format(
            m=min_genecnt, M=max_genecnt, N=num_was-num_now, NOW=num_now))
    return get_b2aset(go2genes_prun), gos_rm
599790
def missing_pdb_structures(self):
        """"""
        return [x.id for x in self.genes if not self.genes_with_experimental_structures.has_id(x.id)]
138082
def get_alt(self, individual=0, nucleotides_only=True):
        """"""
        #not i.startswith(',') is put in to handle cases like <DEL:ME:ALU> where we have no alternate allele
        #but some reference
        if isinstance(individual, str):
            individual = self.individuals[individual]
        if nucleotides_only:
            return [self.alt[i-1].replace('.', '') for i in self.genotype[individual] if i > 0 and not self.alt[i-1].startswith('<')]
        else:
            return [self.alt[i-1].replace('.', '') for i in self.genotype[individual] if i > 0]
685830
def delete_gene(self, *gene_ids):
        """"""
        self.gene_ids = [gene_id for gene_id in self.gene_ids
                         if gene_id not in gene_ids]
745574
def locus(args):
    
    from jcvi.formats.sam import get_minibam
    # See `Format-lobSTR-database.ipynb` for a list of TREDs for validation
    INCLUDE = ["HD", "SBMA", "SCA1", "SCA2", "SCA8", "SCA17", "DM1", "DM2",
               "FXTAS"]
    db_choices = ("hg38", "hg19")

    p = OptionParser(locus.__doc__)
    p.add_option("--tred", choices=INCLUDE,
                 help="TRED name")
    p.add_option("--ref", choices=db_choices, default="hg38",
                 help="Reference genome")
    p.set_home("lobstr")
    opts, args = p.parse_args(args)

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

    bamfile, = args
    ref = opts.ref
    lhome = opts.lobstr_home
    tred = opts.tred

    tredsfile = datafile("TREDs.meta.csv")
    tf = pd.read_csv(tredsfile, index_col=0)
    row = tf.ix[tred]
    tag = "repeat_location"
    ldb = "TREDs"
    if ref == "hg19":
        tag += "." + ref
        ldb += "-" + ref
    seqid, start_end = row[tag].split(":")

    PAD = 1000
    start, end = start_end.split('-')
    start, end = int(start) - PAD, int(end) + PAD
    region = "{}:{}-{}".format(seqid, start, end)

    minibamfile = get_minibam(bamfile, region)
    c = seqid.replace("chr", "")
    cmd, vcf = allelotype_on_chr(minibamfile, c, lhome, ldb)
    sh(cmd)

    parser = LobSTRvcf(columnidsfile=None)
    parser.parse(vcf, filtered=False)
    items = parser.items()
    if not items:
        print("No entry found!", file=sys.stderr)
        return

    k, v = parser.items()[0]
    print("{} => {}".format(tred, v.replace(',', '/')), file=sys.stderr)
222479
def generate_complementary_strand(strand1):
        """"""
        rise_adjust = (
            strand1.rise_per_nucleotide * strand1.axis.unit_tangent) * 2
        strand2 = NucleicAcidStrand.from_start_and_end(
            strand1.helix_end - rise_adjust, strand1.helix_start - rise_adjust,
            generate_antisense_sequence(strand1.base_sequence),
            phos_3_prime=strand1.phos_3_prime)
        ad_ang = dihedral(strand1[0]["C1'"]._vector, strand1.axis.start,
                          strand2.axis.start + rise_adjust,
                          strand2[-1]["C1'"]._vector)
        strand2.rotate(
            225.0 + ad_ang, strand2.axis.unit_tangent,
            point=strand2.helix_start)  # 225 is the base adjust
        return strand2