TORMES pipeline를 이용한 bacterial WGS analysis

일루미나 플랫폼에서 생산된 미생물의 whole genome sequencing(WGS) analysis를 편리하게 도와주는 파이프라인인 TORMES를 소개한다. TORMES는 일루미나 시퀀싱 자료를 메타데이터 파일과 함께 제공하면 read에 대한 QC, 조립, reference 서열에 대한 순서 결정, MLST, annotation, 항생제 내성 유전자 예측 등을 자동적으로 실시하며, -g/–genera 옵션으로 genus를 특정하면(Escherichia와 Salmonella만 가능) plasmid replicon screening과 serotyping 등 상세한 분석을 추가적으로 실시한다. 결과 보고서는 RMarkdown 코드 파일로 제공된다.

TORMES pipeline을 실행하기에 앞서서 각 파일(fastq)의 시퀀싱 방향과 샘플 정보를 수록한 tab-delimited text file인 메타데이터 파일을 준비해야 한다. 다음의 명령어(‘tormes example-metadata’)를 실행하여 samples_metadata.txt를 얻은 뒤 두 번째 줄부터 실제 샘플에 맞게 수정하도록 한다. 네 번째 컬럼부터 샘플에 대한 설명을 담게 되는데, 필요하다면 컬럼의 수를 더 늘려도 된다. 최소 4개의 컬럼이 있어야 에러 없이 제대로 분석 파이프라인이 실행된다. 5번째 컬럼부터는 필수 항목이 아니므로 없어도 프로그램 실행에는 문제가 없다. 실수로 메타데이터 파일 맨 끝에 빈 칸을 남겨두면 오류가 발생하므로 주의하기 바란다. 또한 라인의 맨 앞에 ‘#’을 삽입한다고 해서 이를 코멘트로 생각하고 무시한다고 생각하면 안 된다. 또한 Description 컬럼에 작은 따옴표를 넣으려면 “’”와 같이 큰 따옴표로 둘러싸야 한다.

$ tormes example-metadata
$ cat samples_metadata.txt
Samples	Read1	Read2	Description	Use_as_many_descrpition_columns_as_wanted
Sample1	Directory1/S1_R1.fastq.gz	Directory1/S1_R2.fastq.gz	E.coli isolated in 2018	Colistin resistant

다음은 실제 reference genome sequence가 있을 경우의 파이프라인 실행 사례이다. Read에 대한 QC 후 quality filtering을 적용하게 되는데 기본 조건은 –min_len 125이다. 따라서 101 nt cycle로 생산한 raw data를 처리하려면 –min_len 옵션의 값을 101 bp보다 작게 설정해야만 한다.

$ tormes --metadata salmonella_metadata.txt --output Salmonella_TORMES_2018 --reference S_enterica-CT02021853.fasta --threads 32 --genera Salmonella

TORMES에서 기본적으로 쓰이는 SPAdes assembler가 항상 좋은 결과를 만들어내지는 못한다. 따라서 CLC Genomics Workbench 등에서 조립한 뒤 average coverage가 미흡한 것을 제거하여 정리한 contig 서열 파일을 갖고 있다면, 파이프라인 내부 assembler를 건너뛰면서 외부 제공 assembly를 활용할 수 있는 수정된 스크립트인 tormes-hyjeong을 사용하면 된다. 사전에 준비된 contig 서열 파일(.fasta)은 genomes 디렉토리에 모은 뒤 ‘tormes-hyjeong --assembler external’ 명령으로 실행한다. 나머지 옵션은 위에서 설명한 방법을 그대로 따른다. 이 방법을 사용하면 NCBI에서 입수한 유전체 서열을 포함하여 분석하는 것도 가능하다. 단, fake read를 적절히 만들어서 다른 read와 함께 입력물로 넣어 주어야 한다.

RMarkdown으로 작성된 리포트 파일은 Rscript 명령을 실행하여 html 포맷으로 전환할 수 있다. R은 tormes-1.0 환경에는 설치하지 않았으므로 base 환경 혹은 conda 외부에서 실시한다.

$ Rscript -e 'library(rmarkdown); rmarkdown::render("tormes_report.Rmd", "html_document", encoding="UTF-8")'

분석 결과 리포트의 사례는 https://nmquijada.github.io/tormes/files/에서 참조할 수 있다. 리포트에 포함된 섹션 중에서는 필요하지 않은 것도 있을 것이다. 이를 제거한 깔끔한 리포트를 얻으려면, 결과 파일 중에서 report_files.tgz를 압축 해제하여 나오는 파일 중 tormes_repord.Rmd 파일을 적절히 편집하라. 그 다음에 같은 위치에 있는 render_report.sh 스크립트를 실행하면 새롭게 만들어진 tormes_report.html 파일을 얻게 될 것이다. 동일 종에 속하지 않은 여러 균주의 Illumina sequencing 데이터를 점검하기 위하여 TORMES를 이용하는 것도 가능하다. 이 때에는 –no_pangenome과 –no_mlst을 주는 것이 상식적으로 옳다. TOMES 1.1부터는 raw data가 아니라 이미 조립된 contig 서열을 분석 대상에 포함시킬 수 있게 되었다. Contig 서열을 사용하려면 metadata file의 두 번째 컬럼(Read1)을 ‘GENOME’으로 표기하고 세 번째 컬럼(Read2)에 Directory/FASTA_file을 넣는다. 2020년 10월 공개된 TORMES 1.2부터는 KRAKEN2와 RDP Classifier를 채용하는 등 성능이 확장되었다. TORMES 1.3에서는 설치 방법이 개선되었다.