세균의 종을 구성하는 유전체는 모든 균주에게 공통적으로 존재하는 유전자 집합인 core genome(orthologs)과 각 균주에 특이적으로 존재하는 유전자의 모임인 accessory genome의 두 부분으로 나누어 정의할 수 있다. Core genome과 accessory genome의 합을 pan-genome이라 부르며, 종 수준에서 세균 유전체의 진화가 일어날 수 있는 실질적인 경계로 여겨진다. 물론 미생물 세계에서는 종 단위를 뛰어넘은 유전자의 수평적 이동이 빈번하게 일어나는 것으로 알려져 있다.

Pan-genome analysis는 비교유전체학 분석의 대표적인 방법으로서 core gene sequence를 연결하여 모든 균주의 진화적 상관관계를 보여주는 계통수를 작성하거나, accessory gene list를 탐색하여 균주 특이적인 형질을 결정하는 유전자를 찾을 수도 있다. Pan genome analysis의 기본적인 개념은 Torsten Seemann의 [SlideShare] Comparing bacterial isolates 자료를 참고하도록 한다. Pan-genome 분석 도구에 대한 비교적 최근의 리뷰는 NCBI Books의 The Pangenome: Diversity, Dynamics and Evolution of Genomes가 좋은 참고가 될 것이다.

기존에 작성한 문서: Roary를 이용한 미생물 유전체의 비교 분석 - 겹치는 내용이 많음!

기본적인 pan-genome analysis는 Roary 파이프라인을 통하여 진행할 수 있다. Roary는 주석화된 유전체 서열 정보로부터 단백질 서열을 추출하여 CD-HIT로 사전 클러스터링을 실시한 뒤, BLASTP로 all-against-all comparison을 거친 MCL algorithm(Markov Cluster Algorithm)으로 클러스터를 만드는 과정으로 이어진다. '-e [--mafft]' 옵션을 주면 core gene에 대한 alignment까지 실시하므로 작업 시간이 오래 소요된다. 이러한 기본 흐름 외에도 여러 부속 유틸리티를 포함한다.

Roary는 동일 종에 속하는 미생물의 유전체 집합을 대상으로 하므로, 종 범주를 벗어나는 유전체가 섞여 있는 경우 core genome을 제대로 생성하지 못한다. 따라서 dRep 등으로 동일 종에 속하는 유전체를 미리 선별하여 사용하는 것이 바람직하다.

필요한 입력물은 유전체 및 유전자의 염기서열이 뒷부분에 포함된 GFF3 파일이다. Prokka로 만든 GFF3은 roary에 그대로 사용하면 된다. 그러나 일반적인 경우 GFF3 파일에는 염기서열 정보가 첨부되지 않으므로, NCBI에서 다운로드한 GenBank flat file을 BioPerl의 부속 스크립트인 bp_genbank2gff3.pl로 처리하여 새롭게 만드는 것이 더욱 바람직하다. BioPerl의 최신 패키지에서는 bp_genbank2gff3라는 이름으로 바뀌어 있을 것이다. Prokka로 만든 GenBank 파일을 bp_genbank2gff3로 처리하여 GFF3 파일을 만든 다음 roary에 사용하면 'Could not extract any protein sequences from … Does the file contain the assembly as well as the annotation?'에러가 발생한다. 로컬 환경에서 PGAP을 실행하였거나 NCBI에서 주석화하여 아직 공식 release가 되기 전의 GenBank 파일은 ACCESSION과 VERSION 정보가 비어 있는데, 이를 그대로 GFF3로 전환하면 roary가 제대로 작동하지 않는다. 따라서 ACCESSION은 GenBank 파일 첫 줄의 LOCUS와 같은 값으로, VERSION은 LOCUS.1과 같이 적당히 채워 넣도록 한다. 여러 서열로 구성된 경우에는 다소 번거로울 것이다.

genbank 디렉토리에 존재하는 모든 GenBank flat file을 GFF로 전환해 보자. INPUT.gbk에 대하여 INPUT.gbk.gff라는 결과 파일이 생긴다. 이를 INPUT.gff로 바꾸지 않아도 roary 실행에는 문제가 없다.

$ bp_genbank2gff3 -d genbank
$ bp_genbank2gff3 genbank/*gbk # 결과는 같다.

'-f outdir'로 출력 디렉토리를 따로 지정하지 않으면 현 디렉토리에 결과를 저장한다. '-o STR' 옵션은 출력 디렉토리가 아니라 clustered_protein 결과 파일을 다른 이름으로 지정하기 위한 것이다. Roary의 실행 시간을 줄이려면 다중 서열 정렬에 의한 core gene alignment('-e' or '-e --mafft') 를 생략하면 된다. 그러나 이렇게 실행하면 core_gene_alignment.aln 및 pan_genome_reference.fa 파일이 생기지 않아서 유전자 서열 관련 작업을 하기가 어렵다. 다음은 8개의 쓰레드를 사용하는 사례이다. -e 옵션을 주면 prank에 의한 codon-aware alignment를 실시하므로 mafft보다는 시간이 좀 더 걸린다.

$ roary -p 8 -f outdir /path/to/*.gff  # core gene alignment 생략
$ roary -p 8 -f outdir -e --mafft /path/to/*gff  # mafft에 의한 core gene alignment
$ roary -p 8 strain1.gff strain2.gff strain3.gff  # 일부 균주에 대해서만 분석

현 작업 디렉토리에 모든 GFF 파일이 존재하지만 이중에서 일부에 대해서만 분석을 실시하고 싶다면 위에 보인 명령행 사례 중 가장 마지막 방법처럼 대상 GFF 파일을 전부 공백으로 분리하여 roary의 인수로 제공하면 된다. 그러나 그 개수가 많아지면 이를 명령행에 일일이 타이핑하기가 곤란해 진다. Roary에 사용할 GFF 파일의 목록을 list.txt 파일 내에 갖고 있다면(한 줄에 하나의 GFF 파일), 다음과 같이 paste 명령어를 활용하면 된다.

$ roary -p 24 $(cat list.txt | paste -sd' ')
# $(paste –sd' ' < list.txt)라 해도 결과는 같다.
# $(paste –sd' ' list.txt)라 해도 결과는 같다.

Rory 실행의 첫 단계에서는 GFF 파일을 체크하여 수정된 파일을 fixed_input_files 디렉토리에 저장한다. 이때 나오는 메시지는 다음과 같다. 단, prokka로 직접 annotation하여 만든 GFF 파일은 특별히 수정되지 않는다.

2023/06/26 16:11:25 Input file contains duplicate gene IDs, attempting to fix by adding a unique suffix, new GFF in the fixed_input_files directory: /media/sf_Shared_Folder/ani_r_exercise/gff/GCF_000146875.3.gbk.gff 
2023/06/26 16:11:26 Input file contains duplicate gene IDs, attempting to fix by adding a unique suffix, new GFF in the fixed_input_files directory: /media/sf_Shared_Folder/ani_r_exercise/gff/GCF_000463565.1.gbk.gff 

수정 전후의 GFF 파일을 비교해 보자. 하나의 gene에 부속된 mRNA, CDS 및 exon feature가 서로 다른 ID를 갖게 되는 것이 흥미롭다. Roary는 이들 feature 중에서 CDS에 해당하는 것만을 취하여 분석 작업에 사용한다.

# 수정 전
NC_014483       GenBank gene    440     1786    .       +       1       ID=PPE_RS00005;Name=dnaA;locus_tag=PPE_RS00005;old_locus_tag=PPE_00001
NC_014483       GenBank mRNA    440     1786    .       +       1       ID=PPE_RS00005.t01;Parent=PPE_RS00005
NC_014483       GenBank CDS     440     1786    .       +       1       ID=PPE_RS00005.p01;Parent=PPE_RS00005.t01;gO_function=GO:0003677 - DNA binding [Evidence IEA],GO:0003688 - DNA replication origin binding [Evidence IEA],GO:0005524 - ATP binding [Evidence IEA];gO_process=GO:0006270 - DNA replication initiation [Evidence IEA],GO:0006275 - regulation of DNA replication [Evidence IEA];Name=dnaA;Note=Derived by automated computational analysis using gene prediction method: Protein Homology.;codon_start=1;inference=COORDINATES: similar to AA sequence:RefSeq:WP_014279007.1;locus_tag=PPE_RS00005;old_locus_tag=PPE_00001;product=chromosomal replication initiator protein DnaA;protein_id=WP_023986414.1;transl_table=11;translation=length.448
# 수정 후 - 각 genome을 뜻하는 긴 문자열이 생성되었다. 보기에 썩 아름답지는 않다.
NC_014483       GenBank gene    440     1786    .       +       1       ID=f15c2d87eb71223f47ba3d741ffa81b7_1;Name=dnaA;locus_tag=PPE_RS00005;old_locus_tag=PPE_00001
NC_014483       GenBank mRNA    440     1786    .       +       1       ID=f15c2d87eb71223f47ba3d741ffa81b7_2;Parent=PPE_RS00005
NC_014483       GenBank CDS     440     1786    .       +       1       ID=f15c2d87eb71223f47ba3d741ffa81b7_3;Parent=PPE_RS00005.t01;gO_function=GO:0003677 - DNA binding [Evidence IEA],GO:0003688 - DNA replication origin binding [Evidence IEA],GO:0005524 - ATP binding [Evidence IEA];gO_process=GO:0006270 - DNA replication initiation [Evidence IEA],GO:0006275 - regulation of DNA replication [Evidence IEA];Name=dnaA;Note=Derived by automated computational analysis using gene prediction method: Protein Homology.;codon_start=1;inference=COORDINATES: similar to AA sequence:RefSeq:WP_014279007.1;locus_tag=PPE_RS00005;old_locus_tag=PPE_00001;product=chromosomal replication initiator protein DnaA;protein_id=WP_023986414.1;transl_table=11;translation=length.44

이후에 작업에서는 모두 수정된 유전자 ID를 기준으로 정보가 다루어지므로 주의가 필요하다. 특히 roary가 끝난 뒤 query_pan_genome 스크립트를 이용할 경우, 인수로 제공해야 할 GFF 파일은 수정이 이루어진 GFF가 되어야 한다. 따라서 fixed_input_files에 존재하는 수정된 GFF 파일을 현 디렉토리로 모두 복사하여 기존의 것을 덮어 쓴 다음에 query_pan_genome 스크립트를 실행해야 한다.

Roary가 생성하는 파일은 Roary 웹사이트의 Output files에서 상세히 설명하였다. * 표시가 있는 3개의 파일은 '-e [--mafft]' 옵션과 함께 roary를 실행할 때에만 생성되는 것이다.

accessory.header.embl             *core_gene_alignment.aln
accessory.tab                     fixed_input_files(디렉토리)
accessory_binary_genes.fa         gene_presence_absence.Rtab
accessory_binary_genes.fa.newick  gene_presence_absence.csv
accessory_graph.dot               number_of_conserved_genes.Rtab
blast_identity_frequency.Rtab     number_of_genes_in_pan_genome.Rtab
clustered_proteins                number_of_new_genes.Rtab
core_accessory.header.embl        number_of_unique_genes.Rtab
core_accessory.tab                *pan_genome_reference.fa
core_accessory_graph.dot          summary_statistics.txt
*core_alignment_header.embl

정상적으로 프로그램이 종료되면 pan genome에 대한 수치 정보 요약은 다음과 같은 형식의 summary_statistics.txt 파일에 실린다. Core gene의 수가 형편없이 적다면 뭔가 잘못한 것이다.

Core genes	(99% <= strains <= 100%)	3018
Soft core genes	(95% <= strains < 99%)	424
Shell genes	(15% <= strains < 95%)	2315
Cloud genes	(0% <= strains < 15%)	7524
Total genes	(0% <= strains <= 100%)	13281

clustered_proteins 파일은 클러스터링된 모든 core와 accessory gene을 수록한 tab-delimited file이다. 첫 번째 컬럼이 클러스터 식별자이고 나머지 컬럼은 클러스터를 구성하는 모든 멤버 서열의 ID(corrected)이며, 컬럼1과 컬럼2 사이의 구분자는 스페이스이다. 실제 계산된 결과로부터 clustered_proteins 파일의 한 줄을 인용해 본다. Gene ID는 roary의 시작 단계에서 매우 건조한 문자열로 바뀌었다.

dnaA: a1a82563115b75beda45b943a4982aba_3      c89560a585264e8dfaed7ccb9edd3b24_19403  b468e412484598827c4d804e35e27cd8_7      ...(생략)

clustered_proteins 파일의 총 라인 수는 summary_statistics.txt 파일의 'Total genes' 숫자와 동일하다. 클러스터 식별자(ID)는 genes_presence_absence.csv(or .Rtab)의 'Gene' 컬럼 및 pan_genome_reference.fa의 description field에서도 공통적으로 쓰인다. 각 클러스터를 구성하는 서열에서 가장 많이 나타나는 gene name을 가져와서 클러스터 식별자로 쓰는 것이 원칙이지만 여의치 않은 경우 cluster_XXX의 일반명을 사용한다. 클러스터는 단일 멤버 유전자도 전부 포함하기 때문에 summary_statistics.txt 파일의 total genes와 그 숫자가 일치한다. Clustered_proteins 파일은 나중에 설명할 query_pan_genome 명령에도 GFF(fixed) 파일과 함께 쓰이게 된다.

Prank를 이용하여 느리지만 정확한 core gene의 codon-aware alignment를 하려면 -e 옵션을, mafft를 사용한 빠른 alignment를 하려면 -e --mafft 옵션을 사용한다. Core gene alignment를 생략하면 전체 계산 과정이 매우 빠르게 끝나지만 pan_genome_reference.fa 파일이 생성되지 않음에 유의하라. pan_genome_reference.fa는 pan genome(core 및 accessory) 클러스터 각각에 대한 대표 서열(nucleotide)을 수록하고 있다. clustered_proteins 파일에서 가장 앞에 위치한 유전자가 각 클러스터의 대표 서열에 해당한다. 염기서열의 ID는 수정된 것을 따르며, 바로 뒤의 description field에는 클러스터 번호정보가 나온다.

Roary가 생성하는 결과 파일을 이용하면 diagnostic plot 작성, pan-genome에 대한 그룹 단위의 연산, 유전차 서열 추출 등 다양한 추가 결과를 얻을 수 있으므로 상세한 설명은 웹사이트의 정보를 잘 참조하기 바란다. Core gene alignment를 이용한 계통수(phylogenetic tree)를 작성하려면 FastTree나 PhyML, https://cme.h-its.org/exelixis/web/software/raxml/}RAxML 등의 프로그램을 이용한다.

간단한 활용 사례부터 알아본다. query_pan_genome이 요구하는 인수는 clustered_proteins와 GFF 파일명이다. fixed_input_files 디렉토리 아래의 수정된 GFF 파일을 원래의 위치에 덮어 씌워야 함은 위에서 이미 설명하였다. 시각화 스크립트인 roary_plots.py와 roary_2svg.py는 자동으로 설치되지 않으므로 Roary GitHub 사이트의 contrib 디렉토리에서 가져와야 한다.

# 특정 유전자의 아미노산 서열 출력하기(pan_genome_results.GENE.fa 파일이 3개 생긴다)
$ query_pan_genome -a gene_multifasta -n dnaA,recF,gyrA *gff
# Core gene 정보 얻기(pan_genome_results 파일이 생긴다. FASTA file이 아니라 clustered_proteins의 부분집합임.)
$ query_pan_genome -a intersection *gff
# Accessory gene 정보 얻기
$ query_pan_genome -a complement *gff
$ roary_plots.py name_of_your_newick_tree_file.tre gene_presence_absence.csv
$ roary2svg.pl genes_presence_absence.csv > pan_genome.svg
# core gene alignment를 이용한 계통수 만들기
$ fasttree –nt –gtr core_gene_alignment.aln > core_gene_alignment.aln.newick

다음으로는 roary 결과물에서 여러 가지 유용한 정보를 추출하는 방법에 대하여 구체적인 사례와 함께 알아보도록 한다.

'query_pan_genome -a intersection *gff' 명령을 사용하면 균주의 ≥99%가 공유하는 core gene의 cluster ID를 알아낼 수 있다(결과 파일의 첫 컬럼을 취하여 콜론을 제거). GFF를 특별히 선택하지 않고 전체 유전체에 대하여 core gene 정보를 추출하려면 awk와 clustered_proteins 파일을 이용하여 같은 목적의 작업을 할 수 있다. 만약 5개의 GFF 파일을 이용하여 roary를 실행했다고 가정하자. 그러면 clustered_proteins 파일에서 컬럼의 수가 5 + 1 = 6인 라인을 추출하여 적절히 조작하면 된다. 공백과 탭문자가 섞여 있지만 awk는 이를 전부 컬럼 구분자로 알아서 인식한다.

$ awk 'NF==6{print $1}' clustered_proteins | sed 's/://' > final.txt
# 다음과 같이 명령해도 결과는 같다. 콜론(':')을 언제 제거하느냐의 문제이다.
$ awk 'NF==6{sub(":", "",$1); print $1}' clustered_proteins > final.txt

똑같은 목적의 작업을 아주 어렵게 하는 방법을 알아보자. roary 분석의 모든 결과는 gene_presence_absence.csv(또는 .Rtab)에 고스란히 들어있다. 따라서 clustered_proteins와 개별 GFF 파일에 의존하는 query_pan_genome을 사용하지 않고 대신 이 결과 파일을 엑셀이나 R에서 이를 읽어 들인 후 적절히 조작하여 조건을 만족하는 row의 첫 번째 컬럼(명칭은 'Gene'; 각 클러스터를 구성하는 서열에서 가장 흔하게 나타나는 gene name을 가져다 쓰며, 적절한 것이 없으면 grpup_XXX 형식의 식별자를 붙임)을 추출하면 된다. 다음의 사례는 R에서 gene_presenece_absence.Rtab 파일을 불러들여서 (1) core gene의 ID(클러스터 식별자)을 final.txt 파일에 출력하고, (2) 균주 특이적 유전자와 이들이 유래한 균주명(정확하게는 GFF 파일명에서 확장자를 제거한 것)을 unique_genes.txt 파일에 출력하는 것이다. R에서 데이터프레임을 다루는 좋은 연습이 되므로 상세하게 설명하였다. 이 R 스크립트에서는 모든 균주 유전체에 공통적으로 존재하는 유전자(즉 100% presence)를 core로 간주하였다.

> df = read.table("gene_presence_absence.Rtab",sep="\t",header=T,check=F)
# check=F는 행과 열의 라벨이 원본 파일에 있는 그대로 쓰이게 만든다. 
# 그렇지 않으면 라벨 내의 공백을 밑줄로 바꾸거나 숫자로 시작하는 라벨을 수정하게 된다.
> total.strains = ncol(df) - 1 # 전체 균주 수 확인
# unique gene에 대한 작업만 하려면 total.strains 값을 구하지 않아도 된다.
# core 또는 unique gene에 해당하는 row만을 골라서 별도의 데이터프레임에 저장하는 것이 최종 목표이다.
# 다음 명령은 core gene의 목록(Gene 컬럼)을 벡터로 추출한다. 맨 끝의 '1'을 생략하면 df 데이터프레임에서 조건을 만족하는 row를 반환하게 된다.
> core = df[apply(df[,c(2:ncol(df))],1,sum) == total.strains,1]
> core.row.df = df[apply(df[,c(2:ncol(df))],1,sum) == total.strains,]
# 다음 명령도 같은 결과를 반환한다.
# core.1 = df[which(rowSums(df[,c(2:ncol(df))])==total.strains),1]
# subset() 함수를 쓸 수도 있다. 컬럼을 하나로 제한하여도 반환되는 자료형은 데이터프레임이 된다. 
# core.df.1 = subset(df, select=Gene, subset=(rowSums(df[,c(2:ncol(df))])==total.strains))
# core.df.2 = subset(df, subset=(rowSums(df[,c(2:ncol(df))])==total.strains))
# Gene 목록을 파일로 저장한다.
> write.table(core,"final.txt",quote=F,row.names=F,col.names=F)
# df에서 균주 특이적(unique) 유전자 자료를 뽑아 uniq 데이터프레임에 저장하자.
# Gene ID(클러스터 번호)만으로는 이 유전자가 어느 균주의 것인지 파악하기 어려우므로
# unique 유전자에 해당하는 row를 추출하여 데이터프레임 형태를 유지하는 것이 좋다.
# 이번에는 subset() 함수를 사용해 본다.
> uniq = subset(df, subset=(rowSums(df[,c(2:ncol(df))])==1))
# 다음의 두 명령행도 결과는 같다.
# uniq = df[which(rowSums(df[,c(2:ncol(df))])==1),]
# uniq = df[apply(df[,c(2:ncol(df))],1,sum) == 1,]
# 각 유전자에 균주 번호를 부가하여 y 데이터프레임에 저장한다.
# 첫 번째 방법은 for 반복문을 쓰는 번거로움이 있지만, 각 균주별로 unique gene을
# 모아서 보여준다는 장접이 있다.
> y = data.frame(Gene = factor(), source = factor())
> strains = colnames(uniq[c(2:ncol(uniq))])
> for (i in c(1:length(strains))) {
> x=cbind(subset(uniq,select=Gene,subset=(uniq[[i+1]]==1)), source=strains[[i]])
> y=rbind(y,x)
> }
# 두 번째 방법에서는 균주 정보를 뽑아내어 uniq 데이터프레임의 새로운 컬럼(source)에
# 채워 넣은 뒤 Gene과 source 컬럼만을 추출하여 y에 저장하는 것이다.
> uniq$source = names(uniq)[apply(uniq, 1, function(i) which(i == 1))]
> y.1 = df.2[,c("Gene","source")]
# 세 번째 방법은 두 번째 방법과 유사하다.
> uniq$source = apply(uniq[,c(2:ncol(uniq))], 1, function(i) names(i[i == 1]))
> y.2 = uniq[,c("Gene","source")]
# 확인 후 파일로 저장
> View(y)
> write.table(y,"unique_genes.txt",sep="\t",quote=F,row.names=F,col.names=F)

이상의 방법으로 얻어낸 gene ID(=cluster ID)를 사용하여 pan_genome_reference.fa 파일로부터 각 클러스터의 대표 서열을 추출하려면 어떻게 하면 좋을까? Gene ID는 다름아닌 클러스터 번호이지만 pan_genome_reference.fa의 각 유전자 서열은 '>(개별_sequence_ID) (클러스터_번호)'의 헤더 형식을 갖고 있다. FASTA 파일의 서열 ID 자체가 아니라 description field를 참조하여 염기서열을 뽑아내야 하므로 일반적인 서열 조작용 유틸리티를 적용하기가 불편하다. 따라서 pan_genome_reference.fa를 tab-delimited file로 전환한 다음 description field에 해당하는 컬럼 값을 이용하여 원하는 join을 사용하여 서열을 추출하는 트릭을 이용한다. Gene ID(클러스터 번호)가 수록된 파일은 results.txt라는 이름으로 준비하도록 한다. 마지막 단계에서는 EMBOSS 패키지의 seqret 프로그램이 필요하다.탭문자를 sort 명령어의 필드 구분자로 제공하기 위해 백슬래시를 사용하여 -t $'\t' 형식으로 표현하는 것에 주의하라. 후속 분석 작업에서 유전자 서열이 직접적으로 필요하다면 일반적인 경우 pan_genome_reference.fa로 충분할 것이다.

$ awk '/^>/ {gsub(">", "", $1); printf("%s%s\t%s\t",(N>0?"\n":""),$1,$2); N++; next;} {printf("%s",$0);} END {printf("\n");}' pan_genome_reference.fa | \
sort -t $'\t' -k2,2 > 01_pan_genome_reference_sorted.tab
$ sort -t $'\t' -k1,1 results.txt > 02_resulted.txt.sorted
$ join -t $'\t' -1 2 -2 1 -o 1.1,1.2,1.3 01_pan_genome_reference_sorted.tab 02_resulted.txt.sorted > 03_extracted_sequences.tab
$ awk '{printf(">%s %s\n%s\n", $1, $2, $3)}' 03_extracted_sequences.tab | seqret fasta::stdin fasta:04_final.fa
# 서열 ID와 desc(클러스터 ID)의 순서를 바꾸고 싶다면 위 명령에서 $1과 $2의 순서를 서로 바꾼다.
$ awk '{printf(">%s %s\n%s\n", $2, $1, $3)}' 03_extracted_sequences.tab | seqret fasta::stdin fasta:04_final.fa

Strain 특이적 유전자를 추출하려면 query_pan_genome 명령을 사용할 수 있다. 1.gff, 2.gff,.. 5.gff를 이용하여 roary를 실행했다고 가정하자. 만약 1.gff만이 갖고 있는 유전자를 추출하고 싶다면 다음의 명령어를 실행하면 된다. 결과물은 두 세트 공통(set_difference_common), 첫번째 세트 특이적(set_difference_unique_set_one) 및 두번째 세트 특이적(set_difference_unique_set_two) 유전자로 구분되어 얻어진다. 세트 특이적 유전자라 해도 세트의 모든 멤버가 전부 공유하는 것으로만 이루어지는 것은 아님에 유의해야 한다. Query_pan_genome 실행 시에는 roary에 사용한 원본 GFF 파일이 아니라 ID가 수정되어 fixed_input_files 디렉토리에 보관된 파일을 가져다가 사용해야 함을 잊지 말자.

$ query_pan_genome -a difference --input_set_one 1,gff --input_set_two 2.gff,3.gff,4.gff,5.gff

Accessory gene을 뽑는 방법은 이미 설명하였다.

oary를 이용하여 특정 균주의 유전체로부터 core와 accessory gene을 서로 구분해 냈다고 가정하자. 이어지는 질문은 accessory gene에 어떠한 생물학적 기능 혹은 경로가 더 많은지를 알아보는 것이다. 또는 분석에 이용한 모든 균주의 pan genome으로부터 특정 균주 그룹만이 공유하는 유전자 세트를 발굴하였다고 가정하자. 이러한 유전자 서브셋이 갖고 있는 enriched function or pathway에는 무엇이 있을까? Over-Representation Analysis(ORA)는 가장 오랫동안 쓰여온 pathway analysis 도구로서, GO나 KEGG Pathway 등 기능 정보가 부가된 전체 유전자 집합에 대하여 이로부터 유래한 유전자 서브셋이 어떠한 기능을 더 많이 갖고 있는지를 통계적 테스트와 함께 추출해 준다. 여기에서는 WebGestalt(WEB-based GEne SeT AnaLysis Toolkit)의 사용 방법을 알아보고자 한다. 총 다섯 개의 균주 유전체를 비교하는 것으로 가정한다. 목표는 LMT16-62가 보유한 전체 유전자에 대하여 특이적인 유전자에 더욱 많이 존재하는 기능을 찾는 것이다.