Contig circularity를 점검하는 스크립트, check_circularity.pl

박테리아 유전체를 long read 기술로 시퀀싱하여 최신 assembler로 조립하면 원형 구조의 replicon 구조를 쉽게 얻는다. Contig가 원형으로 닫힌 구조인지 어떻게 확인하는가? Contig 앞뒤에 나오는 서열이 겹치면 그렇게 판단할 수 있다. 자잘한 contig를 정리하고, 부족한 부분은 채워서 circular contig로 만들고, 마지막으로 dnaA 유전자를 찾아서 염색체 시작 위치를 조정해 주는('fixstart') 종합 도구로는 circlator만한 것이 없다. 아쉽지만 이 프로그램이 앞으로 업데이트될 가능성은 별로 없다고 보는 것이 옳다. roary도 마찬가지일 것이다. 두 프로그램은 모두 2015년에 영국의 Sanger Institute에서 개발하였다.

요즘 나오는 long read용 조립 프로그램(flye나 unicycler가 대중적임)은 대부분 자체적으로 circularization을 해 준다. 단, canu는 예외이다.

Circlator와 같은 종합 도구를 실행하지는 않고 단지 주어진 염기서열의 시작과 끝에 같은 방향으로 반복하여 나타나는 서열이 있는지를 점검하고 싶다고 가정하자. gepard를 써서 시각적으로 확인할 수도 있고, contig 앞부분을 조금 잘라서 전체 염기서열에 대해 blastn을 해도 되지만 너무 귀찮다. 좀 편하게 확인할 길이 없을까?

Circlator 논문을 읽다가 sprai assembler 패키지의 check_circularity.pl이라는 스크립트를 개선하여 additional file 3에 추가하였다는 내용을 보았다. 만약 작은 플라스미드처럼 길이가 read length보다 짧은 경우 tandem repeat 형태의 contig가 나올 수 있는데, 오리지널 check_circularity.pl 스크립트는 이를 1 반복 단위로 줄이지는 못한다. Circlator 패키지에는 당연히 이 기능이 들어 있으며, 친절한 개발자들은 이미 알려진 다른 프로그램의 부속 스크립트도 손을 보아서 논문에 추가해 놓았다. 논문의 'Optimizing existing method' 섹션 일부를 인용해 본다.

출처: Genome Biology(2015)

SPRAI(Single Pass Read Accuracy Improver)는 long reads 시퀀싱 기법이 나온지 얼마 되지 않던 시절 그 이름을 들어 본 일이 있다. Sprai는 지금은 유지되지 않는 웹사이트에만 프로그램이 공개되었고 그 소프트웨어 자체를 발표한 논문은 나온 일이 없다. 나는 sprai가 제공하는 check_circularity.pl 스크립트의 단순함이 마음에 들어서 이를 구해 보기로 하였다. 검색을 거듭한 결과 anaconda package 웹사이트에서 *.bz2 파일의 링크를 겨우 구할 수 있었다. 파이썬 2.7을 쓰는 오래 된 패키지이니 conda environment를 만들어서 설치할 필요까지는 없고, 압축 파일을 받아서 푼 뒤 내가 원하는 스크립트만을 찾아서 shebang line을 수정하였다.

자, 그러면 시작과 끝 부분에 중복된 서열이 나오는 contig를 만들어서 시험해 보아야 한다. 앞서도 이야기했지만 요즘 인기 있는 assembler는 circularization을 다 해서 제공을 하니 check_circularity.pl을 이용할 여지가 없다. 그러나 canu는 다르다. Canu 공식 웹문서에서 제공하는 대장균의 25x PacBio read를 가져다가 canu v2.2로 조립을 하였다.

>tig00000001 len=4665109 reads=10220 class=contig suggestRepeat=no suggestBubble=no suggestCircular=yes trim=20235-4660825

참으로 친절한 설명이다. 이제 check_circularity.pl을 실행해 본다.

$ check_circularity.pl ecoli.contigs.fasta tempdir
tig00000001	circular	[1-24522] matches [4640588-4665109] with 99.833% identity.

서열의 match 여부는 blastn으로 점검한다. 그런데 canu contig의 서열 ID를 열어보다가 예전에는 없었던 'trim=20235-4660825'이라는 정보가 눈에 확 들어온다. 여기에서 제시하는 coordinate만 취하면 서열 앞과 뒤에 겹쳐서 나타나는 영역을 제거할 수 있음이 아니겠는가? 허 참, 세상은 넓고 새로운 것은 없구나. 새로워 보이는 것은 이미 존재하는 것을 뒤늦게 발견한 것에 불과할 뿐...