정렬된 펩타이드 서열로부터 아미노산 교체 오즈비 (substitution odd-ratio) 구하는 기능

1. 개요

 우리는 이미 잘 계산된 blosum62과 같은 아미노산 간의 교체되는 것이 기대치에 비해 얼만큼 빈번히 일어나는지를 나타내는 오즈비 (odds-ratio) 매트릭스를 사용하고 있다. 그런데, 스스로 이것을 계산하고 싶을 때가 있다. 이것은 그 때를 위한 파이썬 스크립트이다.

 

2. 내용

 blosum62는 오즈비 (odds-ratio)를 구한 후에 2*log2(odds-ratio)로 변환된 값이기 때문에 아래에서 계산된 값에서 한 번 더 처리가 필요하다. 다만, 알아둬야할 것은 아래 기능에서는 negative inifinite가 나올 경우 계산할 방법이 없다는 것이다.

정렬된 펩타이드 서열 (aligned-read)가 너무 적은 경우에 교체가 일어나지 않는 아미노산들도 있기에 pseudo-count를 넣는 것도 생각해봤지만, 값이 심하게 왜곡되는 것 같은 인상을 받았기에 넣지 않았다.

따라서, 아래의 파이썬 기능을 통해 산출되는 0은 2가지 의미가 있다.

 

1) 오즈비를 산출하는 아미노산 어느 하나에서라도 출현한 횟수가 0인 경우 0 ÷ 

0 으로 오류가 발생해서 0을 집어넣은 경우2) 아미노산 1과 아미노산 2가 정렬된 펩타이드 서열에서 동시에 존재 (co-occurence) 한 경우가 없는 경우
(i.e P(A ∩ B) 

= 0)

from itertools import combinations,permutations
from collections import defaultdict

def calculate_substitution_odds_ratio(aligned_peptides):
    amino_acids = 'ARNDCQEGHILKMFPSTWYV'
    substitution_counts = defaultdict(int)
    # Count substitutions
    for peptide_pair in zip(*aligned_peptides):
        for aa1, aa2 in combinations(peptide_pair, 2):
            substitution_counts[(aa1, aa2)] += 1
            substitution_counts[(aa2, aa1)] += 1
    # Adjust count for identical pairs (ex. a->a, d->d) by dividing 2
    for aa in amino_acids:
        substitution_counts[(aa,aa)]/=2
    # Count aas
    amino_acid_count=dict(zip(list(amino_acids),[0]*20))
    for aa in amino_acids:
        keys=[sorted(key) for key in substitution_counts.keys() if aa in key]
        keys = [list(x) for x in set(tuple(x) for x in keys)]
        counts=sum([substitution_counts[(key[0],key[1])] for key in keys])
        amino_acid_count[aa]=counts
    # Calculate observed probability of occurence (qij of each amino acid pair i,j)
    total_pairs=sum(amino_acid_count.values())
    # Calculate observed probability of amino acid
    odds_ratios = {}
    for aa1 in amino_acids:
        for aa2 in amino_acids:
            observed = substitution_counts[(aa1, aa2)]
            expected = ((amino_acid_count[aa1])) * (amino_acid_count[aa2]) / total_pairs
            try:
                odds_ratios[(aa1, aa2)] = observed / expected
            except:
                odds_ratios[(aa1, aa2)] = 0 # the only case that evokes error is dividing 0 by 0
    return odds_ratios