BLAST+(최신 버전)가 완전히 정착된 것이 언제인데 나는 아직도 'legacy' BLAST(v2.2.26)를 즐겨 사용한다. blastall 커맨드 라인 사용법에 너무나 길들여져 있기 때문이다. BLAST+의 옵션 설명서를 인쇄하여 책상 앞에다 붙여놓아야 겨우 익숙해질까?
예전, 하지만 그렇게 오래 전도 아니던 시절, 미생물 유전체의 functional annotation을 하려면 Swiss-Prot이나 NCBI NR 데이터베이스를 로컬 컴퓨터에 받은 다음 수천 개의 query를 blastp로 날리는 일을 수시로 했었다. 요즘은 RAST server에 올려 버리거나, 로컬 머신을 쓴다 하여도 Prokka를 쓰면 단 몇 분에 해결이 되이 큰 용량의 공개 단백질 서열 DB에 대해서 BLAST를 직접 돌릴 일이 많지 않다. 또 그럴 일이 있으면 유료 서비스에 가입한 Blast2GO에서 검색을 수행하면 된다.
그래도 간혹 연구 목적을 위해 특별한 목적으로 만든 서열 데이터베이스에 대해서 BLAST를 실행할 일이 있다. 그 결과물을 파싱하는 프로그램을 직접 구현하는 것이 예전에는 매우 보편적인 생명정보학 실습 과제이기도 했다. 누가 그랬더라? 이제 제발 그런 거 하지 말라고...
사람의 눈으로만 보기 편한 BLAST 결과물을 컴퓨터에게 읽히려고 하니 프로그래밍으로 처리하기에 얼마나 불편하겠는가? 사실 내가 tabular output format(blastall -m 8 또는 -m 9; BLAST+에서는 -outfmt 6)에 관심을 가진 것도 얼마 되지 않았다. 이를 엑셀로 불러들이면 완벽하지는 않더라도 대량의 결과를 처리하기에 좋으니까 말이다. 엑셀의 귀재라면 각 query에 대한 hit 중에서 best를 골라내는 방법을 함수로 구현할 수 있을지도 모른다.
R을 사용하면 BLAST의 tabular format output으로부터 각 query에 대한 best hit을 쉽게 뽑을 수 있을 것이라는 생각을 오랫동안 하고 있었다. 사실 R 환경 자체에서 BLAST를 실행할 수 있는 rBLAST라는 패키지가 있으니 이 환경 안에서 best hit을 찾을 수 있을 것이고 또 더 검색을 해 보면 다른 도구가 있을 것이다. 범위를 좀 더 넓힌다면 reciprocal best hit을 찾아주는 유틸리티도 있을 것이다. 하지만 나는 좀 더 '하등한' 수준의 해결 방법을 알고 싶었다.
[SEQanswers] How to output only best-hit result using standalone blast <= 심지어 sort 명령어를 사용한 방법도 있다. one-liner를 즐기는 사람에게는 좋은 소일거리가 될 것이다.
그러면 본론으로 들어가서, blastall -m 9(tabular with comment lines; #으로 시작) 명령으로 작성한 결과물을 R에서 읽어들여서 best hit을 찾는 방법을 공부해 보자. 예전에 마이크로어레이 자료 분석을 공부하면서 control spot 때문에 한 유전자에 대한 형광 측정값이 여럿 존재하는 경우 이것의 median 값을 취하도록 R 코드를 짰던 것 같은데 도저히 기억이 나질 않는다. 아마 apply() 혹은 by() 계열의 함수를 쓰지 않았었나 생각한다.
그런데 구글링을 계속 해 보니 base R만을 이용한 아주 원초적인 방법을 제시한 글이 있어서 소개하고자 한다. 이 안에는 외부 패키지와 함수를 이용한 복잡한 방법도 포함되어 있다. 공부를 위해서 나머지 방법도 따라서 해 보아야 한다.
[Stack Overflow] Extract the maximum value within each group in a dataframe
blastall -m 9 옵션을 사용하여 출력 파일(tblastn.m9,.txt)을 만든 뒤 각 query에 대한 best hit을 뽑아보자.
> d = read.table("tblastn.m9.txt",sep="\t",comment.char="#",header=F)V1(query)의 값이 동일한 hit group에 대해서 V12(bit score)가 최대인 것을 출력하는 것이다. 어떻게 이보다 더 간단명료할 수 있겠는가?
> d
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12
1 BfmR bfmR 100.00 238 0 0 1 238 1 714 8.0e-156 423.0
2 BfmR bfmS 30.56 36 25 0 21 56 1437 1330 2.3e-01 22.3
3 BfmS bfmS 99.82 549 1 0 1 549 1 1647 0.0e+00 1083.0
4 BfmS bfmS 27.78 54 39 1 133 186 658 533 7.1e+00 18.9
5 BfmS nfu 38.46 13 8 0 117 129 442 480 8.1e+00 18.5
6 CsuA csuA 100.00 192 0 0 1 192 1 576 2.0e-129 353.0
7 CsuA csuE 43.75 16 9 0 177 192 970 1017 4.5e-01 20.8
8 CsuB csuB 100.00 172 0 0 1 172 1 516 4.0e-114 312.0
...중간 생략
> do.call(rbind, lapply(split(d,d$V1), function(x) {return(x[which.max(x$V12),])}))
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12
BfmR BfmR bfmR 100.00 238 0 0 1 238 1 714 8e-156 423
BfmS BfmS bfmS 99.82 549 1 0 1 549 1 1647 0e+00 1083
CsuA CsuA csuA 100.00 192 0 0 1 192 1 576 2e-129 353
CsuB CsuB csuB 100.00 172 0 0 1 172 1 516 4e-114 312
CsuC CsuC csuC 100.00 277 0 0 1 277 1 831 0e+00 568
CsuD CsuD csuD 100.00 832 0 0 1 832 1 2496 0e+00 1650
CsuE CsuE csuE 100.00 339 0 0 1 339 1 1017 0e+00 666
EntA EntA entA 100.00 256 0 0 1 256 1 768 0e+00 502
EpsA EpsA epsA 99.73 366 1 0 1 366 1 1098 0e+00 727
NfuA NfuA nfu 100.00 212 0 0 1 212 1 636 2e-162 438
OmpA_partial OmpA_partial ompA 100.00 315 0 0 1 315 1 945 0e+00 592
pgaA pgaA pgaA 100.00 747 0 0 1 747 1 2241 0e+00 1544
PgaB PgaB pgaB 99.80 510 1 0 1 510 1 1530 0e+00 980
PgaC PgaC pgaC 100.00 392 0 0 1 392 1 1176 0e+00 805
PgaD PgaD pgaD 100.00 154 0 0 1 154 1 462 6e-101 278
PlcD PlcD plcD 100.00 487 0 0 1 487 1 1461 0e+00 1018
Ptk Ptk ptk 100.00 727 0 0 1 727 1 2181 0e+00 1372
댓글 없음:
댓글 쓰기