[하루에 한 R] reshape 패키지를 이용하여 쌍(pairwise) 자료를 매트릭스로 전환하기

dRep(논문, GitHub)은 수십, 혹은 그 이상의 미생물 유전체를 신속하게 비교하여 중복을 제거하는 도구이다. 메타게놈 데이터의 처리에도 유용하지만 다량의 유전체 자료를 비교하여 유사도(gANI) 기준으로 샘플을 묶는 데에도 편리하게 쓰일 수 있다. 그러나 pyani 처럼 매트릭스 형태의 자료를 제시하지 않는다. ANI 값들이 일단 매트릭스로 주어진다면 분석자의 입맛에 맞는 가공과 시각화가 가능하지 않겠는가?

결과 디렉토리 안에 활용 가능한 자료가 있을 것이라 생각하고 data_table 디렉토리를 열어보았다. 다음과 같은 네 개의 csv 파일이 존재한다.

Bdb.csv  Cdb.csv  Mdb.csv  Ndb.csv

하나씩 열어보면 무슨 자료를 담고 있는지 누구나 파악할 수 있다. 이 중에서 Ndb.csv가 모든 유전체를 1:1로 pairwise 비교를 하여 계산한 수치를 담은 파일이다. 파일의 형식은 다음과 같이 첫 줄에 나온다. 왜 query가 아니고 querry인가? 내가 잘못 타이핑한 것이 아니라 실제 결과 파일을 복사해서 여기에 가져다 놓았을 뿐이다.

alignment_coverage,ani,querry,reference,MASH_cluster

굵은 글씨로 표현된 컬럼만 뽑아내면 된다. 이는 awk를 적절히 사용하는 것으로 가능하다.

$ awk -F"," -v OFS="\t" '!/^alignment/{print $3, $4, $2}' Ndb.csv > matrix.txt

query와 reference 컬럼은 실제 dRep 계산에 사용한 파일 이름 그대로이다. 나의 경우는 다음과 같다.

Lactobacillus_rhamnosus_116_GCF_000801045.1.fna

너무 길어서 보기가 좋지 않다. Lacto..에서 GCF까지를 non-greedy하게 제거하고, 다음으로 .fna 확장자를 제거하는 것이 좋을 것이다. GCF_까지 없애면 나중에 R로 읽어들였을 때 숫자로 인식을 하여 약간 번거로워진다. non greedy하게 매치된 문자열을 없애려면 vim의 last line mode에서 다음과 같이 입력한다.

:1,$s/Lacto.\{-}_GCF//g
:1,$s/.fna//g

그 다음으로는 R에서 데이터프레임으로 읽어들여서 처리하면 된다. Ndb.csv는 두 샘플이 각각 query와 reference일 때의 gANI 값을 두 개의 라인으로 기록해 놓은 것이다. 따라서 매트릭스로 전환한 다음에는 대각선에 대하여 대칭이 아니다. 이는 t() 함수로 역행렬을 만든 것을 자기 자신에게 더한 뒤 둘로 나눔으로써 해결한다.

R 코드는 다음과 같다. reshape 패키지(링크)가 큰 도움이 되었으며, 힌트는 여기에서 얻었다. 이 패키지는 R 개발자로 널리 알려진 해들리 위캠의 작품이다.

> install.packages("reshape")
> library(reshape)
> d = read.table(file="matrix.txt",sep="\t",header=F)
> d2 = cast(d,V1~V2)
> rownames(d2) = d2[,1]
> d2 = d2[,-1]
> d3 = (d2 + t(d2))/2
# 확인용
> dim(d3)
> View(d3)

[업데이트]

위에서 보인 R 코드는 언제나 적용 가능한 것이 아니다. 자세한 것은 [하루에 한 R] pair 형태의 데이터를 matrix로 전환할 때 주의할 점을 참고하라.

[2021년 8월 26일 업데이트]

무려 3년이 지난 상황에 위 R code의 오류를 발견하였다. cast() 함수로 읽어들인 데이터는 매트릭스의 대각선을 기준으로 절반에 대한 것이다. 나머지 절반은 같은 값으로 채워야 한다. 따라서 d3 = (d2 + t(d2))/2가 아니라 d3 = (d2 + t(d2))가 되어야 한다. 으이그...