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)을 업로드해 두었다.

Heamap 데이터셋에서는 각 노드에 연결된 복수의 수치 자료(필드, row 형태로 나열)가 색상으로 표현된다. ANI 혹은 BSR 등이 heatmap으로 나타내기에 매우 적합한 자료 유형이다.

수치는 annotation file(template)에 정의된 color gradient(COLOR_MIN, COLOR_MAX)에 매핑된 결과에 의하여 색으로 나타난다. 각 필드에 해당하는 라벨은 annotation file 내에서 FIELD_LABEL로 정의되는데, FIELD_TREE에서 그 순서에 대한 정보를 별도로 제공할 수 있다. FIELD_TREE 라인이 없거나 주석화되어 있으면 FIELD_LABEL과 실제 데이터에서 각 필드의 순서를 주의 깊게 결정해야 한다.

앞 과정에서 만들었던 ANIm_modified_labels_sorted.csv를 heatmap으로 표현하는 방법을 알아보자. 컬럼/로 라벨은 전부 assembly accession(예: GCF_001573585.1)으로 바꾸어 놓은 상태이므로, heatmap에서 윗부분에 표현될 FIELD_LABEL은 실제 균주명 정보를 활용하여 “Clostridium botulinum CDC42961” 등의 형태로 표현되게 만들어야 한다. 그러기 위해서는 전 단계에서 만든 labels.txt 파일을 이용하여 컬럼 라벨을 되돌려 놓는 것이 핵심이다. 역시 R 환경에서 named vector를 이용한 데이터 조작 작업을 해야 한다. Labels.txt 파일은 실제 DATA 라인만을 추출하여 accession-label.txt 파일에 수록해 놓되, 구분자는 쉼표라고 가정하자.

# 이전 작업 내용을 지운다.
> rm(list=ls())

> ani = read.table("ANIm_modified_labels_sorted.csv",sep=",",header=T,row.names=1)
> data = read.table("accession-label.txt",sep=",",stringsAsFactors=F)
> key = data[,1]
> names(key) = data[,2]

# 만약 일부의 컬럼만 취하고 싶다면 해당되는 accession number를 수록한 벡터(key2)를 만들고
# 컬럼 라벨을 바꾸기 전에 다음을 실행한다.
> key2 = c("CBAW1310","GCF_000063585.1","GCF_000017025.1","GCF_001276985.1","GCF_001020205.1","GCF_001889325.1","GCF_000020285.1","GCF_000020165.1","GCF_000204565.1","GCF_000768655.1","GCF_000171095.1")

> ani.subset = ani[,which(colnames(ani) %in% key2)]

# key2 벡터에 들어갈 값이 너무 많아서 별도의 텍스트 파일(data.txt)에 line-by-line 
# 형태로 들어 있다면, 이를 읽어들여서 처리하면 된다.
# data.txt 파일이 multi-column file이라면 $V2, $V3..을 붙여넣어 지정한다.
> key2 = read.table(“data.txt”,sep=”\t”, stringsAsFactors)$V1
# 꼭 which(colnames() %in% ) 구문을 쓰지 않아도 된다. 다음과 같이 하면
# key2 벡터에 출현하는 순서대로 컬럼 이름이 정렬되는 효과를 갖는다.
> ani.subset = ani[,key2]

# key2의 순서에 맞게 컬럼을 정렬한다.
> ani.subset = ani.subset[,key2]

# column name substitution using named vector("key") and for loop
> for (i in names(ani)) {
  temp = names(key)[key==i]
  names(ani)[names(ani)==i] = temp
  } 

> for (i in names(ani.subset)) {
  temp = names(key)[key==i]
  names(ani.subset)[names(ani.subset)==i] = temp
  }

> write.table(ani,"ANIm_with_strain_names.csv",sep=",",quote=F)
> write.table(ani.subset,"ANIm_subset.csv",sep=",",quote=F)

https://itol.embl.de/help/dataset_heatmap_template.txthttps://itol.embl.de/help/dataset_heatmap_template.txtheatmap.txt 템플레이트 파일에서 SEPARATOR를 COMMA로 선언한 뒤 ANIm_with_strain_names.csv의 내용을 적절히 재배치하고 수정하면 된다. 첫줄은 “FIELD_LABELS,” 바로 뒤에 한 줄로 넣도록 하고, 나머지 라인 전체는 “DATA” 라인 바로 밑에 넣으면 된다. 아울러 FIELD_TREE 라인은 주석 처리한다.

이상의 방법에서는 ANI 매트릭스를 미리 정렬해 두고 iTOL에서도 그대로 표현되게 하느라 상당히 공을 들였다. 행 단위의 정렬은 iTOL에 최초로 업로드한 트리 자료를 따르게 되므로, 여기에서 설명한 정렬이란 열(컬럼) 단위, 즉 수치화된 필드의 정렬을 의미한다. 그런데 열 단위의 정렬 정보를 별도의 트리 파일로 보유하고 있다면 이를 heatmap.txt 데이터 파일에 넣는 것으로 작업이 단순해진다(FIELD_TREE 항목). 단, 콤마가 섞인 newick format의 자료가 데이터 파일에 포함되므로 field separator는 콤마가 아닌 것으로 바꾸어야 하고, leaf label에는 괄호 등이 들어가지 않게 유의해야 한다. 또한 newick file의 leaf label은 field label과 같아야 한다. 다시 말해서 GCA_ 또는 GCF_(assembly accession)로 시작하는 파일명을 지닌 유전체 염기서열을 이용하여 ANI 매트릭스를 만들었고 이것으로부터 hierarchical clustering을 하여 newick file을 만들게 되면 당연히 assembly accession이 leaf label이 될 것이다. 그러나 iTOL annotated tree에서 heatmap 데이터의 필드를 assembly accession으로 표현하고 싶지는 않을 것이다. 따라서 accession-label.txt 파일을 참조하여 assembly accession을 균주명 등의 정보로 치환하는 것이 바람직하다. 그러나 괄호나 콜론 등 균주명에서는 허용될지라도 newick 파일에서는 특수한 용도로 쓰이는 문자는 다른 것으로 치환해야 하니 세심한 작업이 필요하다.

LS-BSR에서 얻어진 마커 유전자에 대한 blast score ratio 매트릭스도 마찬가지 방법으로 정리하여 iTOL에서 heatmap 데이터셋으로 표현할 수 있다. LS-BSR의 설치 및 일반적인 사용 방법은 본 매뉴얼 내의 다른 곳에서 설명하였다.

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

> bsr = read.table("RUN_toxins_bsr_matrix.txt",sep="\t",row.names=1,header=T,check=F)
> bsr.t = t(bsr)
# 자료 구조 확인용
> dim(bsr.t)
> View(bsr.t)
> row.names(bsr.t) = gsub("^.*_GCF_","GCF_",row.names(bsr.t))
> write.table(bsr.t,file="bsr_matrix_modified.csv",sep=",",quote=F)

컬럼으로 전환된 마커 유전자의 순서는 markers.fasta(or .pep) 파일에 있던 것과 달라질 수 있다. 아마도 LS-BSR 프로그램 내부에서 FASTA 파일의 서열 ID를 추출하여 문자열 기준으로 다시 정렬을 하는 것으로 보인다. 원래 상태로의 재정렬이 필요하다면 R 환경에서 하는 것이 좋을 것이다. 바로 위의 “Heatmap 자료의 표현” 항목에서 컬럼의 일부만을 추출하는 방법을 참고하기 바란다.