2021년 5월 16일 일요일

[QIIME 2] Nanopore로 생산한 full-length 16S rRNA gene sequencing data 다루기 - 중간 처리가 끝난 데이터에서 각 샘플별 서열을 꺼내는 방법

인터넷에서 구할 수 있는 QIIME 1과 2의 학습용 자료는 매우 많다. 이를 이용하여 2년이 넘게 여러 차례에 걸쳐서 연습을 해 보았으나 실제 내 손에 주어진 데이터가 없으면 근본적인 이해도 잘 되지 않고 진도도 잘 나가지 않는다.

QIIME은 버전 1과 2의 사용 방법과 철학에 현격한 차이가 있어서 두 가지 전부에 익숙해지기가 쉽지 않다. 만약 QIIME을 처음으로 배우는 사람이라면 당연히 버전 2를 권할 것이다. 그런데 최근 새로 주어진 미션에서는 오래전에 454로 생산한 16S microbiome raw data를 분석해야 한다. 아마도 이것이 QIIME 1을 써야 할 마지막 이유가 될 것이다.

나노포어로 시퀀싱하여 demultiplexing까지 끝난 24개의 샘플을 들고 약 일주일 넘게 씨름을 하는 중이다. 미생물 profiling 연구를 하는데 nanpore sequencing 기술이 정말 유용할 것인가? 지금까지 표준 방법으로 널리 쓰이는 일루미나 시퀀싱 자료는 미생물 16S rRNA gene의 V3-V4 region이라는 제한된 영역을 다루지만, nanopore는 1 kb가 훨씬 넘는 16S rRNA 서열 거의 전체(심지어 16S-23S-5S 전체!)를 다룰 수 있으니 species 수준의 해상도는 훨씬 높을 것이 자명하다. 2020년 공개된 논문인 A preliminary study on the potential of Nanopore MinION and Illumina MiSeq 16S rRNA gene sequencing to characterize building-dust microbiomes(full text)의 초록은 이러한 문장으로 결론을 내리고 있다.

At the level of family and above, there was no significant difference between the microbial compositions as revealed by the two platforms. However, at the genus, and particularly at the species level, the ONT MinION reported greater taxonomic resolution than Illumina MiSeq.

EPI2ME에 전체 데이터를 업로드하고 벌써 열흘 정도가 지났는데 아직도 분석이 진행 중이라 한다. 무료 서비스의 한계인지도 모르겠다. 다른 분석 도구로는 무엇을 쓰면 좋을까? 2020년 Nucleic Acids Research에 공개된 NanoCLUST(a species-level analysis of 16S rRNA nanopore sequencing data)라는 도구를 쓰니 별로 어려울 것이 없이 각 샘플에 대하여 10~20개의 species 할당 정보를 손에 쥘 수 가 있었다. NanoCLUST에서 사용하는 클러스터링 방법은 매우 독특하다. 서열 identity를 기준으로 하는 것이 아니라 Uniform Manifold Approximation and Projection(UMAP)와 Hierarchical Density-Based Spatial Clustering of Applications with Noise(HDBSCAN)이라는 방법을 사용하여 'unsupervised read clustering'을 실시한다는 것이다. 이 두 가지의 방법은 캐나다의 Leland McInnes라는 사람이 개발한 것이다.

https://twitter.com/leland_mcinnes
일단 NanoCLUST를 사용하여 일찌감치 전체 데이터 처리를 끝내 놓은 다음, 다른 '전통적'인 방법이 없는지 알아보았다. Deni Ribicic이라는 사람이 q2ONT라는 QIIME 2용 bash pipeline을 만들어 놓았기에 이를 응용해 보기로 했다. GitHub 웹페이지에 나온 QIIME 2 설치 설명에 약간의 실수가 보이길래 이를 바로잡고자 한다.

yml 파일은 qiime2-2019-10을 위한 것인데 conda 환경명을 2019.7로 만드는 것으로 잘못 지시하였다. q2ONT.sh 스크립트 내에서는 2019.10 환경을 activate하는 것으로 짜여져 있으니 이를 고치지 않으면 에러가 발생한다.

q2ONT에서는 OTU clustering의 기준값을 85%로 잡으라고 하였다. 97%가 아닌 85%를 쓰는 이유는 물론 nanopore read의 에러 때문에 그러한 것이다. Deni Ribic이 이런 기준을 쓰게 된 것은 다음의 논문을 근거로 한다. 
Curren, Emily, et al. "Rapid profiling of tropical marine cyanobacterial communities." Regional Studies in Marine Science 25 (2019): 100485.
그래서 나도 q2ONT 스크립트에 24개의 demux read를 넣어 보았다. 이미 24개의 파일로 분리가 되어 있어서 스크립트를 약간 고쳐야 했다. 그러나 그 결과는 처참했다. 약 5천개의 chimera까지 제거된 100만개의 read는 85% OTU clustering을 거쳐도 여전히 60만개 수준이었던 것이다!  Silva 99 database(16S only)의 서열 수 전체가 369,953개라는 것을 생각하면 클러스터링이 거의 되지 않은 것이나 다를바가 없다.

같은 나노포어 데이터를  QIIME 1.9.1의 스크립트로 처리하여 de novo uclust clustering을 해 보았다. 이번에도 threshold를 85%로 잡았지만 q2ONT의 결과와 별로 다르지 않았다. 이쯤에서 중간 결론을 내려 보자. 내가 이번에 매만진 나노포어 데이터가 아주 유별난 것이 아니라면, 85% 서열 기준의 clustring-first method는 적합하지 않아 보인다.

그렇다면 서열 데이터베이스에 대한 검색을 통한 read의 분류를 먼저 이용하는 방법이 나을 수도 있다. 국내에서 나노포어 데이터를 이용하여 microbiome을 분석한 선구자라고 할 수 있는 KAIST 조병관 교수 연구실의 논문을 두 개 찾아서 읽어 보았다. 나노포어 데이터를 사전에 클러스터링하는 것이 아니라 LAST aligner를 써서 16S reference database에 대하여 검색을 한 것으로 나온다. Nanopore의 데이터 특성에 맞는 align parameter가 몇 개의 논문에 소개되어 있었지만, LAST Cookbook에 의하면 last-train을 이용하여 high-indel-error long read를 'genome'에 정렬하는 방법을 소개하였다. Reference DB가 genome이 아닌 16S rRNA 서열이라면 이러한 방법을 그대로 적용하는 것이 옳은지는 아직 알 수가 없다.

LAST website: http://last.cbrc.jp/에서 https://gitlab.com/mcfrith/last로 이전함

LAST를 트레이닝과 더불어서 실행하거나, 또는 mat-convert를 쓰는 경우 상당히 시간이 많이 걸렸다. 실제 데이터를 투입하여 분석을 한 다음 각 query의 best hit을 찾고 이것에 연결된 taxonomy를 뽑으려면 KOBIC의 도움을 받는 것이 좋을 것 같다.

편의성 측면에서는 NanoCLUST가 제일이다. 그런데 전통적인 서열 기반이 아닌 k-mer frequency vector를 이용한 클러스터링을 사람들이 얼마나 받아들이려 할까? 기준을 제시하지 않고도 데이터가 '스스로 뭉치게' 만드는 것이 요즘 데이터 과학의 대세이니 생소한 방법이라고 무시할 수는 없다.

NanoCLUST나 LAST를 쓰기 전에 나노포어 데이터의 전처리를 해 두는 것이 옳지 않을까? QIIME(2) 과정의 일부를 활용하는 것이 필요할 수 있다. 즉 quality와 길이 기준으로 read를 잘라내고, chimera까지 제거해 두는 것이다. 여기까지 진행한 QIIME table과 data로부터 각 샘플별로 FASTA file을 뽑아내면 NanoCLUST 또는 LAST를 쓰기에 매우 좋을 것이다. 물론 프로그램에 따라서는 아무런 전처리를 하지 않은 raw data를 그대로 적용하는 것을 원할 수도 있다. 이는 좀 더 테스트를 해 볼 필요가 있지만.

q2ONT를 돌리고 나면 6.1_uchime-ref-out 디렉토리에 다음의 두 결과 파일이 생긴다.

  • table-nonchimeric-wo-borderline.qza
  • rep-seqs-nonchimeric-wo-borderline.qza
이로부터 특정 샘플에 해당하는 서열을 FASTA file로 뽑아내고 싶다. 어떻게 하면 좋을까? QIIME 2 documentation 웹사이트에서는 이를 직접적으로 설명하지 않기에 Q&A를 뒤져서 겨우 방법을 알아냈다.  우선 추출할 샘플의 ID를 수록한 tsv 파일을 만든다. table-nonchimeric-wo-borderline.qza 파일을 'qiime feature-table summarize'로 처리하여 visualization file을 만든 뒤 컬럼 이름('Sample ID')와 추출할 실제 샘플의 명칭이 정확히 어떻게 되어 있는지를 확인하는 것이 바람직하다. 실수로 공백이 하나 포함되어 있다거나 한다면 모든 샘플이 전부 필터링 아웃되어 제거될 수도 있기 때문이다.

$ cat samples-to-keep.tsv 
Sample ID
0200731-0847-barcode01

다음에 입력할 명령어는 이러하다. qiime feature-table filter-samples 명령어를 사용하여 테이블에 대한 필터링을 수행하고, qiime feature-table filter-seqs 명령에 이를 제공하여 데이터 파일에서 조건을 만족하는 서열을 추출하는 것이다. 오류가 없는지(조건을 잘못 기재하였거나 타이핑 실수 등) 중간에 틈틈이 visualization file을 만들어 확인하는 것이 바람직할 것이다.

$ qiime feature-table filter-samples \
  --i-table 6.1_uchime-ref-out/table-nonchimeric-wo-borderline.qza \
  --m-metadata-file samples-to-keep.tsv \
  --o-filtered-table id-filtered-table.qza

$ qiime feature-table summarize \
  --i-table id-filtered-table.qza \
  --o-visualization temp.qzv

$ qiime tools view temp.qzv # 38435 sequences

$ qiime feature-table filter-seqs \
  --i-data 6.1_uchime-ref-out/rep-seqs-nonchimeric-wo-borderline.qza \
  --i-table id-filtered-table.qza \
  --o-filtered-data id-filtered-seqs.qza

$ qiime tools export \
  --input-path id-filtered-seqs.qza \
  --output-path export_test

$ grep -c '>' export_test/dna-sequences.fasta
38435 # OK!

 '0200731-0847-barcode01' 샘플에 존재하는 38,435개의 서열이 들어 있는 파일인 export_test/dna-sequences.fasta를 얻게 되었다. 이 과정을 24개의 샘플 전체에 대하여 반복을 하면 된다. 좀 귀찮은 일이기는 하다.

'qiime feature-table filter-samples' 명령을 실행할 때 메타데이터 파일에 정의된 샘플의 속성을 기반으로 추출하는 것도 가능하다. 다음의 예시를 보자.

$ qiime feature-table filter-samples \
   --i-table 6.1_uchime-ref-out/table-nonchimeric-wo-borderline.qza \
   --m-metadata-file metadata_all_samples.tsv \
   --p-where "[Sample ID]='0200731-0847-barcode01'" \
   --o-filtered-table subject-1-filtered-table.qza

'--p-where' 뒤에 입력할 조건은 SQLite의 where syntax를 참조하여 만들면 된다.

QIIME 2 문서를 뒤적거리면서 한숨이 나왔다. 기본 개념을 잡는 것이 왜 이렇게 어려울까? 내가 공부가 부족한 것인가? 다들 잘 이해하고 있는데 나만 모르고 있었나? 그저 RFTM(Read the fuxxing manual)을 하지 못한 내 탓으로나 여기고 Amanda Birmingham의 QIIME 2 튜토리얼을 인쇄하여 집에 들고 가련다. 사실 이것도 최신 자료는 아니다. 그리고 오늘은 일요일이다...

댓글 없음: