phyloFlash를 수동으로 설치해 보니...

편리함으로 인하여 얻는 것이 있다면, 당연히 잃는 것도 있다. 편리함이라는 껍데기 속에 숨어있는 기초적인 원리를 알 수 없게 된다. phyloFlash는 whole-metagenome shotgun sequencing data로부터 full-length 16S rRNA 염기서열을 재구성해주는 Perl script 모음이다. Conda를 이용하여 쉽게 설치할 때에는 전혀 모르고 있던 16S rRNA reference DB의 설치 과정('phyloFlash_makedb.pl --remote')이 생각보다 간단하지 않았다. vsearch, bowtie, BBMap(bbmask.sh & bbmap.sh) 등 추가적으로 필요한 프로그램이 몇 개 더 있었다. 이를 다 확인한 다음에는 NCBI에서 UniVec 데이터베이스를 다운로드하고, 마지막으로 가장 중요한 단계인 SILVA rRNA DB를 가져다가 재가공을 거쳐 설치한다. DB의 설치 위치는 phyloFlash.pl 메인 스크립트 실행 때 -dbhome 파라미터 혹은 PHYLOFLASH_DBHOME 변수로 지정한다.

SILVA DB의 다운로드 및 재가공 과정은 다음과 같다.

1. SIVA DB ftp site에서 /current/Exports/*_SSURef_N?99_tax_silva_trunc.fasta.gz 파일을 다운로드하여 압축 해제(오늘 기준으로 SILVA_138.1_SSURef_NR99). 라이선스 정보(Creative Commons Attribution 4.0, CC-BY 4.0)를 읽어 볼것.
2. barrnap_HGV를 이용하여 SSU에 섞인 LSU contamination을 검색하여 제거. 120개의 염기서열이 사라졌다.
3. bbmask를 사용하여 low entropy region을 마스킹.
4. bbmap을 사용하여 UniVec contamination을 제거.
5. 'vsearch --makeudb_usearch'를 사용하여 UDB 파일 생성. Robert Edgar(drive5.com)의 설명에 의하면 UDB database file이란 서열과 그것의 k-mer index를 수록한 파일이라 한다(링크). 바이너리 파일이므로 .udb 파일을 텍스트 편집기로 열어 보는 것은 의미가 없다.
6. 'vsearch --cluster_fasta --id 0.96' 조건으로 사용하여 클러스터링.
7. bowtie로 인덱싱(bowtie2도 아니고...)

상당히 신경을 써서 후처리를 거치고 있다. Sequence identity threshold를 96%로 맞추어서 클러스터링을 한다는 점이 중요하다. 단계 6번이 실제로 어떻게 돌아가는지를 알기 위하여 화면에 표시되는 메시지를 확인해 보았다.

[09:50:44] running subcommand:
	   /usr/bin/vsearch   --cluster_fast
	   ./138.1/SILVA_SSU.noLSU.masked.trimmed.NR99.fasta --id 0.96
	   --centroids ./138.1/SILVA_SSU.noLSU.masked.trimmed.NR96.fasta
	   --notrunclabels --threads 32 
vsearch v2.14.1_linux_x86_64, 125.7GB RAM, 32 cores
https://github.com/torognes/vsearch

Reading file ./138.1/SILVA_SSU.noLSU.masked.trimmed.NR99.fasta 100%  
609727540 nt in 421417 seqs, min 803, max 3706, avg 1447
Masking 100% 
Sorting by length 100%
Counting k-mers 100% 
Clustering 100%  
Sorting clusters 100%
Writing clusters 100% 
Clusters: 156798 Size min 1, max 1184, avg 2.7
Singletons: 107548, 25.5% of seqs, 68.6% of clusters

클러스터링을 거쳐 최종적으로 만들어지는 파일은 다음과 같다. SILVA_SSU 뒤에 새로 붙은 문자열로부터 어떤 가공을 거쳤는지를 알 수 있다.

./138.1/SILVA_SSU.noLSU.masked.trimmed.NR96.fixed.fasta

Metagenome assembly를 일상적으로 실시하게 된 요즘, 동일한 데이터로부터 16S rRNA 염기서열을 찾아내는 일이 과연 무슨 의미가 있겠느냐고 할 수도 있겠으나, 꼭 그렇지만도 않다는 것이 나의 생각이다.