[하루에 한 R] ANI 매트릭스로 덴드로그램 그리기

오늘 작성하는 글에서 다루는 내용은 기술적으로 정확하지 않을 수도 있음을 사전에 밝히는 바이다.

다양성 측면에서 미생물 유전체를 다루는 사람이라면 서열 정보를 가지고서 두 스트레인이 동일한 종에 속하는지 혹은 그렇지 않은지를 판단해야 하는 상황에 종종 직면하게 된다. 복잡한 이론은 다 잊어버리고, 다음과 같이 4줄로 요약해 보겠다.

• 16S rRNA gene의 similarity가 97% 미만이면, 두 균주는 서로 다른 종에 속한다. 역은 성립하지 않는다(무슨 말이냐 하면 97% 혹은 그 이상일 경우 같은 종일 수도 있고 아닐 수도 있다는 말이다).
• 97% 16S rRNA gene similarity는 70% DDH(DNA-DNA Hybridization)에 해당한다.
• 70% DDH에 해당하는 ANI(Average Nucleotide Identity) 값은 95~96%이다.
• 일반적인 경우라면 동일 종에 속하는 균주는 TETRA > 0.99를 만족시켜야 한다.

종 구분을 위한 도구 중에 specI라는 것도 있다. 이는 40개의 phylogenetic marker gene을 사용하는 도구이다. annotation이 필요하지만 genome 수준의 ANI 계산이 없어서 빠르다. Nature Methods(2013)에 실린 논문은 너무 함축적이라서 읽기가 힘들다. 웬만한 것은 supplementary material을 보라고 해 놓으니...(국립국어원: '웬만해서'가 옳은 표기입니다).

ANI 값을 계산하여 행렬 형태로 만들어 주는 대단히 편리한 소프트웨어로서 JSpecies라는 것이 있다. 이 프로그램을 돌리면 다음 그림과 같은 매트릭스가 나온다. 대각선은 자기 자신에 대한 계산이므로 '---'표시되는데, 이를 100(단위는 퍼센트)으로 치환하여 보았다. 컬럼 라벨과 행(row) 라벨이 동일하다. 이 자료를 가지고 덴드로그램을 그리는 것이 오늘의 학습 목표이다.

덴드로그램(dendrogram)이란 계층적 클러스터링 결과를 가지가 달린 나무와 같은 그림으로 표현하는 것이다. 

위 그림에서 보인 데이터가 data.txt하는 파일에 저장되어 있다고 가정하자. 각 값의 구분자는 탭이다. 아무런 생각 없이 그림을 그리려면 다음과 같이 하면 된다.

> data = read.table("data.txt", sep="\t", header=TRUE, row.names=1)

> plot(hclust(dist(data)))

설정을 전혀 변경하지 않았으므로, complete agglomeration 방법이 사용되었다. hclust() 함수 안에서 method="average"라고 지정하면 UPGMA(Unweighted Pair Group with Arithmetic Mean) 방법으로 클러스터링이 된다. 이상의 경우는 대칭 매트릭스라서 상당히 간편하게 클러스터링이 되었다(엄밀히 말하면 대각선을 중심으로 각 쌍의 값이 완벽하게 같지는 않다. blast 혹은 mummer 계산을 하면서 어느 것을 query로 하느냐에 따라서 결과 수치가 매우 근소하게 차이가 나기 때문이다). 여기서 조심할 것은 각 데이터는 row를 기준으로 클러스터링을 했다는 것이다. 만약 expression data와 같이 각 데이터가 컬럼 단위로 차곡차곡 쌓인 경우라는 t() 함수로 트랜스포즈를 해야 된다. 덴드로그램은 'column dendrogram'이지만 사용되는 데이터는 row 단위임을 명심하자. mtcars라는 내장 데이터셋을 가지고 시험해 보라.

제목과 설명이 마음에 들지 않으면 plot() 함수의 매개변수를 적절히 조절하면 된다.

> plot(hclust(dist(data), method="average"), main="Title", xlab="X-label", ylab="Y-label", sub="subtitle")

UPGMA에 의한 클러스터링을 그림으로 쉽게 설명한 사이트가 있어서 링크를 소개한다. 같은 사이트 내에 distance matrix를 이용한 클러스터링 서비스 페이지가 있다. 혹은 좀 더 심오해 보이는 웹 서비스를 이용하고 싶다면 DendroUPGMA를 방문해 보라.

자, 여기서 궁금증이 생겼다. 모든 궁금증이 오늘의 포스팅 내에서 해결되지는 않을 것이다.

1. 세로축 "Height"의 단위는 무엇인가? 공부해!
2. 이미 대칭인 매트릭스에 dist() 함수를 적용하는 것은 중복이 아닐까? 즉, row1 = (a, b, c...)와 같은 일반적인 형태의 "비대칭" 데이터가 더 적합한 것이 아닐까?
3. data 매트릭스를 이루는 수치는 percent identity이다. 만약 100 - data를 하면 이는 자연스럽게 dissimilarity matrix와 비슷한 것이 된다. 이것을 가지고 직접 클러스터링하면 안되나?
4. 세로축의 단위를 실제의 difference로 나타나게 할 수는 없을까?

3을 해결해 보자. 그저 R에서 100 - data라고 하여 새로운 오브젝트에 지정하면 된다. R은 이럴 때 참으로 편리하다.

논문을 찾아보면 100 - (ANI)를 한 매트릭스를 가지고서 덴드로그램을 그렸다는 글귀가 보인다. 그렇다고 해서 무턱대로 plot(hclust(data.2))를 하면 안된다. hclust() 함수가 받는 인자는 아무 매트릭스나 되는 것이 아니라 "a dissimilarity structure as produced by 'dist'"라고 명시되어 있기 때문이다. 만약 R에서 dist() 함수를 쓴다면 인자로 data를 주든 100-data를 주든 똑같은 트리가 나온다.  백문이 불여일견이다.. 플롯 2개를 동시에 그려보자.

> par(mfrow=c(1,2))

> plot(hclust(dist(data)))

> plot(hclust(dist(data.2)))

data.2 그 자체는 dissimilarity matrix(또는 distance matrix)가 아니다. dist()함수로 넘길 것이라면 data(=ANI)면 어떻고 data.2(100-ANI)면 어떤가?

그러면 구체적인 사례로서 논문을 직접 찾아보자.

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4154481/ 

100% - ANI를 ANI divergence로 정의하고 이를 hcluster Python package를 사용한 스크립트에 투입하였다고 한다. 이 논문의 Supplementary Data 2에 데이터와 파이썬 코드가 있다. 다음은 코드 전체이다. 결과 그림은 여기에 있다. 생각보다 그림이 예쁘다. 아마도 identity level에 따라 색깔이 입혀지도록 만든 것 같다.

# coding: utf-8

import pylab

import matplotlib.pyplot as plt

from pylab import loadtxt

from hcluster import *

from pylab import xlabel

X = loadtxt('SD2_ani_data.txt')

lab = open('SD2_ani_names.txt').read().splitlines()

Y = pdist(X)

Z = complete(Y)

dendrogram(Z,orientation='right',labels=lab,distance_sort='ascending',leaf_font_size=8,color_threshold=0)

xlabel('Euclidean Distance')

plt.savefig('ani_plot.svg')

입력물을 distance matrix로 전환한 다음 클러스터링하였다. 그렇다면 100 - ANI를 하지 않아도 상관이 없었을 것이다. 도대체 무엇이 정답일까? 구글신에게 뭐라고 여쭈어야 내가 원하는 대답이 나올까?

그렇다. Use own Distance Matrix in HCLUST? 완벽한 질문이다. 이 한 문장에 답이 들어있다. "Try to cast you distance matrix before using it in hclust()"

> data.3 = as.dist(data.2)

> par(mfrow=c(1,2))

> plot(hclust(dist(data.2)),main="Using dist() function")

> plot(hclust(data.3),main="Using as.dist() function")

100% 확신할 수는 없지만, 이제 height는 data.3 매트릭스에 실린 숫자와 동등한 것으로 보면 될 것이다. 마지막으로 원본 매트릭스와 dist()로 만들어진 distance matrix의 모양이 어떻게 다른지 확인해 보자. 아래 그림에서 data.2는 9x9 크기의 원본 매트릭스이고, data.3은 8x8로 줄어든 distance matrix이다. 원본에서 대각선에 위치한 자료는 distance = 0이므로 이를 제거한 조작이 가해진 것이다.

오늘도 예상외로 많은 지식을 쌓았다. 휴!