[하루에 한 R] Expression data로 heatmap 그리기

오래전의 microarray 실험 데이터를 가지고 뒤늦은 논문을 쓰고 있다. 대단한 발견이라고 할만한 것은 없고, 너무 오래되어서 더 이상 시의성이 없는 데이터로 추락(?)하기 전에 일단락하여 발표하기 위함이다. 사실 과학에서는 과거의 연구성과가 상당한 시간이 흐른 뒤에 발굴되어 재평가되면서 과학계의 역사를 새로 쓰는 일이 다반사이다. 따라서 변신에 변신을 꾀하면서 지나치게 유행을 좇아 연구 분야 변신을 하는 것은 그렇게 바람직하지는 않다고 본다. '시대적 요청'이라는 것이 있기는 하지만 말이다.

heatmap 또는 heat map이란 데이터를 2차원 형태로 늘어놓고 각각의 값을 색으로 표현한 데이터 시각화 기법의 하나이다. 데이터의 배열을 행 또는 열에 따라 적절히 조절(clustering)함으로써 눈에 확 뜨이는 두드러지는 패턴을 찾아낼 수 있다.

21C 미생물프론티어 사업단 시절을 풍미했던 바로 그 대장균 microarray를 이용한 데이터를 이용해 보고자 한다. 한 개의 slide에 같은 프로브 셋이 두 개의 블록으로 나뉘어 찍혀 있으므로, 칩 하나로부터 두 개의 샘플 데이터가 얻어진다. 엄밀히 따지면 이렇게 나온 데이터는 technical replicate에 해당한다. 그러나 편의상 마치 biological replicate인 것처럼 다루었다. 이를 보고 '말도 안돼!'라고 외치실 분이 있을지도 모른다^^ 이번 포스팅에서는 R을 이용하여 QC, normalization, DEG 추출 등과 같은 expression analysis를 하는 방법을 다루고자 함이 아니다. 단지 hierarchical clustering을 동반한 heatmap을 그리는 작업을 설명하기 위해 고전적인 microarray data를 이용하는 것일 뿐이다.

• BioProject PRJNA206430
• GEO platform GPL7395
• GEO accession GSE47589
• Eighteen samples: GSM1152845-GSM1152862

처음에는 CLC Genomics Workbench의 Transcriptomics Analysis(과거에는 "Expression Analysis")로 heatmap을 그려 보았는데 도무지 마음에 드는 그림이 나오지 않았었다. CLC에서는 GEO 데이터를 어떻게 임포트하는가? GSE####_series_matrix.txt.gz 형식의 series matrix 파일을 다운로드하였다면 압축을 풀고 Import -> Standard Import 메뉴에서 파일 유형을 Expression array/data: GEO SOFT format series file (.txt/.text)로 설정하여 해당 파일을 읽어들이면 된다.

그러면 R을 사용해서 좀 더 예쁜 heatmap을 그려보자. R에 기본으로 내장된 heatmap() 함수는 그다지 예쁘지 않은데다가 사용자 입맛대로 조정을 할 여지가 별로 많지 않다. 따라서 gplots 패키지가 제공하는 heatmap.2() 함수가 인기를 끄는 모양이다. HeatPlus 라는 패키지도 있으니 참고하면 되겠다. 몇 가지 유용한 링크부터 먼저 소개하고 진행하겠다.

• gplots 매뉴얼(pdf 파일)
• R: heatmatp with gplots
• A short tutorial for decent heat maps in R
• Making heatmaps with R for microbiome analysis
• Heatmaps in R, two ways

입력물로는 GEO에서 다운로드한 series matrix 파일을 이용하겠다. !로 시작하는 코멘트 라인은 전부 제거하여 헤더와 데이터 부분만을 남겼다. 이 파일은 이미 normalization과 log2 transformation이 된 상태이다. 파일명은 data.txt로 수정하여 접근하기 쉬운 폴더(C:\R_work)에 복사해 두었다. 첫번째 시도는 우선 되는대로 막 그리는 것.

> install.packages("gplots")
> library("gplots")
> setwd("C:R_work")
> data = read.table("data.txt", header=TRUE, sep="\t", row.names="ID_REF")
> jpeg("plot1.jpg")
> heatmap.2(as.matrix(data), col=greenred(10), trace="none")
> dev.off()

jpeg으로 출력한 그림을 보자. 어떤 문제점이 존재하는가? microarray data의 heatmap은 각 샘플이 얼마나 유사한지를 보여주기도 하므로 일종의 QC 개념으로 쓸 수도 있다. 동일 조건에서 나온 각 6개의 샘플들이 한데 묶이므로 일단 큰 문제는 없다. 단 하나의 슬라이드에서 나온 두 블록이 항상 가장 가까이 묶이지는 않았음을 기억은 해 두자. 샘플의 이름이 마진에 비해 너무 길어서 잘렸고, 4721개나 되는 probe의 수는 하나의 heatmap으로 표현하기에는 너무 많다. probe ID와 dendrogram이 제대로 보이질 않으니 말이다. 그리고 대부분의 데이터가 0 근처에 있어서 너무 밋밋한 그림이 되었다.

Plot의 마진을 설정하는 방법은 완전히 다른 주제라서 여기에서는 다루지 않겠다. 이제 그림을 좀 더 돋보이게 하는 방법을 알아보겠다. 

모든 컬럼에서 특정 조건(값의 범위)을 충족시키는 row만 남기기(미완?)

Two-color microarray experiment에서 얻어진 데이터이므로 log2(Fold_change)가 0 근처에 있는 것은 그림을 복잡하게만 만들뿐 heatmap에 남겨둘 필요가 없다. 특정 column의 값을 기준으로 하여 전체 row를 삭제하는 방법은 인터넷 검색에서 흔히 찾을 수 있다. 예를 들어 첫번째 컬럼의 값이 -0.5보다 크고 0.5보다 작은 것을 데이터에서 제외해 보자. 이를 "dropping"이라 한다.

> data.2 = data[-which(data[,1] > -0.5 & data[,1] < 0.5), ]

그렇다면 18개의 모든 컬럼에 대해서 -0.5 < 데이터값 < 0.5를 만족시키는 row를 일괄적으로 삭제하려면 어떻게 하면 좋을까? 어느 한 컬럼이라도 -0.5 .. 0.5의 범주에 들면 그 row를 삭제한다는 것이 아니고(이렇게 되면 너무 많은 row가 제거될 것이다), 모든 컬럼이 다 그럴 때에만 제거하지는 것이다. 데이터가 카운트 값이면 row 단위로 합하여 그것을 기준으로 하면 되겠지만 이번의 사례에서는 '범위'가 판단의 기준이라서 쉽지가 않다.

아주 단순하게는(그리고 미련하게는) which() 함수 내부의 조건에 해당하는 부분을 모든 컬럼으로 확장하면 된다. (data[,1] > -0.5 & data[,1] < 0.5) & (data[,2] > -0.5 & data[,2] < 0.5) & (data[,3] > -0.5 & data[,3] < 0.5) ... (data[,18] > -0.5 & data[,18] < 0.5) 오타가 안나도록 주의깊게 타이프를 치면 된다. 그런데 이 방법은 너무 우습지 않은가? 만약 dropping을 마친 row가 아직도 너무 많아서 cutoff 값을 바꾸고 싶다면? 구글링을 열심히 해 보았지만 이 경우 딱 맞는 솔루션은 보이질 않는다. R에서 별로 바람직하게 여겨지지 않는 for 문을 사용하여 어찌어찌 만들어 보았는데 최종적으로 선별된 프로브의 수가 너무 많아서 역시 구별이 가능한 수준의 그림이 만들어지질 않는다.

row를 값의 범위에 의해 선별하여 버리는 방법을 알아내는 것이 문제가 아니라, 남은 row 자체가 너무 많다는 것이 문제다. 결국은 논문 작업을 위해서 뽑아둔 DEG 목록을 활용하기로 하였다. 몇 그룹의 조건으로부터 총 813개의 DEG를 확보하여 all_DEG.txt 파일에 저장하였다. 한 줄에 하나씩의 probe ID가 들어있고 중복을 없애는 작업이 필요하다. scan() 함수를 사용하여 매트릭스 형태가 아닌 데이터 파일을 입력하는 것은 조금 복잡하다. 이 링크를 참조하여 추후에 조금 더 공부하도록 하자.

> deg = scan(file="all_DEGs.txt", what=list(probe=character()))
> deg.nr = unique(deg$probe)

여기까지 왔으면 heatmap 작성 대상 유전자(probe ID)의 목록이 deg.nr이라는 벡터에 들어있는 상태이다.

목록에 존재하는 row name에 해당하는 레코드만을 데이터프레임에서 추출하기

추출할 probe의 목록을 deg.nr이라는 벡터에 담았으니, 이를 참조하여 data 데이터프레임의 일부를 꺼내면 된다. 참조할 곳은 특정 컬럼이 아니라 row name에 해당한다.

> data.final = data[which(row.names(data) %in% deg.nr), ]
> heatmap.2(as.matrix(data.final), col=greenred(10), trace="none")

probe가 457개로 줄어들어서 한층 보기가 수월해졌다.

다음의 숙제는 plot의 마진을 설정하는 일이다. 아마도 par() 함수를 쓰게 되지 않을까 싶다.