Multifasta에서 빠르게 일부 레코드만 추출해오기

이런 작업은 사실 일상적으로 하는 작업인 관계로 Perl이나 Python 과 같은 스크립트 언어를 조금이라도 아는 사람이라면 그닥 어렵지 않게 할 수 있는 일이긴 하지만, 이런 작업도 수GB 이상의 파일이라면 그리고 불러올 파일의 갯수가 늘어나면 적지 않은 시간이 소요된다. 좀 더 빠르게 이런 작업을 마칠 수 없을까?

Profile HMM 에 의해 도메인 및 단백질 검색을 하는 유명한 소프트웨어인 HMMER 3.0 에 딸려있는 유틸리티인 esl-sfetch를 이용하면 아주 빠르게 이런 작업을 할 수 있다.

HMMER3 을 설치하고 hmmer3 의 실행파일에 PATH를 걸어둔 다음

esl-sfetch -h

하면

$ esl-sfetch -h
# esl-sfetch :: retrieve sequence(s) from a file
# Easel h3.0 (March 2010)
# Copyright (C) 2010 Howard Hughes Medical Institute.
# Freely distributed under the Janelia Farm Software License.
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Usage: esl-sfetch [options]          (one seq named )
Usage: esl-sfetch [options] -f   (all seqs in )
Usage: esl-sfetch [options] --index        (index )

 where general options are:
  -h     : help; show brief info on version and usage
  ....중략.

가령 지금 나는 vibrio.fasta라는 274MB짜리 fasta 파일이 있고, 이 파일은

$tail vibrio.fasta
SDGATWLDRRLIHGETADLAPTGFGQQVREAMDQRREHHIEQGDATRSRDSRVFYRRNLLAILREREVAG
VGSDMALSKGLPFRAATDGESVSGKFTGTVHLSSGKFAVVEKSHEFTLVPWRPIIDRQLGREVMGIVQGG
SVSWQLGRQRGLER
>gi|424025505|ref|ZP_17765144.1| drug resistance transporter, Bcr/CflA subfamily protein [Vibrio cholerae HC-62B1]
MTTTRPAWAYTLPAALLLMAPFDILASLAMDIYLPVVPAMPGILNTTPAMIQLTLSLYMVMLGVGQVIFG
PLSDRIGRRPILLAGATAFVIASLGAAWSSTAPAFVAFRLLQAVGASAMLVATFATVRDVYANRPEGVVI
YGLFSSMLAFVPALGPIAGALIGEFLGWQAIFITLAILAMLALLNAGFRWHETRPLDQVKTRRSVLPIFA
SPAFWVYTVGFSAGMGTFFVFFSTAPRVLIGQAEYSEIGFSFAFATVALVMIVTTRFAKSFVARWGIAGC
VARGMALLVCGAVLLGIGELYGSPSFLTFILPMWVVAVGIVFTVSVTANGALAEFDDIAGSAVAFYFCVQ
SLIVSIVGTLAVALLNGDTAWPVICYATAMAVLVSLGLVLLRLRGAATEKSPVV

요런 식으로 되어 있는 multifasta 파일이다. 여기서 gi|424025505|ref|ZP_17765144.1| 와 같은 accession 으로 파일을 찾고싶다고 치자.

gi|260362785|ref|ZP_05775654.1|
gi|260880072|ref|ZP_05892427.1|
gi|28901226|ref|NP_800881.1|
gi|153835952|ref|ZP_01988619.1|
gi|308095182|ref|ZP_07663192.1|
gi|254224936|ref|ZP_04918551.1|
gi|229529188|ref|ZP_04418578.1|
gi|254285269|ref|ZP_04960234.1|
gi|421351525|ref|ZP_15801890.1|
gi|229515173|ref|ZP_04404633.1|

요런 내용을 list.txt로 저장했다고 치자
제일 먼저 할 일은 fasta 파일을 인덱싱하여 빨리 검색을 할 수 있게 하는거.

$esl-sfetch --index vibrio.fasta
Creating SSI index for vibrio.fasta...    done.
Indexed 665311 sequences (665311 names).
SSI index written to file vibrio.fasta.ssi

이러면 파일이름.ssi 식의 인덱스 파일이 생기고 이제 검색을 할 수 있다.

$esl-sfetch -f vibrio.fasta list.txt

하면 바로 쫙 stdout 으로 파일이 추출된다. ㅋ

사실 이런 작업은 hmmer 등을 수행해서 나온 결과파일을 가지고 해당하는 fasta 파일을 데이터베이스에서 추출할때 유용한데, hmmer를 수행해서 테이블 형식으로 다음과 같은 출력파일이 나왔다고 하자.
Screen Shot 2012-12-21 at 12.59.44 AM

이 출력물을 awk 등으로 파싱해서 원하는 gi 값만 esl-sfetch에 전달해 주면 곧바로 hmmer hit 에 상응하는 시퀀스 파일을 뽑을 수 있다.

awk '{if (substr($0,0,1)!="#") {print $1}}' VopL.fasta.hmmer | esl-sfetch -f vibrio.fasta - > hmmerresults.fasta

참고로 samtools 에도 faidx 라는 동일한 기능 (fasta 를 indexing한후 추출) 이 있으나 테스트해 본 결과 esl-sftech가 더 빠르다.
약 247MB의 fasta에서 10개 정도의 fasta를 추출할때 시간을 재보면

time awk '{if (substr($0,0,1)!="#") {print $1}}' VopL.fasta.hmmer | xargs samtools faidx vibrio.fasta
real	0m1.033s
user	0m0.921s
sys	0m0.110s

time awk '{if (substr($0,0,1)!="#") {print $1}}' VopL.fasta.hmmer | esl-sfetch -f vibrio.fasta
real	0m0.007s
user	0m0.005s
sys	0m0.006s

수십 GB 정도의 파일이라면 좀 더 차이날 것 같다.;;

-C 옵션을 이용하면 multifasta 내의 subsequence을 잘라오는 것도 가능하다. (BED file 등을 파싱해서 빠르게 시퀀스를 얻는것도 가능..아 그건 bedtools쓰면 되나..여튼.

어쩄든 NGS 시대가 되면서 처리해야 할 데이터의 양이 부쩍 늘어난 상황에서 파싱과 같은 간단한 작업의 경우에도 가능하면 unix 기반의 텍스트 툴및 파이프, 리다이렉션을 이용하는 것이 수행속도 면에서도 훨씬 이점이 있는 것 같다는.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s