[하루에 한 R] ANI 계산 결과 다루기 - column 혹은 row 이름을 통한 조작

R을 이용해서 '4차 산업혁명' '빅데이터' 시대에 걸맞는 거창한 일 - 말하자면 방대한 데이터 속에서 새로운 지혜를 발굴하는 - 을 해야 하는데, 아직까지도 기초 수준의 자료 작성 기법을 익히는데 머물고 있으니 어찌 한심하지 않겠는가!

기본적인 자료 조작을 할 줄 알아야 통계분석, 시각화 등등의 고급스런 일을 비로소 할 수 있다. 앞으로 더 멋진 일을 하기 위한 기초 소양을 쌓아나간다고 생각하고 오늘도 부끄럽지만 [하루에 한 R]을 기록해 나간다.

매트릭스 형태의 자료 파일(csv, txt, tab-delimited...)를 R에서 데이터프레임으로 읽어들인 다음에 어떤 일을 할 수 있을까? 그동안은 각 cell을 채운 자료 자체에 집중을 하는 편이었지만 요즘은 컬럼과 로우(한글로는 '로'인가 '로우'인가? '윈도'가 표준 표기법이니 '로'라고 하는 것이 맞겠지만 어색함을 떨치기 어렵다, '로우'라고 쓰자.)의 이름을 활용하는 것이 더 흥미가 있다.

• key,value1(value2,value3...) 형태의 파일 여러개를 key 값 기준으로 합치고 싶다(key 값은 동일할 수도 있고 그렇지 않을 수도 있다).
• 컬럼이나 로우의 이름을 기준으로 재배열하고 싶다.
• 컬럼이나 로우의 이름을 일괄적으로 바꾸고 싶다(스트링 조작 혹은 외부 데이터 사용).
• 컬럼이나 로우의 이름을 이용하여 데이터 프레임의 일부를 덜어내고 싶다(slicing).
• 최대값(혹은 최저값)에 해당하는 컬럼 혹은 로우의 이름을 알아내고 싶다.
• 특정 로우 그룹의 평균값, 최대값(혹은 최저값)을 추출하고 싶다.

내가 종사하는 생명과학 분야, 특히 유전체학과 관련한 분야에서는 microarray data, RNA-seq (expression value) data, OTL table, BLAST result(tabular format) 등 다양한 종류의 매트릭스 형태 데이터 파일을 접하게 된다. 대개 하나의 로우는 하나의 레코드에 해당하지만 생명과학 분야 데이터에서는 레코드가 컬럼 단위로 기록되는 경우가 많다. 따라서 hclust() 함수를 쓰거나 할 때 트랜스포즈를 미리 하는 등의 변환이 필요함을 염두에 두어야 한다.

아마 과거에 작성한 글에 이미 이에 대한 방법을 일부 기록해 놓은 것이 있을지도 모른다. 오늘 쓰는 글은 최종 정리판이라고 보면 된다. 사용할 데이터 파일은 pyani에서 만든 ANIb_percentage_identity.tab(다운로드 링크)이다. 어떻게 생겼는지 살짝 열어보자. OpenOffice를 사용해서 command line에서 이런 종류의 파일을 한 번에 스프레드시트로 열 수 있다(ooffice --calc ANIb_percentage_identity.tab). 

R에서 읽어들이려면 다음과 같이 하면 된다. 컬럼 및 로우 네임을 정확히 지정하는 것, 그리고 대각선을 중심으로 대칭이 되게 만드는 것이 관건이다. 이 경우에는 컬럼과 로우의 네임 순서가 정확히 같다는 것을 가정하고 있다.

> d = read.table("ANIb_percentage_identity.tab",sep="\t",header=T,row.names=1)
> d2 = t(d)
> d3 = (d + d2)/2

만약 컬럼과 로우의 이름 순서가 일치하지 않는다는 의심이 든다면 다음과 같이 하면 된다. 전처리가 끝난 데이터프레임의 이름은 d3임을 기억해 두자. 자료 유형이 데이터프레임인 경우 colnames()나 names()나 차이가 없다.

> d3 = d3[, order(names(d3))]
> d3 = d3[order(row.names(d3)), ]

다음으로는 컬럼/로우 이름이 너무 길다는 것이 불만이다. 전부 동일 종에서 유래한 서로 다른 스트레인의 유전체 정보를 사용한 것이니 종 이름을 남겨둘 이유가 없다. 이는 사용자 함수를 만들어서 적용하면 된다.

> a = names(d3)
> trim.species = function (x) sub("^Lactobacillus_rhamnosus_","",x)
> names(d3) = trim.species(a)
> rownames(d3) = trim.species(a)

row.names()와 rownames() 함수는 무엇이 다른가? 이것은 천천히 알아보자. colnames() 함수는 있는데 col.names() 함수는 존재하지 않는다(?col.names를 입력하면 아무런 문서도 나오지 않으므로 그렇게 생각하는 것이다). 데이터프레임의 외형을 다듬는 일은 전부 끝났다. 이번에는 특정한 이름을 갖는 컬럼 하나를 뽑아서 어떤 데이터가 들어있는지 보고싶어졌다. 방법은 다음의 두 가지가 가능하다.

> d3[, which(names(d3) == "DSM_20021_GCF_001435405.1")]
> d3[, which(names(d3) == "DSM_20021_GCF_001435405.1"),drop=F]

> d4 = d3[, which(names(d3) == "DSM_20021_GCF_001435405.1"),drop=F]

drop=F 파라미터를 하나 추가하는 것이 훨씬 바람직한 결과를 낳는다. d4의 구조를 살펴보라. 컬럼과 로우의 이름이 전부 지정된 상태이다.

컬럼의 이름을 완벽하게 지정할 필요도 없다. DSM_20021로 시작하는 이름을 갖는 컬럼을 골라내도록 하자.

> d3[, grep("^DSM_20021",names(d3)),drop=F]

컬럼이 여러 개 선택될 가능성도 염두에 두어야 한다. 다음으로는 뽑아낼 컬럼이 여러개인 경우를 생각해 보자. 해당되는 컬럼의 이름은 벡터에 미리 지정해 둔다. 정규식 패턴을 벡터에 저장하여 적용하는 것은 아마 어려울 것이다.

> vec = c("Pen_GCF_00207655.1","ATCC_8530_GCF_000233755.1","4B15_GCF_002158925.1")
> d3[, which(names(d3) %in% vec)]

다음으로는 특정 컬럼에 대해여 최대값이 들어있는 cell의 로우 네임을 뽑아보는 일을 해 보겠다. 쉽게 해서 특정 유전체에 대하여 ANI가 가장 높게 나오는 것은 누구인지를 찾아내자는 것이다. 그런데 조심할 것은 당연히 자기 자신에 대한 ANI가 가장 높으므로 이것이 출력되는 것을 배제해야 한다. JSpecies에서는 자기 자신에 대해서는 값을 계산하지 않지만 pyani는 정직하게 1을 산출해 놓는다. 이를 전부 NA로 바꾸어 놓자.

> d3[d3 == 1] = NA

각 컬럼에서 최대값을 갖는 것들의 row name을 찾아보자.

> names(d3)[apply(d3,1,which.max)]

이보다 미려한 출력을 만들고 싶다면 좀 더 연구를 해야 한다. 데이터 프레임의 컬럼 혹은 로우 네임을 named vector를 사용하여 일괄적으로 바꾸려면 http://genoglobe.kr/kribb/data_visualization 를 참고하도록 한다.