Fetch simple gene-term assocaitions from Golr using bioentity document type, one line per gene.
547269
def main():
    """"""
    import argparse
    prs = argparse.ArgumentParser(__doc__,
                                  formatter_class=argparse.ArgumentDefaultsHelpFormatter)

    prs.add_argument('--taxon_id', type=str,
                     help='NCBI taxon ID, must match exact species/strain used by GO Central, e.g. 4896 for S Pombe')
    prs.add_argument('--golr_url', default='http://golr.geneontology.org/solr/', type=str,
                     help='NCBI taxon ID, must match exact species/strain used by GO Central, e.g. 4896 for S Pombe')
    prs.add_argument('-o', default=None, type=str,
                     help="Specifies the name of the output file")
    prs.add_argument('--max_rows', default=100000, type=int,
                     help="maximum rows to be fetched")

    args = prs.parse_args()

    solr = pysolr.Solr(args.golr_url, timeout=30)

    sys.stderr.write("TAX:"+args.taxon_id+"\n")
    results = solr.search(q='document_category:"bioentity" AND taxon:"NCBITaxon:'+args.taxon_id+'"',
                          fl='bioentity_label,annotation_class_list', rows=args.max_rows)
    sys.stderr.write("NUM GENES:"+str(len(results))+"\n")
    if (len(results) ==0):
        sys.stderr.write("NO RESULTS")
        exit(1)
    if (len(results) == args.max_rows):
        sys.stderr.write("max_rows set too low")
        exit(1)

    file_out = sys.stdout if args.o is None else open(args.o, 'w')
    for r in results:
        gene_symbol = r['bioentity_label']
        sys.stderr.write(gene_symbol+"\n")
        if 'annotation_class_list' in r:
            file_out.write(r['bioentity_label']+"\t" + ';'.join(r['annotation_class_list'])+"\n")
        else:
            sys.stderr.write("no annotations for "+gene_symbol+"\n")
            

    if args.o is not None:
        file_out.close()
        sys.stdout.write("  WROTE: {}\n".format(args.o))
647307
def ensure_annotations(resources, data):
    
    transcript_gff = tz.get_in(["rnaseq", "transcripts"], resources)
    if transcript_gff and utils.file_exists(transcript_gff):
        out_dir = os.path.join(tz.get_in(["dirs", "work"], data),
                               "inputs", "data", "annotations")
        resources["rnaseq"]["gene_bed"] = gtf.gtf_to_bed(transcript_gff, out_dir)
    return resources
31071
def from_seqfeature(s, **kwargs):
    
    source = s.qualifiers.get('source', '.')[0]
    score = s.qualifiers.get('score', '.')[0]
    seqid = s.qualifiers.get('seqid', '.')[0]
    frame = s.qualifiers.get('frame', '.')[0]
    strand = _feature_strand[s.strand]

    # BioPython parses 1-based GenBank positions into 0-based for use within
    # Python.  We need to convert back to 1-based GFF format here.
    start = s.location.start.position + 1
    stop = s.location.end.position
    featuretype = s.type
    id = s.id
    attributes = dict(s.qualifiers)
    attributes.pop('source', '.')
    attributes.pop('score', '.')
    attributes.pop('seqid', '.')
    attributes.pop('frame', '.')
    return Feature(seqid, source, featuretype, start, stop, score, strand,
                   frame, attributes, id=id, **kwargs)
384218
def get_gene_name(cls, entry):
        
        gene_name = entry.find("./gene/name[@type='primary']")

        return gene_name.text if gene_name is not None and gene_name.text.strip() else None
605932
def get_gene_leaves(graph) -> Iterable[BaseEntity]:
    
    for node in get_nodes_by_function(graph, GENE):
        if graph.in_degree(node) != 0:
            continue

        if graph.out_degree(node) != 1:
            continue

        _, _, d = list(graph.out_edges(node, data=True))[0]

        if d[RELATION] == TRANSCRIBED_TO:
            yield node
396515
def df_kegg_metadata(self):
        """"""
        kegg_pre_df = []
        df_cols = ['gene', 'kegg', 'refseq', 'uniprot', 'num_pdbs', 'pdbs', 'seq_len', 'sequence_file', 'metadata_file']

        for g in self.genes:
            kegg_mappings = g.protein.filter_sequences(KEGGProp)
            for kegg_prop in kegg_mappings:
                kegg_dict = kegg_prop.get_dict(df_format=True, only_attributes=df_cols)
                kegg_dict['gene'] = g.id
                kegg_pre_df.append(kegg_dict)

        # Save a dataframe of the file mapping info
        df = pd.DataFrame.from_records(kegg_pre_df, columns=df_cols).set_index('gene')
        if df.empty:
            log.warning('Empty dataframe')
            return df
        else:
            return ssbio.utils.clean_df(df)
402292
def gene_one_per_line(in_file, data):
    
    if in_file:
        # Report all duplicated annotations one-per-line
        one_per_line_file = "%s-opl.bed" % utils.splitext_plus(in_file)[0]
        if not utils.file_uptodate(one_per_line_file, in_file):
            with file_transaction(data, one_per_line_file) as tx_out_file:
                with open(tx_out_file, 'w') as out:
                    for r in pybedtools.BedTool(in_file):
                        for g in r.name.split(','):
                            out.write('\t'.join(map(str, [r.chrom, r.start, r.end, g])) + '\n')
        return one_per_line_file
684424
def get_gene_ids(fusion_bed):
    
    with open(fusion_bed, 'r') as f:
        gene_to_id = {}
        regex = re.compile(r'(?P<gene>ENSG\d*)')
        for line in f:
            line = line.split('\t')
            transcript, gene_bit, name = line[3].split(';')
            m = regex.search(gene_bit)
            if m:
                gene_to_id[name] = m.group('gene')
    return gene_to_id
507483
def parse_collection(obj: dict) -> BioCCollection:
    """"""
    collection = BioCCollection()
    collection.source = obj['source']
    collection.date = obj['date']
    collection.key = obj['key']
    collection.infons = obj['infons']
    for doc in obj['documents']:
        collection.add_document(parse_doc(doc))
    return collection
1167150
def _init_objaartall(self):
        """"""
        kws = {
            'sortgo':lambda nt: [nt.NS, nt.dcnt],
            # fmtgo=('{p_fdr_bh:8.2e} {GO} '
            # Formatting for GO terms in grouped GO list
            'fmtgo':('{hdr1usr01:2} {NS} {GO} {s_fdr_bh:8} '
                     '{dcnt:5} {childcnt:3} R{reldepth:02} '
                     '{D1:5} {GO_name} ({study_count} study genes)\n'),
            # Formatting for GO terms listed under each gene
            'fmtgo2':('{hdr1usr01:2} {NS} {GO} {s_fdr_bh:8} '
                      '{dcnt:5} R{reldepth:02} '
                      '{GO_name} ({study_count} study genes)\n'),
            # itemid2name=ensmusg2symbol}
            }
        return AArtGeneProductSetsAll(self.grprdflt, self.hdrobj, **kws)
176247
def load_sequences_to_reference(self, sc=None, force_rerun=False):
        """"""
        log.info('Loading sequences to reference GEM-PRO...')
        from random import shuffle
        g_ids = [g.id for g in self.reference_gempro.functional_genes]
        shuffle(g_ids)

        def _load_sequences_to_reference_gene_sc(g_id, outdir=self.sequences_by_gene_dir,
                                                 g_to_pickle=self.gene_protein_pickles,
                                                 strain_infodict=self.strain_infodict,
                                                 orth_matrix=self.df_orthology_matrix, force_rerun=force_rerun):
            """Load orthologous strain sequences to reference Protein object, save as new pickle"""
            import ssbio.utils
            import ssbio.io
            from Bio import SeqIO
            import os.path as op

            protein_seqs_pickle_path = op.join(outdir, '{}_protein_withseqs.pckl'.format(g_id))

            if ssbio.utils.force_rerun(flag=force_rerun, outfile=protein_seqs_pickle_path):
                protein_pickle_path = g_to_pickle[g_id]
                protein_pickle = ssbio.io.load_pickle(protein_pickle_path)

                for strain, info in strain_infodict.items():
                    strain_sequences = SeqIO.index(info['genome_path'], 'fasta')
                    strain_gene_functional = info['functional_genes'][g_id]
                    if strain_gene_functional:
                        # Pull the gene ID of the strain from the orthology matrix
                        strain_gene_key = orth_matrix.at[g_id, strain]
                        new_id = '{}_{}'.format(g_id, strain)
                        if protein_pickle.sequences.has_id(new_id):
                            continue
                        protein_pickle.load_manual_sequence(seq=strain_sequences[strain_gene_key],
                                                            ident=new_id,
                                                            set_as_representative=False)
                protein_pickle.save_pickle(outfile=protein_seqs_pickle_path)

            return g_id, protein_seqs_pickle_path

        if sc:
            genes_rdd = sc.parallelize(g_ids)
            result = genes_rdd.map(_load_sequences_to_reference_gene_sc).collect()
        else:
            result = []
            for g in tqdm(g_ids):
                result.append(self._load_sequences_to_reference_gene(g, force_rerun))

        log.info('Storing paths to new Protein objects in self.gene_protein_pickles...')
        updated = []
        for g_id, protein_pickle in result:
            self.gene_protein_pickles[g_id] = protein_pickle
            updated.append(g_id)
        not_updated = set(list(self.gene_protein_pickles.keys())).difference(updated)
        log.info('No change to {} genes, removing from gene_protein_pickles'.format(len(not_updated)))
        log.debug(not_updated)
        for rem in not_updated:
            del self.gene_protein_pickles[rem]
151257
def _get_tcntobj(goids, go2obj, **kws):
        """"""
        # kws: gaf (gene2go taxid)
        if 'gaf' in kws or 'gene2go' in kws:
            # Get a reduced go2obj set for TermCounts
            _gosubdag = GoSubDag(goids, go2obj, rcntobj=False, prt=None)
            return get_tcntobj(_gosubdag.go2obj, **kws)
1304015
def to_list(self):
        
        src = self._source or ''
        coll = self._collection or ''
        desc = self._description or ''

        l = [self._id, src, coll, self._name,
             ','.join(sorted(self._genes)), desc]
        return l
1064170
def load_bioassembly_info_from_file(self, biomol_num):
        

        # current functionality is to take in a pre-assembled bioassembly file, parse it's MODELs and get info from that
        pass
        from Bio.PDB import PDBParser
        p = PDBParser(PERMISSIVE=True, QUIET=True)

        structure = StructureIO(self.structure_path, self.file_type)

        structure = p.get_structure('tmp', infile)
        bioass_to_chain_stoich = defaultdict(int)
        for model in structure:
            for chain in model:
                bioass_to_chain_stoich[chain.id] += 1
        return dict(bioass_to_chain_stoich)
231986
def epcrparse(self):
        
        from Bio.Blast.Applications import NcbiblastnCommandline
        while True:
            sample, record, line = self.epcrparsequeue.get()
            # Split the data on tabs
            gene, chromosome, strand, start, end, m_match, gaps, act_len_exp_len = line.split('\t')
            # Extract the gene sequence from the contigs
            # The record dictionary has the contig name, and the sequence. Splice out the data using the start and
            # end coordinates specified by ePCR
            genesequence = record[chromosome][int(start) - 1:int(end)]
            # Set up BLASTn using blastn-short, as the probe sequences tend to be very short
            blastn = NcbiblastnCommandline(db=sample[self.analysistype].probes.split('.')[0],
                                           num_threads=12,
                                           task='blastn-short',
                                           num_alignments=1,
                                           outfmt="'6 qseqid sseqid positive mismatch gaps "
                                                  "evalue bitscore slen length'")
            # Run the BLASTn, with the gene sequence as stdin
            out, err = blastn(stdin=genesequence)
            # Split the output string on tabs
            results = out.rstrip().split('\t')
            # Populate the raw blast results
            sample[self.analysistype].rawblastresults[gene] = results
            # Create named variables from the list
            positives = float(results[2])
            mismatches = float(results[3])
            gaps = float(results[4])
            subjectlength = float(results[7])
            # Calculate the percent identity
            percentidentity = float('{:0.2f}'.format((positives - gaps) / subjectlength * 100))
            # Create a dictionary with the desired values to store in the metadata object
            resultdict = {
                'matches': positives,
                'mismatches': mismatches,
                'gaps': gaps,
                'subject_length': subjectlength,
                'percent_identity': percentidentity,
                'match_length': results[8].split('\n')[0]
            }
            # Populate the metadata object with the dictionary
            sample[self.analysistype].blastresults[gene] = resultdict
            self.epcrparsequeue.task_done()
699016
def download_annotation(genomebuild, gene_file): 
    
    pred_bin = "genePredToBed"
    pred = find_executable(pred_bin)
    if not pred:
        sys.stderr.write("{} not found in path!\n".format(pred_bin))
        sys.exit(1)

    tmp = NamedTemporaryFile(delete=False, suffix=".gz")

    anno = []
    f = urlopen(UCSC_GENE_URL.format(genomebuild))
    p = re.compile(r'\w+.Gene.txt.gz')
    for line in f.readlines():
        m = p.search(line.decode())
        if m:
            anno.append(m.group(0))

    sys.stderr.write("Retrieving gene annotation for {}\n".format(genomebuild))
    url = ""
    for a in ANNOS:
        if a in anno:
            url = UCSC_GENE_URL.format(genomebuild) + a
            break
    if url:
        sys.stderr.write("Using {}\n".format(url))
        urlretrieve(
                url,
                tmp.name
                )
         
        with gzip.open(tmp.name) as f:
            cols = f.readline().decode(errors='ignore').split("\t")

        start_col = 1
        for i,col in enumerate(cols):
            if col == "+" or col == "-":
                start_col = i - 1
                break
        end_col = start_col + 10
       
        cmd = "zcat {} | cut -f{}-{} | {} /dev/stdin {}"
        print(cmd.format(tmp.name, start_col, end_col, pred, gene_file))
        sp.call(cmd.format(
            tmp.name, start_col, end_col, pred, gene_file), 
            shell=True)

    else:
        sys.stderr.write("No annotation found!")
944690
def _get_chrbands(self, limit, taxon):
        

        if limit is None:
            limit = sys.maxsize   # practical limit anyway
        model = Model(self.graph)
        line_counter = 0
        myfile = '/'.join((self.rawdir, self.files[taxon]['file']))
        LOG.info("Processing Chr bands from FILE: %s", myfile)
        geno = Genotype(self.graph)
        monochrom = Monochrom(self.graph_type, self.are_bnodes_skized)

        # used to hold band definitions for a chr
        # in order to compute extent of encompasing bands

        mybands = {}
        # build the organism's genome from the taxon
        genome_label = self.files[taxon]['genome_label']
        taxon_id = 'NCBITaxon:' + taxon

        # add the taxon as a class.  adding the class label elsewhere
        model.addClassToGraph(taxon_id, None)
        model.addSynonym(taxon_id, genome_label)

        geno.addGenome(taxon_id, genome_label)

        # add the build and the taxon it's in
        build_num = self.files[taxon]['build_num']
        build_id = 'UCSC:' + build_num
        geno.addReferenceGenome(build_id, build_num, taxon_id)

        # process the bands
        col = ['scaffold', 'start', 'stop', 'band_num', 'rtype']
        with gzip.open(myfile, 'rb') as f:
            for line in f:
                line_counter += 1
                # skip comments
                line = line.decode().strip()
                if line[0] == '#' or line_counter > limit:
                    continue
                # chr13	4500000	10000000	p12	stalk
                row = line.split('\t')
                scaffold = row[col.index('scaffold')]
                start = row[col.index('start')]
                stop = row[col.index('stop')]
                band_num = row[col.index('band_num')].strip()
                rtype = row[col.index('rtype')]

                # NOTE some less-finished genomes have
                # placed and unplaced scaffolds
                # * Placed scaffolds:
                #       the scaffolds have been placed within a chromosome.
                # * Unlocalized scaffolds:
                #   although the chromosome within which the scaffold occurs
                #   is known, the scaffold's position or orientation
                #   is not known.
                # * Unplaced scaffolds:
                #   it is not known which chromosome the scaffold belongs to
                #
                # find out if the thing is a full on chromosome, or a scaffold:
                # ex: unlocalized scaffold: chr10_KL568008v1_random
                # ex: unplaced scaffold: chrUn_AABR07022428v1
                placed_scaffold_pattern = r'(chr(?:\d+|X|Y|Z|W|M))'
                unlocalized_scaffold_pattern = placed_scaffold_pattern+r'_(\w+)_random'
                unplaced_scaffold_pattern = r'chr(Un(?:_\w+)?)'

                mch = re.match(placed_scaffold_pattern + r'$', scaffold)
                if mch is not None and len(mch.groups()) == 1:
                    # the chromosome is the first match of the pattern
                    chrom_num = mch.group(1)
                else:
                    # skip over anything that isn't a placed_scaffold
                    # at the class level
                    LOG.info("Found non-placed chromosome %s", scaffold)
                    chrom_num = None

                m_chr_unloc = re.match(unlocalized_scaffold_pattern, scaffold)
                m_chr_unplaced = re.match(unplaced_scaffold_pattern, scaffold)

                scaffold_num = None
                if mch:
                    pass
                elif m_chr_unloc is not None and len(m_chr_unloc.groups()) == 2:
                    chrom_num = m_chr_unloc.group(1)
                    scaffold_num = chrom_num + '_' + m_chr_unloc.group(2)
                elif m_chr_unplaced is not None and len(m_chr_unplaced.groups()) == 1:
                    scaffold_num = m_chr_unplaced.group(1)
                else:
                    LOG.error(
                        "There's a chr pattern that we aren't matching: %s", scaffold)

                if chrom_num is not None:
                    # the chrom class (generic) id
                    chrom_class_id = makeChromID(chrom_num, taxon, 'CHR')

                    # first, add the chromosome class (in the taxon)
                    geno.addChromosomeClass(
                        chrom_num, taxon_id, self.files[taxon]['genome_label'])

                    # then, add the chromosome instance (from the given build)
                    geno.addChromosomeInstance(
                        chrom_num, build_id, build_num, chrom_class_id)

                    # add the chr to the hashmap of coordinates for this build
                    # the chromosome coordinate space is itself
                    if chrom_num not in mybands.keys():
                        mybands[chrom_num] = {
                            'min': 0,
                            'max': int(stop),
                            'chr': chrom_num,
                            'ref': build_id,
                            'parent': None,
                            'stain': None,
                            'type': self.globaltt['chromosome']}

                if scaffold_num is not None:
                    # this will put the coordinates of the scaffold
                    # in the scaffold-space and make sure that the scaffold
                    # is part of the correct parent.
                    # if chrom_num is None,
                    # then it will attach it to the genome,
                    # just like a reg chrom
                    mybands[scaffold_num] = {
                        'min': start,
                        'max': stop,
                        'chr': scaffold_num,
                        'ref': build_id,
                        'parent': chrom_num,
                        'stain': None,
                        'type': self.globaltt['assembly_component'],
                        'synonym': scaffold}

                parents = list()
                if band_num is not None and band_num != '':
                    # add the specific band
                    mybands[chrom_num+band_num] = {
                        'min': start,
                        'max': stop,
                        'chr': chrom_num,
                        'ref': build_id,
                        'parent': None,
                        'stain': None,
                        'type': None}

                    # add the staining intensity of the band
                    if re.match(r'g(neg|pos|var)', rtype):
                        mybands[chrom_num+band_num]['stain'] = self.resolve(rtype)

                    # get the parent bands, and make them unique
                    parents = list(
                        monochrom.make_parent_bands(band_num, set()))
                    # alphabetical sort will put them in smallest to biggest,
                    # so we reverse
                    parents.sort(reverse=True)
                    # print('parents of',chrom,band,':',parents)

                    if len(parents) > 0:
                        mybands[chrom_num + band_num]['parent'] = chrom_num + parents[0]

                # loop through the parents and add them to the hash
                # add the parents to the graph, in hierarchical order
                # TODO PYLINT Consider using enumerate
                # instead of iterating with range and len
                for i in range(len(parents)):
                    rti = getChrPartTypeByNotation(parents[i])

                    pnum = chrom_num+parents[i]
                    sta = int(start)
                    sto = int(stop)
                    if pnum not in mybands.keys():
                        # add the parental band to the hash
                        bnd = {
                            'min': min(sta, sto),
                            'max': max(sta, sto),
                            'chr': chrom_num,
                            'ref': build_id,
                            'parent': None,
                            'stain': None,
                            'type': rti}
                        mybands[pnum] = bnd
                    else:
                        # band already in the hash means it's a grouping band
                        # need to update the min/max coords
                        bnd = mybands.get(pnum)
                        bnd['min'] = min(sta, sto, bnd['min'])
                        bnd['max'] = max(sta, sto, bnd['max'])
                        mybands[pnum] = bnd

                        # also, set the max for the chrom
                        chrom = mybands.get(chrom_num)
                        chrom['max'] = max(sta, sto, chrom['max'])
                        mybands[chrom_num] = chrom

                    # add the parent relationships to each
                    if i < len(parents) - 1:
                        mybands[pnum]['parent'] = chrom_num+parents[i+1]
                    else:
                        # add the last one (p or q usually)
                        # as attached to the chromosome
                        mybands[pnum]['parent'] = chrom_num

        f.close()  # end looping through file

        # loop through the hash and add the bands to the graph
        for bnd in mybands.keys():
            myband = mybands.get(bnd)
            band_class_id = makeChromID(bnd, taxon, 'CHR')
            band_class_label = makeChromLabel(bnd, genome_label)
            band_build_id = makeChromID(bnd, build_num, 'MONARCH')
            band_build_label = makeChromLabel(bnd, build_num)
            # the build-specific chrom
            chrom_in_build_id = makeChromID(myband['chr'], build_num, 'MONARCH')
            # if it's != part, then add the class
            if myband['type'] != self.globaltt['assembly_component']:
                model.addClassToGraph(
                    band_class_id, band_class_label, myband['type'])
                bfeature = Feature(
                    self.graph, band_build_id, band_build_label, band_class_id)
            else:
                bfeature = Feature(
                    self.graph, band_build_id, band_build_label, myband['type'])
                if 'synonym' in myband:
                    model.addSynonym(band_build_id, myband['synonym'])

            if myband['parent'] is None:
                if myband['type'] == self.globaltt['assembly_component']:
                    # since we likely don't know the chr,
                    # add it as a part of the build
                    geno.addParts(band_build_id, build_id)
            elif myband['type'] == self.globaltt['assembly_component']:
                # geno.addParts(band_build_id, chrom_in_build_id)
                parent_chrom_in_build = makeChromID(
                    myband['parent'], build_num, 'MONARCH')
                bfeature.addSubsequenceOfFeature(parent_chrom_in_build)

            # add the band as a feature
            # (which also instantiates the owl:Individual)
            bfeature.addFeatureStartLocation(myband['min'], chrom_in_build_id)
            bfeature.addFeatureEndLocation(myband['max'], chrom_in_build_id)
            if 'stain' in myband and myband['stain'] is not None:
                bfeature.addFeatureProperty(
                    self.globaltt['has_sequence_attribute'], myband['stain'])

            # type the band as a faldo:Region directly (add_region=False)
            # bfeature.setNoBNodes(self.nobnodes)
            # to come when we merge in ZFIN.py
            bfeature.addFeatureToGraph(False)

        return