2021년 11월 16일 화요일

Roary 'query_pan_genome'의 고약한 출력물

Roary: the Pan Genome Pipeline은 내가 밥 먹듯이 즐겨 사용하는 프로그램이다. 최신판은 2019년 11월 6일에 릴리즈되었으니 더 이상 개선의 여지도 없는 프로그램이 아닐까 한다. MUMmer와 같이 거의 수정도 거치지 않으면서 오랜 기간 사랑받는.

Roary의 부속 유틸리티인 query_pan_genome은 두 그룹 간의 교집합, 합집합, 또는 차집합에 해당하는 gene cluster를 추출하는 다재다능한 프로그램이다. 그룹 정보는 GFF 파일로 제공하면 된다. 예를 들어 그룹 one에는 있지만 그룹 two에는 없는 유전자를 찾고 싶다면, 다음과 같이 실행하면 된다.

$ query_pan_genome -a difference --input_set_one 1.gff,2.gff --input_set_two 3.gff,4.gff,5.gff

만약 각 그룹을 구성하는 GFF 파일이 많다면 명령행에 이를 일일이 타이핑하기가 곤란할 것이다. 그런 경우에는 fofn(file of file name)을 만든 뒤 paste 명령을 써서 쉼표를 사이에 두고 이어 붙이면 된다.

$ cat one.fofn
1.gff
2.gff
3.gff
$ one=$(paste -sd, one.fofn)

두 그룹에 대하여 이렇게 변수를 만든 다음, '--input_set_one $one'과 같은 형식으로 옵션을 주면 된다. GFF 파일이 현재 디렉토리에 없다면, fofn을 만들 때 find 명령어를 잘 이용하면 된다.

'query_pan_genome -a difference' 명령을 수행하면 'set_difference_'로 시작하는 11개의 파일이 생긴다. 그룹 one에만 있는 유전자에 대한 정보는 다음의 파일에 수록된다.

  • set_difference_unique_set_one
  • set_difference_unique_set_one_reannotated
  • set_difference_unique_set_one_statistics.csv
여기서 한 가지 주의할 것이 있다.  그룹 one에만 있는('unique') 유전자란, 그룹 one의 모든 멤버에 다 있는 것만을 포함하는 것이 아니다. 예를 들어 그룹 one이 10개의 균주(A, B, C ... J)로 이루어졌다고 하자. 그러면 set_difference_unique_set_one에는 A에만 있는 것, A와 B에 있는 것 등 별의별 녀석들이 다 섞여 있다. 만약 그룹 특이적인 유전자를 찾는 것이 목적이라면, A부터 J까지 모든 멤버에 다 있는 것을 파악해야 한다. 그러려면 위에서 나열한 세 개의 결과 파일 중 set_difference_unique_set_one_statistics.csv를 참조해야 한다.

여기에서 문제가 발생한다. 바로 gene ID가 바뀐다는 것이다. set_difference_unique_set_one의 첫 번째 컬럼은 Roary 실행 결과 파일에서 보편적으로 쓰이는 unique cluster ID를 유지하지만(cluster는 accessory gene도 포함, 즉 pan genome을 전부 아우름), set_difference_unique_set_one_reannotated에서는 GFF 파일을 참조하여 일부가 gene name으로 바뀌고, 이것이 set_difference_unique_set_one_statistics.csv로 그대로 전달된다.

실제 사례를 들어 보자. 오늘 아침에 query_pan_genome 스크립트로 Strepcococcus genus에 속하는 두 종(N=301, set one은 232개, set two는 69개)의 roary 결과를 투입하여 13683개의 set one-specific gene을 얻었다. 그런데 gene ID를 sort하여 uniq 명령어를 통과시키면 13553개만 남는다. reannotated gene ID에서 중복이 발생했다는 뜻이다.

set_difference_unique_set_one과 set_difference_unique_set_one_reannotated 파일의 나머지 컬럼에는 실제 이를 구성하는 유전자들의 ID가 있으니 그룹 자체를 건드리지는 않는다. 그러나 pan_genome_reference.fa 파일에서 set one에 특이적인 유전자(대표 서열로서)를 찾으려면 여간 번거로운 일이 아니다.

query_pan_genome 스크립트를 전혀 쓰지 않고 roary의 기본 결과 파일인 gene_presence_absence.csv를 R에서 조작하여 각 그룹 특이적인 유전자를 찾을 수는 있는데 여간 성가신 일이 아니다. 왜 query_pan_genome 스크립트는 다시 annotation을 하여 식별자를 유일하지 않은 것으로 만들어 버리는지 알 수가 없다. QueryRoary.pm를 내가 고칠 수 있는 수준의 일도 아니고...


댓글 없음: