define ribosomal proteins and location of curated databases
191635
def find_databases(databases):
    
    # 16 ribosomal proteins in their expected order
    proteins = ['L15', 'L18', 'L6', 'S8', 'L5', 'L24', 'L14',
            'S17', 'L16', 'S3', 'L22', 'S19', 'L2', 'L4', 'L3', 'S10']
    # curated databases
    protein_databases = {
                'L14': 'rpL14_JGI_MDM.filtered.faa',
                'L15': 'rpL15_JGI_MDM.filtered.faa',
                'L16': 'rpL16_JGI_MDM.filtered.faa',
                'L18': 'rpL18_JGI_MDM.filtered.faa',
                'L22': 'rpL22_JGI_MDM.filtered.faa',
                'L24': 'rpL24_JGI_MDM.filtered.faa',
                'L2': 'rpL2_JGI_MDM.filtered.faa',
                'L3': 'rpL3_JGI_MDM.filtered.faa',
                'L4': 'rpL4_JGI_MDM.filtered.faa',
                'L5': 'rpL5_JGI_MDM.filtered.faa',
                'L6': 'rpL6_JGI_MDM.filtered.faa',
                'S10': 'rpS10_JGI_MDM.filtered.faa',
                'S17': 'rpS17_JGI_MDM.filtered.faa',
                'S19': 'rpS19_JGI_MDM.filtered.faa',
                'S3': 'rpS3_JGI_MDM.filtered.faa',
                'S8': 'rpS8_JGI_MDM.filtered.faa'}
    protein_databases = {key: '%s/%s' % (databases, database) \
            for key, database in list(protein_databases.items())}
    return proteins, protein_databases
1187703
def main():
    
    parser = argparse.ArgumentParser(
        description='monoseq: pretty-printing DNA and protein sequences',
        epilog='If INPUT is in FASTA format, each record is pretty-printed '
        'after printing its name and ANNOTATION (if supplied) is used by '
        'matching chromosome/record name. If INPUT contains a raw sequence, '
        'only the first chromosome in ANNOTATION is used.')
    parser.add_argument(
        'sequence_file', metavar='INPUT', nargs='?', default=sys.stdin,
        type=argparse.FileType('r'), help='file to read sequence(s) from, '
        'can be in FASTA format (default: standard input)')
    parser.add_argument(
        '-b', '--block-length', metavar='LENGTH', dest='block_length',
        type=int, default=10, help='block length in letters (default: 10)')
    parser.add_argument(
        '-l', '--blocks-per-line', metavar='BLOCKS', dest='blocks_per_line',
        type=int, default=6, help='blocks per line (default: 6)')
    parser.add_argument(
        '-a', '--annotation', metavar='POS', dest='annotation', nargs=2,
        action='append', type=int, help='first and last positions of '
        'subsequence to annotate (allowed more than once)')
    parser.add_argument(
        '-e', '--bed', metavar='ANNOTATION', dest='annotation_file',
        type=argparse.FileType('r'), help='file to read annotation from in '
        'BED format')

    args = parser.parse_args()
    pprint(_until_eof(args.sequence_file), annotation=args.annotation,
           annotation_file=args.annotation_file,
           block_length=args.block_length,
           blocks_per_line=args.blocks_per_line)
790806
def retrieve_data_from_rcsb(cls, ligand_code, pdb_id = None, silent = True, cached_dir = None):
        ''''''
        if not silent:
            colortext.printf("Retrieving data from RCSB")
        if cached_dir:
            assert(os.path.exists(cached_dir))

        ligand_info_path, ligand_info, pdb_ligand_info, pdb_ligand_info_path = None, None, None, None
        if cached_dir:
            ligand_info_path = os.path.join(cached_dir, '{0}.cif'.format(ligand_code))
            if os.path.exists(ligand_info_path):
                ligand_info = read_file(ligand_info_path)
        if not ligand_info:
            ligand_info = retrieve_ligand_cif(ligand_code)
            if cached_dir:
                write_file(ligand_info_path, ligand_info)

        # Parse .cif
        l = cls(ligand_code)
        l.parse_cif(ligand_info)
        l.pdb_id = pdb_id or l.pdb_id
        has_pdb_id = l.pdb_id and (len(l.pdb_id) == 4) and (l.pdb_id != '?')  # the last case is unnecessary and will be short-cut but I included it to show possible values

        # Parse PDB XML
        if has_pdb_id:
            if cached_dir:
                pdb_ligand_info_path = os.path.join(cached_dir, '{0}.pdb.ligandinfo'.format(l.pdb_id.lower()))
                if os.path.exists(pdb_ligand_info_path):
                    pdb_ligand_info = read_file(pdb_ligand_info_path)
                else:
                    pdb_ligand_info = retrieve_pdb_ligand_info(l.pdb_id)
                    write_file(pdb_ligand_info_path, pdb_ligand_info)
            else:
                pdb_ligand_info = retrieve_pdb_ligand_info(l.pdb_id)
        if pdb_ligand_info:
            l.parse_pdb_ligand_info(pdb_ligand_info)

        # Retrive the diagram image
        l.get_diagram()

        return l
1330502
def parse(self, limit=None):
        
        if limit is not None:
            LOG.info("Only parsing first %d rows", limit)

        protein_paths = self._get_file_paths(self.tax_ids, 'protein_links')
        col = ['NCBI taxid', 'entrez', 'STRING']
        for taxon in protein_paths:
            ensembl = Ensembl(self.graph_type, self.are_bnodes_skized)
            string_file_path = '/'.join((
                self.rawdir, protein_paths[taxon]['file']))

            with gzip.open(string_file_path, 'rb') as reader:
                dataframe = pd.read_csv(reader, sep=r'\s+')
            p2gene_map = dict()

            if taxon in self.id_map_files:
                LOG.info("Using string provided id_map files")
                map_file = '/'.join((self.rawdir, self.id_map_files[taxon]['file']))

                with gzip.open(map_file, 'rt') as reader:
                    line = next(reader).strip()
                    if line != '# NCBI taxid / entrez / STRING':
                        LOG.error(
                            'Expected Headers:\t%s\nRecived Headers:\t%s\n', col, line)
                        exit(-1)

                    for line in reader.readlines():
                        row = line.rstrip('\n').split('\t')
                        # tax = row[col.index(''NCBI taxid')].strip()
                        gene = row[col.index('entrez')].strip()
                        prot = row[col.index('STRING')].strip()

                        genes = gene.split('|')
                        p2gene_map[prot.replace(taxon + '.', '')] = [
                            "NCBIGene:" + entrez_id for entrez_id in genes]
            else:
                LOG.info("Fetching ensembl proteins for taxon %s", taxon)
                p2gene_map = ensembl.fetch_protein_gene_map(taxon)
                for key in p2gene_map:
                    for phen, gene in enumerate(p2gene_map[key]):
                        p2gene_map[key][phen] = "ENSEMBL:{}".format(gene)

            LOG.info(
                "Finished fetching ENSP ID mappings, fetched %i proteins",
                len(p2gene_map))

            LOG.info(
                "Fetching protein protein interactions for taxon %s", taxon)

            self._process_protein_links(dataframe, p2gene_map, taxon, limit)
285676
def read_sbml(filename):
    
    import libsbml

    if not os.path.exists(filename):
        raise IOError('File {0!r} does not exist'.format(filename))

    reader = libsbml.SBMLReader()
    document = reader.readSBML(filename)

    sbml_model = document.getModel()
    if not sbml_model:
        raise ValueError('Cannot parse SBML model from {0!r}'.format(filename))

    species = sympy.symbols([s.getId() for s in sbml_model.getListOfSpecies()])
    initial_conditions = [s.getInitialConcentration() for s in sbml_model.getListOfSpecies()]
    compartments = sympy.symbols([s.getId() for s in sbml_model.getListOfCompartments()])
    compartment_sizes = [s.getSize() for s in sbml_model.getListOfCompartments()]

    reactions = map(_parse_reaction, sbml_model.getListOfReactions())

    # getListOfParameters is an attribute of the model for SBML Level 1&2
    parameters_with_values = [(sympy.Symbol(p.getId()), p.getValue())
                              for p in sbml_model.getListOfParameters()]
    parameter_values = dict(parameters_with_values)
    parameters = map(lambda x: x[0], parameters_with_values)
    if not parameters:
        track_local_parameters = True
        parameters = set()
        parameter_values = {}
    else:
        track_local_parameters = False


    stoichiometry_matrix = np.zeros((len(species), len(reactions)), dtype=int)
    propensities = []
    for reaction_index, reaction in enumerate(reactions):
        if track_local_parameters:
            for param, value in reaction.parameters:
                parameters.add(param)
                parameter_values[param] = value
        reactants = reaction.reactants
        products = reaction.products
        propensities.append(reaction.propensity)
        for species_index, species_id in enumerate(species):
            net_stoichiometry = products.get(species_id, 0) - reactants.get(species_id, 0)
            stoichiometry_matrix[species_index, reaction_index] = net_stoichiometry


    if track_local_parameters:
        # sympy does not allow sorting its parameter lists by default,
        # explicitly tell to sort by str representation
        sorted_parameters = sorted(parameters, key=str)
    else:
        sorted_parameters = parameters

    parameter_values_list = [parameter_values[p] for p in sorted_parameters]

    # We need to concatenate compartment names and parameters as in our framework we cannot differentiate the two
    compartments_and_parameters = compartments + sorted_parameters
    parameter_values_list = compartment_sizes + parameter_values_list

    model = Model(species, compartments_and_parameters, propensities, stoichiometry_matrix)

    return model, parameter_values_list, initial_conditions
347100
def _biotype_lookup_fn(gtf):
    
    db = get_gtf_db(gtf)
    sources = set([feature.source for feature in db.all_features()])
    gene_biotypes = set([feature.attributes.get("gene_biotype", [None])[0]
                         for feature in db.all_features()])
    biotypes = set([feature.attributes.get("biotype", [None])[0]
                    for feature in db.all_features()])
    if "protein_coding" in sources:
        return lambda feature: feature.source
    elif "protein_coding" in biotypes:
        return lambda feature: feature.attributes.get("biotype", [None])[0]
    elif "protein_coding" in gene_biotypes:
        return lambda feature: feature.attributes.get("gene_biotype", [None])[0]
    else:
        return None
553483
def parse_uniprot_xml_metadata(sr):
    
    # TODO: What about "reviewed" status? and EC number
    xref_dbs_to_keep = ['GO', 'KEGG', 'PDB', 'PROSITE', 'Pfam', 'RefSeq']

    infodict = {}

    infodict['alt_uniprots'] = list(set(sr.annotations['accessions']).difference([sr.id]))
    infodict['gene_name'] = None
    if 'gene_name_primary' in sr.annotations:
        infodict['gene_name'] = sr.annotations['gene_name_primary']
    infodict['description'] = sr.description
    infodict['taxonomy'] = None
    if 'organism' in sr.annotations:
        infodict['taxonomy'] = sr.annotations['organism']
    infodict['seq_version'] = sr.annotations['sequence_version']
    infodict['seq_date'] = sr.annotations['sequence_modified']
    infodict['entry_version'] = sr.annotations['version']
    infodict['entry_date'] = sr.annotations['modified']

    tmp = defaultdict(list)
    for xref in sr.dbxrefs:
        database = xref.split(':', 1)[0]
        xrefs = xref.split(':', 1)[-1]

        if database in xref_dbs_to_keep:
            if database == 'PDB':
                tmp['pdbs'].append(xrefs)
            else:
                tmp[database.lower()].append(xrefs)

    infodict.update(tmp)

    return infodict
188179
def create_psm_lookup(fn, fastafn, mapfn, header, pgdb, unroll=False,
                      specfncol=None, decoy=False,
                      fastadelim=None, genefield=None):
    
    proteins = store_proteins_descriptions(pgdb, fastafn, fn, mapfn, header,
                                           decoy, fastadelim, genefield)
    mzmlmap = pgdb.get_mzmlfile_map()
    sequences = {}
    for psm in tsvreader.generate_tsv_psms(fn, header):
        seq = tsvreader.get_psm_sequence(psm, unroll)
        sequences[seq] = 1
    pgdb.store_pepseqs(((seq,) for seq in sequences))
    pepseqmap = pgdb.get_peptide_seq_map()
    psms = []
    for row, psm in enumerate(tsvreader.generate_tsv_psms(fn, header)):
        specfn, psm_id, scan, seq, score = tsvreader.get_psm(psm, unroll,
                                                             specfncol)
        if len(psms) % DB_STORE_CHUNK == 0:
            pgdb.store_psms(psms)
            psms = []
        psms.append({'rownr': row,
                     'psm_id': psm_id,
                     'seq': pepseqmap[seq],
                     'score': score,
                     'specfn': mzmlmap[specfn],
                     'scannr': scan,
                     'spec_id': '{}_{}'.format(mzmlmap[specfn], scan),
                     })
    pgdb.store_psms(psms)
    pgdb.index_psms()
    store_psm_protein_relations(fn, header, pgdb, proteins)
971976
def load_pdb(filename):
    
    with open(filename) as f:
        numbers = []
        coordinates = []
        occupancies = []
        betas = []
        for line in f:
            if line.startswith("ATOM"):
                symbol = line[76:78].strip()
                numbers.append(periodic[symbol].number)
                coordinates.append([float(line[30:38])*angstrom, float(line[38:46])*angstrom, float(line[46:54])*angstrom])
                occupancies.append(float(line[54:60]))
                betas.append(float(line[60:66]))
    if len(numbers) > 0:
        molecule = Molecule(numbers, coordinates)
        molecule.occupancies = np.array(occupancies)
        molecule.betas = np.array(betas)
        return molecule
    else:
        raise FileFormatError("No molecule found in pdb file %s" % filename)
674996
def build_blast_db(blastdb)
      Dir.glob(File.expand_path(blastdb) + "/*").each do |d|
        if File.directory?(d)
          contigs  = File.join(d, "contigs.fa")
          peptides = File.join(d, "peptides.fa")

          found = false

          if File.exists?(contigs) && File.readable?(contigs)
            execute_makeblastdb("nucl", d, contigs)
            found = true
          end
          if File.exists?(peptides) && File.readable?(peptides)
            execute_makeblastdb("prot", d, peptides)
            found = true
          end

          unless found
            raise "Extracted data not found for #{contigs} or #{peptides}. " <<
            "Make sure you supplied the correct data directory and file names."
          end
        end
      end
    end
845456
def get_atlas_summary_df(self):
        
        all_info = []
        for g in self.reference_gempro.genes_with_a_representative_sequence:
            info = {}
            info['Gene_ID'] = g.id
            info['Gene_name'] = g.name

            # Protein object
            p = g.protein
            info['Protein_sequences'] = len(p.sequences)
            info['Protein_structures'] = len(p.structures)

            # SeqProp
            rseq = p.representative_sequence
            info['RepSeq_ID'] = rseq.id
            info['RepSeq_sequence_length'] = rseq.seq_len
            info['RepSeq_num_sequence_alignments'] = len([x for x in p.sequence_alignments if x.annotations['ssbio_type'] == 'seqalign'])
            info['RepSeq_num_structure_alignments'] = len([x for x in p.sequence_alignments if x.annotations['ssbio_type'] == 'structalign'])

            # SeqRecord annotations (properties calculated that summarize the whole sequence)
            for annotation_name, annotation in rseq.annotations.items():
                info['RepSeq_' + annotation_name] = annotation

            # SeqRecord alignment annotations
            all_num_mutations = []
            all_num_deletions = []
            all_len_deletions = []
            all_num_insertions = []
            all_len_insertions = []
            all_percent_identity = []
            all_percent_similarity = []
            for aln in p.sequence_alignments:
                # Gather the strain speicific stuff
                if '{}_'.format(p.id) not in aln.annotations['b_seq']:
                    continue
                info[aln.annotations['b_seq'].split('{}_'.format(p.id))[1]] = aln.annotations['percent_identity']

                # Gather the percent identities/similarities
                all_percent_identity.append(aln.annotations['percent_identity'])
                all_percent_similarity.append(aln.annotations['percent_similarity'])

                # Gather the number of residues that are mutated (filter for different mutations of same residue)
                num_mutations = len(list(set([x[1] for x in aln.annotations['mutations']])))
                all_num_mutations.append(num_mutations)

                # Gather the number of deletions as well as the length of the deletion
                if not aln.annotations['deletions']:
                    num_deletions = 0
                    len_deletions = [0]
                else:
                    num_deletions = len(aln.annotations['deletions'])
                    len_deletions = [x[1] for x in aln.annotations['deletions']]
                all_num_deletions.append(num_deletions)
                # Get the total length of the deletion for this one strain
                avg_len_deletions = np.sum(len_deletions)
                all_len_deletions.append(avg_len_deletions)

                # Gather the number of insertions as well as the length of the insertion
                if not aln.annotations['insertions']:
                    num_insertions = 0
                    len_insertions = [0]
                else:
                    num_insertions = len(aln.annotations['insertions'])
                    len_insertions = [x[1] for x in aln.annotations['insertions']]
                all_num_insertions.append(num_insertions)
                # Get the total length of insertion for this one strain
                avg_len_insertions = np.sum(len_insertions)
                all_len_insertions.append(avg_len_insertions)

            info['ATLAS_mean_num_mutations'] = np.mean(all_num_mutations)
            info['ATLAS_mean_num_deletions'] = np.mean(all_num_deletions)
            info['ATLAS_mean_len_deletions'] = np.mean(all_len_deletions)
            info['ATLAS_mean_num_insertions'] = np.mean(all_num_insertions)
            info['ATLAS_mean_len_insertions'] = np.mean(all_len_insertions)
            info['ATLAS_mean_percent_identity'] = np.mean(all_percent_identity)
            info['ATLAS_mean_percent_similarity'] = np.mean(all_percent_similarity)

            # Other mutation analysis
            single, fingerprint = p.sequence_mutation_summary()

            # Mutations that show up in more than 10% of strains
            singles = []
            for k, v in single.items():
                k = [str(x) for x in k]
                if len(v) / len(p.sequence_alignments) >= 0.01:
                    singles.append(''.join(k))  # len(v) is the number of strains which have this mutation
            info['ATLAS_popular_mutations'] = ';'.join(singles)

            # Mutation groups that show up in more than 10% of strains
            allfingerprints = []
            for k, v in fingerprint.items():
                if len(v) / len(p.sequence_alignments) >= 0.01:
                    fingerprints = []
                    for m in k:
                        y = [str(x) for x in m]
                        fingerprints.append(''.join(y))
                    allfingerprints.append('-'.join(fingerprints))
            info['ATLAS_popular_mutation_groups'] = ';'.join(allfingerprints)

            # StructProp
            rstruct = p.representative_structure
            if rstruct:
                if rstruct.structure_file:
                    info['RepStruct_ID'] = rstruct.id
                    info['RepStruct_is_experimental'] = rstruct.is_experimental
                    info['RepStruct_description'] = rstruct.description
                    info['RepStruct_repseq_coverage'] = p.representative_chain_seq_coverage

                    # ChainProp
                    rchain = p.representative_chain
                    info['RepChain_ID'] = rchain

                    # ChainProp SeqRecord annotations
                    rchain_sr = rstruct.chains.get_by_id(rchain).seq_record
                    for annotation_name, annotation in rchain_sr.annotations.items():
                        info['RepChain_' + annotation_name] = annotation

            all_info.append(info)

        cols = ['Gene_ID', 'Gene_name', 'Protein_sequences', 'Protein_structures',
                'RepSeq_ID', 'RepSeq_sequence_length',
                'RepSeq_num_sequence_alignments', 'RepSeq_num_structure_alignments',
                'RepStruct_ID', 'RepChain_ID', 'RepStruct_description',
                'RepStruct_is_experimental', 'RepStruct_repseq_coverage',
                'ATLAS_mean_percent_identity', 'ATLAS_mean_percent_similarity', 'ATLAS_mean_num_mutations',
                'ATLAS_popular_mutations', 'ATLAS_popular_mutation_groups', 'ATLAS_mean_num_deletions',
                'ATLAS_mean_num_insertions', 'ATLAS_mean_len_deletions', 'ATLAS_mean_len_insertions',
                'RepSeq_aromaticity', 'RepSeq_instability_index', 'RepSeq_isoelectric_point', 'RepSeq_molecular_weight',
                'RepSeq_monoisotopic', 'RepSeq_num_tm_helix-tmhmm', 'RepSeq_percent_acidic', 'RepSeq_percent_aliphatic',
                'RepSeq_percent_aromatic', 'RepSeq_percent_B-sspro8', 'RepSeq_percent_basic',
                'RepSeq_percent_buried-accpro', 'RepSeq_percent_buried-accpro20', 'RepSeq_percent_C-sspro',
                'RepSeq_percent_C-sspro8', 'RepSeq_percent_charged', 'RepSeq_percent_E-sspro',
                'RepSeq_percent_E-sspro8', 'RepSeq_percent_exposed-accpro', 'RepSeq_percent_exposed-accpro20',
                'RepSeq_percent_G-sspro8', 'RepSeq_percent_H-sspro', 'RepSeq_percent_H-sspro8',
                'RepSeq_percent_helix_naive', 'RepSeq_percent_I-sspro8', 'RepSeq_percent_non-polar',
                'RepSeq_percent_polar', 'RepSeq_percent_S-sspro8', 'RepSeq_percent_small',
                'RepSeq_percent_strand_naive', 'RepSeq_percent_T-sspro8', 'RepSeq_percent_tiny',
                'RepSeq_percent_turn_naive', 'RepChain_percent_B-dssp', 'RepChain_percent_C-dssp',
                'RepChain_percent_E-dssp', 'RepChain_percent_G-dssp', 'RepChain_percent_H-dssp',
                'RepChain_percent_I-dssp', 'RepChain_percent_S-dssp', 'RepChain_percent_T-dssp',
                'RepChain_SSBOND-biopython']
        cols.extend([x.id for x in self.strains])

        df_atlas_summary = pd.DataFrame(all_info, columns=cols)
        # Drop columns that don't have anything in them
        df_atlas_summary.dropna(axis=1, how='all', inplace=True)

        return df_atlas_summary
1313511
def parse_kegg_gene_metadata(infile):
    
    metadata = defaultdict(str)

    with open(infile) as mf:
        kegg_parsed = bs_kegg.parse(mf.read())

    # TODO: additional fields can be parsed

    if 'DBLINKS' in kegg_parsed.keys():
        if 'UniProt' in kegg_parsed['DBLINKS']:
            unis = str(kegg_parsed['DBLINKS']['UniProt']).split(' ')
            # TODO: losing other uniprot ids by doing this
            if isinstance(unis, list):
                metadata['uniprot'] = unis[0]
            else:
                metadata['uniprot'] = unis
        if 'NCBI-ProteinID' in kegg_parsed['DBLINKS']:
            metadata['refseq'] = str(kegg_parsed['DBLINKS']['NCBI-ProteinID'])
    if 'STRUCTURE' in kegg_parsed.keys():
        metadata['pdbs'] = str(kegg_parsed['STRUCTURE']['PDB']).split(' ')
    else:
        metadata['pdbs'] = None
    if 'ORGANISM' in kegg_parsed.keys():
        metadata['taxonomy'] = str(kegg_parsed['ORGANISM'])

    return metadata
550438
def get_yeast_gene_location(gene_name):
    '''

    '''
    from intermine.webservice import Service
    service = Service('http://yeastmine.yeastgenome.org/yeastmine/service')

    # Get a new query on the class (table) you will be querying:
    query = service.new_query('Gene')

    # The view specifies the output columns
    query.add_view('primaryIdentifier', 'secondaryIdentifier', 'symbol',
                   'name', 'organism.shortName',
                   'chromosome.primaryIdentifier',
                   'chromosomeLocation.start', 'chromosomeLocation.end',
                   'chromosomeLocation.strand')

    # Uncomment and edit the line below (the default) to select a custom sort
    # order:
    # query.add_sort_order('Gene.primaryIdentifier', 'ASC')

    # You can edit the constraint values below
    query.add_constraint('organism.shortName', '=', 'S. cerevisiae',
                         code='B')
    query.add_constraint('Gene', 'LOOKUP', gene_name, code='A')

    # Uncomment and edit the code below to specify your own custom logic:
    # query.set_logic('A and B')
    chromosomes = {'chrI': 1,
                   'chrII': 2,
                   'chrIII': 3,
                   'chrIV': 4,
                   'chrV': 5,
                   'chrVI': 6,
                   'chrVII': 7,
                   'chrVIII': 8,
                   'chrIX': 9,
                   'chrX': 10,
                   'chrXI': 11,
                   'chrXII': 12,
                   'chrXIII': 13,
                   'chrXIV': 14,
                   'chrXV': 15,
                   'chrXVI': 16}
    first_result = query.rows().next()

    return [chromosomes[first_result['chromosome.primaryIdentifier']],
            first_result['chromosomeLocation.start'],
            first_result['chromosomeLocation.end'],
            int(first_result['chromosomeLocation.strand'])]
911466
@Override
	public void parseHeader(String header, S sequence) {
		//uniptrot
		// tr|Q0TET7|Q0TET7_ECOL5 Putative uncharacterized protein OS=Escherichia coli O6:K15:H31 (strain 536 / UPEC) GN=ECP_2553 PE=4 SV=1
		sequence.setOriginalHeader(header);
		String[] data = getHeaderValues(header);

		if (data.length == 1) {
			sequence.setAccession(new AccessionID(data[0]));
		} else  if (data[0].equalsIgnoreCase("sp") || data[0].equalsIgnoreCase("tr")) {
			if (data[0].equalsIgnoreCase("sp")) {
				sequence.setAnnotationType(AnnotationType.CURATED);
			} else {
				sequence.setAnnotationType(AnnotationType.PREDICTED);
			}

			sequence.setAccession(new AccessionID(data[1], DataSource.UNIPROT));
			if (data.length > 2) {
				sequence.setDescription(data[2]);
			}

		} else if (data[0].equalsIgnoreCase("gi")) {
			DataSource giSource = DataSource.UNKNOWN;
			if (data.length >= 3) {
				if (data[2].equalsIgnoreCase("gb")) {
					giSource = DataSource.GENBANK;
				} else if (data[2].equalsIgnoreCase("emb")) {
					giSource = DataSource.ENA;
				} else if (data[2].equalsIgnoreCase("dbj")) {
					giSource = DataSource.DDBJ;
				}
				sequence.setAccession(new AccessionID(data[3], giSource));
			} else {
				sequence.setAccession(new AccessionID(header, giSource));
			}
		} else if (data[0].equalsIgnoreCase("pir")) {
			sequence.setAccession(new AccessionID(data[2], DataSource.NBRF));
		} else if (data[0].equalsIgnoreCase("prf")) {
			sequence.setAccession(new AccessionID(data[2], DataSource.PRF));
		} else if (data[0].equalsIgnoreCase("pdb")) {
			sequence.setAccession(new AccessionID(data[1] + ":" + data[2], DataSource.PDB1));
		} else if (data[0].startsWith("PDB")) {
			String[] pdbe = data[0].split(" ");
			String[] pdbaccession = pdbe[0].split(":");
			sequence.setAccession(new AccessionID(pdbaccession[1], DataSource.PDBe));
		} else if (data[0].indexOf(":") != -1 && data.length > 1 && data[1].equals("PDBID")) {
			sequence.setAccession(new AccessionID(data[0], DataSource.PDB2));
		} else if (data[0].equalsIgnoreCase("pat")) {
			sequence.setAccession(new AccessionID(data[2], DataSource.PATENTS));
		} else if (data[0].equalsIgnoreCase("bbs")) {
			sequence.setAccession(new AccessionID(data[1], DataSource.GENINFO));
		} else if (data[0].equalsIgnoreCase("gnl")) {
			sequence.setAccession(new AccessionID(data[2], DataSource.GENERAL));
		} else if (data[0].equalsIgnoreCase("ref")) {
			sequence.setAccession(new AccessionID(data[1], DataSource.NCBI));
		} else if (data[0].equalsIgnoreCase("lcl")) {
			sequence.setAccession(new AccessionID(data[1], DataSource.LOCAL));
		} else {
			sequence.setAccession(new AccessionID(data[0])); // avoid the common problem of picking up all the comments original header in getOriginalHeader
		}


	}
1052698
def df_proteins(self):
        """"""
        pre_df = []
        df_cols = ['gene', 'id', 'sequences', 'num_sequences', 'representative_sequence',
                   'repseq_gene_name', 'repseq_uniprot', 'repseq_description',
                   'num_structures', 'experimental_structures', 'num_experimental_structures',
                   'homology_models', 'num_homology_models',
                   'representative_structure', 'representative_chain', 'representative_chain_seq_coverage',
                   'repstruct_description', 'repstruct_resolution',
                   'num_sequence_alignments', 'num_structure_alignments']

        for g in self.genes:
            # Get per protein DataFrame
            protein_dict = g.protein.protein_statistics
            protein_dict['gene'] = g.id
            pre_df.append(protein_dict)

        df = pd.DataFrame.from_records(pre_df, columns=df_cols).set_index('gene')
        if df.empty:
            log.warning('Empty dataframe')
            return df
        else:
            return ssbio.utils.clean_df(df)
1170579
def process_catalog(self, limit=None):
        
        raw = '/'.join((self.rawdir, self.files['catalog']['file']))
        LOG.info("Processing Data from %s", raw)
        efo_ontology = RDFGraph(False, "EFO")
        LOG.info("Loading EFO ontology in separate rdf graph")
        efo_ontology.parse(self.files['efo']['url'], format='xml')
        efo_ontology.bind_all_namespaces()
        LOG.info("Finished loading EFO ontology")

        so_ontology = RDFGraph(False, "SO")
        LOG.info("Loading SO ontology in separate rdf graph")
        so_ontology.parse(self.files['so']['url'], format='xml')
        so_ontology.bind_all_namespaces()
        LOG.info("Finished loading SO ontology")

        with open(raw, 'r', encoding="iso-8859-1") as csvfile:
            filereader = csv.reader(csvfile, delimiter='\t')
            header = next(filereader, None)  # the header row
            header_len = len(header)
            LOG.info('header length:\t %i', header_len)

            for row in filereader:
                if not row:
                    pass
                else:
                    if header_len != len(row):
                        LOG.error('BadRow: %i has %i columns', filereader.line_num, row)

                    (date_added_to_catalog,
                     pubmed_num,
                     first_author,
                     pub_date,
                     journal,
                     link,
                     study_name,
                     disease_or_trait,
                     initial_sample_description,
                     replicate_sample_description,
                     region,
                     chrom_num,
                     chrom_pos,
                     reported_gene_nums,
                     mapped_gene,
                     upstream_gene_num,
                     downstream_gene_num,
                     snp_gene_nums,
                     upstream_gene_distance,
                     downstream_gene_distance,
                     strongest_snp_risk_allele,
                     snps,
                     merged,
                     snp_id_current,
                     context,
                     intergenic_flag,
                     risk_allele_frequency,
                     pvalue,
                     pvalue_mlog,
                     pvalue_text,
                     or_or_beta,
                     confidence_interval_95,
                     platform_with_snps_passing_qc,
                     cnv_flag,
                     mapped_trait,
                     mapped_trait_uri,
                     study_accession,
                     GENOTYPING_TECHNOLOGY
                    ) = row

                    if self.test_mode:
                        continue

# 06-May-2015	25917933
#   Zai CC	20-Nov-2014	J Psychiatr Res	http://europepmc.org/abstract/MED/25917933
# A genome-wide association study of suicide severity scores in bipolar disorder.
# Suicide in bipolar disorder
# 959 European ancestry individuals	NA
# 10p11.22	10	32704340	C10orf68, CCDC7, ITGB1	CCDC7
# rs7079041-A	rs7079041	0	7079041	intron	0		2E-6	5.698970

                    variant_curie, variant_type = self._get_curie_and_type_from_id(
                        strongest_snp_risk_allele)

                    if strongest_snp_risk_allele.strip() == '':
                        LOG.debug(
                            "No strongest SNP risk allele for %s:\n%s",
                            pubmed_num, str(row))
                        # still consider adding in the EFO terms
                        # for what the study measured?
                        continue

                    if variant_type == 'snp':
                        self._add_snp_to_graph(
                            variant_curie, strongest_snp_risk_allele, chrom_num,
                            chrom_pos, context, risk_allele_frequency)

                        self._add_deprecated_snp(
                            variant_curie, snp_id_current, merged, chrom_num, chrom_pos)

                        self._add_snp_gene_relation(
                            variant_curie, snp_gene_nums, upstream_gene_num,
                            downstream_gene_num)
                    elif variant_type == 'haplotype':
                        self._process_haplotype(
                            variant_curie, strongest_snp_risk_allele, chrom_num,
                            chrom_pos, context, risk_allele_frequency, mapped_gene,
                            so_ontology)
                    elif variant_type is None:
                        LOG.warning(
                            "There's a snp id i can't manage: %s",
                            strongest_snp_risk_allele)
                        continue

                    description = self._make_description(
                        disease_or_trait, initial_sample_description,
                        replicate_sample_description,
                        platform_with_snps_passing_qc, pvalue)

                    self._add_variant_trait_association(
                        variant_curie, mapped_trait_uri, efo_ontology,
                        pubmed_num, description)

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

        # TODO loop through the location hash,
        # and make all snps at that location equivalent
        for l in self.id_location_map:
            snp_ids = self.id_location_map[l]
            if len(snp_ids) > 1:
                LOG.info("%s has >1 snp id: %s", l, str(snp_ids))
        return
56338
def parse_rRNA(insertion, seq, gff):
    
    offset = insertion['offset']
    strand = insertion['strand']
    for rRNA in parse_masked(seq, 0)[0]:
        rRNA = ''.join(rRNA)
        Start = seq[1].find(rRNA) + 1
        End = Start + len(rRNA) - 1
        if strand == '-':
            Start, End = End - 2, Start - 2
        pos = (abs(Start + offset) - 1, abs(End + offset) - 1)
        Start, End = min(pos), max(pos)
        source = insertion['source']
        annot = '%s rRNA' % (source.split('from', 1)[0])
        gff['#seqname'].append(insertion['ID'])
        gff['source'].append(source)
        gff['feature'].append('rRNA')
        gff['start'].append(Start)
        gff['end'].append(End)
        gff['score'].append('.')
        gff['strand'].append(strand)
        gff['frame'].append('.')
        gff['attribute'].append('Name=%s' % (annot))
    return gff