Table of Contents
Data visualization
데이터 시각화의 목적은 화려한 그림으로 독자나 청중을 현혹하거나 기술을 과시하고자 함이 아니라, 어떤 메시지를 효과적으로 전달하기 위함이다. R을 활용함으로서 데이터의 시각화가 가능하다.
Matrix data의 시각화(+ tree 또는 dendrogram 곁들이기)
X1 = (V1, V2, V3, …Vn) 형태의 자료가 여럿(X1, X2, X3… Xn) 존재한다고 가정하자. 각 X들의 관계는 어떻게 될까? 어떤 유의미한 패턴을 보이지는 않을까? 이 데이터를 한데 모으면 결국 엑셀에서 열어볼 수 있는 행렬 형태의 데이터가 된다. Xi의 관계를 표현하는 전통적인 방법은 tree 형태의 그림이다. 혹은 각 Xi 자료의 내부를 구성하는 변수들의 패턴을 보려면 heatmap(+dendrogram)이 매우 유용하다. 여기에서는 LS-BSR(논문 PDF 매뉴얼)에서 소개한 BSR matrix의 시각화 방법을 알아본다.
유용한 링크들
- [정해영의 블로그] 매트릭스 형태의 텍스트 파일 입력과 헤더의 처리 (중요함!) - 첫 줄의 컬럼 수가 나머지보다 하나 적으면 header=T,row.names=1 파라미터가 필요없다는 것을 기억해 두면 편리하다.
- [정해영의 블로그] BLAST score ratio(BSR)을 이용한 유전체 비교 분석 - BSR의 개념 소개
- EvolView: an online phylogenetic tree viewer and customization tool
매트릭스 형태의 텍스트 파일을 R에서 읽어들이기(1)
Marker gene을 제공하지 않고 LS-BSR을 실행했을 때 만들어지는 결과 파일 중 하나인 ls-bsr_out_bsr_matrix.txt의 구조는 다음과 같다. 이 그림은 ooffice –calc <file>로 연 것을 캡쳐한 화면이다. 각 컬럼은 유전체, 행은 48개의 유전체에서 예측한 유전자(단백질로 번역)를 클러스터링한 뒤 추출된 대표 유전자들의 ID이다. 그러니 row의 수가 2만개가 넘는다. 다시말해서 48개 미생물에 대해서 계산을 통해 pangenome을 얻었고, 이를 구성하는 모든 유전자가 row에 자리잡았다는 뜻이다.
R에서 다음과 같이 입력하면 헤더를 정확히 인식하고 첫 컬럼을 row name으로 사용하여 BSR matrix를 데이터프레임으로 읽어들일 수 있다.
> data = read.table("ls-bsr_out_bsr_matrix.txt",sep="\t",header=T, row.names=1)
> View(data) # 그림으로 자료 구조 확인
> dim(data)
[1] 20177 48
이 데이터프레임은 실습용으로 그대로 쓰기에 곤란하다. 왜냐하면 row의 수가 너무 많기 때문이다(이 경우에는 20,177). 따라서 맨 위에 있는 것 100개(100줄)만을 골라내도록 한다. 다음으로는 transpose를 해서 각 유전체 자료가 row 단위로 배열되게 만들어야 한다. 각각의 자료를 column by column으로 배열하는 것은 생물학쪽에서 흔히 쓰는 관습이나(예: microarray or RNA-seq data analysis) 이것이 보편적인 것은 아니기 때문이다. 처리가 끝난 데이터프레임을 data3이라 하자.
> data2 = data[1:100,] > data3 = t(data3) > View(data3) > dim(data3) [1] 48 100
컬럼의 배열 순서 뒤집기
marker gene(g1, g2, g3,… gn)을 사용하여 LS-BSR을 실행한 뒤 OUT_bsr_matrix.txt(OUT은 -x OUT으로 설정)을 열어보니 gn, gn-1, gn-2,… g1 순서로 row가 배열되어 있었다. t() 함수를 써서 이를 컬럼 단위로 배열하면 원래 의도했던 유전자의 역순이라서 불편하다. 만약 컬럼의 배열 순서를 역으로 바꾸고 싶다면 다음과 같이 하라.
> data2_r = data3[,ncol(data2):1]
컬럼의 배열 순서를 임의로 바꾸기
myord라는 벡터를 하나 생성하여 컬럼 이름을 여기에 순서대로 넣는다. 그리고 다음과 같이 하면 그만이다. 이는 row에 대해서도 똑같이 적용할 수 있다. 이렇게 쉬운 것을 모르고 for 구문을 쓸 궁리를 하다니!
> data_new = data3[, myord]
특정 row만을 뽑아내어 후속 분석에 사용하고 싶다면
위에서는 단순히 맨 앞의 100개 row를 뽑아내는 방법을 소개하였다. 만약 어떤 목록에 들어있는 유전자에 대한 row만을 추출하여 subset을 만들고 싶다면 어떻게 하면 좋을까? 뽑아낼 유전자는 genes라는 벡터에 들어있다고 가정하자. 다음과 같이 %in% operator를 사용하면 아주 쉽게 목적한 바를 이룰 수 있다. 추출된 데이터는 transpose를 거쳐서 최종적으로 data4 데이터프레임에 저장하였다. 이 기법은 [정해영의 블로그] 하루에 한 R 시리즈에도 상세하게 소개하였다. Row name이 지정되지 않은 상태라면 data$V1과 같이 variable을 지정하면 된다.
> genes [1] "pmx_cluster_contig_1" "pmx_cluster_contig_10" "pmx_cluster_contig_13" [4] "pmx_cluster_contig_14" "pmx_cluster_contig_15" "pmx_cluster_contig_6" [7] "pmx_cluster_contig_7" "pmx_cluster_contig_8" "pmx_cluster_contig_9" > temp = data[which(row.names(data) %in% genes),] > data4 = t(temp) > View(data4) > dim(data4) [1] 48 9
테이블 병합하기
<key>,<value>를 여러 row에 걸쳐 수록한 텍스트 파일이 N 개 존재한다고 가정하자. 이를 R에서 데이터프레임으로 한꺼번에 읽어들인 뒤 key를 기준으로 병합하는 방법을 알아보자. 최종 결과물은 <key>,<value 1>,<value 2>,<value 3>,…<value N> 형식이 되기를 원하는 것이다. 예를 들자면 여러 FASTA 파일을 phylosift로 분석하여 각 샘플에 대해 동정된 마커의 수를 수록한 marker_summary.txt를 하나로 병합하는 것이다. 이에 대해서는 [정해영의 블로그] PhyloSift 결과 종합하기(tab-delimited CSV 파일의 일괄 입력 및 병합 등)에서 자세히 다루었다.
그런데 모든 텍스트 파일의 key가 동일한 set이 아니라면 위에서 소개한 방법으로는 잘 되지 않을 것이다.
매트릭스 형태의 텍스트 파일을 R에서 읽어들이기(2)
이번에는 pyani를 사용하여 만든 ANI matrix를 가공하여 특정한 컬럼만을 추출하는 방법을 공부해 본다. 원래 이 일은 iTOL에 업로드하여 phylogenetic tree 옆에 ANI heatmap을 그려넣기 위하여 시작한 것이다. 서로 밀접하게 관련된 몇 개의 speies에서 유래한 51개의 스트레인의 유전체 서열을 사용하여 ANI 계산을 한 뒤, 5개의 type strain에 대한 ANI 값만을 뽑아내어 51 x 5 매트릭스를 만들고자 한다. 이 작업은 생각보다 복잡하였다. 컬럼 네임을 이용하여 특정 컬럼만을 뽑아내는 column slice 작업은 간단하지만 컬럼 이름을 바꾸고 user defined order로 정렬하는 것은 쉽지 않았다. 그러나 최종 구현된 코드는 매우 단순하다. 만약 dplyr 패키지를 사용한다면 좀 더 논리적인 접근이 가능할 것이나, 되도록 R base package를 써 보자는 것이 나의 철학이다.
ANIb_percentage_identity.tab 파일의 구조는 다음과 같다.
column과 row의 이름은 assembly accession이다.
> data = read.table("ANIb_percentage_identity.tab",sep="\t",header=T)
> rownames(data) = data[,1]
> data = data[,-1]
# 사실 위의 세 줄은 지극히 초보적인 방법이다. 다음 한 줄이면 족하다.
# data = read.table("ANIb_percentage_identity.tab",sep="\t",header=T,row.names=1)
> #View(data)
# 대각선에 대하여 대칭인 구조이지만 query -> subject 방향에 따라 identity가 똑같지 않다.
# 따라서 transpose를 하여 두 데이터프레임을 더한 뒤 2로 나누는 편법으로 평균을 구한다.
> data2 = t(data)
> data3 = (data + data2) / 2
> key = c("GCF_000146875.3","GCF_000217775.1","GCF_000236805.1","GCF_002240415.1","GCF_000235585.1")
# 다음의 값으로 나중에 데이터프레임의 컬럼 이름을 바꾸어야 한다.
> names(key) = c("E681","P.polymyxa","P.peoriae","P.kribbensis","P.terrae")
# column slicing by column name
> data4 = data3[,which(colnames(data3) %in% key)]
# data4 데이터프레임의 컬럼 이름을 바꾼다. key 벡터의 각 원소값을 이용하여 name을 구하는 것이 관건이다.
> for (i in names(data4)) {
# i = c("GCF_000146875.3","GCF_002240415.1",...)
+ temp = names(key)[key==i]
+ names(data4)[names(data4)==i] = temp
+ }
# user-defined ordering
> data4[, names(key)]
#View(data4) # 가끔 데이터 구조를 확인하라.
> write.table(data4,"ANI.key.txt",sep=",")
이 사례에서는 이름을 바꿀 컬럼이 몇 개 되지 않았으므로 named vector를 사용하여 원하는 목적을 달성하였다. 만약에 목록이 길어진다면 이렇게 할 수는 없고, 별도의 외부 파일에 저장된 값을 불러와서 작업해야 할 것이다. 최종적으로 저장된 데이터 파일은 다음과 같다. 컬럼의 수는 여기에서 보이는 것이 전부이다.
gplots 패키지의 heatmap.2() 함수 이용
> library(gplots)
> heatmap.2(as.matrix(data4),col=greenred(10),trace="none",scale="none",dendrogram=c("row"))
dendrogram=c(“row”)는 column과 row 전부 clustering을 하지만 단지 row의 dendrogram을 보이지 않을 뿐이다. Colv=F라고 하면 경고 메시지가 나오지만 컬럼의 clustering을 하지 않는다. 그림은 잘 그려졌으나 각 축의 라벨이 너무나 길어서 그림 영역을 벗어난다. 이를 조정하여 heatmap의 완성도를 높이려면 R에서 plot의 margin을 건드리는 법을 알아야 한다. 일단은 [정해영의 블로그] heatmap.2 그림을 조정해 보자 또는 R 그래픽스 생기초 - plot area와 margin 등의 이해(매우 중요!)를 참조하라.
다음은 레이아웃을 변경하는 예제이다. 여러가지 방법으로 응용해 보라. 가장 유용한 인수 설정은 density.info=“none”였다. 이를 적용하면 Color Keys and Histogram(두 줄)이 Color Key로 바뀌면서 텍스트 라벨과 컬러 키가 닿는 불편한 모습이 사라진다.
> heatmap.2(as.matrix(d2),col=greenred(10),trace="none",scale="none",dendrogram=c("row"),cexCol=1,cexRow=0.2,margins=c(5,20),lhei=c(1,5),density.info="none")
색상 체계가 마음에 들지 않는다면 다음과 같이 팔레트를 자체 정의할 수 있다.
> my_palette = colorRampPalette(c("turquoise", "yellow", "red"))(n = 100)
> heatmap.2(x, col=my_palette,...)
heatmap.2()로 그려지는 plot의 각 요소(heatmatp, row dendrogram, column dendrogram, key)의 위치와 크기를 조절하는 방법은 다음을 참조하라.
Label에 색깔 입히기
Row label에 쓰인 값들을 사전에 특별한 규칙에 의해서 그룹으로 묶어 두었다고 가정하자. 실제로 heatmap에서 dendrogram으로 묶이는 것과는 어떤 차이가 있을까? Row label에 그룹별로 색깔을 입힐 수 있다면 가능하다. key_info.txt 파일을 준비하여 다음과 같은 정보를 기록하라. 여기에서는 4개 컬럼으로 구성하여 첫번째가 key, 네번째가 color인 것으로 정하였으나 이는 사용자 맘대로이다. 다음을 a 라는데이터프레임에 저장하였다고 가정하자.
key1 info1-a info1-b black key2 info2-a info2-b red key3 info3-a info3-b blue
색깔에 해당하는 컬럼을 cols = a[,4] 명령으로 슬라이스한 뒤 heatmap.2(x,…,colRow=cols,) 함수를 사용하면 된다. 단, key_info.txt와 데이터프레임 x 내부에서 key는 같은 순서로 존재해야 한다. R에서 사용 가능한 컬러에 대해서는 R color cheetsheet를 참조하라. R에는 총 657개의 내장 색상을 갖고 있으며 colors() 명령으로 이를 확인할 수 있다.
Heatmap 그림과 phylogenetic tree의 결합
위에서 보인 그림은 고작 9개의 유전자에 대한 blast score를 가지고 왼쪽의 dendrogram을 그린 것이었다. 각 strain의 실제 phylogenetic tree를 별도의 방법으로 만든 다음 이를 그림의 왼쪽에 배치하면 더욱 바람직한 것으로 생각된다. 본격적인 phylogentic tree는 Roary에서 얻어진 core genome alignment로부터 FastTree로 그려도 되고, 아니면 Harvest suite를 이용하여 그려도 된다. 그러나 유의할 점이 있다. 분석을 위해서 내가 수집한 genome들은 서로 충분히 유사해야 하나의 core genome SNP를 이용한 그림에 들어가게 된다. Harvest suite의 parsnp을 예로 들자면 기본적으로 MUMi distance ≤ 0.01인 것만을 분석에 포함하게 된다. 만약 너무나 거리가 먼 genome이 섞여있다면 제대로 된 tree가 그려질 수 없다.
gplots 패키지의 heatmap.2() 함수를 사용하여 그려진 dendrogram을 뽑아내는 것은 구글을 통해 뒤져보면 방법을 찾을 수 있다. 아래에서 사용한 pyani matrix 파일은 첫 줄이 tab으로 시작한다.
> library(ape)
> d = read.table("matrix_pyani.tab",sep="\t",header=T,row.names=1)
> fit = hclust(dist(as.matrix(d)),method="average")
# heatmap.2() 함수는 complete agglomeration method를 클러스터링에 사용한다. average method는 UPGMA와 같다.
> my_tree=as.phylo(fit)
> write.tree(phy=my_tree,file="test.newick")
How to create a newick file from a cluster in R?
그 반대로, 다른 도구를 통해 만든 tree를 heatmap.2() 함수에 공급할 방법은 없을까? 내가 알기로 그것은 불가능하고, 아래에서 소개할 iTOL을 이용하는 것이 바람직하다. 물론 이것이 유일한 방법은 아니다. 바로 ComplexHeatmap 라이브러리가 있다.
다음은 보다 상세한 heatmap 작성 관련 기술 문서이다.
Heatmap in R: Static and interactive visualization - heatmap(), heatmap.2(), interactive heatmap using 3dheatmap(), complexheatmap…