[LNP] Step7.1_create_top_hexagon

 

 

 

Step 7.1 노트북을 실행한 뒤, 다시 Step 7_Solvate 노트북으로 돌아간다.

 

이 노트북의 목적은 hexagon_coat_solv.gro 를 읽어서 "어떤 분자가 어떤 순서로 몇 개 있는지"를 GROMACS 용 .top 파일로 적는 것이다. 이를 위해 구조 파일 - MDAnalysis 을 읽는 과정과 분자 정의(.itp) - vermouth 를 읽는 과정이 필요하다.

 

 

# 라이브러리 불러오기

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

• MDAnalysis 는 분자 구조/궤적 분석 라이브러리이다. universe.residues 로 잔기를 파일 순서대로 하나씩 훑으면서, 이름이 MC3H, CHOL, WN 인지 등을 확인한다. 순서가 중요한 이유는 .top의 [ molecules ] 순서가 .gro의 순서와 정확히 같아야 하기 때문이다.
• ForceField 는 분자 정의들을 담는 그릇 역할의 클래스이고, read_itp 는 GROMACS .itp 파일을 해석하는 함수이다. Block은 분자 하나의 정의를 의미하는데, 이 정의로 분자별 비드 수를 알 수 있어서 RNA 인식과 검산이 가능해진다.

 

---

# Step 1: Generate topology for the solvated system

Step 1은  물을 채운 시스템의 분자 목록(.top)를 새로 만드는 과정이다.

노트북 이름은 토폴로지 생성이지만, 실제로 하는 일은 .gro 에 있는 분자를 순서대로 세어서 적는 것이다.

 

GROMACS는 .gro - 각 비드의 좌표, .top - 그 비드들이 어떤 분자인지(.itp 정의 + [ molecules ] 목록) 두 파일을 함께 읽는다.

두 파일은 비드 수와 순서가 정확히 일치해야 한다. 

 

insane이 WN를 추가했는데, Step 6 노트북의 .top 에는 WN이 없어서 grompp가 원자 수 불일치 에러를 낸다. 

그래서 물이 포함된 목록이 새로 필요하다.

 

 

# 1. .itp 파일 하나를 읽는 함수

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 파일 하나를 읽어 분자 정의로 바꾸는 함수들을 만든다.

• read_itp 는 값을 돌려주는 대신 넘겨받은 forcefield 를 직접 채운다.
• 결과에서 forcefield.blocks['CHOL'] 처럼 분자 이름으로 정의를 꺼낼 수 있다.

 

 

# 2. .itp 6개를 읽어 all_blocks로 합치기

# 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)

 

RNA, 인지질, 이온화 지질, 스테롤, 용매, 이온 정의를 all_blocks 하나로 합친다.

같은 이름이 두 번 나오면 나중 정의가 앞의 것을 덮어쓰므로 경고를 출력한다.

 

 

 

# 3. 분자별 비드 목록 만들기

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

 

분자마다 (잔기명, 비드명) 목록을 만든다. 이 목록의 길이가 곧 분자의 비드 수이다.

• block은 비드를 노드로 갖는 그래프라서, .nodes(data=True) 로 비드 번호(key)와 속성(values)을 함께 얻는다.
• 속성 중 잔기명과 비드명만 뽑아서 정리한다.
• 이후 셀에서 이 목록의 길이(비드 수) 만 사용한다. RNA 인식과 원자 수 검산에 쓰인다.

 

 

# 4. 구조 파일 불러오기

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

 

hexagon_coat_solv.gro 를 불러온다.

• 물을 채운 구조를 읽어 universe에 저장한다. 
• #!vmd 줄은 원작성자가 구조를 눈으로 볼 때 쓰던 명령이 주석으로 남은 것으로 무시해도 된다.

 

 

# 5. 잔기를 분자 이름으로 바꾸기

# 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)

 

잔기를 처음부터 하나씩 분자 이름으로 바꿔 matches에 저장한다.

맨 위의 주석의 요지는 "완전히 일반적인 방법은 아니니, 분자마다 잔기명이 고유하면 안전하다"를 의미한다.

 

변수는 다음과 같다.

• matches : 결과를 담을 리스트
• in_RNA : 지금까지 누적한 RNA 비드 수 (아니면 False)

 

universe.residues 를 파일 순서대로 돌면서 잔기마다 네 갈래 중 하나로 처리한다.

1. 잔기명이 곧 분자 이름인 경우.
   MC3H, DSPC, CHOL, W, WN 이 여기에 해당한다.
2. 이온인 경우
   이온은 잔기명이 모두 ION 이라 구분이 안된다. 그래서 잔기의 첫 번째 비드 이름(CL)을 분자 이름으로 쓴다.
3. RNA인 경우
   RNA는 잔기 42개가 이어진 한 분자이다. 잔기마다 기록하면 안 되므로 비드 수를 더해 가다가, 357개가 되는 순간 RNA 하나만 기록한다. 길이만 보는 단순한 방식이지만, 이 시스템에서는 RNA 종류가 하나라 충분하다.
4. 모르는 잔기인 경우
   경고를 출력하고 반복을 즉시 멈춘다. 그러면 matches 가 중간에서 끊겨 검산이 실패하고, 파일도 만들어지지 않는다. 결과는 잔기별 이름이 나열된 긴 리스트이다.

 

 

# 6. 연속 구간으로 압축하기

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))

 

같은 분자가 연속되면 .top 에 적을 (이름, 개수) 형태로 압축하는 코드이다.

 

끝에서 CHOL은 아직 저장되지 않는다. 구간은 다음 분자가 나올 때 저장 되므로 마지막 구간 뒤에는 다음 분자가 없기 때문에 저장되지 않는다. 그래서 반복문 밖의 마지막 줄이 필요하다.

 

if previous_match is not None 은 첫 분자일 때 (None, None) 이 저장되는 것을 막는다.

 

 

 

#7. 원자 수 검산

# 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 in correct: {all_fine}')

 

분자 1개의 (비드 수 * 개수) 를 모두 더해 .gro의 실제 원자 수와 비교한다.

all_fine은 다음 셀에서 파일을 쓸지 결정한다. False 이면 .top 목록 파일이 만들어지지 않는다.

 

이 검산은 총 비드 수만 본다.

분자 순서는 .gro 순서대로 읽기 때문에 맞춰지고, 비드 이름 불일치는 이후 grompp가 경고로 알려준다.

 

 

 

# 8. 분자 목록 파일 쓰기

# 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.'

 

검산을 통과했을 때만 분자 목록을 파일로 쓴다.

Hexagon_coat_solv_mol_counting.top 을 작성한다.

 

 

 

# 9. 메인 토폴로지 파일 쓰기

topol = open('hexagon_coat_solv.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_solv_mol_counting.top"\n')
topol.close()

 

.itp 둘을 #include 하고, 위 목록을 불러오는 메인 파일 작성 단계이다.

hexagon_coat_solv.top 을 만드는 셀로, topol.write 로 고정된 내용을 차례로 적고 topol.close() 로 저장한다.

• martini_v3.0.0.itp : 힘장 기본 규칙([defaults]), 비드 타입, 비드 간 비결합 상호 작용
• martini_v3.0_ffbonded.itp : 결합/각도 파라미터를 b_C1_C1 같은 이름으로 정의한 매크로
• RNA. itp : RNA 분자 정의
• martini_v3.0.0_phospholipids_v1.itp : DSPC 등 인지질 정의
• ILs.itp : MC3H 등 이온화 지질 정의
• martini_v3.0_sterols_v1.0.itp : CHOL 정의
• martini_v3.0.0_solvents_v1.itp : W, WN 등 용매 정의
• martini_v3.0.0_ions_v1.itp : NA, CL 등 이온 정의

 

순서가 중요하다. 비드 타입과 매크로를 분자 정의들이 가져다 쓰므로, 힘장 파일 두 개가 먼저 나와야 한다.

[ molecules ] 의 순서와 개수가 .gro 와 일치해야 한다.