[하루에 한 R] POCP 매트릭스를 이용하여 미생물 균주를 genus 별로 군집화하기

POCP(percentage of conserved proteins)에 대한 간략한 소개는 여기를 참조하라. 결론만 이야기한다면 POCP 50%를 기준으로 하여 세균이 같은 속(genus)에 속하는지를 판별하는 것을 제안한 것이다. 50%라는 기준은 평균값을 통해서 제안된 것이라 실제의 분포는 꽤 넓다.

POCP 매트릭스를 만들면 R에서 계층적 군집화(hierarchical clustering) 분석을 통해서 같은 genus에 속하는 균주를 구분하는데 도움이 될 것으로 기대를 하였다. heatmap을 그리면 균주들 사이의 대략적인 관계를 파악할 수 있지만, 실제로 군집이 어떻게 이루어지고 그 멤버는 누구인지를 결정해 주지는 않기 때문이다.

군집 분석에 대하여 알기 쉽게 설명한 글이 있어서 먼저 소개해 본다. 거리를 계산하는 방법과 군집 알고리즘의 차이에 대해서 특히 유념하여 읽도록 한다.

군집분석(Cluster Analysis) - Amazon AWS
[R 군집분석 (Cluster Analysis)] 군집분석의 개념 및 유형

그러면 34개의 유전체 서열을 이용하여 만든 POCP 매트릭스로 작업을 해 보자. 이 매트릭스는 d3이라는 데이터 프레임에 들어 있는데, 아직 완성되지 않은 일이라서 genus 정보만 남기고 균주 정보를 살짝 가리기 위해 약간의 트릭을 썼다. 수정된 데이터 프레임은 d4로 저장하였다. paste(vector1, vector2, sep=" ")을 실행하면 vector1의 모든 원소는 같은 위치에 해당하는 vector2의 원소와 공백을 경계로 하여 연결된 상태로 반환된다. POCP 매트릭스를 dist() 함수로 처리하여 euclidean distance를 구하여 이를 기반으로 계층적 군집화를 실시한다. hclust() 함수에 넘겨주는 method를 average로 지정하면 바로 우리에게 친숙한 UPGMA(Unweighted Pair Group Method with Arithmetic mean)이 된다.

출처: R Friend

실제 코드와 실행 결과는 다음과 같다. 요즘 몰두하고 있는 미생물은 Agathobaculumn, Butyricicoccus, 그리고 Eubacterium의 일부이다.

d4 = d3
k = row.names(d3)
km = gsub(" .*$", " sp. ", k)
n = 1:34
row.names(d4) = paste(km, n, sep="")
x = hclust(dist(as.matrix(d4)),method="average")
plot(x)

그런데 이렇게 그린 그림에서 세로축의 height는 우리에게 별다른 느낌을 주지 않는다. POCP 값(퍼센트)을 100에서 뺀 수치를 나만의 거리로 사용하면 안 되는 것일까? ANI도 마찬가지 성격의 지표가 된다. 이렇게 자체적으로 만든 매트릭스를 그대로 디스턴스로 삼으려면 as.distance() 함수를 쓰면 된다.

mydist = 100 - d4
x = hclust(as.dist(as.matrix(mydist)),method="average")
plot(x)
rect.hclust(x, border="red",h=50)

덴드로그램의 구조는 약간 달라진다. 그러나 여기에서는 (100 - POCP)가 height로 나타나므로, 일정 기준을 만족시키는 클러스터 주변에 사각형 경계를 씌울 수 있다. 동일 클러스터를 이루는 균주들은 POCP가 50을 넘는 것들이니 같은 genus에 속하는 것으로 볼 수 있다.

각 클러스터를 이루는 멤버들의 이름을 확인하고 싶다면? 트리 형태의 데이터를 데이터 그룹으로 잘라내려면 cutree() 함수를 쓰면 된다. cuttree가 아님에 유의하라. height = 50을 기준으로 잘라내는 코드와 결과를 같이 살펴보도록 하자. 클러스터의 수를 지정하여 잘라낼 수도 있다(k = 원하는 수).

> groups = cutree(x2,h=50)
> table(groups)
groups
 1  2  3  4 
31  1  1  1 
> names(groups[groups==1])
 [1] "Butyricicoccus sp. 1"     "Agathobaculum sp. 2"     
 [3] "Clostridia sp. 5"         "Butyricicoccus sp. 7"    
 [5] "Butyricicoccus sp. 8"     "Agathobaculum sp. 9"     
 [7] "Butyricicoccus sp. 10"    "Agathobaculum sp. 11" 

 ....(후략)

벡터의 원소를 가리키는 이름(names)이 groups 벡터를 다루는데 얼마나 요긴한지를 알 수 있다.