Retrieve and download PFAM results from the HMMER search tool.

    Args:
        seq:
        outpath:
        searchtype:
        force_rerun:

    Returns:

    Todo:
        * Document and test!
483715
def manual_get_pfam_annotations(seq, outpath, searchtype='phmmer', force_rerun=False):
    
    if op.exists(outpath):
        with open(outpath, 'r') as f:
            json_results = json.loads(json.load(f))

    else:
        fseq = '>Seq\n' + seq
        if searchtype == 'phmmer':
            parameters = {'seqdb': 'pdb', 'seq': fseq}
        if searchtype == 'hmmscan':
            parameters = {'hmmdb': 'pfam', 'seq': fseq}
        enc_params = urllib.urlencode(parameters).encode('utf-8')
        request = urllib2.Request('http://www.ebi.ac.uk/Tools/hmmer/search/{}'.format(searchtype), enc_params)
        url = (urllib2.urlopen(request).geturl() + '?output=json')
        request = str(url)
        request_read = urlopen(request).read().decode("utf-8")

        with open(outpath, 'w') as f:
            json.dump(request_read, f)

        json_results = json.loads(request_read)

    return json_results['results']['hits']
653681
def annotate_orfs(fa, seqs, threads, threshold = 1e-6):
    
    db = '%s/pfam/Pfam-A.hmm' % os.environ['databases']
    base = fa.rsplit('.', 1)[0]
    out = '%s.pfam.hmmscan' % (base)
    tblout = '%s.pfam.hmmscan-tblout' % (base)
    orf2pfam = {} # orf2pfam[scaffold id][orf id] = pfam
    if os.path.exists(tblout) is False:
        p = subprocess.Popen('hmmscan --cpu %s --tblout %s -o %s %s %s'\
                % (threads, tblout, out, db, fa), shell = True)
        p.communicate()
    for line in open(tblout):
        if line.startswith('#'):
            continue
        line = line.strip().split()
        hit, id, e, inc = line[0], line[2], float(line[4]), int(line[17])
        contig = id.split('_', 1)[0]
        if e > threshold:
            continue
        if contig not in orf2pfam:
            orf2pfam[contig] = {}
        if id not in orf2pfam[contig]:
            orf2pfam[contig][id] = hit
    for seq, info in list(seqs.items()):
        for insertion in info[2]:
            annotations = {}
            for orf in insertion[4]:
                orf  = orf[0].split('>')[1].split()[0]
                contig = orf.split('_', 1)[0]
                if contig in orf2pfam and orf in orf2pfam[contig]:
                    annotations[orf] = orf2pfam[contig][orf]
                else:
                    annotations[orf] = 'n/a'
            insertion.append(annotations)
    return seqs
1252211
def seqs_to_stream(seqs, ih):
    
    if ih == '_input_as_multiline_string':
        recs = FastaFinder(seqs.split('\n'))
    elif ih == '_input_as_string':
        recs = FastaFinder(open(seqs))
    elif ih == '_input_as_seqs':
        recs = [['>'+str(i), s] for i, s in enumerate(seqs)]
    elif ih == '_input_as_lines':
        recs = FastaFinder(seqs)
    else:
        raise TypeError, "Unknown input handler %s" % ih
    return recs
1273575
def write_sequences_to_fasta(path, seqs):
    
    from Bio import SeqIO
    from Bio.Seq import Seq
    from Bio.SeqRecord import SeqRecord

    path = Path(path)
    records = []

    for id, seq in seqs.items():
        record = SeqRecord(Seq(seq), id=id, description='')
        records.append(record)

    SeqIO.write(records, str(path), 'fasta')
78286
def get_protein_seq_for_transcript(self, transcript_id):
        
        
        headers = {"content-type": "text/plain"}
        
        self.attempt = 0
        ext = "/sequence/id/{}?type=protein".format(transcript_id)
        
        return self.ensembl_request(ext, headers)
60111
def library_sequencing_results(self):
        
        sres_ids = self.sequencing_result_ids
        res = {}
        for i in sres_ids:
            sres = SequencingResult(i)
            res[sres.library_id] = sres
        return res
1292129
def seq(self):
        """"""

        if self.sequence_file:
            file_to_load = copy(self.sequence_path)
            log.debug('{}: reading sequence from sequence file {}'.format(self.id, file_to_load))
            tmp_sr = SeqIO.read(file_to_load, 'fasta')
            return tmp_sr.seq

        else:
            if not self._seq:
                log.debug('{}: no sequence stored in memory'.format(self.id))
            else:
                log.debug('{}: reading sequence from memory'.format(self.id))

            return self._seq
1232741
def extract_hit_list_from_hmmsearch_results(self):

        '''
        
        '''

        combined_list_of_hits = []

        for result in self.hmm_result_list:
            fullpath = self.hhsearch_result_folder + result
            

            se = SearchIO.read(fullpath, 'hmmer3-text')

            sublist = []
            for hit in se:
                combined_list_of_hits.append(hit)
                sublist.append(hit.id)
            print 'extracted', str(len(sublist)), 'hits for', result
            
        return combined_list_of_hits
90740
def get_records(self, sequence_id, gt=None, gte=None, lt=None, lte=None, limit=None,
                    query_ascending=True, results_ascending=True):
        
654724
def get_sequence(queries,fap=None,fmt='fasta',
            organism_taxid=9606,
                 test=False):
    
    url = 'http://www.ebi.ac.uk/Tools/dbfetch/dbfetch'
    params = {
    'id':' '.join(queries),
    'db':'uniprotkb',
    'organism':organism_taxid,    
    'format':fmt,
    'style':'raw',
    'Retrieve':'Retrieve',
    }
    response = requests.get(url, params=params)
    if test:
        print(response.url)
    if response.ok:
        if not fap is None:
            with open(fap,'w') as f:
                f.write(response.text)
            return fap
        else:
            return response.text            
    else:
        print('Something went wrong ', response.status_code)
560910
def get_next_rngruns(self):
        
        available_runs = [result['params']['RngRun'] for result in
                          self.get_results()]
        yield from DatabaseManager.get_next_values(available_runs)
150800
def _get_SEQRES_sequences(self):
        ''''''

        pdb_id = self.get_pdb_id()
        SEQRES_lines = self.parsed_lines["SEQRES"]

        modified_residue_mapping_3 = self.modified_residue_mapping_3
        # I commented this out since we do not need it for my current test cases
        #for k, v in self.modified_residues.iteritems():
        #    assert(v['modified_residue'] not in modified_residues)
        #    modified_residues[v['modified_residue']] = v['original_residue_3']

        for x in range(0, len(SEQRES_lines)):
            assert(SEQRES_lines[x][7:10].strip().isdigit())

        if not SEQRES_lines:
            #colortext.warning("WARNING: No SEQRES records were found. Kyle is trying to handle this gracefully, but Shane may need to fix it")
            return

        seqres_chain_order = []
        SEQRES_lines = [line[11:].rstrip() for line in SEQRES_lines] # we cannot strip the left characters as some cases e.g. 2MBP are missing chain identifiers

        # Collect all residues for all chains, remembering the chain order
        chain_tokens = {}
        for line in SEQRES_lines:
            chainID = line[0]
            if missing_chain_ids.get(self.pdb_id):
                chainID = missing_chain_ids[self.pdb_id]
            if chainID not in seqres_chain_order:
                seqres_chain_order.append(chainID)
            chain_tokens[chainID] = chain_tokens.get(chainID, [])
            chain_tokens[chainID].extend(line[6:].strip().split())

        sequences = {}
        self.chain_types = {}

        canonical_acid_types = set(residue_type_3to1_map.keys())
        canonical_acid_types.remove('UNK')
        assert(len(canonical_acid_types) == 20)

        for chain_id, tokens in chain_tokens.iteritems():

            # Determine whether the chain is DNA, RNA, or a protein chain
            # 1H38 is a good test for this - it contains DNA (chains E and G and repeated by H, K, N, J, M, P), RNA (chain F, repeated by I, L, O) and protein (chain D, repeated by A,B,C) sequences
            # 1ZC8 is similar but also has examples of DU
            # 4IHY has examples of DI (I is inosine)
            # 2GRB has RNA examples of I and U
            # 1LRP has protein chains with only CA atoms
            # This will throw an exception when a non-canonical is found which is not listed in basics.py. In that case, the list in basics.py should be updated.

            chain_type = None
            set_of_tokens = set(tokens)
            if (set_of_tokens.union(all_recognized_dna) == all_recognized_dna):# or (len(set_of_tokens) <= 5 and len(set_of_tokens.union(dna_nucleotides)) == len(set_of_tokens) + 1): # allow one unknown DNA residue
                chain_type = 'DNA'
            elif (set_of_tokens.union(all_recognized_rna) == all_recognized_rna):# or (len(set_of_tokens) <= 5 and len(set_of_tokens.union(dna_nucleotides)) == len(set_of_tokens) + 1): # allow one unknown DNA residue
                chain_type = 'RNA'
            elif len(set_of_tokens) == 1 and 'UNK' in set_of_tokens:
                chain_type = 'Unknown'
            elif not(set_of_tokens.intersection(canonical_acid_types)):
                # Zero canonical residues may imply a ligand or a heterogen chain
                chain_type = PDB._determine_heterogen_chain_type(set_of_tokens)
            else:
                assert(len(set(tokens).intersection(dna_nucleotides)) == 0)
                assert(len(set(tokens).intersection(rna_nucleotides)) == 0)
                assert(len(set(tokens).intersection(rna_nucleotides)) == 0)
                chain_type = 'Protein'
                if not self.chain_atoms.get(chain_id):
                    # possible for biological unit files
                    continue
                if self.chain_atoms[chain_id] == set(['CA']):
                    chain_type = 'Protein skeleton'

            # Get the sequence, mapping non-canonicals to the appropriate letter
            self.chain_types[chain_id] = chain_type

            sequence = []
            if chain_type == 'DNA':
                for r in tokens:
                    if dna_nucleotides_2to1_map.get(r):
                        sequence.append(dna_nucleotides_2to1_map[r])
                    else:
                        if non_canonical_dna.get(r):
                            sequence.append(non_canonical_dna[r])
                        else:
                            raise Exception("Unknown DNA residue %s." % r)
            elif chain_type == 'RNA':
                for r in tokens:
                    if r in rna_nucleotides:
                        sequence.append(r)
                    else:
                        if non_canonical_rna.get(r):
                            sequence.append(non_canonical_rna[r])
                        else:
                            raise Exception("Unknown RNA residue %s." % r)
            else:
                token_counter = 0
                for r in tokens:
                    token_counter += 1
                    if residue_type_3to1_map.get(r):
                        sequence.append(residue_type_3to1_map[r])
                    else:

                        if self.modified_residue_mapping_3.get(r):
                            sequence.append(residue_type_3to1_map[self.modified_residue_mapping_3.get(r)])
                        elif non_canonical_amino_acids.get(r):
                            #print('Mapping non-canonical residue %s to %s.' % (r, non_canonical_amino_acids[r]))
                            #print(SEQRES_lines)
                            #print(line)
                            sequence.append(non_canonical_amino_acids[r])
                        elif r == 'UNK':
                            continue
                        # Skip these residues
                        elif r == 'ACE' and token_counter == 1:
                            # Always allow ACE as the first residue of a chain
                            sequence.append('X')
                        elif r == 'ACE' and pdb_id in cases_with_ACE_residues_we_can_ignore:
                            sequence.append('X')
                            #continue
                        # End of skipped residues
                        else:
                            #print(modified_residue_mapping_3)
                            if modified_residue_mapping_3.get(r):
                                if modified_residue_mapping_3[r] == 'UNK':
                                    sequence.append('X')
                                else:
                                    assert(modified_residue_mapping_3[r] in residue_types_3)
                                    sequence.append(residue_type_3to1_map[modified_residue_mapping_3[r]])
                            else:
                                raise Exception("Unknown protein residue %s in chain %s." % (r, chain_id))
            sequences[chain_id] = "".join(sequence)

        self.seqres_chain_order = seqres_chain_order

        # Create Sequence objects for the SEQRES sequences
        for chain_id, sequence in sequences.iteritems():
            self.seqres_sequences[chain_id] = Sequence.from_sequence(chain_id, sequence, self.chain_types[chain_id])
680179
def sequences_from_fasta(path):
    
    from Bio import SeqIO
    return {x.description: x.seq for x in SeqIO.parse(path, 'fasta')}
981282
def query_sequence(self):
    
    if not self.entries.seq: return None
    if self.check_flag(0x10): return rc(self.entries.seq)
    return self.entries.seq
218926
def get_gene_seqs(database_path, gene):
    
    gene_path = database_path + "/" + gene + ".fsa"
    gene_seq = ""
    # Open fasta file
    with open(gene_path) as gene_file:
        header = gene_file.readline()
        for line in gene_file:
            seq = line.strip()
            gene_seq += seq
    return gene_seq
1130037
def pp_event(seq):
    """"""

    if isinstance(seq, Event):
        return str(seq)

    # Get the original sequence back if seq is a pretty name already
    rev_curses = dict((v, k) for k, v in CURSES_NAMES.items())
    rev_curtsies = dict((v, k) for k, v in CURTSIES_NAMES.items())
    if seq in rev_curses:
        seq = rev_curses[seq]
    elif seq in rev_curtsies:
        seq = rev_curtsies[seq]

    pretty = curtsies_name(seq)
    if pretty != seq:
        return pretty
    return repr(seq).lstrip('u')[1:-1]
754917
def fetch(self, seq_id, start=None, end=None):
        
        rec = self._db.execute("""select * from seqinfo where seq_id = ? order by added desc""", [seq_id]).fetchone()

        if rec is None:
            raise KeyError(seq_id)

        if self._writing and self._writing["relpath"] == rec["relpath"]:
            logger.warning("""Fetching from file opened for writing;
            closing first ({})""".format(rec["relpath"]))
            self.commit()

        path = os.path.join(self._root_dir, rec["relpath"])
        fabgz = self._open_for_reading(path)
        return fabgz.fetch(seq_id, start, end)