16S rRNA gene amplification이 아니라 whole-genome shotgun sequencing으로 얻은 metagenome read로부터 샘플 집단의 미생물 구성을 추정하고 싶다면 어떻게하는 것이 가장 적당한 방법일까? Reference mapping을 통해서 어떤 taxon에 속하는 read가 몇%라고 단순히 분류하는 것이 가장 단순해 보이지만(read classification), 여기에는 genome size에 의한 바이어스가 필연적으로 발생하게 된다. 이에 대한 설명은 metaphyler 논문에 포함된 그림 1을 참조하라. 따라서 single copy로 존재하는 marker를 활용하는 것이 더욱 바람직하다.

다시 한 번 강조하자면, read classification과 abundance estimation은 서로 다르다는 것이다.

• [논문] Accurate and fast estimation of taxonomic profiles from metagenomic shotgun sequences. BMC Genomics 12(Supple 2)LS4
• [웹사이트] http://metaphyler.cbcb.umd.edu/

MetaPhyler에서 선정된 31개의 marker gene은 Genome Biol 2008에서 사용된 것으로, 대부분 genome 내에 single copy로 존재하면 정보 처리(replication, transcription, and translation) 또는 central metabolism에 관여하는 것이라 HGT와는 비교적 무관한 것들이다.

• dnaG, frr, infC, nusA, pgk, pyrG, rplA, rplB, rplC, rplD, rplE, rplF, rplK, rplL, rplM, rplN, rplP, rplS, rplT, rpmA, rpoB, rpsB, rpsC, rpsE, rpsI, rpsJ, rpsK, rpsM, rpsS, smpB, and tsf

웹사이트에 게시된 V1.25(05/23/2012)보다는 short Illumina read에 특화된 MetaPhylerSRV0.115.tar.gz를 쓰는 것이 좋다. Fasta file만을 입력물로 받기 때문에 fastq 파일은 fastq2fasta.pl, fq2fa(idba에 포함), seqtk(seqtk seq -a in.fq.gz > out.fa) 등을 사용하면 된다.

./metaphyler.pl <number of threads> <in.fasta> <output prefix>

결과물은 prefix.classify.tab 파일과 prefix.genus|family|order|class|phylum.tab의 두 가지 부류로 만들어진다. 실제 .family.tab 파일의 사례는 metaphyler_example여기를 클릭해 보라. 결과파일에 나타나는 depth of coverage는 31개 phylogenetic marker gene의 coverage에 대한 median 값이며, 이를 근거로 하여 abundance가 추정되는 것이다.

"prefix".classify.tab
    column 1: query sequence id
    column 2: phylogenetic marker gene name
    column 3: best reference gene hit
    column 4: % similarity with best hit
    column 5: classification rule
    column 6-10: taxonomic label at genus,family,order,class,phylum level

"prefix".genus|family|order|class|phylum.tab
    column 1: taxonomic clade name
    column 2: % relative abundances
    column 3: depth of coverage of genomes
    column 4: number of sequences binned to this clade
    column 5: similarity with reference genes (only available at the genus level)
    

metaphyler는 taxonomic rank에 따라서 서로 다른 classification threshold를 적용하므로, taxon을 할당할 충분한 증거가 없으면 하위의 taxonomic group으로 내려보내지 않는다. 예를 들어 어떤 read들이 family Enterobacteriaceae로 분류되었으나 그 이하의 속(genus)으로 내려가기 어려우면, 이러한 novel species는 Enterobacteriaceae{family}라고 표기된다.

예전에 블로그에 포스팅한 Metaphyler 결과물로 krona chart 그리기에 잘 설명되어 있다.

• [논문] PhyloSift: phylogenetic analysis of genomes and metagenomes. PeerJ. 2014; 2: e243.
• [웹사이트] https://phylosift.wordpress.com/ Intro Tutorial Input data and command line instructions GitHub
• [일루미나 2 x 100 nt 데이터를 사용한 튜토리얼] 링크
• [tube server 내 설치 위치] /usr/local/apps/phylosift_v1.0.1

이 경우는 fastq file이 아니라 contig sequence를 대상으로 하여 마커 유전자가 온전하게 존재하는지를 테스트하는 것이다. phylosift의 search mode를 사용하면 된다. 참조 사이트

$ phylosift search --unique --besthit --isolate contigs.fa

PS_temp/contigs.fa/blastDir/marker_summary.txt 파일을 열어서 DNGNGWU00001-00040의 37개 마커가 온전히 한 벌씩 있는지 확인한다(PhyloSift reference marker genes. 마커의 번호는 40까지 있지만 4/8/38번은 더 이상 사용되지 않으므로 총 37개이며, 13번 마커는 가끔 복수로 존재하기도 한다. 물론 이 분석 방법이 완벽한 것은 아니다. 출력 디렉토리를 바꾸려면 –output <directory name> 옵션을 준다.

만약 여러 assembly에 대하여 phylosift를 동시에 실행하려면 GNU Parallel을 쓰는 것이 좋다. 작업한 contig sequence file의 목록이 in_file.list에 수록된 경우 다음과 같이 실행한다.

$ parallel -a in_file.list phylosift search --unique --besthit --isolate

각 assembly에서 확인된 특정 마커 서열을 뽑아내어 tree 작성 등 후속 작업을 하고 싶다면 추출된 서열을 이용하여 align mode로 phylosift를 돌려야 한다. phylosift search를 실행한 후에 해야 하며, 출력 디렉토리도 이전 것을 그대로 쓴다는 것을 잊지 말아라! 단, 이전 단계(search)에서 만들어진 결과 디렉토리에 대해서 그대로 작업을 해야 하므로 –force 옵션이 필요하다.

$ phylosift align --force --unique --besthit --isolate contigs.fa
$ (parallel 쓰기) parallel -a in_file.list phylosift align --force --unique --besthit --isolate

phylosift align을 실행하면 PS_temp/contigs.fa/alignDir에 ${markerID}.1.codon.updated.1.fasta(aligned)가 생긴다. 마커 유전자 하나에 다하여 이들을 전부 모아서 하나의 파일(aligned fasta)로 만들려면 post_phylosift.pl 스크립트를 사용하면 된다. 다음 명령의 결과 DNGNGWU00013.aln이 생긴다. 16S rRNA gene sequence marker의 이름은 16s_reps_bac이다. 모든 genome assembly에서 16S rRNA gene marker를 성공적으로 찾는 것은 아니니 조심해야 한다.

$ post_phylosift.pl PS_temp DNGNGWU00013

marker_summary.txt를 전부 종합하여 결과 매트릭스를 만드는 방법에 대해서는 [하루에 한 R] PhyloSift 결과 종합하기(tab-delimited CSV 파일의 일괄 입력 및 병합 등)을 참고하라.

Version 1.0 매뉴얼 프로그램 소스 및 miniKraken DB

가장 표준적인 사용법은 다음과 같다. –preload 옵션이나 kraken-filter 적용은 필수 사항은 아니다. kraken 결과물의 confidence filtering에 대한 상세한 설명은 여기를 참고하라. threshold는 지정하지 않으면 0이다. query file을 표준 입력으로 제공하는 것이 아니라면, fasta와 fastq 포맷 여부를 자동으로 검출하여 작업하게 된다.

$ DB=/data/Utilities/DB/standard_kraken # 각자 상황에 맞게
$ kraken --db $DB seqs.fa # 화면으로 출력
$ kraken --threads 16 --db $DB infile.fastq --preload | \
kraken-filter --db $DB --threshold 0.05 > infile.kraken.filtered
# 파이프로 한번에 실행하지 않고 kraken과 kraken-filter를 별도로 돌려도 무방하다.

각 read에 부여된 taxonomic name이 궁금하면 kraken-translate를 이용한다.

$ kraken-translate --db $DB infile.kraken.filtered > infile.kraken.filtered.label

kraken 결과, 이를 kraken-filter 적용한 것, 그리고 kraken-translate로 이름을 부여한 결과 파일을 서로 비교해 보자. 먼저 kraken만 실시한 결과를 보자. 첫번째 컬럼(C)은 classify되었다는 것, 두번째 컬럼은 read ID, 세번째는 NCBI tax id, 네번째는 read length, 그리고 마지막은 tax id:k-mer number이다. 즉 첫번째 k-mer 7개는 taxid = 0에, 그 다음 k-mer 2개는 taxid = 1491에 매치하였다는 뜻이다. 여기에서 알 수 있듯이 하나의 read가 갖는 k-mer는 여러 taxid에 거쳐서 분포할 수 있고, 이를 종합하여 confidence filtering을 하는 것이다.

# kraken only
C	K00171:228:H7FCKBBXX:5:1101:12591:2668	1491	101	0:7 1491:2 0:62
C	K00171:228:H7FCKBBXX:5:1101:12591:2668	1491	101	0:39 1491:5 0:27
# kraken-filter
C	K00171:236:H7JTCBBXX:3:1101:26433:1543	441771	101	P=0.437	36826:15 441771:31 36826:25
C	K00171:236:H7JTCBBXX:3:1101:26433:1543	36826	101	P=1.000	36826:71
# kraken-lable
00171:228:H7FCKBBXX:5:1101:20679:2650	root;cellular organisms;Bacteria;Terrabacteria group;Firmicutes;Clostridia;Clostridiales;Clostridiaceae;Clostridium;Clostridium botulinum
K00171:228:H7FCKBBXX:5:1101:13839:2686	root;cellular organisms;Bacteria;Terrabacteria group;Firmicutes;Clostridia;Clostridiales;Clostridiaceae;Clostridium;Clostridium botulinum

필터링을 해도 출력 파일의 라인 수는 그대로이다. 필터링을 거치면 classify가 되지 않은 read('U')가 약간 늘어난다. 다시 말하자면 부정확하게 classify되었던 read('C')가 unclassified read로 바뀐다고 보면 된다.

세분화된 각 taxon에 대한 read의 비율이 tap delimited 형식으로 나온다(표준출력). Kraken report라고 하면 바로 이 kraken-report 명령으로 얻은 출력물을 의미한다(나중에 Braken을 사용할 때 직접적으로 필요함).

$ kraken-report --db $DB kraken.output > kraken.output.report

kraken 결과물 중 두번째와 세번째 컬럼만을 떼어내어야 한다. Kraken report를 이용하는 것이 아니다!

$ cut -f2,3 infile.kraken.filtered > infile.kraken.filtered.krona
$ ktImportTaxonomy infile.kraken.filtered.krona -o infile.kraken.filtered.krona.html

MIDAS = Metagenomic Intra-species Diversity Analysis System

• [논문] An integrated metagenomics pipeline for strain profiling reveals novel patterns of bacterial transmission and biogeography. Genome Res (2016) 26(11): 1612
• [GitHub] https://github.com/snayfach/MIDAS

설치한 위치에 주의하라.

$ export PYTHONPATH=/home/hyjeong/.linuxbrew/Cellar/python/2.7.12_1/lib:/data/Utilities/git/MIDAS
$ export PATH=$PATH:/data/Utilities/git/MIDAS/scripts
$ export MIDAS_DB=/data/Utilities/DB/midas_db_v1.2
$ run_midas.py species output_dir -1 file_1.fastq.gz -2 file_2.fastq.gz