iTOL(Interactive Tree Of Life)는 트리라는 기본 데이터 위에 다른 분석 자료를 얹어서 시각화할 수 있는 웹 도구이다. 사용자 등록을 하면 데이터의 관리 및 공유가 가능하며, 논문으로 투고용으로 그대로 쓸 수 있는 고해상도 이미지 파일을 제공한다. 사용법에 대한 상세한 설명은 Help 페이지나 비디오 튜토리얼을 숙독하기 바란다. 2020년 10월부터 무료 구독자에 대한 지원이 크게 줄어든 것은 매우 아쉬운 점으로 남는다.

iTOL에서 가장 핵심이 되는 자료는 트리이며, annotation을 위해 추가되는 데이터셋은 leaf name을 key로 하여 연결된다. 따라서 트리 파일에서는 긴 taxon name보다는 accession 번호 등을 사용하고, 추후에 label annotation(template)을 적용하는 것이 훨씬 바람직하다.

업로드 가능한 트리 자료는 텍스트로 이루어진 Newick, Nexus 또는 PhyloXML 포맷의 파일이다. Roary에서 얻어진 core gene의 다중 서열 정렬 결과물을 FastTree로 처리하여 얻은 트리 파일 또는 ANI 분석물을 가공하여 얻은 트리 파일을 전부 이용할 수 있다. 후자의 경우 선행 분석 프로그램으로부터 쌍 혹은 매트릭스 형태의 ANI 자료를 얻게 되는데, 이로부터 트리 자료를 만들어 내려면 R을 사용해야 한다. 이번 실습에서는 pyani를 사용하여 만든 51개 Paenibacillus strain의 ANI 계산 자료(tab delimited text file)을 사용하자. ani_r_exercise.zip을 다운로드하여 압축을 푼 뒤 ANIb_alignment_coverage.tab 파일을 활용하도록 한다.

> library(ape)
# check=F는 행과 열의 라벨이 원본 파일에 있는 그대로 쓰이게 만든다. 
# 그렇지 않으면 R이 라벨 내의 공백을 밑줄로 바꾸거나 숫자로 시작하는 라벨을
# 수정하는 일이 벌어진다.
ani = read.table("ANIb_percentage_identity.tab",sep="\t",header=T,row.names=1,check=F)

# gsub() 함수를 이용하여 행과 열 라벨에서 불필요한 텍스트를 제거한다.
# 실제 상황에 맞게 함수를 적용해야 한다.
> colnames(ani) = gsub("^.*GCF_", "GCF_", colnames(ani))
> rownames(ani) = gsub("^.*GCF_", "GCF_", rownames(ani))

# 계층적 클러스터링 방법 중 average 메소드는 흔히 UPGMA라 불리는 것이다.
# 결과를 newick 포맷으로 저장한다.
fit = hclust(dist(as.matrix(ani)),method="average")
plot(fit,hang=-1,cex=.8,main="UPGMA Clustering")
my_tree = as.phylo(fit)
write.tree(phy=my_tree,file="ANI_UPGMA.newick")

# average 메소드는 흔히 UPGMA라 불리는 것이다.
fit = hclust(dist(as.matrix(ani)),method="average")
plot(fit,hang=-1,cex=.8,main="UPGMA Clustering")
my_tree = as.phylo(fit)
write.tree(phy=my_tree,file="ANI_UPGMA.newick")

# 데이터를 tree에 나타난 라벨 순서대로, 즉 클러스터링 결과에 맞추어 정렬하기
label_sorted = fit$labels[c(fit$order)]
ani.sorted = ani[,label_sorted]
ani.sorted = ani.sorted[label_sorted,]

# 데이터의 저장 I: 라벨만 수정한 데이터프레임
write.table(ani,"ANIm_modified_labels.txt",sep="\t",quote=F)

# 데이터의 저장 II: 클러스터링된 순서대로 정렬
write.table(ani.sorted,"ANIm_modified_labels_sorted.csv",sep=",",quote=F)

# 클러스터를 나누기 위해 다시 plot() 함수를 써서 시각화
plot(fit,hang=-1,cex=.8,main="UPGMA Clustering")

# 5개의 클러스터로 나눈다.
rect.hclust(fit,k=5)

# 각 클러스터의 멤버를 추출한다.
groups = cutree(fit,k=5) # height(h)를 기준으로 트리를 자를 수도 있다.
table(groups)
names(groups[groups==1])

이상의 과정에서 얻어진 ANI_UPGMA.newick을 iTOL에 업로드한 뒤 다음과 같은 형태로 준비된 labels.txt 파일을 마우스로 드래그하여 넣으면 제대로 라벨이 표시된다.

LABELS

SEPARATOR COMMA

DATA

GCF_001956155.1,Paenibacillus peoriae FSL J3-0120
GCF_000146875.3,Paenibacillus polymyxa E681
GCF_000217775.1,Paenibacillus polymyxa ATCC 842
GCF_000520755.1,Paenibacillus sp. 1-18
GCF_000520795.1,Paenibacillus polymyxa 1-43
...

labels.txt 파일을 만드는 방법은 비교적 쉽다. 실습용 zip 파일의 압축을 해제하였을 때 나오는 accessions.txt 파일을 ncbi-genome-download('-n' or '--$dry-run' 옵션)에 적용한 뒤, 출력물을 파일로 리다이렉트한 다음 적절히 편집하여 DATA 라인 아래에 채우면 된다. 구분자인 콤마 뒤에 공백을 넣지 않도록 주의한다.

$ ncbi-genome-download -A accessions.txt -n bacteria

실제로 이 명령을 실행해 보면 accessions.txt에 수록된 51개 genome 중 50개에 대한 정보만이 반환될 것이다(2022-08-11 기준). 이는 GCF_001663585.1(Paenibacillus polymyxa ND24)이 'unverified source organism'이라는 이유로 RefSeq에서 suppress되었기 때문이다. 이를 적절히 수정한 labels.txt 파일(zipped)을 업로드해 두었다.

앞에서 설명한 ANI 매트릭스는 column/raw의 이름이 iTOL 트리 자료의 leaf name 그대로이고 대각선에 대하여 대칭 구조이므로 원하는 자료만을 선별하거나 이름을 바꾸는 것 정도의 조작만 가하면 된다. 그러나 BSR 매트릭스는 각 유전체에 대한 BSR 값이 컬럼 단위로 위치하므로 이를 트랜스포즈해야 iTOL의 데이터셋 파일로 쓸 수 있게 된다. R에서 이를 실제로 어떻게 처리하는지 다음의 코드를 통해 알아보도록 한다.