# -*- coding: utf-8 -*-
# Copyright (c) 2016-2017, Zhijiang Yao, Jie Dong and Dongsheng Cao
# All rights reserved.
# This file is part of the PyBioMed.
# The contents are covered by the terms of the BSD license
# which is included in the file license.txt, found at the root
# of the PyBioMed source tree.
"""
Created on Sat Jul 13 11:18:26 2013
This module is used for downloading the PDB file from RCSB PDB web and
extract its amino acid sequence. This module can also download the protein
sequence from the uniprot (http://www.uniprot.org/) website. You can only
need input a protein ID or prepare a file (ID.txt) related to ID. You can
obtain a .txt (ProteinSequence.txt) file saving protein sequence you need.
Authors: Zhijiang Yao and Dongsheng Cao.
Date: 2016.06.04
Email: gadsby@163.com
"""
# Core Library modules
import ftplib
import gzip
import os
import signal
import string
import sys
import threading
import time
[docs]class TimeoutException(Exception):
pass
HOSTNAME = "ftp.wwpdb.org"
DIRECTORY = "/pub/pdb/data/structures/all/pdb/"
PREFIX = "pdb"
SUFFIX = ".ent.gz"
# List of three and one letter amino acid codes
_aa_index = [
("ALA", "A"),
("CYS", "C"),
("ASP", "D"),
("GLU", "E"),
("PHE", "F"),
("GLY", "G"),
("HIS", "H"),
("HSE", "H"),
("HSD", "H"),
("ILE", "I"),
("LYS", "K"),
("LEU", "L"),
("MET", "M"),
("MSE", "M"),
("ASN", "N"),
("PRO", "P"),
("GLN", "Q"),
("ARG", "R"),
("SER", "S"),
("THR", "T"),
("VAL", "V"),
("TRP", "W"),
("TYR", "Y"),
]
# AA3_TO_AA1 = dict(_aa_index)
# AA1_TO_AA3 = dict([(aa[1],aa[0]) for aa in _aa_index])
[docs]class InterruptableThread(threading.Thread):
def __init__(self, func, *args, **kwargs):
threading.Thread.__init__(self)
self._func = func
self._args = args
self._kwargs = kwargs
self._result = None
[docs] def run(self):
self._result = self._func(*self._args, **self._kwargs)
@property
def result(self):
return self._result
[docs]class timelimited(object):
def __init__(self, sec):
self._sec = sec
def __call__(self, f):
def wrapped_f(*args, **kwargs):
it = InterruptableThread(f, *args, **kwargs)
it.start()
it.join(self._sec)
if not it.is_alive():
return it.result
raise TimeoutException("Network timeout, skip!!")
return wrapped_f
[docs]def unZip(some_file, some_output):
"""
Unzip some_file using the gzip library and write to some_output.
"""
f = gzip.open(some_file, "r")
g = open(some_output, "w")
g.writelines(f.readlines())
f.close()
g.close()
os.remove(some_file)
[docs]def pdbDownload(
file_list, hostname=HOSTNAME, directory=DIRECTORY, prefix=PREFIX, suffix=SUFFIX
):
"""
Download all pdb files in file_list and unzip them.
"""
success = True
# Log into server
print("Connecting...")
ftp = ftplib.FTP()
ftp.connect(hostname)
ftp.login()
# Remove .pdb extensions from file_list
for file_index, file in enumerate(file_list):
try:
file_list[file_index] = file[: file.index(".pdb")]
except ValueError:
pass
# Download all files in file_list
to_get = ["%s/%s%s%s" % (directory, prefix, f, suffix) for f in file_list]
to_write = ["%s%s" % (f, suffix) for f in file_list]
for i in range(len(to_get)):
try:
ftp.retrbinary("RETR %s" % to_get[i], open(to_write[i], "wb").write)
final_name = "%s.pdb" % to_write[i][: to_write[i].index(".")]
unZip(to_write[i], final_name)
print("%s retrieved successfully." % final_name)
except ftplib.error_perm:
os.remove(to_write[i])
print("ERROR! %s could not be retrieved!" % file_list[i])
success = False
# Log out
ftp.quit()
if success:
return True
else:
return False
[docs]def GetPDB(pdbidlist=[]):
"""
Download the PDB file from PDB FTP server by providing a list of pdb id.
"""
for i in pdbidlist:
templist = []
templist.append(i)
pdbDownload(templist)
return True
[docs]def pdbSeq(pdb, use_atoms=False):
"""
Parse the SEQRES entries in a pdb file. If this fails, use the ATOM
entries. Return dictionary of sequences keyed to chain and type of
sequence used.
"""
# Try using SEQRES
seq = [l for l in pdb if l[0:6] == "SEQRES"]
if len(seq) != 0 and not use_atoms:
seq_type = "SEQRES"
chain_dict = {l[11]: [] for l in seq}
for c in chain_dict.keys():
chain_seq = [l[19:70].split() for l in seq if l[11] == c]
for x in chain_seq:
chain_dict[c].extend(x)
# Otherwise, use ATOM
else:
seq_type = "ATOM "
# Check to see if there are multiple models. If there are, only look
# at the first model.
models = [i for i, l in enumerate(pdb) if l.startswith("MODEL")]
if len(models) > 1:
pdb = pdb[models[0] : models[1]]
# Grab all CA from ATOM entries, as well as MSE from HETATM
atoms = []
for l in pdb:
if l[0:6] == "ATOM " and l[13:16] == "CA ":
# Check to see if this is a second conformation of the previous
# atom
if len(atoms) != 0:
if atoms[-1][17:26] == l[17:26]:
continue
atoms.append(l)
elif l[0:6] == "HETATM" and l[13:16] == "CA " and l[17:20] == "MSE":
# Check to see if this is a second conformation of the previous
# atom
if len(atoms) != 0:
if atoms[-1][17:26] == l[17:26]:
continue
atoms.append(l)
chain_dict = {l[21]: [] for l in atoms}
for c in chain_dict.keys():
chain_dict[c] = [l[17:20] for l in atoms if l[21] == c]
AA3_TO_AA1 = dict(_aa_index)
tempchain = chain_dict.keys()
tempchain.sort()
seq = ""
for i in tempchain:
for j in chain_dict[i]:
if j in AA3_TO_AA1.keys():
seq = seq + (AA3_TO_AA1[j])
else:
seq = seq + ("X")
return seq, seq_type
try:
# Python 3
from urllib.request import urlopen
except ImportError:
# Python 2
from urllib2 import urlopen
##################################################################################################
[docs]def GetProteinSequence(ProteinID):
"""
#########################################################################################
Get the protein sequence from the uniprot website by ID.
Usage:
result=GetProteinSequence(ProteinID)
Input: ProteinID is a string indicating ID such as "P48039".
Output: result is a protein sequence.
#########################################################################################
"""
ID = str(ProteinID)
localfile = urlopen("http://www.uniprot.org/uniprot/" + ID + ".fasta")
temp = localfile.readlines()
res = ""
for i in range(1, len(temp)):
res = res + temp[i].strip()
return res
##################################################################################################
[docs]def GetProteinSequenceFromTxt(path, openfile, savefile):
"""
#########################################################################################
Get the protein sequence from the uniprot website by the file containing ID.
Usage:
result=GetProteinSequenceFromTxt(path,openfile,savefile)
Input: path is a directory path containing the ID file such as "/home/orient/protein/"
openfile is the ID file such as "proteinID.txt"
savefile is the file saving the obtained protein sequences such as "protein.txt"
#########################################################################################
"""
f1 = file(path + savefile, "wb")
f2 = file(path + openfile, "r")
# res=[]
for index, i in enumerate(f2):
itrim = i.strip()
if itrim == "":
continue
else:
temp = GetProteinSequence(itrim)
print("--------------------------------------------------------")
print("The %d protein sequence has been downloaded!" % (index + 1))
print(temp)
f1.write(temp + "\n")
print("--------------------------------------------------------")
# res.append(temp+'\n')
# f1.writelines(res)
f2.close()
f1.close()
return 0
[docs]def GetSeqFromPDB(pdbfile=""):
"""
Get the amino acids sequences from pdb file.
"""
f1 = file(pdbfile, "r")
res1, res2 = pdbSeq(f1)
f1.close()
return res1
[docs]class Seq:
def __init__(self, name, seq, no):
self.name = name
self.seq = seq.upper()
self.no = no
self.length = len(seq)
def __str__(self):
"""Output seq when 'print' method is called."""
return "%s\tNo:%s\tlength:%s\n%s" % (
self.name,
str(self.no),
str(self.length),
self.seq,
)
[docs]def IsFasta(seq):
"""
#################################################################
Judge the Seq object is in FASTA format.
Two situation:
1. No seq name.
2. Seq name is illegal.
3. No sequence.
:param seq: Seq object.
#################################################################
"""
if not seq.name:
error_info = "Error, sequence " + str(seq.no) + " has no sequence name."
print(seq)
sys.stderr.write(error_info)
return False
if -1 != seq.name.find(">"):
error_info = "Error, sequence " + str(seq.no) + " name has > character."
sys.stderr.write(error_info)
return False
if 0 == seq.length:
error_info = "Error, sequence " + str(seq.no) + " is null."
sys.stderr.write(error_info)
return False
return True
[docs]def ReadFasta(f):
"""
#################################################################
Read a fasta file.
:param f: HANDLE to input. e.g. sys.stdin, or open(<file>)
Return Seq obj list.
#################################################################
"""
name, seq = "", ""
count = 0
seq_list = []
lines = f.readlines()
for line in lines:
if not line:
break
if ">" == line[0]:
if 0 != count or (0 == count and seq != ""):
if IsFasta(Seq(name, seq, count)):
seq_list.append(seq)
else:
sys.exit(0)
seq = ""
name = line[1:].strip()
count += 1
else:
seq += line.strip()
count += 1
if IsFasta(Seq(name, seq, count)):
seq_list.append(seq)
else:
sys.exit(0)
return seq_list
####################################################################
if __name__ == "__main__":
print("-" * 10 + "START" + "-" * 10)
@timelimited(10)
def run_GetPDB():
GetPDB(["1atp", "1efz", "1f88"])
@timelimited(10)
def run_GetSeqFromPDB():
seq = GetSeqFromPDB("1atp.pdb")
print(seq)
@timelimited(10)
def run_GetProteinSequence():
GetProteinSequence("O00560")
print(ReadFasta(open("../test/test_data/protein.fasta")))
run_GetPDB()
print("-" * 25)
run_GetSeqFromPDB()
print("-" * 25)
run_GetProteinSequence()
print("-" * 10 + "END" + "-" * 10)
