Whole-metagenome shotgun(WGS) 시퀀싱 자료의 성공적인 분석을 위하여 host의 유전체 서열을 제거하는 일은 매우 중요하다. 그런데 어떤 방법으로? 즐겨 사용하는 CLC Genomics Workbench에서 인간 유전체 GRCh38.p7을 레퍼런스로 하여 매핑을 한 뒤 unmapped read를 모아서 후속 작업을 할 수도 있겠지만, metagenomic assembly는 외부 프로그램을 써야 하니 이를 export해야 한다. 파일의 수가 많으면 상당히 번거로운 일이 아닐 수 없다. 차라리 human genome sequence screen 등의 전처리를 아예 명령행 환경에서 하는 것이 낫다.
CLC Genomics Workbench의 Microbial Genomics Module을 써서 'de novo assemble metagenome' 기능을 실행할 수도 있지만, 내가 갖고 있는 버전은 너무 오래 되었다.
요즘 공부하고 있는 논문에서는 그저 Bowtie2를 이용하여 sequencing read를 인간 유전체 염기서열에 붙여서 제거하였다고만 간단히 기술하였다. SAM/BAM 파일을 만들어서 unmapped read를 추출하는 작업은 명령어 한 두 줄로 해결되지는 않는다. 좀 간단하게 할 방법은 없을까? 검색을 해 보니 이런 논문이 눈에 뜨인다.
Evaluation of methods for detecting human reads in microbial sequencing datasets. Microbial Genomics 2020; 6 DOI 10.1099/mgen.0.000393
매핑 후 Picard Tools를 써서 BAM 파일을 정제해야 하고, SAMtools view에 복잡한 파라미터를 넣어서 두 단계의 작업을 거치고... 할 일이 너무나 많다. 최적의 검출 방법을 결정하기 위한 연구 논문이니 이렇게 정성을 기울이는 것도 이해는 간다. 보다 심플한 방법이 없을까? 한참을 알아보다가 FastQC로 잘 알려진 Babraham Bioinformatics의 FastQ-Screen을 쓰기로 하였다. 일루미나 PE read는 BBduk과 BBmerge를 거쳐서 어댑터를 제거하고 fwd/rev read를 병합하여 사용하기로 했다.
Babraham Bioinformatics에서 개발하는 프로그램은 버전 번호를 올리는데 상당히 인색하다. FastQC는 2010년 v0.1에서 시작하여 지금은 v0.11.9이고, FastQC는 2011년 v0.1로 시작하여 현재는 0.14.0이다. NGS의 종말에 임박하면 비로소 v1이 되려는 것인지?
FastQ-Screen은 Perl로 짜여진 프로그램이다. Reference genome을 내려받는 옵션이 있어서 실행해 보았지만 도무지 다운로드가 될 기미가 보이지 않아서 일루미나 iGenomes 웹사이트에서 인간 유전체 염기서열 자료만을 다운로드하였다. Configuration file에서 설정 가능한 genome database는 human, mouse, E. coli, PhiX, adapters, vectors의 여섯 가지. 여섯개의 DB를 다 적재한 상태에서 첫 번째 것에만 unique하게 매핑하고 나머지에는 하나도 붙지 않는 read만를 뽑아내려면 다음과 같이 명령어를 날리면 된다.
fastq_screen --tag --filter 100000 reads.fq.gz피
--filter 뒤에 오는 연속적인 숫자는 다음의 규칙에 의해서 만들면 된다.
- 0: Read does not map
- 1: Read maps uniquely
- 2: Read multi-maps
- 3: Read maps (one or more times)
- 4: Passes filter 0 or filter 1
- 5: Passes filter 0 or filter 2
- -: Do not apply filter to this genome
'fastq_screen --tag --filter' 명령과 달리 'fastq_screen'으로 단순히 탐색적인 분석을 실시하면 10만 개의 read만 샘플링하여 사용한다. 매우 합리적인 선택이다. 하긴, CLC Genomics Workbench에서 read mapping 후 paired reads distance distribution을 산출할 때에도 모든 mapped read를 다 쓰지는 않는다.
FastQ-Screen 전단계에서 어댑터 서열을 제거하므로, 본 분석 단계에서는 human genome을 걸러내는 일만 하면 충분할 것이다. 따라서 다른 genome DB를 다운로드하여 설치할 필요는 없다고 본다. 단, PE read를 투입한 경우 어느 한 read라도 human genome에 매핑이 되면 나머지 것도 제거하는 후처리를 거쳐야 하고, 최종적으로 read pair를 수선해야 된다. 이러한 post-process용 스크립트를 짜는 것이 숙제로 남았다. BBMap/repair.sh를 쓰면 될 것이다. 프로그램 개발자인 Brian Bushnell에게 감사를! 더불어 Babraham Bioinformatics의 생명정보학 트레이닝 과정 웹사이트에 빼곡하게 소개된 교육 자료에도 감동을 받았다.
Brian Bushnell의 최근 행적을 찾아보려다가 비교적 최근의 논문을 발견하였다. 어머, 이건 꼭 봐야 해!
DOE JGI Metagenome Workflow. mSystems 2021; DOI 10.1128/mSystems.00804-20
