[하루에 한 R] clustering 결과에서 Newick format의 트리 파일 추출하기

R에서 데이터프레임을 hclust() 함수로 처리하여 얻어진 클러스터링 결과물을 시각적으로 분석하는 일이 흔하다. 여기에서 만들어진 dendrogram을 Newick file로 추출하여 다른 응용프로그램에서 확인하려면 어떻게 하면 좋을까?

내가 구하는 대부분의 질문은 인터넷을 잘 검색만 해도 그 답이 나온다. 이번에는 StackOverflow에서 그 해답을 찾았다.

How to create a newick file form a cluster in R?

재료로 사용한 파일은 pyani가 생성한 ANIb_percentage_identity.tab이다. 이 파일은 헤더 라인의 첫 컬럼이 비어있으므로 row.names=1을 설정하면 아주 쉽게 데이터프레임으로 읽어들일 수 있다. R 코드는 아주 심플하다. as.phylo() 함수에 들어갈 값은 hclust class여야 한다.

> install.package("ape")
> library(ape)
> d = read.table("ANIb_percentage_identity.tab",sep="\t",header=T,row.names=1)
> fit = hclust(dist(as.matrix(d)),method="average")
> my_tree=as.phylo(fit)
> write.tree(phy=my_tree,file="test.newick")

참조했던 링크에서는 heatmap.2() 함수로 만든 자료를 연속적으로 변화시켜 가면서 newick 파일을 추출하는 방법도 소개되어 있는데 약간 번거롭다. Row와 column dendrogram 중 필요한 것을 골라서 변형시키면 된다. 이 방법 역시 ape 패키지를 로드해야 하며, heatmap.2() 함수를 쓰려면 gplots 패키지도 로드해야 한다.

> heat = heatmap.2(as.matrix(data))
> row.dendro = heat$rowDendrogram
> row.hcclust = as.hclust(row.dendro)
> row.phylo = as.phylo(row.hcclust)
> row.newick = write.tree(row.phylo)

마지막 명령은 좀 이상하다. 저런 방법으로 실행하면 R prompt에서 row.newick이라고 쳐서 화면에 newick 파일의 내용을 출력할 수는 있지만 이를 마우스로 긁어서 파일로 저장할 수는 없지 않은가. 다음과 같이 하는 것을 추천한다.

> write.tree(phy=row.phylo,file="test.newick")