Metaphor.pl 스크립트의 문제점 발견

Genome 사이에서 bi-directional best hit(BBH)에 의한 ortholog candidate를 찾아주는 프로그램이 있는지 검색을 하다가 Metaphor라는 것을 발견하였다. 예전에는 inparanoid나 orthoMCL을 대충 쓰고는 했었는데, 사실 이 도구는 달랑 두 개의 유전체를 비교하기에는 썩 적합하지 않다. 양방향 blast parsing script를 직접 작성하려면 또 얼마나 귀찮은가.

Metaphor: Finding Bi-directional Best Hit homology relationships in (meta)genomic datasets. Genomics 104(6): 459-463.

프로그램 다운로드 페이지에는 README 파일이 있다고 하였지만 패키지를 받아서 압축을 풀어보니 펄 스크립트 하나만 들어있다. 입력 파일은 아미노산 서열 파일(반드시 .fasta 확장자를 가져야 함) 또는 GenBank flat file이다. GenBank file을 공급하면 자동으로 translation을 뽑아낸다. 검색 엔진으로는 legacy blast(blastall)을 사용하되 xml 형식으로 만들어진 blast 출력 파일을 파싱하여 최종적으로 core gene을 뽑아낸다.

논문을 훑어본 다음 테스트 러닝을 실행하였다.

-------------- ERROR --------------Cannot find match for the first hit accession identifier.-----------------------------------

이건 도대체 무슨 에러인가? BioPerl module을 쓰는 것도 아니니 Metaphor.pl 스크립트에서 이런 메시지를 출력하는 라인이 있을 것으로 생각하여 찾아보았다.

 642 ## Extract first hit
 643 if ($blast_output_hit =~ m/1<\/Hit_num>/){
 644      #if($blast_output_hit =~ m/(\w+).*<\/Hit_def>/){                                              
 645      if($blast_output_hit =~ m/([^\s]+)<\/Hit_def>/gi){                                                  
 646          $hit_accession = $1;
 647      }else{
 648          error("Cannot find match for the first hit accession identifier.");
 649      }
 650 }

645번 라인에서 ... 사이에 존재하는 문자열을 제대로 추출하지 못한다는 뜻이다. xml 파일을 찾아보니 다음과 같은 형식이다.

MAE_RS00010 MAE_RS00010 hypothetical protein 969:1514 reverse MW:20589

패턴 매칭을 위한 정규식이 ([^\s]+)<\/Hit_def>이므로, $hit_accession이 제대로 값을 얻어내려면 결과 파일에는 MAE_RS00010라고 기록되어야 한다. 즉 빨강색으로 표시된 부분이 존재해서는 안된다는 뜻이다. 이렇게 하려면  (1) query file의 서열 ID 부분에는 오직 >seqID라고만 되어 있어야 함을 의미한다. 아니면 (2) 코멘트 처리된 644번 라인을 살리고 645번 라인을 무력화하면 query file을 건드리지 않고서 그대로 사용할 수 있다. 실제로 (1)과 (2)의 두 가지 방법을 전부 사용하여 성공적으로 Metaphor.pl을 실행하였다.

Perl에 까막눈이 아니라서 정말 다행이다. 그러나 이보다 더 심오한 에러가 발생한다면 수정할 자신은 별로 없다.