단백질 서열로부터 BLOSUM (Blocks substitution matrix)을 만들기 위한 내용들
TTSR2024. 4. 1. 11:17
728x90
반응형
1. 개요
어떤 목적이 있어서 정렬된 펩타이드들로부터 blosum을 만들어야할 필요가 있기도 한다. 그런데, 그런 과정을 하나하나 정리한 글은 찾기가 어렵다. 이 블로그 글은 이러한 목적을 가진 사람들을 위해 작성된 것으로 원본 데이터의 출처와 어떻게 만들어가는 지를 정리했다. 다만, blosum을 만들 때 쓰인 서열이 정리된 파일은 없고 frequency테이블로부터 출발하므로 여기서는 그 테이블에서 만드는 것은 다루지 않는다. 언어는 파이썬3 (python3)로 작성됐다.
block의 각 column별로 match/mismatch의 수를 구한다. 예를들어, 블락 내에 첫 번 째 컬럼이 A가 9개이고 S가 1개라면 여기에는 9AA matches와 1AS mismatch가 존재한다.
이 계산 과정을 전체 블락에 대해 전체 컬럼에 대해 진행하여 하나의 테이블로 만든다. 예를들어 위의 A와 S 예제에서는 AA의 조합은 8+7+...+1 = 36가지의 가능한 AA pairs가 나오고 AS 혹은 SA는 9가지의 가능한 조합이 나오며 SS pair는 없다.
위의 과정에서 산출된 column별 pair값을 하나로 합친다. 그래서, 한 개의 block이 종류가 w개인 아미노산으로 이루어진 s의 길이인 block의 경우 ws(s-1)/2 의 기여도를 갖느다. 위의 A와 S 예제는 45개이다.
이 과정을 통해 frequency table은 20+19+...+1=210개의 amino acid pairs에 대한 값으로 이루어진다. 이 부분의 과정을 파이썬으로 짜면 아래와 같다.
def get_frequency_table_from_block(peptides):
'''
# parameters
peptides (list) : 정렬된 펩타이드 리스트로 모두 동일한 길이를 가져야 한다.
# returns
frequency (dictionary) : 아미노산 짝끼리 얼만큼의 빈도로 변환되는지 담아 놓은 딕셔너리
형태는 frequency[('A','A')]=int 이다.
'''
from itertools import combinations,permutations
from collections import defaultdict
amino_acids = 'ARNDCQEGHILKMFPSTWYV'
substitution_counts = defaultdict(int)
# 교체된 숫자 세기
for peptide_pair in zip(*peptides):
for aa1, aa2 in combinations(peptide_pair, 2):
substitution_counts[(aa1, aa2)] += 1
substitution_counts[(aa2, aa1)] += 1
# 동일한 아미노산 짝에 대한 값을 2로 나눠서 보정하기 (ex. a->a, d->d)
for aa in amino_acids:
substitution_counts[(aa,aa)]/=2
return substitution_counts
위의 과정을 거치면 아래와 같이 예시데이터로 나온다
peptides=['A']*9+['S']*1
get_frequency_table_from_block(peptides)
"""
defaultdict(int,
{('A','A'):36,
('A', 'S'): 9,
('S', 'A'): 9,
('S', 'S'): 0.0,
('V', 'V'): 0.0,
...생략...)
"""
# AS와 SA 그리고 SS는 논문의 예시와 동일한 숫자가 나오는 것을 확인할 수 있다.
다만, 후술하겠지만 이것은 클러스터 사이즈로 조정해줘야 한다. 해당 기능 (c파트)은 아래에 있다.
b. 로그오즈비를 계산하는 과정 (Computing a logarithm of odds (Lod) matrix)
이제 위의 과정을 통해 계산된 아미노산 i,j 짝들 (1 =< j =< i =< 20)의 빈도수 (frequency)를 f_ij라고 하자. 이것을 이제 관측된 발생 확률값 (the observed probability of occurrence for each i, j pair)로 변경해준다. 이는 아래의 수식1로 한다.
수식 1
위의 펩타이드 예시에서는 f_AA = 36이고 f_AS = 9이므로 q_AA = 36/45 = 0.8이 되며 q_AS = 9/45 = 0.2이다.
def get_qij_from_fij(fij):
# Get total number of pairs
total_number=sum(fij.values())
# Remove the redundant number
# 예를들어 SA와 AS와 같이 중복되는 짝들의 숫자가 합쳐진 것을 보정해야한다.
amino_acids = 'ARNDCQEGHILKMFPSTWYV'
pairs_without_0=[list(set(key)) for key in fij.keys() if fij[key]>0]
pairs_without_0=[','.join(key) for key in pairs_without_0]
for p in set(pairs_without_0):
n=pairs_without_0.count(p)
if n == 2:
p1=tuple(p.split(','))
total_number-=fij[p1]
# Dividing the numbers by total number
qij={key: fij[key]/total_number for key in fij.keys()}
return qij
peptides=['A']*9+['S']*1
fij=get_frequency_table_from_block(peptides)
get_qij_from_fij(fij=fij)
"""
{('A', 'S'): 0.2,
('S', 'A'): 0.2,
('S', 'S'): 0.0,
('A', 'A'): 0.8,
...중략...}
"""
이제 각각의 i,j 짝들이 발생할 기대 확률값을 계산해야한다. 여기서 관찰된 쌍의 빈도는 집단의 실제 빈도를 나타낸다고 가정합니다. 예를들어, 36개의 짝은 A를 짝들의 양쪽에 가지고 있고 9개는 A를 한쪽에만 가지고 있습니다. 그래서, A의 기대 확률값은 [36+(9/2)]/45 = 0.9가 되고 S는 (9/2)/45 = 0.1이 됩니다. 일반적으로 i번 째 아미노산이 i,j 짝에서 발생할 확률은 수식 2를 따릅니다.
수식 2
따라서, 각 i,j 짝마다 기대 발생 확률 (the expected probability of occurrence, e_ij)는 아래와 같습니다.
i=j일 때 p_i*p_j
i≠j일 때 p_i*p_j + p_j*p_i = 2p_i*p_j
이것은 아래의 기능으로 구현됩니다.
def get_pi_from_qij(qij):
# Get p_i
amino_acids = 'ARNDCQEGHILKMFPSTWYV'
p={}
for i in amino_acids:
p_i=0
for j in amino_acids:
try:
if i!=j:
p_i+=qij[(i,j)]/2
else:
p_i+=qij[(i,i)]
except:
p_i+=0
p[i]=p_i
return p
peptides=['A']*9+['S']*1
fij=get_frequency_table_from_block(peptides)
qij=get_qij_from_fij(fij=fij)
get_pi_from_qij(qij)
"""
{'A': 0.9,
'R': 0.0,
'N': 0.0,
...중략...
'S': 0.1,
...중략...
'V': 0.0}
"""
위의 예제에서는 AA의 기대값은 0.9*0.9=0.81이고 AS+SA는 2*(0.9*0.1) = 0.18입니다. 그리고 SS는 0.1*0.1=0.01입니다.
아래의 기능으로 이것은 구현됩니다. 부동소수점 문제로 아주 작은 오차가 발생합니다.
def get_eij_from_pi(pi):
# Get e_ij
amino_acids = 'ARNDCQEGHILKMFPSTWYV'
eij={}
for i in amino_acids:
for j in amino_acids:
if i==j:
eij[(i,j)]=pi[i]*pi[j]
else:
eij[(i,j)]=2*pi[i]*pi[j]
return eij
peptides=['A']*9+['S']*1
fij=get_frequency_table_from_block(peptides)
qij=get_qij_from_fij(fij=fij)
pi=get_pi_from_qij(qij=qij)
get_eij_from_pi(pi=pi)
"""
{('A', 'A'): 0.81,
('A', 'R'): 0.0,
...중략...
('A', 'S'): 0.18000000000000002,
('A', 'T'): 0.0,
...중략...
('S', 'S'): 0.01000000000000002
}
"""
이제 오즈비 (odds-ratio)는 q_ij/e_ij로 계산하면 됩니다. 이것은 log2로 변환해서 bit unit으로 계산됩니다.
즉, s_ij = log2(q_ij/e_ij)입니다. 만일 관측 빈도들이 기대치와 동일하면 s_ij=0이고 빈도수보다 작으면 s_ij < 0이며 많이 관측되면 s_ij > 0이 됩니다.
def get_sij_from_qij_and_eij(qij,eij):
import numpy as np
# Get s_ij
amino_acids = 'ARNDCQEGHILKMFPSTWYV'
sij={}
for i in amino_acids:
for j in amino_acids:
pair=(i,j)
try:
sij[pair]=np.log2(qij[pair]/eij[pair])
except:
sij[pair]=np.nan
return sij
peptides=['A']*9+['S']*1
fij=get_frequency_table_from_block(peptides)
qij=get_qij_from_fij(fij=fij)
pi=get_pi_from_qij(qij=qij)
eij=get_eij_from_pi(pi=pi)
get_sij_from_qij_and_eij(qij=qij,eij=eij)
{('A', 'A'): -0.017921907997262457,
('A', 'R'): nan,
...중략...
('A', 'P'): nan,
('A', 'S'): 0.15200309344504975,
('A', 'T'): nan,
...중략...
('S', 'S'):-inf
('V', 'V'): nan}
이값은 스케일링 팩터인 2를 곱한 후 가장 가까운 정수로 변환해서 half-bit score로 저장합니다. 개별적인 교체 행렬들에 대해 아미노산 짝 별로 평균적인 상호적인 정보량 (average mutual information per amino acid pair H also called relative entropy)와 기대되는 점수 E (the expected score E in bit units)를 계산합니다.
def get_scores(sij,qij,pi,lambda_=2):
'''
lambda_ : scaling factor
'''
amino_acids = 'ARNDCQEGHILKMFPSTWYV'
# multiply sij by scale factor
sij2={key: 2*sij[key] for key in sij.keys()}
# get nearest intenger
def roundup(v):
import numpy as np
sign=-1 if v<0 else 1
try:
if (float(abs(v))-int(abs(v)))>=0.5:
return sign*(abs(int(v))+1)
else:
return sign*(abs(int(v)))
except:
return np.nan
sij2={key: roundup(sij2[key]) for key in sij2.keys()}
# Get relative entropy H
H=0
for idx,i in enumerate(amino_acids):
for j in amino_acids[:idx]:
try:
H+=qij[(i,j)]*sij[(i,j)]
except:
continue # 에러가 발생하는 경우는 sij==np.nan임.
# Get expected score E
E=0
for idx,i in enumerate(amino_acids):
for j in amino_acids[:idx]:
try:
E+=pi[i]*pi[j]*sij[(i,j)]
except:
continue
return sij2,H,E
peptides=['A']*9+['S']*1
fij=get_frequency_table_from_block(peptides)
qij=get_qij_from_fij(fij=fij)
pi=get_pi_from_qij(qij=qij)
eij=get_eij_from_pi(pi=pi)
sij=get_sij_from_qij_and_eij(qij=qij,eij=eij)
get_scores(sij=sij,qij=qij,pi=pi,lambda_=2)
"""
({('A', 'A'): 0,
('A', 'R'): nan,
...중략...
('A', 'P'): nan,
('A', 'S'): 0,
('A', 'T'): nan,
('S', 'S'): -inf
('V', 'V'): nan},
0.03040061868900995
0.013680278410054479
"""
c. 블락들 내에서 구획들을 클러스터링하기 (Clustering segments within blocks)
사실은 위의 과정들보다 먼저 나와야하는게 중복된 서열들 (redundant sequences)를 처리해주는 것이다.
이것은 특정 서열이 너무 많이 나와서 생기는 문제를 없애주기 위한 것이다. 서열을 정렬한 후에 동일성 (ientity)가 일정 퍼센트 (예. 62%) 이상이면 하나로 처리하는 식이다. 예를들어, A구획과 B구획이 80% 이상 정렬된 서열에서 일치된다면 A와 B는 함께 클러스터링되고 이것들의 기여도는 평균으로 전체 계산에 기여된다. 만일 C가 A나 B 둘 중 어느 한쪽에 80% 이상 겹친다면 A, B, C 모두 평균으로 처리된다. 비록 C가 A와 B 양쪽에 80% 이상 정렬되는 것과 동일하지 않더라도 말이다. (If C is identical to either A or B at >= 80% of aligned positions, it is also clustered with them and the contributions of A, B, and C are averaged, even though C might not be identical to both A and B at >=80% of aligned positions.)
예를들어 위의 예제에서 A를 가진 9개의 서열 중 8개가 클러스터링된다면 (9A-1S column are clustered) 이 컬럼의 빈도 테이블에서의 기여도는 2A-1S가 되며 2AS의 짝을 기여한다.
def clustering_segments_within_blocks(sequences, threshold=0.62):
"""
정렬된 펩타이드를 pairwise identity로 클러스터링 하는 기능.
Args:
sequences: 펩타이드 리스트
threshold: 클러스터링을 위한 최소한의 일치도 (default: 0.62).
Returns:
클러스터 리스트.
개별 클러스터는 기준치 이상 (greater or equal than threshold)의 서열들의 묶음으로 구성되어 있음.
"""
import numpy as np
from sklearn import cluster
# pair-wise 거리 계산
def calculate_pairwise_identity(seq1, seq2):
"""
Calculates the pairwise identity between two sequences.
서열 1과 2는 동일한 길이라고 가정함.
"""
matches = 0
for i in range(len(seq1)):
if seq1[i] == seq2[i]:
matches += 1
return matches / max(len(seq1), len(seq2))
distances = np.zeros((len(sequences), len(sequences)))
for i in range(len(sequences)):
for j in range(i + 1, len(sequences)):
identity = calculate_pairwise_identity(sequences[i], sequences[j])
distances[i, j] = distances[j, i] = 1 - identity
clustering = cluster.AgglomerativeClustering(
n_clusters=None,
linkage="single", # 논문에 따르면 이것은 min distance임.
distance_threshold=1 - threshold, # Convert identity to distance
metric='precomputed'
)
clustering.fit(distances)
clusters = []
for cluster_id in np.unique(clustering.labels_):
cluster_members = [sequences[i] for i in np.where(clustering.labels_ == cluster_id)[0]]
clusters.append(cluster_members)
return clusters
peptides=['AVVV']*3+['AAVV']*1+['ASSV']+['ABBB']+['SDDD']
clustering_segments_within_blocks(sequences=peptides,threshold=0.62)
"""
[['AVVV', 'AVVV', 'AVVV', 'AAVV'],
['SDDD'],
['ABBB'],
['ASSV']]
# 총 4가지의 클러스터가 생성된다.
# 이것을 활용해서 가장 처음에 나온 frequency table 부분에 넣어줘야한다.
"""
이제 이 클러스터를 활용해서 표준화 상수 (normalization factor)를 넣어서 frequency table을 생성하자.
내가 이해하기로는 단순히 펩타이드의 클러스터의 크기가 4라면 4를 나눠주는 식으로 기여도를 만드는 것 같은데
이 부분은 명확하지 않을 수도 있다. 예를들어 9A-1S에서 8개의 A서열이 하나의 클러스터가 된다면 2A-1S로 처리된다고 했으므로 8A의 비중이 1이 되게 하는 것이므로 클러스터 사이즈인 8로 나눠주면 된다고 생각한다.
다른 것으로는 펩타이드 서열 A와 펩타이드 서열 B가 클러스터 사이즈 4와 3에서 나왔다면 12로 나눠주는 것인가 싶다.
def get_frequency_table_from_block_from_clsuter(cluster):
from itertools import combinations,permutations
from collections import defaultdict
amino_acids = 'ARNDCQEGHILKMFPSTWYV'
substitution_counts = defaultdict(int)
# get normalization factor for peptides
norm_factor=sum([[1/len(c)]*len(c) for c in cluster],[]) # 각 펩타이드별 클러스터의 크기를 구함.
peptides=sum(cluster,[]) # 클러스터의 펩타이드를 풀어줌.
# 만일 어떤 클러스터에 펩타이드가 4개가 들어가 있다면 해당 펩타이드에서 나온 기여도는 1/4로 나옴.
# 펩타이드 A와 펩타이드 B가 4개와 5개의 클러스터에서 유래했다면 1/4 1/5로 하여 1/20의 기여도를 중화시켜준다.
idx_comb=list(combinations(norm_factor,2))
for peptide_pair in zip(*peptides): # 이 부분이 달라짐.
for idx, aa in enumerate(combinations(peptide_pair, 2)):
substitution_counts[(aa[0], aa[1])] += 1*idx_comb[idx][0]*idx_comb[idx][1]
substitution_counts[(aa[1], aa[0])] += 1*idx_comb[idx][0]*idx_comb[idx][1]
# Adjust count for identical pairs (ex. a->a, d->d) by dividing 2
for aa in amino_acids:
substitution_counts[(aa,aa)]/=2
return substitution_counts
peptides=['AVVV']*3+['AAVV']*1+['ASSV']+['ABBB']+['SDDD'] # 여기선 펩타이드의 클러스터링을 위해 변경했다.
cluster=clustering_segments_within_blocks(sequences=peptides,threshold=0.62)
get_frequency_table_from_block_from_clsuter(cluster=cluster)
"""
defaultdict(int,
{('A', 'A'): 3.375,
('A', 'S'): 2.5,
('S', 'A'): 4.0,
('V', 'V'): 1.9375,
('V', 'A'): 0.375,
('V', 'D'): 5.5,
('V', 'B'): 5.5,
('V', 'S'): 3.5,
('A', 'D'): 0.5,
('A', 'B'): 0.5,
('D', 'B'): 6.0,
('D', 'S'): 4.0,
('B', 'S'): 4.0,
('D', 'V'): 2.0,
('B', 'V'): 2.0,
...중략...
('Y', 'Y'): 0.0})
"""
