이 페이지에서는 서로 관련성이 높은 다수 미생물의 미생물 유전체 시퀀싱 데이터를 가지고 수행하는 전형적인 생명정보학적 분석 방법을 다루고자 한다. 관련성이 높다고 함은 동일한 종에 속하는 균주(strain)의 유전체를 분석함을 뜻힌다. 실제 상황에서는 동일 종으로 알려져 있지만 유전체 서열을 이용하여 ANI 혹은 이와 유사한 분석(JSpecies, pyani, MUMi, etc)을 하면 다른 종으로 보아야 하는 것들을 종종 만나게 된다.

가장 기초적인 분석 목표는 각 유전체의 관계를 phylogenetic tree의 형태로 그려내는 것이다. 만약 assembly를 거쳐서 gene sequence를 알아낼 정도가 되었다면 pan genome analysis를 하는 것이 가능하다. 이 과정에서 core genome alignment를 얻어내어 더욱 정확한 계통수를 그려도 되고, SNP만 추출하여 계통수 분석을 할 수도 있다. 동일 소스에서 시간차를 두고 미생물 균주를 분리하여 시퀀싱을 한 longitudinal data가 있다면 molecular clock estimation도 가능하다.

Reference genome sequence에 매핑하여 SNP를 발굴한 뒤 계통수 분석을 하는 것이 가장 일반적이다.

Draft level이든 complete level이든 조립된 상태의 유전체 서열이 있다면 gene prediction이 가능하므로 좀 더 폭넓은 비교 분석이 가능하다. 세균 유전체의 대용량 비교 분석에서는 pan-genome analysis가 가장 대표적인 작업이 될 것이다. Pan-genome은 2005년 Tettelin 등이 주장하여(논문 링크) 이제는 아주 보편적인 개념이 되었다. 어떤 clade(주로 prokaryote species)를 이루는 전체 유전체(pan-genome 또는 supra-genome)은 모든 균주가 공통적으로 보유한 core genome과, 각 균주가 개별적으로 갖고 있는 accessory genome의 총합이라는 것이다.

Pan genome 분석의 필수 요소는 orthologous gene을 찾아내는 방법이라고 할 수 있다. 이에 대해서 깊이 들어갈 수는 없지만 대표적인 자료를 아주 조금 소개하고자 한다. Ortholog 분석 도구는 functional annotation의 한 방법으로 쓰이기도 한다.

• NCBI - Clusters of Orthologous Group (COG)
• The quest for orthologs: finding the corresponding gene across genome. Trends Genet (2008)
• http://inparanoid.sbc.su.se/
• http://eggnogdb.embl.de/#/app/home
• https://questfororthologs.org/

2015년에는 pangenomic를 위한 소프트웨어 도구를 리뷰한 논문이 나오기도 하였다(링크. 물론 그 이후에 나온 소프트웨어도 있다. Core genome이 산출되면 이를 바탕으로 유전자 염기서열의 aliginment를 하여 SNP 발굴 및 phylogenetic tree 작성도 가능하다.

• [논문] dRep: a tool for fast and accurate genomic comparisons that enables improved genome recovery from metagenomes through de-replication. ISME J. (2017)
• GitHub Official documentation

dRep은 원래 대용량 메타게놈 자료의 분석을 위해 만들어진 것이다. time-series metagenome dataset이 있다고 가정하자. 만약 전체 데이터를 한꺼번에 조립한다면(co-assembly) 데이터 분량도 많고 복잡도도 크게 증가하여 오히려 완성도가 떨어진다. 따라서 각 샘플을 별도로 조립한 뒤 dRep을 거치면 근본적으로 동일한 유전체를 찾아내어 각 세트로부터 최적의 유전체를 확보할 수 있다는 것이다. Genome de-replication과 genome comparison으로 이루어지는 dRep의 작동에 대한 좀 더 상세한 설명은 공식문서를 참조하라. 이 페이지에서 dRep을 소개하는 것은 동일 종의 여러 strain에 대한 유전체 조립 결과를 갖고 있을 때 dRep을 이용하여 연관성이 떨어지는 것을 찾아내는 일종의 전처리 용도로 쓸 수 있을 것이라는 기대 때문이다. dereplicate_wf와 compare_wf(wf = workflow) 중에서 일반적인 유전체의 고속 비교 용도에는 compare_wf가 더 유용할 것이다.

$ pyenv shell 3.5.1
$ pyenv versions
$ python test_suite.py # drep_설치_디렉토리/test로 이동하여 실행
$ dRep dereplicate_wf drep_out -g ./*.fna
$ dRep compare_wf comp_out -g ./*.fna

compare_wf를 이용한 명령어 사례를 소개해 본다. Secondary clustering(–S_algorithm)의 기본값은 ANIm이다.

$ dRep compare_wf WORK_DIRECTORY -p 20 --S_algorithm gANI -sa 0.95 -g ./*fna

• gANI를 사용하려면 JGI의 ANIcalculator과 prodigal이 $PATH에 있어야 한다.
• 동일 species 내의 유전체는 ANI >=96.5%가 일반적이다(Varghese 2015). -sa (minimum secondary ANI)의 default value는 99%이다. 이 값은 사실은 종 동정이 아니라 동일한 genome을 골라내기 위한 것이었다.
• conda py35 환경에서는 명령어 문법이 약간 다르다. 예를 들자면 dRep compare_wf가 아니고 dRep compare이다.

• [웹사이트] http://www.sanger.ac.uk/science/tools/roary
• [논문] Roary: rapid large-scale prokaryote pan genome analysis. Bioinformatics (2015) Supplementary material
• GitHub (tutorial 포함)

일반적으로 pan genome을 계산하기 위하여 all-against-all BLAST를 거치므로 대단히 많은 시간이 소요되지만, Roary는 CD-HIT을 사용하여 partial sequence를 사전에 제거하고 pre-clustering을 실시하여 후속 all-against-all 검색에 사용할 서열집합의 수를 크게 줄이므로 다른 도구에 비하여 훨씬 빠르게 분석을 수행하게 된다. 계산이 끝난 뒤 union/intersection/difference 산출을 도와주는 명령어(query_pan_genome) 등 다양한 유틸리티를 제공한다는 점이 특징이다.

설치가 약간 까다로운 편이지만(설치 가이드 요즘은 Virtual Machine과 Docker image를 제공하고 있다. Roary의 기본적인 사용법을 알아보자. 필요한 입력물은 GFF3 파일이 유일하다. 유전체 서열만 준비된 경우라면 Prokka를 이용하여 genome annotation을 하면 되고, NCBI에서 다운로드한 GenBank 파일이 있따면 bp_genbank2gff3.pl(BioPerl에 포함)를 이용하여 GFF3 file을 만들면 된다. GenBank file에서 /pseudo로 표시된 유전자는 초기 단계에서 아예 단백질 서열이 추출되지 않으므로 계산에서 제외된다.

$ roary -a # 모든 소프트웨어가 설치되었는지 점검
$ roary *.gff # core gene의 multiple sequence alignment를 하지 않고 pan genome 분석(기본 동작)
$ roary -e --mafft -p 8 -i 90 -f output_dir -z *.gff # 본문 설명 참조

세 번째 명령에서는 core gene의 alignment를 mafft로 실행하고(-e –mafft) 8 개의 thread를 사용하며(-p 8) blast percent identity cutoff를 90%(기본은 95%)로 한다는 것(-i 90)이다. -e 옵션만 주면 PRANK에 의한 core gene alignment를 실시한다(좀 느리다). -o output_dir을 설정하지 않으면 현재 디렉토리에 결과 파일을 작성한다. -e라고만 하면 prank를 사용하여 core gene multiFASTA alignment를 하게 되는데 codon-aware alignment를 하므로 mafft보다는 정확하지만 느리다. -z는 중간 결과 파일(alignment)을 output_dir/pan_genome_sequences 디렉토리에 남겨둔다는 뜻이다. _로 시작하는 10여 개 중간 결과 파일도 지우지 않고 남기며, 현 디렉토리에는 각 유전체에서 추출한 단백질 서열 파일을 담은 임시 디렉토리도 남기게 된다.

NCBI에서 공식적으로 release한 유전체 정보(GenBank)가 아니라 PGAP만 진행된 정보를 사용하는 경우 주의를 요한다. 왜냐하면 GenBank 파일을 이루는 서열(complete가 아닌 경우 서열들)에 대한 Accession이 아직 부여되기 전이라서(아래 그림 참조), GFF 파일로 전환하는 경우 서열의 ID가 'unknown'으로 기록된다. 이렇게 만들어진 GFF 파일을 Roary에 사용하면 에러가 발생한다.

그러므로 PGAP annotation file은 그대로 사용하지 말고 fillAccession.pl으로 처리하여 accession을 채운 뒤 GFF로 전환해야 한다.

또한 RAST server에서 export한 GFF3 파일도 Roary에서 그대로 쓰일 수가 없다. 왜냐하면 염기서열이 뒷부분에 있지 않기 때문이다. 뿐만 아니라 gene 없이 cds feature만 있다는 것도 문제가 된다. 몇 가지를 테스트해 본 경험으로 가장 바람직한 것은, RAST에서 export한 GFF3 file에서 CDS feature만 골라낸 것 + '##FASTA' 라인 + contig sequence 파일(유전자 염기서열 파일이 아님!)을 합쳐서 새로 만든 GFF3 파일을 사용하는 것이 좋다. 'Name='도 'Product='로 바꾸는 것을 강력 권장한다. 왜냐하면 이것이 gene id로 쓰이게 되면 중간에 공백이 들어 있어서 나중에 매우 불편해지기 때문이다.

Roary를 실행하면 십중팔구는 다음과 같은 메시지와 함께 GFF 파일을 수정하여 fixed_input_files 디렉토리에 복사하게 된다. 나중에 query_pan_genome 스크립트로 GFF 파일을 대상으로 하는 작업을 할 때에는 수정된 것을 써야 한다. 2021/06/30 15:53:28 Input file contains duplicate gene IDs, attempting to fix by adding a unique suffix, new GFF in the fixed_input_files directory:

• gene_presence_absence.csv: Spreadsheet containing statistics on each group and presence and absence of genes in each isolate
• core_gene_alignment.aln: MultiFASTA alignment of the core conserved genes
• clustered_proteins: File with 1 cluster per line and sequence identifiers. gene_presence_absence.*의 첫번째 컬럼(Gene)과 이 파일의 첫번째 컬럼에 쓰인 식별자가 동일하다. 또한 이 파일에 수록된 모든 cluster의 수(=라인 수)는 summary_statistics.txt 파일 맨 아랫줄의 total genes와 같다.
• *.Rtab: Tab Files for producing graphs in R
• *embl: EMBL files for visualizing presence and absence of genes
• accessory_binary_genes.fa.newick: Newick tree based on presence and absence of genes in accessory genome
• pan_genome_reference.fa: A FASTA file with a single REPRESENTATIVE nucleotide sequence for every cluster in the pan genome (core and accessory)
• pangenome_sequences (directory): MultiFASTA alignment files for each cluster

Roary에서는 accessory gene의 여부를 이용한 binary tree만을 만들어 주므로, core gene alignment에서 phylogenetic tree를 작성하려면 FastTree 등의 프로그램을 직접 돌려야 한다.계산 시간을 줄이려면 snp_site(논문)로 필터링을 하여 사용하라. 이 과정에서 recombination은 고려하지 않음에 유의한다.

$ FastTree –nt –gtr core_gene_alignment.aln > my_tree.newick

Roary의 기본 실행 조건에서는 모든 paralog를 분리하여 core gene은 각 유전체에 대하여 하나씩만 존재하게 만든다. 그러나 -s 옵션을 주면 이를 분리하지 않는다. 즉 gene_presence_absence.csv 테이블의 Avg sequences per isolate 컬럼에 1보다 큰 값이 들어갈 수 있는 것이다.

원본 GFF 파일이 아니라 roary가 수정한 것(fixed_input_files/ 디렉토리에 있는 것)을 사용해야 한다. 인수로 주어질 수 있는 것은 *gff 또는 단일 GFF 파일이다.

$ query_pan_genome -h # for help
$ query_pan_genome -a union *.gff # 결과물: pan_genome_results
$ query_pan_genome -a intersection *.gff # 결과물: pan_genome_results
# difference
# complement

-g clustered_proteins 옵션은 항상 필요한 것 같지는 않다. 다음으로 특정 유전자 클러스터의 multiFASTA file을 추출하려면 gene_presence_absence.csv의 첫번째 컬럼에서 유전자 이름을 확인한 뒤 다음과 같이 -n 옵션 뒤에 나열하여라.

$ query_pan_genome -a gene_multifasta -g clustered_proteins -n Rv0498,MT_RS07820,tcrA ../*gff

이렇게 하면 pan_genome_results_tcrA.fa, pan_genome_results_MT_RS07820.fa, 및 pan_genome_results_Rv0498.fa 아미노산 서열 파일이 생긴다. 염기서열을 뽑을 수도 있다고 하지만 어떻게 하는지를 모르겠다. 만약 유전자 염기서열이 필요하다면 pan_genome_sequences 디렉토리 아래에 각 유전자별로 정돈된 aligned fasta file(.aln 확장자인데 clustal 프로그램의 .aln 형식이 아님!)을 사용하는 것이 나을 것이다.

roary2svg.pl을 이용하면 gene_presence_absence.csv 파일을 사용한 그림 파일을 만들어 준다. 웹브라우저에서 읽으면 된다. 다음은 한 사례이다. input file 중 두 건이 잘못되어 유전자 정보가 하나도 없었고, 그로 인하여 core genome이 산출되지 않은 것이다. 이처럼 roary 실행의 진단 목적으로 쓸 수도 있는 것이다.

R에서 gene_presence_absence.Rtab과 gene_presence_absence.csv 두 파일을 다루면 된다. 다음의 사례에서는 Lb_1-46 균주에서 특이적인 유전자의 id를 추출하는 사례를 보여준다. Lb_1-46은 데이터프레임으로 읽어들이면 Lb_1.46으로 바뀌는 것에 유의해라. 구글링을 잘 하면 이를 원래 이름 그대로 유지하는 방법이 있다(하루에 한 R 관련글).

> dat = read.table("gene_presence_absence.Rtab",sep="\t",header=T,row.names=1) 
> dat$Lb_1.46
> dat.s = subset(dat, rowSums(as.matrix(dat))==1)
> Lb.s = dat.s[which(dat.s$Lb_1.46==1),]
> Lb.s.genes = rownames(Lb.s)
> dat.2 = read.table("gene_presence_absence.csv",sep=",",header=T,row.names=1)
> dat.2$Lb_1.46
> dat.f = dat.2[which(rownames(dat.2) %in% Lb.s.genes),]
> dat.f$Lb_1.46
> write.table(dat.f$Lb_1.46,"out.txt",sep="\t",quote=F,col.names=F)

결과 파일을 열어보면 일부 단백질의 ID가 변형되어 쓰인 것을 알 수 있다. 즉 원본 annotation file에서 MT_RS20470라는 locus tag을 갖는 유전자가 “MT_RS20470.p01_16540”으로 바뀐 것이다. .p01은 유전자에서 단백질로 번역됨을 나타내는 것이겠지만 “_16540”은 무엇인가? (원래 '_'가 3회 반복된 것이나 DokuWiki에서 그렇게 타이핑하면 긴 가로선이 나오기 때문에 부득이하게 하나만 타이핑하였다)

Prokka에서 만든 gff3 파일을 사용하였더니 tRNA gene을 결과물 중에 포함시키는 현상이 가끔 관찰된다.

• BPGA- an ultra-fast pan-genome analysis pipeline Scientific Reports (2016)

• [논문] The Harvest suite for rapid core-genome alignment and visualization of thousand of intraspecific microbial genomes. Genome Biol. 2014; 15(11)l 524.
• GitHub Official documentation

Roary는 단백질 암호화 유전자를 단위로 core genome을 산출하지만 Harvest suite는 유전체 서열을 그대로 이용한다는 것이 가장 큰 차이점이다.

설치와 사용은 매우 쉽다. 압축파일을 풀고 그 디렉토리에 들어가서 parsnp, gingr, harvesttools를 실행하면 그만이다. 실행 전에 PATH를 설정한다.

$ PATH=/data/apps/Harvest-Linux64-v1.1.2:$PATH

$ ./parsnp -r reference.fna -d <genome_dir> -p <threads> 

• Reference가 GenBank file이라면 -g reference.gbk라고 입력한다.
• -c를 지정하면 genome_dir에 있는 모든 유전체 서열을 전부 사용한다. 이를 지정하지 않으면 MUMi를 계산하여 지정된 MUMi distance threshold 내에 있는 유전체만 골라서 full alignment를 실시한다. 평소에는 이를 지정하지 않는 것이 바람직할 것이다. 왜냐하면 reference에 비하여 지나치게 먼 유전체 서열이 포함되면 full alignment의 quality가 현저히 떨어질 것이 자명하기 때문이다.
• -o <output_directory>를 지정하지 않으면 P_2017_11_13_130215334092(./P_CURRDATE_CURRTIME)와 같은 포맷으로 자동 생성된다.

출력물 디렉토리에는 1) newick format core genome SNP tree(.tree), SNP 정보(.vcf), gingr용 바이너리 파일(.ggr), 그리고 xmfa 포맷의 MSA 파일(.xmfa)이 생긴다.

Gingr은 대용량의 phylogeny와 multiple sequence alignment를 시각화하는 도구이다. 실행법은 매우 간단하다.

$ path_to_gingr/gingr ./parsnp.ggr

Gingr 대화형 사용법을 설명한 파워포인트 파일(Version 0.8, 2017-11-13 작성)

• [논문] The large-scale blast score ratio (LS-BSR) pipeline: a method to rapidly compare genetic content between bacterial genomes. Peer J (2014)
• [GitHub] 매뉴얼 PDF
• Python >2.7 and ⇐3.5

LS-BSR은 엄밀히 말하자면 SNP를 생성하는 도구는 아니다. BSR은 blast score ratio의 약자로서, 어떤 query protein이 다른 단백질과의 사이에서 얻은 blast score를 자기 자신에 대한 blast score로 나눈 것을 의미한다([정해영의 블로그]에서 소개한 글). 유전자의 상동성 여부를 판별하는 매우 간단한 방법으로서 여러 유전체를 대상으로 비교를 하면 BSR matrix가 얻어지고, 이를 다른 도구를 사용하여 시각화할 수 있다. 비교 대상 단백질 서열은 외부에서 제공하거나(gene screen method) 혹은 gene prediction을 이용하여 자체적으로 만들면 된다. 일종의 pan genome 분석 도구라고 보아도 좋다.

LS-BSR은 Roary와 비슷하면서도 다르다. gene presence or absence matrix를 제공한다는 것은 두 프로그램이 유사하다.

• Gene prediction: LS-BSR만 제공(Roary에서는 gene sequence가 하단에 포함된 GFF file을 입력해야 함)
• Roary는 비교 대상 유전체가 갖고있는 유전자만을 대상으로 계산 작업을 하지만, LS-BSR은 미리 준비된 peptide sequence set를 공급하여 계산을 할 수 있다.
• Roary는 core gene alignment를 제공하지만 LS-BSR에서는 blast score ratio matrix를 제공

Phylogenetic tree를 그린다면 무엇이 유리한가? Core-gene SNP를 사용하여 tree를 그리고자 한다면 Roary가 적합할 것이고(tree file은 alignment로부터 FastTree 등을 돌려야 하고, gene presence absence에 의한 binary tree는 기본적으로 제공한다), heat map을 그리고자 한다면 LS-BSR이 더 나을 것이다.

LS-BSR의 간단한 사용법을 알아보자. genome sequence 파일은 확장자가 .fasta가 아니면 작동을 하지 않는다. markers.fasta는 항상 유전자 염기서열 파일이어야 한다. 검색용 프로그램은 기본이 tblastn이며(단백질로 번역을 먼저 함), 만약 blastn을 쓰고 싶으면 -b blastn을 더해야 한다. -x test라고 지정하면 test라는 디렉토리를 임시로 생성하여 중간 계산 파일을 넣은 뒤, 최종 결과 파일은 'test_'로 시작하는 접두사를 갖고 만들어진다. 분석이 끝나면 test 디렉토리는 없어진다. /usr/local/apps/LS-BSR/ls_bsr.py을 사용해라. 비록 conda를 이용하여 LS-BSR을 설치했다 하여도 이는 environment만 구성할 뿐, 실제 작동 스크립트 등은 git로 깔았기 때문이다. — Haeyoung Jeong 2019/06/27 11:22

# gene screen method (peptides in interest are ready)
$ python path_to_LS-BSR/ls_bsr.py -d genomes -g markers.fasta -x test
# de novo gene prediction method
$ python path_to_LS-BSR/ls_bsr.py -d genomes -c usearch -x test

계산이 끝나면 현 디렉토리에 markers.fasta가 번역된 genes.pep 파일이 생긴다. gene prediction method에서는 prodigal로 유전자를 예측한 뒤 단백질로 번역한 다음 usearch로 클러스터링(default: -i 0.9)하여 사용한다. 클러스터링 방법은 usearch/vsearch/cd-hit이 가능하다($PATH에 있어야 함).

별도의 내부 페이지 data_visualization에서 상세히 다룬다.