QIIME 1을 공부하면서 느낀 점

QIIME 2가 2018년에 공개되면서 QIIME 1에 대한 공식 지원은 중단되었다. 16S rRNA 염기서열에 의한 미생물 군집 분석을 처음 공부하는 입장이라면 당연히 QIIME 2를 가지고 시작하는 것이 바람직할 것이다. 그러나 QIIME 1에 관해서는 워낙 많은 공부자료가 인터넷에 널려있고, SRA에서 raw data를 받아서 실제로 따라서 해 보니 의외로 16S rRNA 염기서열에 기반한 microbial(요즘 유행되는 용어로는 microbiome) taxonomic profiling의 분석 과정에 대해서 많은 것을 알게 되었다. 뒤에서도 언급하겠지만 QIIME 2는 편의성과 일관된 인터페이스(물론 명령행 기반이지만) 안에 기본 원리가 살짝 숨겨져 있다는 느낌이 든다. 또한 부속 스크립트(전체 목록)가 매우 풍부하여 다른 용도로도 사용하기에도 좋다. conda environment가 아니라 CentOS 7에 원래 깔려있는 파이썬 환경에 설치하기를 잘했다고 느끼고 있다. QIIME 1을 사용할 때에는 FastTree와 USEARCH를 최신 버전이 아닌 것으로 잠시 되돌리는 트릭을 써야 하는 불편함이 있지만 말이다.

며칠 동안 QIIME 1을 다루면서 어렵다고 느낀 점은 다음의 두 가지였다. QIIME 2에서는 전부 편리하게 개선이 되었을지도 모르겠다.

1. 입력 파일에 대한 철저한 이해가 필요하다.
2. OTU picking과 chimera 제거의 과정에서 선택의 폭이 넓다. 특히 chimera 제거를 언제 하는 것이 좋은지에 대해서는 좀 더 고민이 필요하다.

지금까지 나의 시퀀싱 경험은 전적으로 미생물 단일 유전체에 대한 것이다. 커다란 유리판에 손으로 부어서 만드는 시퀀싱용 polyacylamide gel에서부터 ABI 3700을 거쳐 Roche/454 Illumina, Ion Torrent, PacBio, 그리고 약간의 MinION까지... 1991년부터 시퀀싱을 해 왔지만 많아봐야 일루미나 장비 기준으로 하나의 레인에서 생성된 최대 12개 정도 샘플에 대한 paired read - 이미 demultiplexing이 되어 샘플 수 x 2에 해당하는 - fastq 파일을 가지고 분석을 시작하는 것이 전부였다. 전형적인 파일 이름의 형태는 L2S357_15_L001_R1_001.fastq.gz와 같다. QIIME 2에서 쓰는 용어로 말하자면 Casava 1.8 paired-end demultiplexed fastq에 해당한다(QIIME 2 documentation - Importing data 링크). 동시에 처리했던 데이터 중 최대 규모였던 것은 150개 정도의 박테리아 균주 유전체에 대한 일루미나 결과였다. 인간 유전체를 다루는 사람에게는 명함도 내밀지 못할 수준이 아닌가.

QIIME '2'의 튜토리얼을 따라가며 공부를 하다가 제일 먼저 놀랐던 것은 바코드를 읽은 fastq 파일이 따로 존재한다는 것이었다. 그러면 SRA에 raw data를 제출할 때에는 어떻게 하지? 실제로 천랩에서 sequencing only service를 요청하면 어떤 결과를 제공할까? 다음주에 곽 박사에게 물어봐야 되겠다. QIIME 1 공식 문서의 Preparing raw Illumina data in different formats for use with QIIME에도 포맷에 따른 일루미나 파일 처리에 대한 유용한 정보가 많다. Single-end read와 barcode read로 구성된 한 쌍의 파일은 EMP(Earth Microbiome Project) protocol의 multiplexed single-end fastq라 부른다. 이것을 처리하는 QIIME 1 script는 split_libraries_fastq.py이다. Paired-end fastq는 join_paired_ends.py로 처리한다. 341F와 805R 프라이머로 16S rRNA의 V3-V4 region을 증폭한 산물의 크기는 중간값이 421 bp이다(링크). Amplicon을 만들어서 MiSeq으로 2 x 250 cycle pair end sequencing을 하면 중간을 연결하기에 매우 적당하다. Join이 얼마나 성공적이었는지 확인하고, 부속되는 바코드 read 파일을 업데이트하는 일이 필요하다. Input read 전부가 성공적으로 연결되지는 않으므로, 후속 분석 작업을 하려면 원본 바코드 read로부터 join된 read의 것만을 남겨야 하기 때문이다. 바코드 입력 데이터 파일 처리 방법은 QIIME 2가 좀 더 일관적이며 합리적인 것으로 느껴진다. QIIME 1에서는 다음과 같이 파일 타입(일루미나)에 따라 다른 방법을 취해야 한다. 

1. EMP multiplexed SE fastq: split_libraries_fastq.py(quality filtering 가능: -q 19)
2. EMP multiplexed PE fastq: join_paired_ends.py + split_libraries_fastq.py
3. Casava 1.8 PE demultiplexed fasqt: join_paired_ends.py + add_qiime_labels.py 이것은 이해하기가 조금 까다롭다. Edamame course의 교육 자료를 보면 도움이 될 것이다. EDAMAME(Exploration in Data Analysis for Metagenomic Advances in Microbial Ecology)에서 개최하는 교육 코스이다.

2018년 Science에 실린 Gopalakrishnan 등의 논문 'Gut microbiome modulates response to anti–PD-1 immunotherapy in melanoma patients'에서 ENA에 등록한 FMT 실험 마우스의 16S sequencing data(study accession PRJEB22895)를 다운로드하여 보았다.

서열 ID(Wargo.R2d7_45029)를 보면 이미 demultiplexing이 된 상태임을 알 수 있다. 아마도 split_libraries_fastq.py를 --store_qual_scores 인수와 함께 실행하여 FASTA 및 .qual 파일을 함께 만든 다음, convert_fastaqual_fastq.py을 실행하여 fastq로 전환한 것 같다. 외견상으로는 single end fastq이지만, 논문에서는 paired end sequencing 결과를 join하였다고 밝히고 있다. 시퀀싱 결과물을 제출하기 편리하게 만들기 위해 전환을 한 것임에 틀림없다. 저러한 형태의 파일은 절대로 일루미나 장비에서 생산된 직후의 모습이 아니다. Wargo는 이 논문의 교신저자 이름이니 말이다.

Chimera 제거는 반드시 필요하다고 생각한다. 천랩의 자료에 의하면 16S rRNA의 454 sequencing 데이터에서 약 15-20%가 chimera였다는 논문을 인용하였다. 그러나 Illumina Overview Tutorial에서는 이 과정이 생략되어 있다는 것이 의외였다. Chimera 제거는 OTU picking 전에 모든 read(seqs.fna)에 대해서 할 수도 있고, 그 후에 할 수도 있다. 공식적으로는 identify_chimeric_seqs.py를 이용하여 제거하지만(ChimericSlayer 또는 USEARCH 6.1 사용; 설명 링크) QC 과정의 일종인 usearch quality filter(pick_otus.py -m usearch로 작동; 설명 링크)를 사용할 수도 있다. Reference sequence가 필요한 방법도 있고 de novo로 검출하는 방법도 있다. OTU picking을 마친 뒤 각 클러스터에서 뽑은 대표 서열을 가지고 chimera를 검출한 경우 phylogenetic tree와 alignment 등에서 이를 제거해야 한다. 개인적인 견해로는 UCHIME을 사용하는 usearch quality filter로 모든 read를 먼저 처리한 다음 이것을 가지고 OTU picking, 대표서열 추출, taxonomic assignment 및 tree building을 순차적으로 실행하는 것이 바람직한 것 같다. pick_open_reference_otus.py를 실행하면 단 하나의 스크립트로 이 모든 것을 다 진행할 수 있다.

OTU picking의 여러 방법, biom 파일의 구조 등 공부할 것이 너무나 많지만 워낙 양이 많아서 여기에서 다 기술하기는 어렵다. QIIME 2에서는 개념이 너무나 많이 바뀌었다. QIIME 1에서는 454 데이터에 대해서만 쓰이는 denoise라는 용어가 QIIME 2로 와서는 아주 일반적으로 쓰이는 것도 같고, OTU cluster는 feature라는 말로 바뀐 것 같으며, 모든 데이터가 .qza라는 'artifact' 파일로 다루어진다는 것이 외견상 깔끔해 보이기는 해도 마치 두꺼운 포장재로 둘러싸인 상자 속 물건을 다루는 것 같아서 아직은 좀 어색하다.

Denoising이란 오류를 교정하는 일이라고 한다. 그런데 genome assembly를 위한 일루미나 데이터의 전처리에서는 error correction이라고는 불러도 denoising이란 말은 잘 쓰지 않는다. 근본적으로 동인할 일을 서로 다른 분야에서 약간 다르게 부르는 것인지, 혹은 16S rRNA amplicon에서 오류 정정을 하는 방법이 독특하기 때문에 별도의 용어로 칭하는 것인지는 잘 모르겠다. 가령 SPAdes assembler의 error correction 모듈을 16S amplicon data에 적용해도 될까?

Alpha 및 beta diversity에 대한 간단한 설명은 metagenomics wiki에 수록되어 있다. OTU table을 만든 다음에 해야 할 일에 해당한다.

일주일 정도 더 QIIME 1을 주물러 봐야 일정 수준의 이해도에 오를 것 같다. 그 다음에는 설치만 해 놓고 방치했던 QIIME 2를 다시 시작해 봐야 한다.

QIIME 2의 특징(논문에서 인용)

• Interactive spatial and temporal analysis and visualization tools
• Support for metabolomics and shotgun metagenomics analysis
• Automated data provenance tracking to ensure reproducible, transparent microbiome data science

참고 자료

• QIIME 1 공식 웹사이트
• QIIME 2 공식 웹사이트
• [GitBook] An introduction to QIIME 1.9.1 [Read] 버튼을 클릭할 것!
• Baochen Shi(UCLA)의 Metagenomics Analysis 워크샵 자료(2017)
• EDAMAME course: Amplicon analysis 자세히 들여다보면 유용한 정보가 많다. 예를 들어 여러 pair의 fastq 파일을 join할 경우 Merge_Reads_Script.sh를 사용하면 한번에 작업을 할 수 있다. 물론 파일 이름(base)의 목록이 필요하다. 5년 동안 진행된 EDAMAME workshop의 성과에 대한 논문('Strategies for Building Computer Skills To Support Microbiome Analysis: a Five-Year Perspective from the EDAMAME Workshop', mSystems 2019)이 나올 정도이니 꽤 유명한 교육 강좌인 것 같다. merge_paired_fastq.sh는 내가 Merge_Reads_Script.sh를 고친 것이다. list.txt 파일을 만들지 않아도 되고, 병합한 fastq 파일을 남겨서 다른 용도로 쓸 수 있게 하였다. 병합한 fastq file은 vsearch 파이프라인에서 쓰려고 한다.
• EZBioCloud(천랩) tutorial - Discovery biomarkers facilliating anti-PD-L1 efficacy
• [The mothur blog] mothur and QIIME