[구버전] 人Co BLOG :: 나만 몰랐던 BLAST 꿀팁

19세기의 자전거란 현재와 달라서 금속의 바퀴로 이루어져 쉽게 사고로 이어지곤 했답니다.
이에 영국의 수학자 던롭은 사랑하는 외동아들 조니를 위해 고무 타이어를 발명하게 되고 공기타이어회사 CEO로 큰 부자가 되게 됩니다.
저 역시 필요 때문에 (더 빨리, 더 쉽게 처리하고 칼퇴하고자) 다음의 방법들을 발견하게 되었으니 평행이론에 따라 언젠가는 저도 큰 부자가 될..... 수 있을까요?

발단은 몹시 사소했습니다.
아래 그림처럼 결과를 만들면 매우 있어 보일 것 같았거든요.

<출처 : GraPhlAn 프로그램 위키 페이지>

그런데 내 손에 있는 건 단백질 서열 단 하나, 계통수를 그리기 위한 DB가 아직 구축되어 있지 않은 상황이었어요.
하지만 걱정 없죠. 백종원 대표님께 만능 간장이 있다면 우리에겐 만능 BLAST가 있으니깐요.
BLAST를 통해 맵핑되는 단백질들을 모은 후 이를 이용하여 계통수를 그려보기로 합니다.
BLAST에 대한 기본 설명이나 분석 방법 등은 위키 (Wiki) 기반의 커뮤니티 형성을 통한 생물정보 분야의 집단 지성 창출을 목적으로 운영되는 지식 커뮤니티인 人CoDOM을 참고해주세요.

그런데 여기서 문제 발생!!!
nr에 대한 BLAST를 다 끝내고 보니 (default setting의 tabular format으로 출력) 결과 파일에서 종명을 확인할 수 없었어요.
종명을 확인 못 하면 계통수를 그려도 계통별로 잘 묶였는지 확인도 어렵고 예쁜 색을 입혀줄 수도 없어요.

MH00089;   gi|761546247|ref|YP_009122458.1|   99.06   530   5       0   1   530   1   530   0.0   1045
MH00089;   gi|1314948409|ref|YP_009444547.1|   94.70   528   28      0   3   530   1   528   0.0   995
MH00089;   gi|1079486692|ref|YP_009307015.1|   94.89   528   27      0   3   530   1   528   0.0   991

< BLAST 수행 결과 예시>

종명 확인을 위해서는 NCBI에 GI number (또는 accession)로 검색해보는 수밖에 없는 듯 하여 매우 절망스러웠지요.
실제로 하나 검색에만 최소 클릭 5~6번이 소요되고 10개 넘어가면 웹 크롤링을 해야하는 건지 고민하게 됩니다.
이때 잘 읽은 메뉴얼 하나가 사람을 살립니다.
역시 오랜 역사를 자랑하는 생물정보 산증인 BLAST는 이미 해답을 제시하고 있었던 거죠.

아래와 같이 2단계를 순서대로 수행해 주시면 BLAST 결과에서 종명, taxid (중요), kingdom (계) 정보를 바로 확인할 수 있습니다.

1단계) taxonomy DB 세팅
먼저 nr로부터 계통 정보를 가져올 수 있도록 NCBI에서 제공하는 taxonomy DB를 세팅해 주어야 합니다.
아래와 같이 최신 버전으로 다운로드 후 환경변수에 추가해 주세요.
(이참에 nr DB도 최신 버전으로 변경해주고, 하는 김에 BLAST도 최신화해주는 게 어떨까요?)
참고로 제가 테스트했던 버전은 BLAST+ 2.2.31입니다.

$ wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/taxdb.tar.gz
$ tar zxvf taxdb.tar.gz
$ export BLASTDB=[download 위치]

2단계) BLAST 실행

BLAST 수행 시 결과 파일의 포맷을 6번, 즉 tabular로 지정하고 format specifiers에 staxids(species taxon id)와 sscinames (species scientific name), sskingdoms (species super kindoms)이 포함되도록 적어줍니다.

$ blastp -query query.faa -db nr -outfmt '6 qseqid sseqid pident evalue staxids sscinames scomnames sskingdoms stitle' -num_threads 20 -out query_vs_nr.table

짜잔~ 커맨드라인 하나로 손목 수명이 일주일 연장되었습니다 (감격).
하는 김에 stitle (species name)도 추가하였더니 이제 정말 NCBI 웹사이트 들락날락할 일 없어졌어요.

MH00089;        gi|761546247|ref|YP_009122458.1|        99.06   0.0     1610689 Sarocladium implicatum  Sarocladium implicatum  Eukaryota       cytochrome oxidase subunit I (mitochondrion) [Sarocladium implicatum]
MH00089;        gi|1314948409|ref|YP_009444547.1|       94.70   0.0     29910   Tolypocladium inflatum  Tolypocladium inflatum  Eukaryota       cytochrome c oxidase subunit 1 (mitochondrion) [Tolypocladium inflatum]
MH00089;        gi|1079486692|ref|YP_009307015.1|       94.89   0.0     71617   Tolypocladium ophioglossoides   Tolypocladium ophioglossoides   Eukaryota       cytochrome oxidase subunit 1 (mitochondrion) [Tolypocladium ophioglossoides]

< BLAST 결과 예시>

한계점)
하지만 결과를 보면 박테리아인지 진핵인지와 같이 계 (kingdom)에 대한 정보만 제시하고 있어요.
근연종의 단백질에 잘 맵핑이 된 것인지 전체 계통 정보를 확인하고 싶고 진핵생물 내에서도 딱정벌레목인지 총채벌레목인지 좀 더 세분화하고 싶은데 이 상태로는 너무 부족합니다.

제가 찾은 방법은 NCBI에서 제공하는 텍스트 포맷의 관련 DB에서 taxid 또는 학명(scientific name)으로 검색하는 것입니다.
참고로 최근 1~2년 사이 NCBI taxonomy DB가 업데이트되면서부터 taxid 별 full lineage 정보를 제공하고 있으니 꼭 최신 버전을 받아주세요.
new_taxdump.tar.gz 파일을 다운로드 후 그 안에 있는 fullnamelineage.dmp 파일을 확인하면 정보를 얻을 수 있습니다.

taxid | scientific name | full lineage
1610689 |       Sarocladium implicatum  |       cellular organisms; Eukaryota; Opisthokonta; Fungi; Dikarya; Ascomycota; saccharomyceta; Pezizomycotina; leotiomyceta; sordariomyceta; Sordariomycetes; Hypocreomycetidae; Hypocreales; Hypocreales incertae sedis; Sarocladium;        |

추가 팁 하나 더!)
BLAST 결과 출력시 파일 포맷을 여러 개로 하고 싶을 때가 있어요.

XML은 기본이니깐 꼭 있어야 할 것 같고 파싱하기 귀찮으니 tabular로도 해야 할 것 같고 또 alignment된 거 확인하고 싶으니 pairwise로도 남겨두고 싶을 때 어떻게 해야하는 거죠??

BLAST를 3번 하면 돼요. 하지만 오래 걸린다는 단점이 있어요.
이때 blast_formatter를 사용하시면 됩니다 (작업시간이 1/3로 줄어드는 매직)!!!
몰랐는데 BLAST 설치 디렉토리에 이미 blastp, blastn이랑 같이 자리잡고 있더라구요.
(역시 메뉴얼은 읽으라고 있는 거였어요.)
중요한 점은 처음 BLAST할 때 반드시 asn 포맷으로 출력해 주어야 한다는 것입니다.

$ blastn -db [nt] -query [query] -outfmt 11 -out [output].asn
$ blast_formatter -archive [output].asn -outfmt 5 -out [output].asn.xml
$ blast_formatter -archive [output].asn -outfmt 6 -out [output].asn.tabular
$ blast_formatter -archive [output].asn -outfmt 0 -out [output].asn.pairwise

출력 포맷은 아래를 참고하셔서 원하는 번호를 기재해 주시면 됩니다.

 0 = pairwise,
 1 = query-anchored showing identities,
 2 = query-anchored no identities,
 3 = flat query-anchored, show identities,
 4 = flat query-anchored, no identities,
 5 = XML Blast output,
 6 = tabular,
 7 = tabular with comment lines,
 8 = Text ASN.1,
 9 = Binary ASN.1,
10 = Comma-separated values,
11 = BLAST archive format (ASN.1),
12 = JSON Seqalign output,
13 = JSON Blast output,
14 = XML2 Blast output

<출력 가능한 BLAST 결과 포맷 목록>

이상 새롭지 않을 수 있을, 그리고 정말 저만 알았다면 너무 슬프고 민망할 것 같은 BLAST 꿀팁 소개를 마칩니다.
이 방법 외에도 오조오억 개의 다른 길이 있을 겁니다. 아시는 분은 제보 바랍니다.
당신은 우리와 함께 가시면 안 될까요??? (=스카우트하고 싶어요).

<출처 : 'Mnet 쇼미더머니8 화면캡처' 후 편집>

작성 : RDC 정명희 선임 연구원

Posted by 人Co

나만 몰랐던 BLAST 꿀팁

Categories

Recent Posts

Archives

Tag Cloud

Site Stats