NCBI 웹사이트 이용도 언젠가는 유료로 전환될까? 허, 생각만 해도 끔찍하다.
KEGG 웹사이트에 제공하는 여러 서비스 중에서 KEGG Mapper는 직관적인 시각화 정보를 제공하기에 매우 유용하다. 유전체 주석화가 끝난 뒤, 각 유전자에 부가된 KO 정보를 KEGG Mapper - Reconstruction Pathway에 입력하면 된다. 입력할 정보(gene list)의 샘플은 여기를 클릭하면 나온다. 유전자 ID와 KO 번호의 두 컬럼으로 이루어진 모습이다. 만약 하나의 유전자가 여러 KO 번호를 갖는다면 여러 줄로 나누어서 입력하면 되는 것으로 알고 있다.
이런 목적으로 사용할 KO annotation 정보를 어디에서 얻을까? 최소한도의 데이터베이스를 사용하는 Prokka에서는 불가능하고, 최소한 Blast2GO를 쓰거나 eggNOG mapper를 써야 한다. Blast2GO는 자체적인 KEGG Pathway Map 생성 기능이 있다. InterPro 검색을 통해서 EC 번호를 먼저 구하고, 이를 기반으로 KEGG Pathway에 매핑을 하는 것이다. 그러나 무료 버전인 Blast2GO Basic에서는 검색 작업에 너무 시간이 많이 걸린다. 따라서 eggNOG mapper를 사용하는 것이 가장 현실적이다. 로컬 서버에서 상식적인 시간에 작업이 끝나기 때문이다.
eggNOG mapper(v2) 검색을 통해 만들어지는 .annotations 파일에서 아홉번째 컬럼이 KEGG KO 정보에 해당한다. 실제 사례는 다음과 같다. 한 라인이 하나의 유전자에 해당한다.
ko:K02796,ko:K19509
ko:K02795,ko:K02796
ko:K00001,ko:K00005,ko:K08317
ko:K00096
하나의 유전자가 복수의 KO 번호를 갖는 경우, 이를 분리하여 유전자 ID와 함께 다음 줄에 인쇄되도록 하고 싶다면 어떻게 하면 좋을까? 예전에는 무조건 vim을 열고 Perl 코드를 짜는 것으로 문제를 해결했었는데, 요즘은 awk나 기타 유틸리티를 bash에서 직접 짜는 것으로 점점 게을러지고 있다. 자, 그러면 꼼수의 나래를 펼쳐보자. 정답(이라고 믿는 것)은 다음의 한 줄이다. 실수는 없었을 것으로 기대한다.
$ awk -F"\t" -vOFS="\t" '$9~/^ko:/{id = "\n"$1"\t"; sub("ko:","",$9); gsub(",ko:",id,$9); print $1, $9}' eggNOG.emapper.annotations
줄괄호 안에서 세미콜론으로 구분되는 awk 명령어 네 도막이 전부이다. $9~/^ko:/는 설명할 필요가 없으니 중괄호 안으로 들어가서 설명을 하겠다.
- id = "\n"$1"\t" - id라는 변수를 하나 만들어서 줄바꿈 + 유전자 ID + 탭으로 내용을 채운다.
- sub("ko:","",$9) - 첫번째 나타나는 KO 번호에서 ko:를 제거한다.
- gsub(",ko:",id,$9) - gsub() 함수는 모든 매치에 대해서 치환을 실시한다. 첫번째 ko:는 이미 윗단계에서 처리되었고 나머지 것들(콤마를 포함 ",ko:") 전체는 id 변수에 들어간 값으로 치환된다.
- print $1, $9 - 결과를 프린트한다. $9가 어떻게 변했는지 생각해 보자. 첫번째 KO는 그 모습을 유지하고 있지만, 나머지는 줄바꿈과 유전자 ID를 앞에 붙인 형태가 되어 새 줄에 인쇄된다.
출력물을 만들어서 원본과 대조를 해 보니 실수는 없는 것으로 판단된다. 실제 유산균 유전체 하나의 eggNOG mapper 결과를 이 awk 코드로 처리하여 2641개 유전자에 대한 KO 번호를 얻은 다음 KEGG mapper에 업로드해 보았다. Global map에서 central carbohydrate metabolite만 하이라이트하였다.
나이를 먹으면 신체 기능의 퇴화를 막기 위해서 자꾸 뭔가를 놀려야, 즉 작동시켜야 한다. 손가락을 놀리기 위하여 납땜질을 하고, 뇌를 굴리기 위해 간단한 one-liner를 계속 만들어 블로그에 소개하는 것이 나의 실천 방안이 되겠다.
댓글 없음:
댓글 쓰기