[LNP] Step6_create_top_hexagon

 

 

 

고지가 얼마 남지 않았다....

힘들어도 킵고잉...!!!

 

 

이번 Step6_create_top_hexagon.ipynb에서는 Step5 노트북에서 코팅까지 끝낸 LNP 구조(coated_test.pdb)로 GROMACS topology(.top) 파일을 만드는 과정이다.

Step4 노트북과 이름이 똑같아서 헷갈릴 수 있지만, 그건 코팅 전 hexagon 용이라면 이번에는 코팅 후 버전의 토폴로지 파일을 만드는 것이다.

 

 

간단히 흐름을 요약하면 다음과 같다.

먼저, ITP를 불러오고 구조(coated_test.pdb)를 불러온다. 그 다음, residue를 순서대로 돌면서 어떤 분자인지 매칭한다. 

같은 분자가 연달아 나오는 구간을 묶어 (분자, 개수) 목록으로 만든다. GROMACS [ molecules ] 는 순서가 중요하기 때문이다.

 

다음으로 원자 수를 검증하고 분자별 개수 목록 (Hexagon_coat_mol_counting.top) 과 최종 topology (hexagon_coat.top) 를 출력한다.

 

 

---

# 라이브러리 불러오기

import MDAnalysis as mda
import numpy as np
import re
from vermouth.forcefield import ForceField
from vermouth.gmx.itp_read import read_itp
import networkx as nx

# Some useful lines for vermouth
#custom_ff.blocks['RNA'].make_edges_from_interaction_type('bonds')
#custom_ff.blocks['RNA'].make_edges_from_interaction_type('constraints')
#print(custom_ff.blocks['RNA'].nodes(data=True)[0])

 

아래 주석은 vermouth 디버깅에 쓸 수 있는 코드를 메모해 둔 것이다.

 

 

1 - import 부분

• MDAnalysis : PDB/GRO 같은 구조 파일을 읽고, 원자/residue 단위로 다루는 라이브러리
• numpy : 수치 배열 계산
• re : 정규표현식
• ForceField (vermouth.forcefield) : vermouth의 force field 객체, 분자 정의(block)를 담는 그릇 역할을 함.
• read_itp (vermouth.gmx.itp_read) : GROMACS .itp 파일의 줄 목록을 파싱해서 ForceField 안에 분자로 채워 넣는 함수
• networkx : 그래프 라이브러리로 vermouth의 Block이 networkx 그래프를 기반으로 만들어져 있다.

 

2 - vermouth 구조 배경

주석을 이해하려면 vermouth가 분자를 어떻게 저장하는지 알아야 한다.

• ForceField.blocks 는 {분자이름: Block} 형태의 dict 이다. ITP의 [ moleculetype ] 하나가 Block 하나가 된다.
• Block은 그래프다. node는 원자(bead) - resname, atomname, charge 같은 속성을 가진다., interactions는 ITP의 bonds, constraints, angles 같은 상호작용 목록이다. edge는 원자 사이의 연결로 interactions 와는 따로 관리된다.

 

3 - 주석 처리된 세 줄

• RNA Block의 bonds 와 constraints 를 보고 그래프에 edge를 만들어 준다.
• 결합이 bond로 정의된 곳도 있고 constraint로 정의된 곳도 있어서, 둘 다 넣어야 분자 연결이 끊김 없이 완성된다.
• 연결 그래프가 있으면 networkx로 연결 성분 개수를 세거나 이웃 원자를 찾는 분석을 할 수 있다.
• RNA Block에서 0번 node의 속성 dict를 출력한다.
• ITP가 제대로 읽혔는지 확인하는 용도이다.

 

 

---

# Generate new topology

ITP 파일들을 읽어서, 분자 이름마다 어떤 원자들로 이루어져 있는지 알려주는 사전을 만드는 과정이다.

 

- itp_to_forcefield() 함수 정의

def itp_to_forcefield(itp):
    """
    Returns a forcefield with all the molecule in the itp added.

    The itp should specify the path to the ITP file.
    """
    path = itp
    with open(path, 'r') as f:
        lines = f.readlines()
    forcefield = ForceField(name='custom')
    read_itp(lines, forcefield)
    return forcefield

 

이 셀의 역할은 ITP 파일 경로 하나를 받아서, 그 파일에 정의된 분자들이 들어 있는 vermouth ForceField 객체를 돌려준다.

 

 

 

---

- ITP 6개를 읽고 하나의 사전으로 합치기

# Import the itps which are sub dependencies
ff_rna = itp_to_forcefield('ITPs/RNA.itp')
ff_lipids = itp_to_forcefield('ITPs/martini_v3.0.0_phospholipids_v1.itp')
ff_lisbeth = itp_to_forcefield('ITPs/ILs.itp')
ff_sterols = itp_to_forcefield('ITPs/martini_v3.0_sterols_v1.0.itp')
ff_solvents = itp_to_forcefield('ITPs/martini_v3.0.0_solvents_v1_BMHB.itp')
ff_ions = itp_to_forcefield('ITPs/martini_v3.0.0_ions_v1.itp')
# Combine all the blocks in one big dictionary (a block is a molecule in our case)
all_blocks = {}
for forcefield in [ff_rna, ff_lipids, ff_lisbeth, ff_sterols, ff_solvents, ff_ions]:
    for name, block in forcefield.blocks.items():
        if name in all_blocks.keys():
            print(f'A double definition was found for {name}, this often results in issues!')
        all_blocks[name] = block
#all_blocks['RNA'].nodes(data=True)

 

먼저 ITP를 읽고 all_blocks로 합친다.

• 6개 force field에 있는 모든 Block을 {분자이름:Block} 형태의 dict 하나로 모은다.
• 같은 분자 이름이 두 파일에 다 정의되어 있으면 경로를 출력한다. 그러나 경고만 하고 멈추지는 않으며, 나중에 읽은 정의가 앞의 정의를 덮어쓴다. 그래서 리스트 순서(RNA, 지질, ILs, 스테롤, 용매, 이온)가 결과에 영향을 줄 수 있다.
• 이 dict의 key(분자 이름)는 나중에 residue.resname in all_blocks.keys()로 쓰인다. 즉, ITP의 moleculetype 이름과 PDB의 resname이 같다는 가정이 깔려있다. 이온의 경우에는 resname이 ION 이라서 맞지 않으므로 나중에 이온만 따로 처리하는 분기가 있다.

 

 

---

- molecule_identifiers 만들기

molecule_identifiers = {}
for molecule in all_blocks.keys():
    molecule_identifiers[molecule] = []
    for key, values in all_blocks[molecule].nodes(data=True):
        molecule_identifiers[molecule].append((values['resname'], values['atomname']))
#molecule_identifiers

 

이 셀의 역할은 분자마다 구성 원자(bead)의 (resname, atomname) 목록을 ITP 순서대로 만드는 것이다.

 

 

 

---

- 구조 파일 읽기

# Read in the GRO file
#!vmd ../hexagon/whole/whole.pdb
universe = mda.Universe('coated_test.pdb')

 

주석에는 GRO라고 써 있지만 실제로 읽는 건 PDB이다. # ! vmd ... 의 경우에는 VMD로 구조를 열어보는 셀 명령인데, 주석 처리되어 있고 경로도 맞지 않는 걸로 보아 사용하지 않는 걸로 보인다.

 

mda.Universe('coated_test.pdb') :

• 이전 노트북에서 만든 코팅된 LNP 구조를 읽는다.
• 파일이 하나뿐이라 topology 정보(원자 이름, resname, resid)와 좌표를 모두 이 PDB에서 가져온다.
• 이후 universe.residues, universe.atoms로 접근한다.

residue를 나누는 기준 :

MDAnalysis는 PDB를 읽을 때 resid나 resname이 바뀌는 지점에서 residue를 나눈다.

그래서 바로 옆에 붙은 두 분자의 resid와 resname이 우연히 같으면 residue 하나로 합쳐질 수 있다.

 

 

---

- residue 마다 분자 판별해서 matches 만들기

# Start matching the resnames
# This is not a very generic approach, but it can be expanded. The issue
#  is that there is no true solution to this problem. For certain molecules
#  might be subgroups of other molecules, meaning it is to a certain degree
#  an ambigous problem. Here we map small segments over long segments for 
#  simplicity, this might or might not be a problem for you system.
#  This can largely be circumvented by making sure that all resnames are
#  unique for every molecule type.
matches = []
active_molecule = []
in_RNA = False
for residue in universe.residues:
    # Matching by residue name if the molecule
    #  is a residue.
    if residue.resname in all_blocks.keys():
        in_RNA = False
        matches.append(residue.resname)
    # The naming of IONS is weird as they are separated
    #  molecules which have the same resname
    elif residue.resname == 'ION':
        in_RNA = False
        matches.append(residue.atoms.names[0])
    # The RNA or any polymer is an issue for
    #  their residues are not molecules. Here
    #  I wrote a basic detection for RNA. But
    #  this might need to be changed for other
    #  polymers. I only match by length, this
    #  might be lazy but for now I think it 
    #  will work.
    elif residue.resname in ['A', 'U']:
        if in_RNA == False:
            in_RNA = len(residue.atoms)
        else:
            in_RNA += len(residue.atoms)
        if in_RNA == len(molecule_identifiers['RNA']):
            in_RNA = False
            matches.append('RNA')
    else:
        in_RNA = False
        print(f'UNKNOWN RESIDUE {residue.resname}')
        break
#print(matches)

 

맨 위 주석의 의미는 다음과 같다.

어떤 분자가 다른 분자의 일부처럼 보일 수 있어서 residue 정보만으로 정하는 데에는 근본적으로 애매하다.

여기서는 단순한 규칙으로 매칭하는데 시스템에 따라 문제가 될 수도 있다.

분자 종류마다 resname을 겹치지 않게 지어 두면 이 문제는 대부분 피할 수 있다는 내용이다.

 

 

이후, residue 마다 4갈래로 나눠서 처리한다.

1. 일반 분자 ( resname in all_blocks ) :
   지질, 콜레스테롤, 물처럼 residue 하나가 곧 분자 하나인 경우에 해당한다.
   resname이 ITP의 moleculetype 이름과 같으면 그 이름을 그대로 추가한다.
   
   
2. 이온 ( resname.== 'ION' ) : 
   Martini 이온은 Na+든 Cl-든 resname이 모두 ION이다. 그래서 resname 만으로는 종류를 알 수 없다.
   대신 첫 원자의 이름(NA, CL 등)을 쓴다. 이 이름이 ion ITP의 moleculetype 이름과 같다.
3. RNA ( resname in ['A', 'U'] ) : 
   RNA는 뉴클레오타이드 하나하나가 residue라서, residue 하나가 분자 하나가 아니다.
   그래서 A/U residue의 원자 수를 계속 더하다가, 합계가 RNA.itp의 원자 수와 정확히 같아지는 순간, 'RNA' 하나로 기록하고 초기화한다. 따라서 주석에 적힌 대로 길이만으로 판정하는 단순한 방식을 사용한다.
   RNA 서열이 A와 U로만 이루어져 있어서 G나 C가 있으면 4번으로 빠진다.
   
   
4. 모르는 residue : 
   경고를 출력하고 루프를 아예 멈춘다. 그러면 matches에는 그 앞까지만 들어 있다.

 

 

---

- 연속된 같은 분자 묶기 (run-length encoding)

successive_elements = []
previous_match = None
counter = None
for match in matches:
    if match == previous_match:
        counter += 1
    else:
        if previous_match is not None:
            successive_elements.append((previous_match, counter))
        counter = 1
        previous_match = match
successive_elements.append((match, counter))

 

이 셀은 matches 에서 연속으로 반복되는 이름을 (이름, 개수)로 압축한다.

 

전체 개수를 합치지 않고 연속 구간으로 나누는 이유는 GROMACS의 [ molecules ] 섹션은 좌표 파일에 원자가 나오는 순서와 똑같아야 하기 때문이다.

 

동작 포인트는 다음과 같다.

• 루프 안에서는 이름이 바뀌는 순간에만 이전 묶음을 저장한다. 그래서 마지막 묶음은 루프가 끝난 뒤에 따로 한 번 더 append 해 줘야 한다.
• previous_match is not None 검사는 첫 번째 원소일 때 아직 빈 묶음 (None, None) 을 저장하지 않게 막는 역할이다.
• matches 가 비어 있으면 루프가 한 번도 돌지 않는다. 그러면 마지막 줄의 match 가 정의되지 않아서 NameError 가 난다.

 

 

---

- 원자 수 검증

# Quick sanity check to see if the amount of atoms
#  based on the segment counting is indeed the equal
#  to the total amount of atoms in the PDB/GRO.
amount_of_atoms = 0
for element in successive_elements:
    amount_of_atoms += len(molecule_identifiers[element[0]])*element[1]
all_fine = amount_of_atoms == len(universe.atoms)
print(f'The amount of atoms based on the segments is correct: {all_fine}')

 

세어 낸 분자 목록으로 계산한 원자 수가 PDB의 실제 원자 수와 같은지 확인한다.

 

successive_elements 의 각 원소는 (분자이름, 개수) 쌍이다.

합계를 len(universe.atoms) 와 비교한 결과를 all_fine(True/False)에 저장한다.

 

이 검사는 총 원자 수만 비교한다는 한계점이 있다. 분자 순서나 원자 이름이 맞는지는 보지 않는다. 따라서 4 bead 분자 하나가 빠지고 다른 곳에서 4 bead가 더 세어진 경우에는 통과할 수도 있다.

 

 

 

---

- 분자 개수 파일 쓰기

# Writing the mol_counting file which can be copied or imported in the top.
if all_fine:
    #with open('Hexagon_coat_solv_mol_counting.top', 'w') as f:
    with open('Hexagon_coat_mol_counting.top', 'w') as f:
        for successive_element in successive_elements:
            f.write(f'{successive_element[0]}\t{successive_element[1]}\n')
    print('The mol_counting.top as been written succesfully.')
else:
    'No ouput has been generated as the amount of atoms does not match.'

 

검증을 통과했을 대만 [ molecules ] 섹션에 들어갈 내용을 따로 파일로 쓴다.

GROMACS는 공백이든 탭이든 구분자로 받아서 탭을 써도 괜찮다.

 

주석 처리된 부분은 같은 노트북을 solvate 이후 구조에도 쓰려고 파일 이름만 바꿔 쓰던 흔적이 남은 것 같다.

 

 

 

---

- 최종 topology 파일 hexagon_coat.top 쓰기

topol = open('hexagon_coat.top', 'w')
            
topol.write(f'#include "ITPs/martini_v3.0.0.itp"\n')
topol.write(f'#include "ITPs/martini_v3.0_ffbonded.itp"\n')
topol.write(f'#include "ITPs/RNA.itp"\n')
topol.write(f'#include "ITPs/martini_v3.0.0_phospholipids_v1.itp"\n')
topol.write(f'#include "ITPs/ILs.itp"\n')
topol.write(f'#include "ITPs/martini_v3.0_sterols_v1.0.itp"\n')
topol.write(f'#include "ITPs/martini_v3.0.0_solvents_v1.itp"\n')
topol.write(f'#include "ITPs/martini_v3.0.0_ions_v1.itp"\n')
topol.write(f'[ system ]\n')
topol.write(f'LNP\n')
topol.write(f'\n')
topol.write(f'[ molecules ]\n')
topol.write(f'#include "Hexagon_coat_mol_counting.top"\n')
topol.close()

 

GROMACS grompp 에 넘길 메인 .top 파일을 만든다.

 

 

주의할 점은 다음과 같다.

• 여기서는 ..._solvents_v1.itp 를 include 한다. 두 파일의 moleculetype 이름이나 원자 수가 다르면 grompp에서 오류가 나거나 다른 정의가 적용된다.
• all_fine 과 상관없이 실행된다. 검증에 실패해도 top이 써지고, 위에서 말한 예전 개수 파일을 include 할 수 있다.
• include 경로는 상대경로다. .top 파일이 있는 폴더 기준으로 찾기 때문에 주의가 필요하다.

 

 

 

 

 

 

이제 .top 파일이 생성됐다.

노트북 6단계 끝~~!!!