Genome Informatics Laboratory at KRIBB

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== R을 이용하여 ANI matrix 자료 다루기 ======
제목은 ANI matrix(행렬)라고 하였으나 실제로는 **dataframe 자료**를 다루는 기술에 관하여 작성한 글이다. 실습용 파일 {{ :ani_r_exercise.zip}}을 다운로드하여 압축을 해제한다. RefSeq genome 51개를 [[https://github.com/widdowquinn/pyani|pyani]](blast)로 처리하여 얻은 결과 중 일부 파일과 newick tree 및 accession 관련 파일이 포함되어 있다.

  $ unzip ANI_R_exercise.zip
    Archive:  ani_r_exercise.zip
      inflating: ANIb_percentage_identity.tab  
      inflating: ANIb_alignment_coverage.tab
      inflating: accession-strain.csv  
      inflating: ezTree_51_pfam.nwk
      inflating: accessions.txt

[[http://stat405.had.co.nz/r-style.html|해들리 위캠이 권하는 R 스타일 가이드]]에 의하면 할당문에는 '=' 말고 '<-'를 쓰라고 권장하는데 이는 나의 습관과는 배치된다. 취향의 문제니까 별로 중요하지는 않다. 대부분의 프로그래밍 언어에서 '='를 쓴다는 현실을 잊지 말도록 하자.
  Use <-, not =, for assignment.
===== 데이터 파일을 읽어들이고 기본 점검하기 =====
ANI(percentage identity) 계산 결과 파일을 읽어서 ani 데이터프레임을 만든다.
  > ani = read.table("ANIb_percentage_identity.tab", sep="\t", header=T, row.names=1)
매트릭스가 아니라 dRep의 Ndb.csv 파일과 같이 <id1>,<id2>,<value> 형태의 파일이라면 임포트한 뒤 reshape 라이브러리를 사용하여 매트릭스로 전환하도록 한다([[https://blog.genoglobe.com/2018/04/r-reshape-pairwise.html|참조글]]).

pyani 결과물에서 행과 열의 이름은 이미 같은 순서이다. 굳이 확인하고 싶다면 다음과 같이 한다.
  > x = colnames(ani)
  > y = rownames(ani)
  > identical(x, y)
  [1] TRUE
데이터프레임의 대칭 구조가 아니라면 원본과 트랜스포즈한 것을 더한 뒤 2로 나누는 방법을 택하여 평균값을 취하게 만든다. t() 함수를 거친 데이터프레임은 행렬이 되지만, 이를 다시 데이터프레임과 더하면 결과적으로는 데이터프레임이 된다.
  > isSymmetric(as.matrix(ani))
  [1] FALSE
  > ani = (ani + t(ani))/2
  > isSymmetric(as.matrix(ani))
  [1] TRUE
  # if statement를 이용하여 테스트와 대칭 전환을 한 번에 수행
  > if (isSymmetric(as.matrix(ani))) {
        print("Dataframe ani is symmmetric")
    } else {
        ani = (ani + t(ani))/2
        print("Testing symmetry after conversion...")
        isSymmetric(as.matrix(ani))
    }

===== 계층적 클러스터링, 시각화 및 트리 자료의 간단한 조작 =====
클러스터링 자체는 R의 기본 기능에 포함된 hclust() 함수로 구현된다. 그러나 트리 자료를 파일로 입출력하거나 정렬된 label을 벡터로 반환하는 등의 작업을 하려면 **ape** 패키지([[http://ape-package.ird.fr/|홈페이지]] [[https://cran.r-project.org/web/packages/ape/vignettes/DrawingPhylogenies.pdf|활용 설명서]])가 필요하다. 물론 뒤에서 설명할 **phylogram 패키지**([[https://ropensci.org/blog/2018/07/12/phylogram/|소개의 글]])도 있다! ape에서는 표준으로 여겨지는 phylo 객체를 사용하지만 phylogram에서는 dendrogram 객체를 사용한다.
  > library(ape)
  > fit = hclust(dist(as.matrix(ani)),method="average")
  > plot(fit,hang=-1,cex=.7,main="UPGMA Clustering")
  # plot(fit)과 plot(tree)는 트리를 그리는 방향이 다르다. 
  > tree = as.phylo(fit)
  > plot(tree,cex=0.7)
  > class(tree)
  [1] "phylo"
  > write.tree(phy=tree,file="ANIb_UPGMA.newick")
  > tree # 트리 형태에 대한 설명
    
  Phylogenetic tree with 51 tips and 50 internal nodes.
  
  Tip labels:
    GCF_000146875.3, GCF_001956155.1, GCF_000597985.1, GCF_000686825.1, ...
  
  Rooted; includes branch lengths.
  
  > tree$tip.label # Leaf label 표시
  # Newick format의 트리 파일을 읽어들이기
  > my_tree = ape::read.tree("ezTree_51_pfam.nwk")
  > plot(my_tree, cex=0.5)
  > is.rooted(my_tree)
  [1] FALSE
ezTree_51_pfam.nwk은 unrooted tree이다. ape 패키지의 root() 함수를 사용하면 사용자가 지정한 outgroup이나 노드를 기준으로 rerooting을 할 수 있다. 그러나 간단하게 midpoint rooting을 하려면 [[http://blog.phytools.org/|phytools]] 패키지의 midpoint.root()함수를 쓰는 것이 편리하다.
  > library(phytools)
  >  my_tree = midpoint.root(my_tree) 
  > is.rooted(my_tree)
  [1] TRUE
  > write.tree(phy=my_tree,file="ezTree_51_pfam_rooted.nwk")

클러스터링은 행과 열 중 어느 것을 대상으로 이루어진 것일까? ANI dataframe은 행과 열의 라벨이 동일해서 언뜻 쉽게 파악이 되지 않는다. 정답은 __row 단위로 계층적 클러스터링__이 되었다는 것이다. R의 데이터프레임에서 개별 자료는 row 단위로 위치한다는 것을 상기하면 쉽게 이해가 될 것이다. 다음을 실행하여 row label을 변경한 뒤 확인해 보자.
  > rownames(ani) = gsub("^GCF_", "ROW_", rownames(pid))
  > fit = hclust(dist(as.matrix(ani)),method="average")
  > plot(fit,hang=-1,cex=.8,main="UPGMA Clustering with row names changed")
클러스터링을 했다고 하여 원본 데이터프레임의 자료 순서가 바뀌지는 않는다.데이터프레임의 행과 열도 클러스터링 결과에 맞추어 재정렬하는 방법을 알아보자.
  > label_sorted = fit$labels[c(fit$order)]
  > ani_sorted = ani[, label_sorted]
  > ani_sorted = ani_sorted[label_sorted, ]
  > write.table(ani_sorted, "ANIb_percentage_identity_clustered.csv", sep=",", quote=F)

  
==== 균주 이름이 들어간 트리를 만들어 보자(데이터프레임 라벨 수정) ====
현재 각 샘플의 식별자(행과 열의 이름)은 "GCF_000146875.3" 형태의 RefSeq assembly accession이라서 실제로 어떤 균주인지 알기가 어렵다. 원본 자료에서 쓰였던 assembly accession과 균주명을 서로 연결한 텍스트 파일('accession_strain.csv')이 있다면 이를 R에서 읽어들인 후 row label과 바꾼 다음 트리를 만들면 된다. [[https://github.com/kblin/ncbi-genome-download|ncbi-genome-download]]를 '%%--%%dry-run' 파라미터와 함께 실행하면 이러한 파일을 쉽게 만들 수 있다.

  $ cat accessions.txt
  GCF_000146875.3
  GCF_001956155.1
  GCF_000597985.1
  ...
  $ ncbi-genome-download -A accessions.txt bacteria --dry-run | \
      awk -F"\t" '/^GCF/{printf "%s,%s %s\n", $1, $2, $3}' > \
      accession_strain.csv
accession_strain.csv 파일은 아쉽지만 수정을 거쳐야 한다. 왜냐하면 두번째 컬럼(organism name) 값의 끝에 strain이 포함된 것이 꽤 있기 때문이다. 이는 유전체 정보를 제출한 사람의 잘못이라고 할 수 있다. 멋진 스크립트를 짜서 이를 수정할 수도 있지만 몇 줄 되지 않는다면 적당히 편집기에서 고치는 것이 나을 것이다. 수정 완료된 accession_strain.csv 파일은 다음과 같아야 한다.
  ...
  GCF_000146875.3,Paenibacillus polymyxa E681
  GCF_002894905.1,Paenibacillus sp. F4
  GCF_000164985.3,Paenibacillus polymyxa SC2
  ...
  
이미 대칭 상태가 되도록 처리를 마친 ani 데이터프레임의 row 라벨을 균주명으로 치환한 ani_s를 생성한 뒤 클러스터링을 해 보자. 아래에 보인 코드의 for{} 블록 내부에 있는 row.names() 함수를 names()로 바꾸면 컬럼 라벨을 바꾸는 효과가 있다. Named vector를 이용하는 매우 유용한 예제이므로 잘 익혀두기 바란다. 바꾸기 전의 라벨이 key 벡터에 들어있고, 바꾸고 나서의 라벨은 names(key)로 지정함을 잊지 말자.

  > data = read.table("accession-strain.csv",sep=",",stringsAsFactors=F)
  > key = data[,1]         # GCF_000463565.1
  > names(key) = data[,2]  # Paenibacillus polymyxa WLY78
  > key
         Paenibacillus polymyxa WLY78      Paenibacillus polymyxa DSM 365 
                    "GCF_000463565.1"                   "GCF_000714835.1" 
  ...
  > ani_s = ani
  > for (i in row.names(ani_s)) {
       temp = names(key)[key==i]
       row.names(ani_s)[row.names(ani_s)==i] = temp
    }
  > fit_s = hclust(dist(as.matrix(ani_s)),method="average")
  > plot(fit_s,hang=-1,cex=.7,main="UPGMA Clustering with strain names")
글씨가 많아져서 약간 어수선하지만 균주 정보가 트리에 표시된 것을 볼 수 있다.

==== 균주들을 클러스터 단위로 나누기 ====
각 클러스터에 균주들이 어떻게 포함되는지를 알아보고 그 중에서 하나를 골라 이것으로만 구성되는 ANI 데이터 서브셋을 만들어 보자. 먼저 클러스터 결과를 화면에 띄워 놓아야 한다.
  > plot(fit,hang=-1,cex=.8,main="UPGMA Clustering")
  # 5개의 클러스터로 나누어 보자.
  > groups = cutree(fit,k=5) # heigt(h)를 기준으로 자를 수도 있다.
  # 빨간색 상자 5개가 나타난다. 실제로 분할을 해 보자.
  > table(groups)
  groups
   1  2  3  4  5 
  21 20  5  2  3
  # 1번 그룹에 속하는 멤버는?
  > names(groups[groups==1])
   [1] "GCF_000146875.3" "GCF_001956155.1" "GCF_000686825.1" "GCF_001709075.1"
   [5] "GCF_900116035.1" "GCF_001707685.1" "GCF_001709135.1" "GCF_000463565.1"
   [9] "GCF_000520775.1" "GCF_000714835.1" "GCF_001956235.1" "GCF_000507205.3"
  [13] "GCF_001272655.2" "GCF_001719045.1" "GCF_001593085.1" "GCF_001922145.1"
  [17] "GCF_002916985.1" "GCF_002894905.1" "GCF_001955935.1" "GCF_001956115.1"
  [21] "GCF_001956225.1"
1번 그룹 멤버를 벡터로 저장한 다음 ani 데이터프레임에서 이에 해당하는 부분을 뽑아낸 서브셋을 만들어 보자. 데이터가 갖는 조건이 아니라 행 또는 열의 라벨을 이용하여 추출하는 것임에 유의하자. 
  # 추출할 행/열 라벨을 우선 벡터 x에 지정한다.
  > x = names(groups[groups==1])
  # 추출할 균주명이 별도의 텍스트 파일('list.txt')에 있다면 readsLines() 함수로 읽어들인다.
  # list.txt 파일의 마지막 줄 다음에 줄바꿈이 없으면 불평을 할 것이다.
  # readLines()가 밑줄('_')을 자동으로 공백으로 바꾸지 않게 하려면?
  > x = readLines("list.txt")
  # subset() 함수로 컬럼을 먼저 뽑아내고...
  > ani_sub = subset(ani, select=x)
  > dim(ani_sub)
  [1] 51 21
  # 이어서 row를 뽑아낸다.
  > ani_sub = subset(ani_sub, rownames(ani_sub) %in% x)
  > dim(ani_sub)
  [1] 21 21
  > write.table(ani_sub, "ANIb_percentage_identity_subset.csv", sep=",", quote=F)
[[https://www.rdocumentation.org/packages/base/versions/3.6.2/topics/subset|subset()]] 함수는 컬럼에 대해서 더 친절한 것 같다. 데이터프레임에서 이름에 의해 컬럼을 선택하는 방법은 이것 외에도 몇 가지가 더 있다. 2015년에 썼던 글 [[https://blog.genoglobe.com/2015/04/r-row.html|[하루에 한 R] 특정 조건에 맞는 row만 출력하기]]를 돌아보며 잠시 복습을 하자.

===== Heatmap 그리기 =====
다음의 자료를 통해 영감을 얻었다!

[[https://www.datanovia.com/en/lessons/heatmap-in-r-static-and-interactive-visualization/|Heatmap in R: Statistics and interactive visualization]]

한동안 gplots 패키지의 heatmap.2() 함수를 즐겨 사용하였으나 위 문서를 통해 ComplexHeatmap 패키지를 접하고 완전히 이리로 전향하기로 하였다. 기존에 작성한 문서로는 [[data_visualization|Data visualization]]이 있다. 여기에서는 ANI 데이터프레임에 대한 heatmap을 만들면서 column은 자체적으로 클러스터링하되 row는 외부에서 임포트한 자료를 사용하는 방법을 설명한다. 아울러서 heatmap에 표시된 label을 자유롭게 바꾸는 방법도 알아보도록 한다.

==== phylogram으로 외부 트리 읽어들이기 ====
ape 패키지의 read.tree() 함수로 읽어들인 phylo 객체를 as.dendrogram() 함수로 전환해도 된다.
  > library(phylogram)
  > tree.2 =  read.dendrogram(file="ezTree_51_pfam_rooted.nwk")
  # dendrogram class에서 tip label 추출하기
  > labels(tree.2)
  > plot(tree.2)

==== Heatmap 그리기 ====
우선 ani 데이터프레임을 자체적으로 클러스터링하여 heatmap을 만들어 보자. gplots heatmap.2()는 레이아웃을 수정하지 않으면 보기에 좀 불편할 것이다. 자세한 사항은 이전에 쓴 글 [[data_visualization|Data visualization]]을 참고하라. 반면 ComplexHeatmap Heatmap()은 특별히 레이아웃을 건드리지 않아도 라벨이 길다고 하여 경계 바깥으로 잘리거나 하지 않아서 매우 좋다.
  > library(gplots)
  > heatmap.2(as.matrix(ani),col=greenred(10),trace="none",scale="none",dendrogram=c("row"))
  > library(ComplexHeatmap)
  > Heatmap(as.matrix(ani), name="ANI", row_names_gp = gpar(fontsize = 5), column_names_gp = gpar(fontsize = 5))
{{ :heatmap_20230215.png?400 |Heatmap()으로 그렸음.}}
이번에는 row와 column을 서로 다른 방식으로 클러스터링해 보자. row는 임포트한 트리(ezTree_51_pfam_rooted.nwk, dendrogram 객체여야 함)를 기준으로 정렬하되 column은 ani 데이터프레임의 값으로 자체적인 클러스터링을 적용한다. 실제로 그려보니 별로 예쁘질 않다.
  > Heatmap(as.matrix(ani), name="ANI", cluster_rows = tree.2, row_names_gp = gpar(fontsize = 5), column_names_gp = gpar(fontsize = 5))
row와 column 전부 tree.2를 기준으로 정렬해 보자.
  > Heatmap(as.matrix(ani), name="ANI", cluster_rows = tree.2, cluster_columns = tree.2, row_names_gp = gpar(fontsize = 5), column_names_gp = gpar(fontsize = 5))
==== Label을 자유롭게 표시하기 ====
위에서 보인 ani 데이터프레임의 라벨 치환 사례에서는 accession-strain.csv 파일로부터 data 데이터프레임을 만들고, 이로부터 key라는 named vector를 생성하여 사용했었다. 그러나 Heatmap() 함수에서는 각 축의 고유한 라벨은 그대로 두되 그림으로 나타낼 라벨만 원하는 형태로 바꿀 수 있다. 단, data 데이터프레임으로부터 아까의 것과는 반대 구조를 갖는 named vector(key.2)를 만들어야 한다. 편의상 row label만 실제 균주의 full name으로 바꾸어서 heatmap을 그려 보자.
  > key.2 = data[,2]        # Paenibacillus polymyxa WLY78
  > names(key.2) = data[,1] # GCF_000463565.1
  > key.2
                        GCF_000463565.1                       GCF_000714835.1 
         "Paenibacillus polymyxa WLY78"      "Paenibacillus polymyxa DSM 365" 
  ...
  > Heatmap(as.matrix(ani), name = "ANI", row_labels = key.2[rownames(ani)], row_names_gp = gpar(fontsize = 5), column_names_gp = gpar(fontsize = 5))
{{ :heatmap_20230215b.png?400 |Row label을 치환하여 그린 heatmap. 데이터프레임 자체의 row label은 바뀌지 않았다.}}
ComplexHeatmap의 활용법을 자세히 알아보려면 [[https://jokergoo.github.io/ComplexHeatmap-reference/book/index.html|Complete Reference]]를 참고하기 바란다.