[하루에 한 R] 주석 정보를 포함하는 서열 파일 읽어들이기

국내에서 분리된 Klebsiella pneumoniae의 유전체 해독 결과물을 분석하면서 항생제 내성 플라스미드를 서로 비교하는 그림을 그리게 되었다. 예전에 클로렐라의 엽록체 게놈 서열을 비교하는 그림을 그리면서 사용했던 genoPlotR을 사용한 적이 있어서 기억을 더듬어 가면서 그림 한 장을 그리는데 꼬박 이틀이 소요되었다.

genoplotR은 GenBank flat file을 읽어들이는 것에서 시작한다.

pJ = read_dna_seg_from_genbank("CP022926.gbk")

이렇게 자료를 읽어들여도 충분한데, 매뉴얼에서는 에러 처리 등을 위하여 우아하게 try()라는 wrapper를 쓰게 하였다. 즉, pJ = try(read_dna_seq_from_genbank()) 이렇게 입력하라는 것이다.

자료를 읽어들인 다음 아무런 장식을 하지 않고 CDS를 표시하였다. 그랬더니 플라스미드의 시작과 끝부분까지 하나로 이어지는 고약한 feature가 보였다.

어이쿠, 순악질 여사의 일자 눈썹처럼 저게 뭔가.

플라스미드는 원형이므로 이를 GenBank 레코드로 표시하려면 부득이하게 어느 한 부분을 끊어서 +1 nt를 정의해야 한다. 그러는 과정에서 경계 부분에 걸치는 유전자가 발견된 것이다. 서열 제출 전에 아주 주의 깊게 annotation을 했다면 아무런 feature가 없는 곳을 +1 nt로 정의했을 것이다. 그러나 요즘은 대부분 NCBI에 염기서열을 제출하면서 annotation까지 일임해 버리니 이러한 일이 생긴다. NCBI가 어떤 위치에서 유전자를 검출하여 이름을 붙이든 그것은 자유다. 다만 제출자가 보낸 염기서열의 시작 위치를 조정하는 일은 하지 않는다.

플라스미드의 복제에 관여하는 rep 유전자를 맨 앞에 오도록 조정하여 서열을 제출했었는데 도대체 어떤 유전자가 시작 염기 위치에 딱 걸렸을까? 솔직하게 고백하자면 나는 플라스미드 자체의 생물학에는 그렇게 정통하지 못하다.

GenBank 파일을 텍스트 편집기(나는 EditPlus를 좋아한다)로 열어보면 되지만, 기왕이면 R 환경을 나가지 않고 그 안에서 파일을 열어보고 싶어졌다. 적당한 함수가 있을 것 같은데... 검색을 해 보니 file.show()라는 것이 있다. 인수로 주어지는 파일 명을 일일이 다 쓰지 않아도 탭을 이용한 자동 완성 기능도 있다! file.show()에서는 매우 단순하게 파일 내용을 보여주시만 라인 번호 표시나 검색 등의 고급 기능을 기대하지는 말자.

GenBank 파일을 열어서 서열 시작 부분을 찾아보았다. 그러면 그렇지! 얄궂은 join(..) 정보가 보인다.

tap이란 무엇에 쓰는 유전자인가? RepA leader peptide Tap이라는 product 정보가 보인다. repA의 위치는 별로 아름답지 못하게 1..858로 잡혀 있다. 이것은 사실 circlator(a tool to circularize genome assemblies 링크)의 소행이다. Bacterial chromosome을 circlator로 처리하면 dnaA 유전자를 첫번째로 오게 하면서 바로 앞에 위치한 intergenic region 안에서 +1 nt가 오도록 적절히 잘 잘라주는데, 플라스미드에서는 repA 코딩 영역의 첫번째 염기를 +1 nt로 싹둑 자르는 만행을 저지른다.

Tap은 translational activator peptide의 약자이다. copA와 repA 유전자 사이의 영역에 코딩이 되는데, 플라스미드의 복제를 조절하는데 중요한 역할을 하는 모양이다. 1992년 EMBO J에 "Replication control of plasmid R1: RepA synthesis is regulated by CopA RNA through inhibition of leader peptide translation PubMed"라는 제목의 논문이 있다. 여기에서는 Tap 펩타이드의 길이가 24 아미노산이라고 하였는데, 위에서 보인 GenBank 자료에서는 이보다 훨씬 길다.

Prokaryote에서는 기본적으로 splicing을 하지 않으므로, join으로 여러 feature가 연결되는 현상을 만나기는 쉽지 않다. +1 nt의 위치 선정 잘못으로 인해서 유전자가 두 동강이 나는 경우 말고는 말이다. genoplotR에서는 이를 어떻게 인식하고 있는지 궁금하여 View(pJ)를 해 보았다. 역시 우려했던 대로 tap 유전자를 여러 조각으로 나누어서 표현하였다.

나는 tap 유전자의 위치 정보를 다음과 같이 조작하는 단순 회피 전략을 구사하였다. 어차피 플라스미드 상호 비교에서 보여줄 영역은 아니었기 때문이다.

join(257970..258210,1..8)  => join(257970..258210)

genoplotR은 GenBank 자료를 필요로 하지만 심각한 수준의 조작을 하기 위해 만들어진 것은 아니다. 따라서 이러한 목적에 맞게 만들어진 패키지를 익히는 것이 더욱 바람직할 것이다. BioPerl을 사용하여 낮은 수준에서 이 기능을 구현할 수도 있지만 요즘은 R을 쓰는 것이 대세가 아니겠는가. 검색을 해 보니 genbankr이라는 것이 있다. 비교적 최근에

Simply put, the genbankr package parses files in the NCBI's GenBank (gb/gbk) format into R.

설치를 해 볼까? 이런, 내가 쓰는 R 3.4.1에서는 깔리지 않는다. 나중에 시도할 숙제로 남겨 놓자. 다른 패키지로는 biofiles(R interface to GenBank/GenPept files)라는 것이 있다. biofiles의 기능을 제대로 쓰려면 NCBI E-utilities와 R reutils 패키지를 이해해야 되겠다.

숙제만 하나 가득히...