Activity for Phylogenetic Tree

UPGMA는 종이나 유전자 서열 간의 관계를 계통수로 표현하는 방법이다. 이 활동에서는 해당 알고리즘을 Python 코드로 구현하는 방법을 이해해본다. 또한, BioPython 라이브러리의 UPGMA 및 Neighbor-Joining(NJ) 알고리즘을 이용하여 계통수를 생성하는 과정을 실습한다.

1. UPGMA

Unweighted Pair Group Method with Arithmetic Means (UPGMA)는 서열 간 distance를 이용해 phylogenetic tree를 구성하는 가장 기본적인 방법이다. 본 실습에서는 해당 알고리즘을 어떻게 구현할 수 있는지 생각해 보고 실습한다.

Algorithm

1. 서열 쌍 간의 거리(distance)를 계산한다.
2. 가장 가까운 서열 쌍을 찾는다.
3. 해당 두 서열 을 하나로 합친다.
4. 남은 서열이 2개 이상이면 2번 단계부터 반복한다.

Pseudo-code

UPGMA 알고리즘은 (1) 가장 가까운 pair를 찾는 과정과 (2) 이를 하나의 그룹으로 병합하는 과정으로 구성된다. 이 과정에서는 index를 관리하는 것이 다소 복잡하므로, 전체 코드를 한 번에 구현하기보다 문제를 작은 단위로 나누어 단계적으로 해결하는 것이 바람직하다.

아래에서는 이 과정을 단계적으로 수행해본다.

서열 쌍의 distance matrix를 계산하고, 그 값을 dmat에 입력한다.

UPGMA(dmat)
    각 서열을 지닌 group을 만든다. (index_group)
    WHILE index_group의 길이가 2 이상일 경우
        가장 가까운 두 group을 찾는다. (find_closest_pair)
        두 group을 합친다. (merge_pair)

find_closest_pair(dmat, index_group):
    모든 pair 중 가장 가까운 pair를 찾는다.

merge_pair(index_group, closest_pair):
    closest_pair를 합친다.  

Python Code

UPGMA 알고리즘 자체는 간단하지만, program 자체는 약간 복잡하다. 아래 code를 JupyterNotebook에서 실행해보고, code가 어떻게 작동하는지 분석해본다.

• 실습을 위해 JupyterLab에서 phylogenetic_tree 디렉토리 생성
• upgma.ipynb 노트북 생성

Main code (전체 실행 code)

아래는 UPGMA 전체 과정을 정리한 code이다. 다만, find_closest_pair와 같은 함수는 아직 구형되어 있지 않고, 이는 아래에서 구현한다.

from Bio.Phylo.TreeConstruction import DistanceTreeConstructor
from Bio.Phylo.TreeConstruction import DistanceCalculator
from Bio import AlignIO
import Bio.Phylo as Phylo
from copy import deepcopy
import numpy as np

def UPGMA(dmat):
    n = len(dmat)
    index_group = []
    for i in range(n):
        index_group.append([i])    
    while len(index_group) > 1:
        closest_pair = find_closest_pair(dmat, index_group)
        index_group = merge_pair(index_group, closest_pair)
        print(len(index_group), index_group)

이 코드에서는 가장 가까운 서열 (혹은 서열 group)을 찾는 과정은 find_closest_pair, 두 서열 group을 합치는 과정은 merge_pair 함수로 분리해서 진행한다.

가장 가까운 pair 찾기 단계 후 두 pair 합치기 단계 진행

두 함수는 아래와 같이 구현할 수 있다. 그런데, 서열 group의 거리 계산 과정은 두 서열 group에 있는 원소들 간의 모든 거리를 계산해야 한다. 이 과정은 calc_dist_group이라는 함수로 분리해서 진행한다.

def find_closest_pair(dmat, index_group):
    m = len(index_group)
    closest_pair = []
    closest_dist = 1.0
    for i in range(m):
        for j in range(i):
            d = calc_dist_group(dmat, index_group[i], index_group[j])
            if d < closest_dist:
                closest_pair = [i, j]
                closest_dist = d
    return closest_pair

def merge_pair(index_group, closest_pair):
    new_index_group = []
    m = len(index_group)
    for i in range(m):
        if i not in closest_pair:
            new_index_group.append(index_group[i])
    new_index_group.append([index_group[closest_pair[0]], index_group[closest_pair[1]]])
    return new_index_group 

거리 계산 함수

• calc_dist_group 함수는 서열 group간 거리 계산을 위한 함수로 아래와 같이 구성할 수 있다.

def flatten_list(nested_list):
    flat_list = []
    for e in nested_list:
        if isinstance(e, list):
            flat_list.extend(flatten_list(e))
        else:
            flat_list.append(e)
    return flat_list

def calc_dist_group(dmat, group1, group2):
    distances = []
    for e1 in flatten_list(group1):
        for e2 in flatten_list(group2):
            if e1 > e2: 
                distances.append(dmat[e1][e2])
            else:
                distances.append(dmat[e2][e1])
    return np.average(distances)  

UPGMA tree 예제

• 이제 UPGMA tree 구성을 위한 과정이 준비되었다.
• 예제 파일을 이용해 UPGMA tree를 만들어 보자.
• 먼저 아래 파일을 받고, 실습 directory에 복사한다.

• 예제 파일은 MSA 결과로 여러 종류의 amylase의 서열을 포함하고 있다. 이 서열들의 관계를 phylogenetic tree를 이용해 표현한다.

aln = AlignIO.read(open('amylase_aln.txt'), 'fasta')
calculator = DistanceCalculator('identity')
dm = calculator.get_distance(aln)
dmat = dm.matrix

UPGMA(dmat)

BioPython을 활용한 UPGMA tree 생성

UPGMA tree 만들기

• BioPython을 이용해 UPGMA tree를 생성하기 위해서는 먼저 tree를 생성하는 constructor를 만들어 주고, 이를 이용해 upgma 함수를 실행시킨다.

constructor = DistanceTreeConstructor()
upgma_tree = constructor.upgma(dm)

UPGMA tree 출력하기

• UPGMA tree의 구조는 아래와 같이 출력할 수 있다.

print(upgma_tree)

UPGMA dendrogram 생성

• 또한, draw 명령어를 이용해 tree를 그릴 수 있다.

Phylo.draw(upgma_tree, do_show=False)

2. Neighbor-Joining

• Neighbor-Joining algorithm의 구현을 위해서는 Q-matrix를 계산하는 중간 과정이 필요하다.
• 도전하고 싶으신 분들만 coding에 도전해보자.
• 여기서는 BioPython을 이용해 tree를 생성하고, 그 결과를UPGMA 방법의 결과와 비교한다.

실습 파일 생성

• ng.ipynb 노트북 파일을 생성한다.
• 아래 code를 이용하여 필요한 프로그램을 불러온다.

from Bio.Phylo.TreeConstruction import DistanceTreeConstructor
from Bio.Phylo.TreeConstruction import DistanceCalculator
from Bio import AlignIO
import Bio.Phylo as Phylo
from copy import deepcopy
import numpy as np

• 아래 code를 이용하여 distance matrix를 생성한다.

aln = AlignIO.read(open('amylase_aln.txt'), 'fasta')
calculator = DistanceCalculator('identity')
dm = calculator.get_distance(aln)

BioPython을 활용한 UPGMA tree 생성

• UPGMA와 같이 constructor를 이용해 tree를 생성한다. 다만, nj 함수를 이용한다.

constructor = DistanceTreeConstructor()
nj_tree = constructor.nj(dm)
Phylo.draw(nj_tree, do_show=False)

• NJ tree는 기본적으로 root가 없는 tree이다.
• 따라서 dendrogram에서 root에 3개의 branch가 있는 것처럼 표현되는데, 이는 root는 아니고, 가장 마지막에 합쳐진 node들이라고 할 수 있다.
• UPGMA와 또 다른 특징은 조상 node를 공유하는 두 branch의 길이가 서로 다르다는 점이다.

Rooted tree 만들기

• Unrooted tree를 rooted tree로 전환할 수도 있다.
• 아래는 root_at_midpoint 방법을 이용한 예시이다.

nj_tree_rooted = deepcopy(nj_tree)
nj_tree_rooted.root_at_midpoint()
Phylo.draw(nj_tree_rooted, do_show=False)