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
여기에서 문제가 발생한다. 바로 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를 내가 고칠 수 있는 수준의 일도 아니고...
댓글 없음:
댓글 쓰기