FastQ-Screen으로 whole-metagenomic shotgun read에서 host genome sequence 없애기

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

약 50분 분량의 강의 동영상 3개를 만들어 보았다

동영상 편집 자체에 익숙하지 못하므로 되도록이면 중간에 끊지 않고 길게 녹화를 하려고 무던히 애를 썼다. 하지만 특별한 '촬영 철학'이 있는 것이 아닌 상태에서 '원 테이크'로 길게 녹화를 한다는 것이 결코 쉬운 노릇은 아니다. 강연으로 먹고 사는 사람이 아닌지라 어색한 모습이 보임은 당연하다. 

화면 캡쳐에는 OBS Studio를, 촬영 후 편집에는 OpenShot Video Editor를 사용했다. 전부 무료 프로그램이고, 아는 소프트웨어는 그것 뿐이라서... 음성 녹음에는 롤랜드 사운드캔버스 SC-D70에 싸구려 다이나믹 마이크로폰을 연결하여 사용하였다. 음성 기록에는 전혀 손을 대지 않았다. 잡음을 없앤다거나 음량을 조정한다거나 하는 등 후처리가 필요하겠지만, 당장의 실력으로는 엄두가 나지 않는다. 내 얼굴을 찍는 웹캠의 높이와 마이크로폰의 각도 및 모니터 위치 등이 전혀 최적화가 되어 있지 않아서 약 네 시간 가까운 촬영을 끝내고 나니 목덜미가 뻐근하다. 새로운 녹음용 장비를 사야 되겠다는 욕심만 늘어났다. 가장 중요한 것은 마이크 스탠드. 집에 있는 탁상용 마이크 스탠드(Audiotrak AMS-11)는 USB 콘덴서 마이크 전용으로 쓰도록 하고, 자바라 형태의 스탠드를 따로 사도록 하자. 기왕이면 마이크도 하나 더...

주름져서 접을 수 있는 물건의 총칭인 '자바라'는 '뱀의 배'를 뜻하는 일본어에서 유래했다고 한다. 오, 그렇구나!(나무위키 링크)

OBS Studio로 화면을 기록할 때의 해상도는 메인 모니터의 그것과 같은 2560x1440 픽셀이었는데, OpenShot에서 편집을 마치고(편집이라고 해 봐야 끝을 조금 잘라내고 이어 붙인 것이 전부) 기본 프로파일로 '비디오 내보내기'를 하니 1024x576 16:9 PAL이 되었다. 너무 저해상도가 아닐까? 만약 터미널 창에 타이핑하는 모습을 제대로 기록하려면 글꼴을 더 키워야 한다. Windows Terminal에서 글꼴을 키우려면 어떻게 해야 하지? 검색을 통해서 알아보아야 한다.

Windows Terminal tips and tricks

공부할 것이 너무나 많다. 몇 시간짜리 강의를 하기 위해 해당 전문 분야의 지식을 축적해 놓는 것과는 차원이 다른 문제이다. 카메라와 마이크가 나를 향해서 작동 중인데, 어떻게 해야 부드럽게 나의 말을 이어 갈 것인가? 내 모습을 찍어서 계속 문제점을 파악해 보고, 연습에 연습을 거듭하는 것 말고는 방법이 없다. 스크립트를 미리 만들어서 이를 보고 읽으면서 녹화를 하고싶은 유혹은 언제나 있다. 그러나 이는 올바른 방법이 아니다. 정 필요하다면 몇 가지 핵심 단어를 적은 메모지를 앞에 놓고 슬쩍 보면서 머릿속에서 할 말을 즉각적으로 조립해 나가는 것이 진정한 '진행자'의 능력이다. 물론 메모지 넘기겠다고 손가락에 침을 발라서는 안 될 것이다.

컴퓨터를 이용한 녹음 작업에 늘 관심이 많았었는데 이번 일과 같이 구체적인 목표가 생기면서 실행에 옮기게 된 것을 매우 즐겁게 생각한다.

강의 동영상 만들기 준비 - 아파치 스트리밍 설정

무료 프로그램인 OBS Studio를 이용하여 컴퓨터 화면(파워포인트 강의자료 및 Windows Terminal) 녹화와 내 얼굴을 향한 웹캠 촬영을 동시에 실시하는 방법을 익혔다. 모니터는 두 대이므로 메인으로 사용하는 것 전체(2560x1440)을 기록하되 출력은 1280x720으로 맞추었고, 기록 파일은 default인 mkv가 아니라 mp4 포맷으로 만들었다.

연구용으로 사용하는 리눅스 서버(CentOS 7)에서 아파치를 통해서 스트리밍을 하려니 약간의 설정이 필요하다. 이것 역시 검색을 통하여 방법을 알아내었다. h264 모듈을 설치하는 것으로부터 시작한다. 그러려면 apxs(APache eXtenSion tool)이 먼저 깔려 있어야 한다. 'yum install httpd-devel pcre-devel'을 관리자 모드로 실행한다. 그러고 나서 다음과 같이 apache_mod_h264_streaming-2.2.7.tar.gz 패키지를 가져다가 압축을 풀고 설치한다.

wget http://h264.code-shop.com/download/apache_mod_h264_streaming-2.2.7.tar.gz
tar xvzf apache_mod_h264_streaming-2.2.7.tar.gz
cd mod_h264_streaming-2.2.7 
./configure
make
make install
...
Libraries have been installed in:
   /usr/lib64/httpd/modules

아파치 설정 파일 /etc/httpd/conf/httpd.conf를 열어서 다음 두 줄을 추가한다.

LoadModule h264_streaming_module /usr/lib64/httpd/modules/mod_h264_streaming.so
AddHandler h264-streaming.extensions .mp4

마지막으로 관리자 권한으로 아파치를 재시작한다.

systemctl restart httpd.service

허름한 html 문서를 만들어서 간단한 설명과 다운로드 및 스트리밍할 파일의 링크를 삽입한 뒤, 웹브라우저에 접속하여 테스트용으로 만든 mp4 파일을 클릭하여 보았다. 오디오 재생을 포함한 스트리밍 상태가 아주 양호하다. 연구소 방화벽 내에 존재하는 컴퓨터라서 보안 같은 것에 크게 신경을 쓰지 않았다. 이 동영상을 다운로드해 가면 또 어떠랴?

이것은 동영상 링크가 아니고 스크린샷 캡쳐 이미지이므로 클릭하지 마시길. 녹음은 Camac SM-221 + 롤랜드 사운드캔버스 SC-D70으로 하였다. Camac 다이나믹 마이크는 아주 오래 전에 이마트에서 구입한 것으로 가정용 노래방에나 어울리는 물건이다. 구글에서 모델명으로 검색하면 내가 쓴 글만 나온다(링크). 

같은 교재를 이용하여 오프라인 강의를 해 본 경험으로는 약 세 시간이 걸렸다. 이것을 실수 없이 녹화하여 몇 회 분량의 강좌를 만들려니 상당히 부담이 된다. 교직에 있는 사람들이라면 비대면 수업을 진행하면서 이러한 방식에 많이 익숙해졌겠지만, 강의를 본업으로 하지 않는 나에게는 쉬운 노릇이 아니다. 혼자서 컴퓨터 화면을 보면서 녹화용으로 강의를 한다는 것이 얼마나 어색한 일인지는 당연하지 않는가?

녹화 후 편집 작업은 일절 하지 않을 예정이므로(과연 그렇게 될지는 알 수 없음...자막을 넣을 일은 없느니까!) 일단 녹화 자체에만 집중해 보자. 녹화가 잘 되면 타이틀만 달아서 내 유튜브 채널에 올려 버릴까? 

취미용으로 갖추어 놓은 웹캠과 녹음용 장비(마이크 및 오디오 인터페이스 등)을 이런 용도로 쓰게 되리라고는 미처 예상하지 못했다. 당장 녹음을 하기에 불편함이 없을 수준의 장비를 갖추어 놓았음에도 불구하고 조금 더 좋은 품질의 다이나믹 마이크와 소형 믹싱 콘솔을 갖고 싶은 것은 도대체 무슨 욕심인지... 실은 MXL Temp USB 콘덴서 마이크로폰이 있으니 더 필요한 것은 없음에도 불구하고 믹싱 콘솔은 이상하게도 갖고 싶은 욕심이 사그러들지를 않는다.

GNU parallel을 readstats.py에 적용하기, 그리고 우연히 알게된 AWS의 RONIN

내가 즐겨사용하는 readstats.py는 FASTA/FASTQ 파일을 읽어들여서 총 sequence 수와 bp를 출력하는 스크립트로서 khmer package에 포함되어 있다. 염기서열 파일을 다루는 많은 유틸리티가 있지만 관성이란 참으로 무서운 것이라서 초기에 접하여 손에 익은 프로그램은 계속 사용하게 된다. 기본적인 사용법은 다음과 같다.

$ readstats.py -o readstats.txt --csv *fastq

|| This is the script readstats.py in khmer.
|| You are running khmer version 2.1.1
|| You are also using screed version 1.0
||
|| If you use this script in a publication, please cite EACH of the following:
||
||   * MR Crusoe et al., 2015. http://dx.doi.org/10.12688/f1000research.6924.1
||
|| Please see http://khmer.readthedocs.io/en/latest/citations.html for details.

... SRR13060952.fastq 0
... found 37750 bps / 250 seqs; 151.0 average length -- SRR13060952.fastq
... SRR13060953.fastq 0
... found 37750 bps / 250 seqs; 151.0 average length -- SRR13060953.fastq
... SRR13060954.fastq 0
...중간 생략...
$ head -n 5 readstats.txt 
bp,seqs,avg_len,filename
37750,250,151.0,SRR13060952.fastq
37750,250,151.0,SRR13060953.fastq
37750,250,151.0,SRR13060954.fastq
37750,250,151.0,SRR13060955.fastq

이틀 전에 소개한 HRGM(Human Reference Gut Microbiome, 블로그 링크)을 공부하기 위하여 NCBI Sequence Read Archive study SRP292575로 등록된 raw sequencing read를 다운로드하였다. 논문에서는 90개의 분변 시료라 하였는데 SRR accession으로는 106개였다.  

SRA에서 파일을 다운로드할 때에도 parallel을 이용할 수 있다. 106개의 모든 SRR accession에 대하여 fasta-dump 명령어를 만들어서 다음과 같은 run.sh이라는 텍스트 파일을 만든다.

fastq-dump --split-files SRR13060942
fastq-dump --split-files SRR13060943
fastq-dump --split-files SRR13060944

이를 parallel에 공급하면 된다. Parallel을 사용하는 아주 기본적인 기법에 들어간다.

$ parallel --jobs 16 < run.sh &
[1] 17592

모든 파일을 압축하고 싶다면 'parallel gzip ::: *fastq'를 실행하면 된다. 

FASTQ 파일을 다운로드하였으니 각 파일이 몇 bp나 되는지를 헤아려 볼 순간이 되었다. 여러 파일의 다운로드 또는 압축 작업은 parallel을 쓰기에 아주 좋다. 왜냐하면 각 파일에 대해서 이루어진 작업에 대하여 출력되는 메시지 따위에 신경을 쓸 필요가 전혀 없기 때문이다. 그러나 readstats.py는 약간 다르다. 인수로 주어지는 FASTA/FASTQ 파일에 대한 집계 수치를 표준 출력/표준에러/'-o filename'으로 지정된 파일에 기록하게 되므로, parallel로 작업을 하려면 세심하게 명령어를 날려야 한다. 예를 들자면 Parallel tutorial 문서를 참조하여 다음과 같이 실행할 수 있다.

$ parallel --results outdir readstats.py ::: *fastq

그러면 새로 만들어진 outdir 디렉토리 아래에 각 fastq 파일에 대한 별도의 디렉토리가 생기고, 그 아래에 seq stderr stdout이라는 파일이 생긴다. 솔직히 말하자면 이것을 정리하는 일이 더욱 귀찮다. 따라서 입력 파일 하나에 대해서 '-o filename'으로 출력 파일을 하나씩 만든 후(이때 parallel을 사용), 가장 마지막 단계에서 모든 출력파일을 집계하여 최종 결과를 만드는 것이 더욱 바람직하다. Parallel 명령어는 다음과 같이 만들어서 실행하라.

$ parallel -j 16 'readstats.py {} > {.}.stat' ::: *fastq

{}는 input line을, {.}는 input line에서 파일 확장자를 제거했음을 의미한다. 매뉴얼을 잘 읽어보면 의미를 알 수 있지만 사용법을 빨리 익히려면 실제 사례를 보는 것이 낫다. 오늘 익힌 이 사례는 Parallelising jobs with GNU Parallel이라는 글에서 참고하였다. 이 글을 쓴 사람은 유전체학 분야에서 일을 하는 것으로 여겨진다. 이 글이 실린 RONIN - our guide to cloud computing's tricky bits이라는 곳이 무슨 목적으로 만들어졌는지도 관심을 가져 보아야 되겠다. RONIN이란, AWS(아마존웹서비스) 위에서 돌아가는 연구용 클라우드 컴퓨팅 서비스 정도로 이해하면 될 것 같다.

• Announcing research computing with RONIN on AWS(2021년 3월)
• RONIN home

AWS를 1년 정도 써 본 경험에 의하면, 생각보다 요금이 많이 들었다는 기억만 남는다. BioStarts에 소개되었던 보다 차원 높은 parallel 사용법을 소개하는 것으로 끝을 맺는다.

Tool:GNU Parallel: Parallelize serial command line programs without changing them

ZGA pipeline에 두 쌍의 paired end read FASTQ 파일을 투입한 결과가 좀 이상하다

서로 다른 종에 속하는 미생물 균주로부터 유래한 일루미나 sequencing read를 한꺼번에 조립하면 결과는 엉망이 된다. 특별히 metagenomic assembler를 쓰지 않은 상황이라면 이런 불량한 결과가 나오는 것은 너무나 당연하다. 가끔은 일부러 테스트를 위해 이런 조립을 해 보기도 한다. 

매우 유사하지만 ANI 기준으로 분명히 서로 다른 종에 속하는 두 균주에서 생산한 일루미나 sequencing read가 있다. 전부 오래전에 NCBI SRA에 등록해 놓았는데 그 분량이 너무나 많은 관계로 다운로드 후 각각 150x 분량으로 서브샘플링을 하였다. 

• 균주 A: A_1.fastq, A_2.fastq
• 균주 B: B_1.fastq, B_2.fastq

SPAdes로 조립을 하려면 다음과 같이 입력해야 한다.

spades.py --threads 16 -o outdir \
--pe-1 1 A_1.fastq --pe-2 1 A_2.fastq \
--pe-1 2 B_1.fastq --pe-2 2 B_2.fastq

혼동하기 쉽게 생겼다. 예전 방식(deprecated syntax)으로는 '--pe1-1 A_1.fastq --pe1-2 A_2.fastq --pe2-1 B_1.fastq --pe2-2 B_2.fastq'의 형태로 입력하란다. Library ID number에 해당하는 숫자(빨간색)의 위치가 미묘하게 다르다. 

이렇게 조립을 하면 당연히 결과는 엉망이 된다. Contig의 수는 많고, 길이는 매우 짧으며, total length도 크게 늘어난다. 그러면 ZGA pipeline을 쓸 때에는 어떻게 명령어를 조합해야 하는가? Assembly engine으로는 SPAdes를 사용한다. 

ZGA 웹사이트의 설명을 보자.

zga -1 Lib1.R1.fq.gz Lib2.R1.fq -2 Lib1.R2.fq Lib2.R2.fq combination of reads from two sequencing libraries

SPAdes보다는 명령어 라인 구성이 간결해서 좋다. 라이브러리의 순서를 지켜서 forward read file은 -1 뒤에, reverse read file은 -2 뒤에 나열하면 된다. 그런데 결과가 너무나 좋다! A 유래 샘플만 사용한 것과 거의 다르지 않은 것이다! 이건 말이 되지 않는다. 그래서 다음과 같이 두 파일을 합쳐서 하나로 만든 뒤 이를 이용하여 ZGA를 실행하였다.

cat A_1.fastq B_1.fastq > C_1.fastq
cat A_2.fastq B_2.fastq > C_2.fastq
zga -o output -1 C_1.fastq -2 C_2.fastq

심하게 조각난 assembly가 산출되었다. 이렇게 나오는 것이 정상이다. ZGA 소스 코드에 약간의 버그가 있는 것 같다. 그것이 아니라면 나의 착각 또는 실수?

SPAdes version 3.14.0부터 도입된 '--isolate' 옵션을 쓰는 이유나 철저하게 이해하도록 하자. metaSPAdes를 돌리는 것이 아닌데도 '--isolate' 옵션을 따로 준다는 것이 잘 받아들여지지는 않지만 말이다.

Special --isolate option for assembly of standard datasets with good coverage (>100x). If you have high-coverage data for bacterial/viral isolate or multi-cell organism, we highly recommend to use --isolate option.

ART sequencing read simulator의 '-m' 파라미터는 정확히 무엇을 가리키는가?

ART sequencing simulator의 여러 파라미터 중 라이브러리의 크기와 직접 관련이 있는 것은 '-m'(the mean size of DNA/RNA fragments for paired-end simulations)이다. 길이 분포는 표정규분포를 가정하되 표준편차는 '-s' 파라미터로 설정한다. 일루미나 시퀀싱 라이브러리의 insert size와 fragment size는 보통 다음 그림과 같이 서로 다르게 정의된다.

출처: BioStars

정리하면 다음과 같다(참조: BioStars).

• insert size = the sequence between adapters
• fragment size = insert size + adapters

Fragment length는 양 끝의 어댑터를 포함하는 길이에 해당한다. CLC Genomics Workbench의 mapping 또는 de novo assembly report에서 만들어지는 paired reads distance 분포는 위 그림에서 inner distance에 해당할 것으로 생각하기 쉬운데, CLC Genomics Workbench 매뉴얼에 의하면 'The paired read distance includes the full read sequences"라고 명시되어 있다. 따라서 위에서 보인 그림의 insert size에 해당한다. 만약 이렇게 하지 않으면 forward와 reverse read가 서로 겹치는 경우 inner distance는 음수의 값을 갖게 될 것이다.

art_illumina의 '-m' 파라미터는 fragment size라고 하였으니 양 끝 어댑터의 길이를 합친 값을 제공해야 할 것처럼 보인다. 그런데 read simulator의 입장에서 어댑터의 크기를 고려할 일이 있을까? 몇 번의 시행착오를 거친 끝에 '-m' 뒤에는 insert size(CLC Genomics Workbench의 paired read distance 대푯값)에 해당하는 숫자를 넣으면 된다는 것을 알았다. 오랫동안 art_illumina를 사용해 왔지만 simulation으로 만들어진 read의 실제 분포를 측정해 본 일은 없었다. 대충 '-m 400' 정도로 놓고 read를 만들어서 사용하기만 했었다. 

다음은 실제 데이터인 SRR1144835의 자료를 샘플링하여 CLC Genomics Workbench에서 de novo assembly를 한 결과 리포트에서 딴 것이다. 분포가 예쁘지는 않지만 대략 470 bp 정도를 목표로 하였다.

다음은 '-p -l 101 -m 470 -s 30 -f 150' 파라미터로 시뮬레이션한 read의 조립 후 분포이다. 너무나도 정직한 정규분포의 모습을 따른다. Fragment(length)와 관계된 것이니 '-f' 파라미터라고 착각하면 안 된다. '-f'는 fold of read coverage를 지정하는 파라미터이다. 나도 처음에는 여러 차례 혼동했었다. 

Torsten Seemann의 블로그 "The Genome Factory"에 paired-end read의 길이에 대한 명쾌한 설명이 나온다. 이미 내 블로그에서 이 링크를 과거에 한 번 소개했었던 것 같다. 오늘 발견한 문제는 ART나 CLC에서 다루는 길이 관련 정의가 일반적으로 받아들여지는 것과는 조금 다르다는 것이다.

한국인 유래 샘플을 포함하는 장내 마이크로바이옴 카탈로그, HRGM(Human Reference Gut Microbiome)

메디톡스에서 파견 근무를 하던 2020년 4월에 생화학분자생물학회에서 발간하는 웹진의 TiBMB 섹션에 인체 장내 마이크로바이옴 참조 유전체 데이터베이스의 개발 현황이라는 글을 실은 적이 있다. TiBMB는 아마도 'Trends in Biochemistry and Molecular Biology'의 약자가 아닌가 싶은데, 정작 학회 웹사이트에서는 이것이 정확히 무엇을 의미하는지 찾기가 어렵다. 처음에는 TiBMB 자체가 학회에서 발간하는 여러 온라인 간행물의 하나라고 생각했었다.

TiBMB 화면 갈무리(링크)

이때 소개했던 Unified Human Gastrointestinal Genome(UHGG) 및 Unified Human Gastrointestinal Protein(UHGP) 관련 논문은 프리프린트 서버인 bioRxiv에만 올라온 상태였었다. 공개된 날짜는 2019년 9월 19일이었다. Nature Biotechnology에 정식으로 출판된 것은 2020년 7월 20일이다.

A unified catalog of 204,938 reference genomes from the human gut microbiome (open access)  

공개된 유전체 자료에서 인체 장에서 유래한 것을 어떻게 모으고 조립하였는지, 일정 identity 기준 이내에 들어오는 것을 어떻게 합쳤는지 등 방대한 자료 처리 과정에 대해서는 정말 공부할 것이 많다. UHGG를 만드는 데에는 isolate genome뿐만 아니라 MAG(metagenome-assembled genome)도 포함되어 있다.

논문 소개 자료를 인쇄해서 늘 책상 위에 놓고 있었는데 이제 비로소 쓸모를 발견하게 되었다.

장내 마이크로바이옴 관련 데이터베이스는 이것이 전부가 아니다. 예를 들어 같은 해에 Microbiome 저널에 실린 HumGut라는 것도 있다.

HumGut: a comprehensive human gut prokaryotic genomes collection filtered by metagenome data (open access)

초록을 조금만 인용해 보자. UHGG과 어떤 점이 다른지를 어렵지 않게 파악할 수 있을 것이다.

In this work, we aimed to create a collection of the most prevalent healthy human gut prokaryotic genomes, to be used as a reference database, including both MAGs from the human gut and ordinary RefSeq genomes. 

We screened > 5,700 healthy human gut metagenomes for the containment of > 490,000 publicly available prokaryotic genomes sourced from RefSeq and the recently announced UHGG collection.

이상의 데이터베이스에서는 한국인 유래 마이크로바이옴 자료는 찾아볼 수 없다. UHGG에는 동아시아인 유래 자료로서 중국인의 것이 꽤 많이 포함되었을 뿐이다. 국내에서도 장내 마이크로바이옴 연구가 점차 활발해지고 있는데, 한국인 분변에서 추출하여 만든 MAG 관련 연구 성과는 과연 무엇이 있을까? 검색을 해 보니 연세대학교 이인석 교수 연구팀에서 2021년에 Genome Medicine에 발표한 HRGM(Human Reference Gut Microbiome) 논문이 나왔다. 

Human reference gut microbiome catalog including newly assembled genomes from under-represented Asian metagenomes (open access)

이미 BRIC에도 소개된 중요 연구 성과였는데(링크) 모르고 지나쳤던 것 같다. 개인적으로는 별로 좋아하지 않는 '유전체 지도 구축'이라는 제목이 쓰였다. KOBIC의 대용량 컴퓨팅 설비가 쓰였다는 점도 매우 반갑다. 실제로 KOBIC에서 어느 부분을 담당했는지는 실무자를 수소문하여 물어 봐야 되겠다. Human genome에 대한 스크리닝인지, 조립 과정인지, 주석화 과정인지?

한국인(90 샘플), 일본인(645 샘플), 그리고 인도인(110 샘플)의 분변 유래 whole-metagenomic shotgun sequencing(WMS) 자료를 모아서 29,082개의 genome을 얻은 뒤, 이를 반복적으로 클러스터링하여 2199개의 species cluster를 얻었다. ㅇ32,098개의 non-redundant genome을 구성하고, 이를 UHGG의 5414 species cluster와 합쳐서 다시 후처리를 계속하여 최종적으로 5414개의 species cluster로 구성된 HRGM을 확보하게 된 것이다. 특히 한국인 샘플로부터는 ultra-deep WMS(>30 Gbp or > 100 million reads pairs)를 생산하여 low-abundance species를 찾는데 만전을 기하였다. 모든 결과물은 https://www.mbiomenet.org/HRGM/에 공개된 상태이다.

HRGM 웹사이트의 화면 갈무리.

요즘 나는 한국인 분변에서 분리한 미생물(메타게놈이 아님)의 특성을 분석하는 일을 하고 있다. 이 균주가 속하는 종이 실제 한국인 분변에 얼마나 존재하는지를 점검하는데 이인석 교수 연구팀에서 생산한 SRA 자료(NCBI Sequence Read Archive SRP292575)가 매우 유용하게 쓰일 것 같다. 각 샘플에 존재하는 특정 taxon의 relative abundance는 KRAKEN 2를 이용하여 판별한 것으로 보인다. 내가 쓰는 컴퓨터에서는 전체 시퀀싱 raw data를 내려받는 것만으로도 벅찰 것만 같다.