Biopython으로 PDB 탐방하기

PDB!! PDB!! 

---

들어가기 전에

PDB는 protein data bank의 약어이기도 하고, 거기서 제공하는 파일 형식이기도 하다. 여기서는 그냥 PDB라고 하면 데이터뱅크, PDB '파일'이라고 하면 PDB 파일이다. 

그리고 이새기들 쿡북쓰기 귀찮았는지 모듈 불러와야 하는 거 빼먹더라... 니들도 일하기 싫었구나 

 

파일 읽기

쓰기도 있긴 한데 그건 생략. 읽는것도 하난가 두갠가 오류나서 안됐다. 이 섹션에서 읽을 파일은 

1) mmCIF
2) MMTF
3) PDB파일
4) PQR

이다. 

 

mmCIF

from Bio.PDB.MMCIFParser import MMCIFParser
parser = MMCIFParser()
structure = parser.get_structure("7f0l", "/home/koreanraichu/7f0l.cif")
print(structure)

from Bio.PDB.MMCIF2Dict import MMCIF2Dict
mmcif_dict = MMCIF2Dict("/home/koreanraichu/7f0l.cif")
print(mmcif_dict)
# FileNotFoundError: [Errno 2] No such file or directory: '1FAT.cif'
# 위 에러가 뜨는 경우, 해당하는 확장자의 파일이 내 PC에 있어야 한다.

둘 중 편한 방식으로 불러오면 된다. 참고로 해당 파일이 컴퓨터에 있어야 불러올 수 있다. 다운받는거는 밑에 나옴. 

 

MMTF

from Bio.PDB.mmtf import MMTFParser
structure = MMTFParser.get_structure_from_url("4CUP")
print(structure)
# Bio.MissingPythonDependencyError: Install mmtf to use Bio.PDB.mmtf (e.g. pip install mmtf-python)
# 뭔데 이거

from Bio.PDB.mmtf import MMTFParser
structure = MMTFParser.get_structure_from_url("4CUP")
print(structure)
# Bio.MissingPythonDependencyError: Install mmtf to use Bio.PDB.mmtf (e.g. pip install mmtf-python)

개발자 나와봐 나랑 면담 좀 하자. 저건 뭔 거지같은 오류냐. 

 

PDB파일

from Bio.PDB.PDBParser import PDBParser
parser = PDBParser(PERMISSIVE=1)
structure_id = "2b10"
filename = "/home/koreanraichu/2b10.pdb"
structure = parser.get_structure(structure_id, filename)
print(structure)

from Bio.PDB import parse_pdb_header
with open(filename, "r") as handle:
    header_dict = parse_pdb_header(handle)
    print(header_dict)

일단 둘 다 되긴 되는데 뭔 Warning이 이렇게 뜨냐... 이럴거면 PDB 가서 봤지... 

 

PQR

from Bio.PDB.PDBParser import PDBParser
pqr_parser = PDBParser(PERMISSIVE=1, is_pqr=True)
structure_id = "7lfv"
filename = "/home/koreanraichu/7lfv.pqr"
structure = parser.get_structure(structure_id, filename, is_pqr=True)
print(structure)
# NameError: name 'parser' is not defined

예제 코드는 이렇게 되어 있는데, 저렇게 쓰면 오류난다. 위에는 pqr_parser로 불러놓고 아래에서 parser쓰면 컴퓨터는 이게 뭔데요 하지 찰떡같이 알아듣질 않음. 에이 이건 알파고도 못알아먹어... 

from Bio.PDB.PDBParser import PDBParser
pqr_parser = PDBParser(PERMISSIVE=1, is_pqr=True)
structure_id = "7lfv"
filename = "/home/koreanraichu/7lfv.pqr"
structure = pqr_parser.get_structure(structure_id, filename)
print(structure)

이렇게 하면 경고는 많이 뜨지만 일단 된다. 

 

쓰기

그냥 이렇게 쓰는구나만 보고 가자. 

 

mmCIF

io = MMCIFIO()
io.set_structure(s)
io.save("out.cif")

io = MMCIFIO()
io.set_dict(d)
io.save("out.cif")

 

MMTF

from Bio.PDB.mmtf import MMTFIO
io = MMTFIO()
io.set_structure(s)
io.save("out.mmtf")

 

PDB

io = PDBIO()
io.set_structure(s)
io.save("out.pdb")

class GlySelect(Select):
    def accept_residue(self, residue):
        if residue.get_name()=="GLY":
            return True
        else:
            return False
io = PDBIO()
io.set_structure(s)
io.save("gly_only.pdb", GlySelect())

밑에껀 글라이신만 매긴건가? ㄷㄷ 

 

PQR

io = PDBIO(is_pqr=True)
io.set_structure(s)
io.save("out.pdb")

 

Data representation

PDB에서 제공하는 데이터는 카테고리가 있다. CSS의 클래스 차일드 이런것처럼. 상위 개념부터 순서대로 

 

Structure>Model>Chain>Residue>Atom

 

으로 간다. 보통 model의 경우 NMR로 규명한거는 겁나 많고, X-ray crystallography로 규명한거는 하나. 

 

Model

Structure는 꺼내고 자시고 할 것도 없으니 모델부터 봅시다. 

 

from Bio.PDB.PDBParser import PDBParser
parser = PDBParser(PERMISSIVE=1)
structure_id = "5zi1"
filename = "/home/koreanraichu/5zi1.pdb"
structure = parser.get_structure(structure_id, filename)
# 뭘 불러와야 하죠

print(structure[0])

<Model id=0>

모델은 ID가 int 형태라 저렇게 불러오면 된다. 하위 카테고리고? 그건 또 따로 불러야된다. 

 

Chain

모델의 하위 카테고리. 얘는 ID가 A, B, C, D 이런 식으로 string 형태이다. 

 

model=structure[0]
print(model["A"])

<Chain id=A>

여기서 짐작하신 분들도 계시겠지만, 리스트 전체 뽑을거면 iteration 들어가야 한다. 그것도 밑에서 서술할 예정. 

 

Residue

얘부터 이제 뭔가 겁나 많아질 예정이다. 일단 Residue 자체는 튜플이다. 낙장불입 구성 요소는 

-Hetero-field-W(물분자)/H(다른 HET 영역에 기재된 분자들-예: Heme)/공백(standard molecule)
-Sequence identifier: 해당 residue의 시퀀스 내 위치 정보를 담은 정수
-Insertion code: 잘은 모르겠지만 insertion mutant 관련된 코드인 듯 하다.

이렇게 세 개. 

 

model=structure[0] # Model
chain=model["A"] # Chain
print(chain[(" ", 100, " ")])

<Residue PHE het=  resseq=100 icode= >

보통 저렇게 치기도 하고 숏컷으로 chain[100] 이런 식으로 치기도 한다. 

 

print(residue.get_resname())
print(residue.is_disordered())
print(residue.get_segid())
print(residue.has_id(10))

PHE
0
    
False

Residue에서 뽑아낼 수 있는 속성은 이름/Disordered 여부(1이면 Disordered)/segment ID/얘 아이디 있냐? 정도. 참고로 이 밑에 가면 정신없어서 내가 주석으로 달아놨다. 

 

Atom

# 이제 atom으로 가보자...
atom=residue["CA"]
print(atom.get_name()) # 이름!
print(atom.get_fullname()) # 이름! (풀네임)
print(atom.get_id()) # ID
print(atom.get_coord()) # 좌투더표
print(atom.get_vector()) # 좌표를 벡터 객체로 받음
print(atom.get_bfactor()) # Isotropic B factor
print(atom.get_occupancy()) # Occupancy
print(atom.get_altloc()) # alternative location specifier
print(atom.get_sigatm()) # 원자 매개변수의 표준편차
print(atom.get_siguij()) # Isotropic B factor의 표준편차
print(atom.get_anisou()) # anisotropic B factor

CA
 CA 
CA
[ 20.149 -36.25   32.947]
<Vector 20.15, -36.25, 32.95>
46.05
1.0
 
None
None
None

Process finished with exit code 0

뭔가 많으니 주석을 참고할 것. 

 

# 아 일일이 빼기 귀찮으시죠? 그럼 한줄로 ㄱㄱ하세요.
atom = structure[0]["A"][100]["CA"]
print(atom.get_name())

이런 식으로 한줄로 빼버릴수도 있다. 복잡해서 글치. 

 

Disorder

Point mutation에 대한 정보라는데 조회 안되던데? 

 

Disordered Atoms

from Bio.PDB.PDBParser import PDBParser
parser = PDBParser(PERMISSIVE=1)
structure_id = "4lqm"
filename = "/home/koreanraichu/4lqm.pdb"
structure = parser.get_structure(structure_id, filename)
# EGFR L858R을 불러보자
model=structure[0]
chain=model['A']
print(chain[858])

<Residue ARG het=  resseq=858 icode= >

일단 시범조교 L858R을 불러보자. (아르기닌 한글자 약어가 R이다)

 

atom.disordered_select("A")
print(atom.get_altloc())

근데 이거 써도 오류남... 

 

Disordered Residue

residue = chain[10]
residue.disordered_select("CYS")

왜 Attribute error가 뜨는지 1000자 이내로 서술좀 해주실까. 

 

Iteration 활용하기

다운로드? 그거 맨 밑에요. 

 

전체 보기

from Bio.PDB.PDBParser import PDBParser
parser = PDBParser()
structure = parser.get_structure("2b10", "/home/koreanraichu/2b10.pdb")
model = structure[0]
chain = model["A"]
residue=chain[10]
atom=residue["CB"]
print(atom)
# 이것이 각개전투다 절망편(...)

이것이 기존 방법이다. 근데 이렇게 보면 어때요? 귀찮아요. 그리고 파일을 까 보지 않는 이상 100번째 아미노산이 뭔지 모른다. 

 

from Bio.PDB.PDBParser import PDBParser
p = PDBParser()
structure = p.get_structure("2m3g", "/home/koreanraichu/2m3g.pdb")
for model in structure:
    for chain in structure:
        for residue in structure:
            for atom in structure:
                print(atom)
# iteration 가즈아!!!

일단 가즈아!!! 

 

structure2 = p.get_structure("2b10", "/home/koreanraichu/2b10.pdb")
atoms = structure2.get_atoms()
for atom in atoms:
    print(atom)
# 단축 가즈아!!!

위에 코드도 길고 번거롭다면 이런 방법도 있다. (아예 구조에서 다이렉트로 아톰 뽑아서 가져오기)

 

atoms = chain.get_atoms()
for atom in atoms:
    print(atom)
# chain에서 가져올 수도 있는 듯 하다.

이런 식으로 상위 카테고리에서 get_ 걸어도 된다. 

 

residues = model.get_residues()
for residue in residues:
    print(residue)

chains=structure2.get_chains()
for chain in chains:
    print(chain)
# chain도 된다.

models=structure2.get_models()
for model in models:
    print(model)
# 모델은 이걸로도 된다

structure단에서 get_뒷부분 바꿔서 가져와도 된다. 

 

취사선택하기

아 사람이 큰 그림만 볼 수는 없지... 또 조건부로 볼 때가 있어요 또. 

 

from Bio.PDB.PDBParser import PDBParser
p = PDBParser()
structure = p.get_structure("2b10", "/home/koreanraichu/2b10.pdb")
# Cytochrome C peroxidase

일단 시범조교 앞으로. 

 

for model in structure.get_list():
    for chain in model.get_list():
        for residue in chain.get_list():
            if residue.has_id("CA"):
                ca = residue["CA"]
                if ca.get_bfactor() > 50.0:
                    print(ca.get_coord())
# List 내 알파탄소 중 B factor가 50 이상인 것의 '좌표'

Residue 중 알파 탄소(아미노산의 몸통)의 B factor가 50 이상인 것의 좌표를 출력한다. 

 

for residue in chain.get_list():
    residue_id = residue.get_id()
    hetfield = residue_id[0]
    if hetfield[0] == "H":
        print(residue_id)
# HET residue를 전체 출력하는 코드

HET residue를 출력한다. 응? 저 ID? 위에 튜플로 된 그거 말하는거다. 튜플도 수정만 안 될 뿐이지 일단 인덱싱은 된다. 순서대로 W or H or 공백/sequence identifier/insertion code니까 가장 왼쪽(0번 인덱스)에 H가 들어간 것을 찾으면 HET residue. W는 물분자고 공백이면 단백질을 구성하는 아미노산이다. 

 

for model in structure.get_list():
    for chain in model.get_list():
        for residue in chain.get_list():
            if residue.is_disordered():
                resseq = residue.get_id()[1]
                resname = residue.get_resname()
                model_id = model.get_id()
                chain_id = chain.get_id()
                print(model_id, chain_id, resname, resseq)
# Disordered atom을 포함하는 모든 residue들을 출력한다. 참고로 여기에는 없다.

Disordered Atom을 포함하는 모든 residue들을 출력하는 코드. 

 

for model in structure.get_list():
    for chain in model.get_list():
        for residue in chain.get_list():
            if residue.get_resname() == "LEU":
                print(residue.get_id())
# 응용: Residue 중 류신 뽑아줘(Leu)

Residue 중 류신의 ID를 뽑는다. (그니까 몇번 몇번이 류신이냐 이 얘기)

 

polypeptide 관련

이거 하려면 뭐 또 불러와야되는데 이놈들 그거 기재 안했음... 

 

from Bio.PDB.PDBParser import PDBParser
from Bio.PDB.Polypeptide import PPBuilder
ppb=PPBuilder()
# 예제에는 안나왔는데 이거 두 줄 추가해야된다...
p = PDBParser()
structure = p.get_structure("7kz1", "/home/koreanraichu/7kz1.pdb")
model_nr = 1
polypeptide_list = ppb.build_peptides(structure, model_nr)
for polypeptide in polypeptide_list:
    print(polypeptide.get_sequence())

KKWTPPRSPFNLVQETLFHDPWKLLIATIFLNRTSGKMAIPVLWKFLEKYPSAEVARTADWRDVSELLKPLGLYDLRAKTIVKFSDEYLTKQWKYPIELHGIGKYGNDSYRIFCVNEWKQVHPEDHKLNKYHDWLWENHEKLS

이거 하려면 모듈을 불러야 하는데 쿡북에 빠졌다. 이녀석들 보게. 

 

Analyze

이거 항목은 많은데 거리 각도 이면각만 보고 끝냅시다... 뭐 깔라는데 안깔려요 그게... 

아, 그리고 각도랑 이면각은 

from Bio.PDB.vectors import calc_angle
from Bio.PDB.vectors import calc_dihedral

이거 있어야된다. 이새기들 이것도 빼먹음. 

 

거리

from Bio.PDB.PDBParser import PDBParser
p = PDBParser()
structure = p.get_structure("5ngo", "/home/koreanraichu/5ngo.pdb")
# 아무튼 불러와보겠습니다.
residues=structure.get_residues()
for residue in residues:
    print(residue)
# 나도 이거 처음보는거라 목록부터 봐야됨...

일단 불러놓고 생각하는 타입. 

residue1=structure[0]["A"][273]["CA"]
residue2=structure[0]["A"][278]["CA"]
print(residue1-residue2)
# Residue간 거리

거리는 뭐 재고 할 것도 없이 간단하다. 

atom1=structure[0]["A"][423]["CA"]
atom2=structure[0]["A"][423]["CB"]
atom3=structure[0]["A"][423]["N"]
vector1=atom1.get_vector()
vector2=atom2.get_vector()
vector3=atom3.get_vector()
angle = calc_angle(vector1, vector2, vector3)
print(angle)
# 각(별)도
# 0.5872530070961592

atom1=structure[0]["A"][415]["CA"]
atom2=structure[0]["A"][423]["CA"]
atom3=structure[0]["A"][431]["CA"]
atom4=structure[0]["A"][439]["CA"]
vector1 = atom1.get_vector()
vector2 = atom2.get_vector()
vector3 = atom3.get_vector()
vector4 = atom4.get_vector()
torsion = calc_dihedral(vector1,vector2,vector3,vector4)
print(torsion)
# Torsion angle 계산 
# 1.28386400091194

각도와 torsion angle은 벡터로 받아서 계산해야 한다. 

 

Internal coordinates for standard residues

적절한 번역 아시는 분 제보 바랍니다. 

 

model=structure[0]
model.atom_to_internal_coordinates()
for r in model.get_residues():
    if r.internal_coord:
        print(
            r,
            r.internal_coord.get_angle("psi"),
            r.internal_coord.get_angle("phi"),
            r.internal_coord.get_angle("omega"),  # or "omg"
            r.internal_coord.get_angle("chi2"),
            r.internal_coord.get_angle("CB:CA:C"),
            (r.internal_coord.get_length("-1C:0N")  # i-1 to i peptide bond
              if r.internal_coord.rprev else "None")
        )

<Residue SER het=  resseq=451 icode= > None 85.88349606675206 167.36084855239682 None 116.39783860086264 None

 

PDB 접속하기

얘네 이것도 모듈 불러야되는데 안썼더라. 

 

from Bio.PDB.PDBList import PDBList # 이거 불러야된다
pdbl = PDBList()
pdbl.retrieve_pdb_file("6WO1")

이렇게 받으면

이렇게 저장된다. (기본 형식은 mmCIF)

 

from Bio.PDB.PDBList import PDBList
pdbl = PDBList()
pdbl.retrieve_pdb_file("6WO1",file_format="pdb")
# 확장자를 따로 입력하지 않으면 CIF파일로 다운로드 된다.
# file_format="확장자"로 입력하면 특정 파일 형식으로 받을 수 있다.
# pdir="경로"를 입력하면 다운로드 경로도 정할 수 있다.

이렇게 file_format으로 지정해주면 된다. 

 

뭐야 PDB파일 돌려줘요 

혹시 경로지정 되냐고? pdir="경로" 하면 된다.