BLAST를 이용한 PDB 자동 다운로드 및 처리

PDB에서 유사한 구조, 혹은 도메인을 모두 검색해서 구조를 서로 비교해 본다면 어떻게 해야 할까? 1. PDB에 가서 2. 키워드로 검색하고 3. 나오는 PDB 파일을 하나하나 클릭해서 저장하고 4. PyMol 등에 불러와서 5. Structure alignment 를 하고..
물론 몇 개 정도의 단백질이라면 이렇게 해도 된다. 그렇지만 비슷한 구조를 가진 단백질을 다 비교해 보고 싶다면? 가령 Kinase 로 구조가 풀려있는 PDB를 모두 검색해서 이들의 구조를 다 비교해 보고 싶다든지 한다면…이거 수작업으로는 쉽지 않다.

특정 아미노산 서열을 가지고 PDB에서 BLAST 검색해서 PDB 다운로드하고, 이걸 PyMol에 로딩한후 제일 처음 파일 기준으로 Structural Alignment를 해주는 스크립트. 자세한 설명은 나중에 시간나면…안나면말구
단, Pymol 버전이 1.3이상이거나 cealign (http://pymolwiki.org/index.php/Cealign) 이 설치되어 있어야 Structural Alignment 가 됨.

#
#
# 1. Blast search of PDB using RESTful web service form www.rcsb.org
# 2. Retrieve PDB ids and PDB files
# 3. Generate PyMOL script which load all of PDBs and structural alignment using cealign (http://pymolwiki.org/index.php/Cealign)
# Written by MadScientist (https://madscientist.wordpress.com
#
import urllib2, urllib, os
import xml.etree.ElementTree as ET
#
# The first Query is BLAST search using sequence and eValue Cutoff
# and second filter is high resolution Cutoff
# the last filter is for excluding title contains "arp" (To exclude Arp 2/3 complex)
# 
url = 'http://www.rcsb.org/pdb/rest/search'
queryXML = """
<?xml version="1.0" encoding="UTF-8"?>
<orgPdbCompositeQuery version=\"1.0\">
	<queryRefinement>
		<queryRefinementLevel>0</queryRefinementLevel>
			<orgPdbQuery>
				<queryType>org.pdb.query.simple.SequenceQuery</queryType>
				<sequence>{0}</sequence>
				<eCutOff>{1}</eCutOff>
				<searchTool>blast</searchTool>
			</orgPdbQuery>
</queryRefinement>
<queryRefinement>
	<queryRefinementLevel>1</queryRefinementLevel>
		<orgPdbQuery>
			<queryType>org.pdb.query.simple.ResolutionQuery</queryType>
			<description>ResolutionQuery </description>
			<refine.ls_d_res_high.comparator>between</refine.ls_d_res_high.comparator>
			<refine.ls_d_res_high.max>{2}</refine.ls_d_res_high.max>
		</orgPdbQuery>
</queryRefinement>
<queryRefinement>
	<queryRefinementLevel>2</queryRefinementLevel>
		<orgPdbQuery>
			<queryType>org.pdb.query.simple.StructTitleQuery</queryType>
			<struct.title.comparator>!contains</struct.title.comparator>
			<struct.title.value>arp</struct.title.value>
		</orgPdbQuery>
</queryRefinement>
</orgPdbCompositeQuery>
"""
#
# Guess what is this sequence🙂 
#
#
sequence = """
MDDDIAALVVDNGSGMCKAGFAGDDAPRAVFPSIVGRPRHQGVMVGMGQKDSYVGDEAQS
KRGILTLKYPIEHGIVTNWDDMEKIWHHTFYNELRVAPEEHPVLLTEAPLNPKANREKMT
QIMFETFNTPAMYVAIQAVLSLYASGRTTGIVMDSGDGVTHTVPIYEGYALPHAILRLDL
AGRDLTDYLMKILTERGYSFTTTAEREIVRDIKEKLCYVALDFEQEMATAASSSSLEKSY
ELPDGQVITIGNERFRCPEALFQPSFLGMESCGIHETTFNSIMKCDVDIRKDLYANTVLS
GGTTMYPGIADRMQKEITALAPSTMKIKIIAPPERKYSVWIGGSILASLSTFQQMWISKQ
EYDESGPSIVHRKCF
"""
resolution="2.8"
#
# if you want to search remote homolog, increase E value cutoff. In this case, this protein is very highly conserved protein, so I used very stringent threshold.
#
eCutOff="1e-30"
biological_unit = True
reprot = """
http://www.rcsb.org/pdb/rest/customReport?pdbids={0}&customReportColumns=structureId,structureTitle,resolution,experimentalTechnique,depositionDate,structureAuthor,classification,structureMolecularWeight&service=wsdisplay&format=xml
"""

queryText = queryXML.format(sequence, eCutOff, resolution)
print "querying PDB...\n"

#
# Fetch pdb list using XML query in queryText
# and pdb list was in pdblist
#

req = urllib2.Request(url, data=queryText)
f = urllib2.urlopen(req)

result = f.read()

if result:	
	pdblist = []
	pdblist = result.split()


	if biological_unit:
		pdburl = "ftp://ftp.wwpdb.org/pub/pdb/data/biounit/coordinates/all/{0}.pdb1.gz"
	else:
		pdburl = "http://www.rcsb.org/pdb/files/{0}.pdb.gz "
	
	f=open('load.pml','w')
	
	for pdb in pdblist:
		print "download {0}".format(pdb)
		urllib.urlretrieve(pdburl.format(pdb.lower()), pdb+".pdb.gz")
		os.system("gunzip {0}.pdb.gz".format(pdb))
		f.write("load {0}/{1}.pdb\n".format(os.getcwd(),pdb))
		
	f.write("hide all\n".format(pdb))
	f.write("show cartoon\n".format(pdb))

	for i in range(1,len(pdblist)):
		f.write("cealign {0}, {1}\n".format(pdblist[0],pdblist[i]))
	f.close()
	
	req = urllib2.Request(reprot.format(",".join(pdblist)))
	f = urllib2.urlopen(req)
	result = f.read()

	if result:
		f=open('list.txt','w')	
		root = ET.fromstring(result)
		for record in root.iter('record'):
			recordline = []
			for field in record:
				recordline.append(field.text)
			f.write("\t".join(recordline)+"\n")
		f.close()


	else:
	    print "Failed to retrieve results" 
			
else:
    print "Failed to retrieve results" 

어쨌든 실행하면

querying PDB...
download 1ATN
download 1C0F
download 1C0G
download 1D4X
download 1DEJ
download 1EQY
download 1ESV
download 1HLU
download 1J6Z
download 1KXP
download 1LOT
download 1MA9
download 1MDU
download 1NLV
download 1NM1
download 1NMD
download 1NWK
download 1P8Z
download 1QZ5
download 1QZ6
download 1RDW
download 1S22
download 1SQK
download 1WUA
download 1YAG
download 1YVN
download 1YXQ
download 2A3Z
download 2A40
download 2A41
download 2A42
download 2A5X
download 2ASM
download 2ASO
download 2ASP
download 2BTF
download 2D1K
download 2FF3
download 2FF6
download 2FXU
download 2GWJ
download 2GWK
download 2HF3
download 2HF4
download 2HMP
download 2OAN
download 2PAV
download 2PBD
download 2Q0R
download 2Q0U
download 2Q1N
download 2Q31
download 2Q36
download 2Q97
download 2V51
download 2V52
download 2VYP
download 3A5L
download 3A5M
download 3A5N
download 3A5O
download 3CHW
download 3CI5
download 3CIP
download 3DAW
download 3EKS
download 3EKU
download 3EL2
download 3HBT
download 3M6G
download 3MMV
download 3MN5
download 3MN6
download 3MN7
download 3MN9
download 3SJH
download 3U4L
download 3U8X
download 3U9D
download 3U9Z
download 3UB5
download 3UE5
download 4EFH

BLAST에서 hit 된 pdb 파일을 쳐묵쳐묵 다운로드하고

이렇게 생성된 load.pml 파일을 더블클릭하거나 PyMol에서 @load.pml 식으로 불러오면

(파일이 수십개가 넘어가면 cealign 단계에서 시간 꽤 걸림. 인내심을 가지고 기다려보삼)

결과물. 헉 단백질이 디글디글. 너무 정신없으니까 좀 엇나간 넘들 좀 끄면..

대충 특정한 단백질 (Actin)에 붙는 바인딩 파트너들이 어떠한 방식으로 붙는지를 개괄적인 양상을 볼 수 있음.

list.txt를 엑셀 등으로 열어보면 다음과 같이 다운로드된 파일의 목록을 볼 수 있음.

위 스크립트를 응용하면 좀 더 다양한 조건을 지정하여 PDB의 목록을 뽑고, 이를 다운로드하고, PDB에 대한 관련정보들을 저장하여 추후분석에 사용할 수 있을 것임.

시간이 나면 CCP4의 ‘superpose’ http://www.ccp4.ac.uk/html/superpose.html 를 이용하여 PyMol 의 cealign 을 거치지 않고 바로 Structure superimpose를 하도록 스크립트를 고칠 예정. 시간나면. ;;;;

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s