A5-miseq 중간 결과물의 재활용

A5-miseq은 일루미나 플랫폼에서 만들어진 미생물의 시퀀싱 결과물을 사용하여 de novo assembly를 실시하는 간편한 파이프라인 형태의 도구이다. PacBio RS II 데이터를 HGAP(또는 Falcon)으로 처리하여 거의 완성 수준에 가깝도록 유전체 서열을 얻는 것이 일상다반사가 된 이상 일루미나 기술도 과거와 같은 인기를 누리지 못할 것이라는 때이른 기대를 할 수도 있다. 그렇지만 아직은 생산성이 매우 높은 기술이기 때문에 대량의 유전체 샘플을 다루거나 메타지놈 시퀀싱을 할 때에는 여전히 일루미나가 유용하다.

Command line interface에서 돌아가는 de novo assembler라고 해서 서열 조립과 관련된 바이너리만 덜렁 주어지는 일은 요즘 거의 없다. 서열의 전-후처리를 위한 보조적인 스크립트 또는 바이너리가 패키지 형태로 같이 배포된다. 조금 더 발전한 형태라면 VelvetOptimiser처럼 조건을 바꾸어가면서 여러차례 assembler를 반복 실행하여 최적화된 결과를 만들어내는 wrapper script가 있다. IDBA_UD도 k-mer의 크기를 20 단위로 바꾸어가면서 가장 좋은 결과를 최종적으로 제공한다. A5-miseq은 이보다 더욱 발전한 형태이다. 즉 서열 조립뿐 아니라 전-후처리에 필요한 다른 프로그램들까지 일괄적으로 실행하는 일종의 파이프라인인 것이다.

최근 논문을 검색해 보면 이보다 더욱 발전한 형태를 볼 수 있다. 설치 및 실행이 편리한 통합 환경 형태의 프로그램이 나오고 있다. 예를 들어 metAMOS는 (meta)genome assembly를 위한 전처리와 조립, 스캐폴드 작성 및 annotation에 이르는 과정을 관장하는 modular framework이다. MyPro는 또 어떠한가? 이것은 원핵생물 유전체의 조립과 주석화를 통합적으로 실시하는 환경으로서 아예 VirtualBox에서 실행 가능한 형태로 배포가 된다. 대신 다운로드할 파일의 사이즈가 꽤 크다.

서론이 좀 길어졌다. 오늘은 A5-miseq 실행 과정에서 만들어지는 중간 결과물, 즉 sga correct가 생성하는 교정된 read를 다른 프로그램에서 쓸 수 있는지를 검토하다가 사소한 문제가 발견되어 이를 해결하고자 포스팅을 작성하게 되었다. A5-miseq은 하나의 error-corrected read file(filebase.ec.fastq.gz)을 생성하는데, 이것을 khmer 패키지의 extract-paired.py로 처리하여 interleaved file과 orphan 파일을 만들려 하였으나 다음과 같은 에러 메시지가 출력되었다.

raise Exception("no paired reads!? check file formats...")
Exception: no paired reads!? check file formats...

이게 뭔 소리지? paired read가 하나도 없다니?  error-corrected file을 열어보았다.

@HWI-ST908R:151:D1M9MACXX:1:1101:1373:2469 1:N:0:CCGTCC

서열 ID 정보가 있는 줄에 존재해야 하는 pair 관련 정보(빨강색)이 하나도 남김없이 제거된 것을 확인하였다. A5-miseq의 어느 단계에서 이 필드가 사라진 것일까? error-corrected read는 혹시 이미 interleaved 형태로 전환된 것은 아닐까? 하지만 내가 훑어본 바에 따르면 A5-miseq 스크립트 안에서는 read file을 쌍으로 다루지 않는 것처럼 보이는데? 이 궁금증을 해소해 보도록 하자.

먼저 error corrected read file이 내부적으로 쌍을 이루고 있는지 알아보았다. 서열 ID를 전부 추출하여 짝이 맞는지 점검하는 스크립트를 짤 수도 있지만, 좀 더 간단하게 일을 하기 위해서 CLC Genomics Workbench에서 이를 임포트하여 Paired status를 "Paired sequences"로 고쳐보았다. 아무런 문제없이 잘 전환된다. 만약 임의로 read 하나를 제거하면 Paired status를 고치지 못한다(Unpaired sequences로 고정).

그러나 CLC Genomics Workbench에서는 pair 상태를 철저하게 확인하지는 못한다. 시험삼아서 인접한 read pair에서 하나씩을 제거하여 짝수 상태를 유지한 채로 Paired sequence를 만들어보니 역시 아무런 오류를 출력하지 않는다. 즉 CLC Genomics Workbench에서는 파일 하나를 읽어들여서 Paired sequence로 전환할 때에는 오직 read의 수가 짝수인지만을 확인하는 것으로 생각된다. 만약 read_1.fastq와 read_2.fastq를 CLC로 임포트할 때 파일 내의 read 순서가 서로 맞지 않는다면? 이것까지 점검을 하는지는 잘 모르겠다.

그렇다면 A5-miseq의 실제 실행 내역을 추적하면서 알아보는 수밖에 없다.

[실제 명령어] a5_pipeline.pl file_1.fastq file_2.fastq a5-out

1. File 2개를 하나로 합침: shuffle_fastq() 함수를 사용한다. 결과 파일은 a5-out.s1/file_1.fastq.both.fastq이다. 
2. trimmomatic을 SE 모드로 실행: a5-out.s1/file_1.fastq.trim.fastq이 생긴다. SE 모드로 실행하였으므로 당연히 결과 파일은 pair 정보 입장에서는 들쑥날쑥이다.
3. sga preprocess --pe-mode=0 실행: outdir.s1filebase_1.fastq.both.pp가 생긴다. 바로 이 파일부터 서열 ID 라인에 pair 관련 정보가 사라졌다.
4. sga index를 위해 outdir.s1/outdir.pp.fastq가 생김
5. sga correct 실행
6. Re-pair('수선'을 의미하는 repair가 아님)를 실시함:  repair_fastq() 함수가 작동.
7. fasta 파일로 전환하여 idba로 de novo assembly 수행

범인은 바로 sga preprocess였다. 이 과정에서는 아직 남아있을 수 있는 N을 제거하는 등의 몇가지 처리를 추가한다. 그러고 나서 pair를 다시 유지시키는 함수를 호출함을 알 수 있었다. 이미 완벽한 interleaved file이 되었으므로 내부에 pair 정보 taq이 있을 필요가 없다. 물론 나중에 이를 다시 두 개의 파일로 분리하려면 프로그램에 따라서는 불편을 토로할 수도 있겠지만 말이다.

trimmomatic을 SE 모드로 실행한다는 것이 오해를 불러일으킨 것이었다. 어차피 후속 과정인 sga preprocess에서 또 read의 pair 상태가 깨질 수 있으니 error correction까지 다 마친 다음에 interleaved paired read를 재생성하는 것이 훨씬 현명하다. idba를 가지고 테스트를 해 본 결과 interleaved read가 아니면 세그멘테이션 오류와 비정상 종료를 한다는 것을 확인하였다.