protein bioinformatics: 1. 단백질 서열 정렬 (NWalign) 및 유사성

 

서열정렬 (sequence alignment)

 

bioinformatics 생명정보학, 생물정보학, 생정보학...

이 분야가 무엇이라고 한마디로 규정하기 어렵지만, 생물과 관련된 정보를 사용하는 분야이고, 대량의 데이터를 다루기에 컴퓨터를 사용하는 일이 많습니다.

대표적으로는 유전체 분석.. 유전체의 서열과 발현을 다루는 일이 있습니다.

사실 단백질도 유전체로부터 나오기에, 단백질 관련 연구도 중요한 주제입니다.

 저는 단백질 구조 및 기능 분야에서 일하고 있고, 당분간 바이오인포메틱스적인 관점에서 단백질을 어떻게 연구하는지, 그리고 어떤 스킬이 필요한지에 대해 이야기하고자 합니다.
 서열 정렬을 어떻게 해야하는지, pdb를 어떻게 파싱해야할지 같은 이야기입니다.

예제코드는 제 github에 있습니다. 

https://github.com/gicsaw/PBI

GitHub - gicsaw/PBI

단백질은 하나 이상의 아미노산 서열로 이루어져있습니다. 그리고 3차원 구조물이고, 고유한 접힘 (folding) 구조 (native structure)를 가지고 있습니다. 단백질이 접히는 단위를 도메인이라 부르는데, 하나의 도메인은 에너지적으로 안정합니다.(즉, 그 부분만큼만 서열을 잘라내도 구조가 틀어지지 않음) 또한 기능의 단위이기도 합니다. 예를들자면, EGFR, VGFR, TGFR 같은 kinase receptor 단백질의 경우, receptor 부분과 kinase 부분이 별개의 도메인으로 나눠져 있습니다. 또한 두개 이상의 도메인이 연관되어 hinge motion을 하기도 합니다. 뿐만 아니라 진화의 단위이기도 합니다. 하나의 도메인 내에서는 공진화 빈도가 높지만, 서로 다른 도메인 사이에선 공진화가 잘 나타나지 않습니다.

 

bioinformatics 적인 측면에서 상당수의 단백질 연구는 homology sequence (상동서열) 탐색에서부터 시작합니다.

 상동 서열 탐색에서 가장 유명한 프로그램으로 blast 가 있습니다. 

 

어떤 프로그램이든 완벽한 방법은 아니기에, 생물 정보를 제대로 활용하려면 프로그램을 돌리는 것만이 아니라 서열 정렬이라는 문제 자체에 대해 이해할 필요가 있습니다. 

일단 아미노산 서열을 기준으로 하겠습니다. 

 인간의 유전체에 코드에 대응하는 아미노산은 기본적으로 20가지가 있습니다. 하지만 아미노산이 20가지만 있는 것은 아니고, 변형된 형태로 존재하기도 합니다. 예를들자면, 두 시스테인의 S가 결합한 형태로 존재하기도 하고, 황 대신 셀레늄이 풍부한 환경에서 발현된 단백질은 시스테인에서 S가 Se로 치환되기도 하고, 단백질의 특정 아미노산이 인산화된 형태로 존재할 수도 있습니다. 인산화는 단백질 신호 전달 과정의 일부인데, 단백질을 인산화 하는 단백질을 kinase 라고 합니다. 
 일단은 20가지 아미노산에 대해서만 이야기하겠습니다. 아미노산을 표현할 때 3문자 를 사용하기도 하고, 1문자 를 사용하기도 합니다.  
이렇게 아미노산의 이름에서 앞의 3문자만 떼어온 것이 3문자 표기입니다. 

'ALA', 'ARG', 'ASN', 'ASP', 'CYS', 'GLN', 'GLU', 'GLY', 'HIS', 'ILE''LEU', 'LYS', 'MET', 'PHE', 'PRO', 'SER', 'THR', 'TRP', 'TYR', 'VAL'

영어 알파벳은 26문자이니, 이 아미노산들을 영문자에 대응시키고도 6개가 남습니다. 
"ARNDCQEGHILKMFPSTWYV"
단백질은 기본적으로 1차원 사슬 형태의 분자이기에, 1문자 표기를 사용하면 문장처럼 쉽게 표기할 수 있기 때문에, 서열 표기할 때는 1문자 표기를 많이 사용합니다. 분자 1차원 문장 표기 방법인 SMILES 같은 것보단 훨씬 쉽죠.

기본적인 서열 정렬 알고리즘부터 이야기해볼까 합니다. 

두 서열은 얼마나 유사할까? 라는 질문을 한다면 어떻게 대답할 수 있을까요?

 유사함의 기준이 무엇이냐? 라고 반문할 수 있습니다. 어떠한 기준 없이는 유사라는 말을 할 수 없습니다.
로봇과 사람이 유사한가, 사람과 고양이가 유사한가 라고 물었는데, 겉모습 이라면 로봇과 사람이 유사하겠지만, 
생물학적인 특성을 묻는다면 사람과 고양이가 유사하겠죠. 

 어떠한 자연 법칙에 의해서 기준이 정해지지 않는 한, 그 기준은 임의적일 수밖에 없습니다. 
 그럼 유사함을 판단하기 위해 어떤 기준을 제시할 수 있을까요?

 사실 거기에 대답하기 위해선 왜 유사한 서열을 찾으려 하는 것인가부터 생각해봐야 합니다만...

 그냥 간단한 예시로 단백질 구조가 유사할 것 같은 서열을 찾고 싶은 것이라 생각해보죠. 

 일단 간단하게 두 문장에서 아미노산이 서로 얼마나 일치하는가를 이야기할 수 있을 것입니다.

 하지만, 서로 길이가 다른 아미노산 서열도 있고, 길이가 같을지라도 같은 위치끼리 정렬했을 때 뭔가 이상하게 보이는 것들이 있을 수 있습니다. 그럴 때는 중간에 갭문자 '-' 을 삽입합니다. 유전적으로는 inserstion 이나 delation 에 해당됩니다. 즉, 두 서열을 단순히 같은 위치끼리 비교만 하는것이 아니라, 어느 문자와 어느 문자를 매치시킬지, 서로 정렬 하는 과정이 필요합니다. 또한, SNP 같은것이 있으니, 두 아미노산 서열의 특정 위치에서, 서로 동일한 문자만이 정렬되게 하는것이 아니라, 서로 다른 문자끼리도 정렬되게 하는 것이 좀 더 좋습니다. 물리적 특성이 유사한 아미노산끼리 정렬된 경우는 높은 점수를 주고, 물리적 특성이 다른 아미노산끼리 정렬된 경우는 패널티를 줄 수 있습니다. 
서열 정령에서 많이 사용되는 아미노산 사이의 치환에 대한 스코어로 BLOSUM62나 PAM250 가 있습니다. 여기에 갭 패널티도 추가로 사용합니다. 

스코어에서 중요한 것중 하나가 갭에 대한 패널티 입니다. 

예를들어 극단적인 경우로,

ABABABABAB 라는 서열과

AAAAA 라는 두 서열이 있습니다. 

이 둘을 

ABABABABAB

A-A-A-A-A- 라고 정렬한다면? 유사성은 매우 높게 나오겠죠 2번째 서열로 노말라이즈 한다면, 서열 유사성이 100%입니다. 서열 유사성은 서로 정렬되어 일치하는 문자의 수를 두 서열 중 하나의 길이로 나눠준 것입니다.  

하지만, 이런 점수는 전혀 의미가 없습니다. 서열의 유사성이 구조의 유사성과 무언가 관계가 있어야 하지만, 저런 정렬은 구조적으론 전혀 의미를 가지지 못합니다. 

이럴거면 차라리

ABABABABAB

AAAAA-----  

라고 정렬한 것이 더 나아보입니다. 아래처럼 정렬되기 위해선 갭이 처음 열릴 때 큰 패널티를 주되, 갭이 이어질 경우에는 패널티가 적게 더해지도록 하면 됩니다. 

 -  하나의 패널티가 -5 라면,

-- 의 패널티는 -5-1 = -6, 

--- 의 패널티는 -5-1-1 = -7

정도로 조절하는 것입니다. 

 

 스코어를 정했으면, 어떻게 정렬해야할까를 생각해볼 수 있습니다. 탐색 방법은 exact solution을 주는 방법과, heuristics 이 있습니다. exact 방법은 확실한 정답을 주지만, 이 방법으로 찾기에 계산 시간이 너무 오래걸리는 경우도 있습니다. heuristics 은 정답을 준다는 보장은 없지만, 빠른 시간안에 정답과 근접한 답을 찾기 위한 방법들입니다. 

 여기서 설명할 dynamic programing 은 exact solution 을 주는 방법입니다. blast 처럼 대량 탐색을 할 경우는빠른 탐색을 위해 heuristics method 를 사용합니다.

 dynamic programming은 하나의 큰 계산이 작은 여러 계산으로 쪼갤 수 있는 경우, 그리고 작은 계산들 중 중복된 계산들이 있는 경우에 사용합니다. 예를 들자면 1!+ 2!+ 3! + 4! + 5! + 6! 이라는 문제를 분다고 생각해봅시다. 
1 + 2*1 + 3*2*1 + 4*3*2*1 + 5*4*3*2*1+ 6*5*4*3*2*1 인데, 이 계산에는 중복된 부분들이 나타납니다.

이 계산을 전부 다 해줄 필요 없이, 이전항의 결과를 메모리에 기록해두고, 항의 번호에 이전항의 결과를 곱해주기만 하면 됩니다.
python으로 한다면,

import numpy as np
factorial = np.zeros(6, dtype=int)
for i in range(0,6):
 if i==0:

  factorial[i] = 1

 else:
  factorial[i] = factorial[i-1]*(i+1)
result = factorial.sum()

print(result)

로 만들면 되겠죠. 
 이렇게 중복되는 계산을 매번 하지 않고 메모리에 저장해두고 사용하여 계산량을 줄이는 것이 dynamic programming 의 장점입니다. 

서열 정렬 문제도, 모든 정렬을 다 시도해보고, 그중에 best 를 찾아야 하는 문제인데, 그렇게 하면 비효율적이니 계산량을 줄여야 합니다. 
서열 정렬 문제에선 보통 Needleman–Wunsch algorithm (NW) 을 많이 사용합니다. 
구체적으로similarity matrix 를 1개만 사용하는 대신 좀 느린 방법과 similarity matrix를 4개 사용하는 대신 더 빠른 방법이 있습니다. (affine?...) 위키피디아의 설명은 전자의 알고리듬이고, 제가 만든 예시 코드는 후자의 알고리듬을 사용했습니다.  
https://en.wikipedia.org/wiki/Needleman%E2%80%93Wunsch_algorithm

위키피디아에 있는 예시는 DNA 서열이지만, DNA나 단백질이나 별로 다르지 않습니다. 스코어링 함수가 조금 다른 정도..? 단백질 서열에 대한 예시는 책에 있는데, 그림 그리기가 너무 귀찮아서 그냥 이걸 퍼왔습니다. 
서열 정렬에 대해서 similarity matrix를 미리 만들고
q = 'GCATGCU'
d = 'GATTACA'

score matrix 를 F라 하면
F[0,j] = d*j, 
F[i,0] = d*i
F[i,j] = max(F[i-1,j-1] + S(q[i], d[j]), F[i,j-1]+d, F[i-1,j] +d)
여기서 d는 갭 패널티 (-1), S는 두 문자 q[i], d[j] 의 similarity score 입니다. 위 에서에선 일치하면 1, 일치하지 않으면 -1 입니다. 
이렇게 표를 채워나가면, 최대 정렬 스코어가 나옵니다. 표에서 제일 오른쪽 아래 값입니다.

 이 값은 3개의 값 중 최대값 입니다. 그 3개의 값은 위에서 왔는지, 대각선 방향에서 왔는지 옆에서 왔는지를 지칭합니다. 이 값을 역추적해나아가면, 최대의 스코어를 주는 정렬을 얻을 수 있습니다. 위에선 경로가 단일하지 않습니다. 즉, 최대 스코어를 주는 정렬이 유일하지 않습니다.


제가 공부할 때 사용한 책은 이미 제 손에 없고, 2013년쯤 공부한거라서 자세히 설명하기도 어렵네요. 
나중에 시간나면 책 구해다가 다시 공부해보면서 잘 적어봐야겠습니다. 
위키피디아에 잘 정리되어 있고, 제가 만든 예시 코드도 있습니다. 

https://github.com/gicsaw/PBI/blob/main/pbi/nwa.py
https://github.com/gicsaw/PBI/blob/main/bin/nwalign.py
다만, python은 좀 느려서 저는 이 코드는 잘 사용하지 않고, 
The Yang Zhang Lab 의 코드를 사용합니다. 
https://zhanglab.dcmb.med.umich.edu/NW-align/

그리고, 최근에 NW score를 단백질 서열의 similarity score 처럼 사용하는 사람들이 있던데, 좋지 않습니다.
이 값은 노말라이즈되지 않았기에, 단백질 서열 길이에도 의존합니다. 

 일반적으로 단백질 서열 유사성은 Seqence identity 를 사용합니다. 정렬된 서열에서 실제로 일치하는 서열의 비율을 봅니다.