Table of Contents
Pan-genome analysis
세균의 종을 구성하는 유전체는 모든 균주에게 공통적으로 존재하는 유전자 집합인 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 활용법
기존에 작성한 문서: Roary를 이용한 미생물 유전체의 비교 분석 ← 겹치는 내용이 많음!
기본적인 pan-genome analysis는 Roary 파이프라인을 통하여 진행할 수 있다. Roary는 주석화된 유전체 서열 정보로부터 단백질 서열을 추출하여 CD-HIT로 사전 클러스터링을 실시한 뒤, BLASTP로 all-against-all comparison을 거친 MCL algorithm(Markov Cluster Algorithm)으로 클러스터를 만드는 과정으로 이어진다. '-e [--mafft]' 옵션을 주면 core gene에 대한 alignment까지 실시하므로 작업 시간이 오래 소요된다. 이러한 기본 흐름 외에도 여러 부속 유틸리티를 포함한다.
입력 파일(GFF) 마련하기
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를 실행할 때에만 생성되는 것이다.
roary_test_outdir_20230627.zip ← 21개의 genome을 이용한 roary 결과 파일 묶음. 분석에 사용한 미생물 유전체의 RefSeq assembly accession은 내부의 00_accessions.txt 파일에서 확인하라.
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) 파일과 함께 쓰이게 된다. genes_presence_absence.csv 파일에는 각 genome(strain) 컬럼에 개별 유전자의 ID가 수록되어 있으나, genes_presence_absence.Rtab은 유전자의 유무에 따라 0 또는 1의 값이으로 바뀌어 있고 유전자 기능 등 보조적인 정보를 담는 컬럼은 전부 삭제되어 있다. 즉 “Gene” 컬럼과 각 균주 컬럼으로만 구성된다.
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 등의 프로그램을 이용한다.
Roary 결과물의 활용법
간단한 활용 사례부터 알아본다. 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 결과물에서 여러 가지 유용한 정보를 추출하는 방법에 대하여 구체적인 사례와 함께 알아보도록 한다.
[1] Core/균주 특이적 gene의 cluster ID 추출
'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)
[2] 각 cluster의 대표 서열 추출
이상의 방법으로 얻어낸 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
[3] 그룹 특이적 유전자 추출
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(pag-genome 중 core gene을 제외한 모든 것)을 뽑는 방법은 이미 설명하였다.
Over-Representatio Analysis(ORA)
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가 보유한 전체 유전자에 대하여 특이적인 유전자에 더욱 많이 존재하는 기능을 찾는 것이다.
Panseq을 이용한 균주 특이적 염기서열 추출
Panseq는 설치 및 사용 방법이 다소 난해하여 이를 사용하는 것을 적극 권장하지는 않는다. 그저 이런 프로그램도 있다는 정도로 이해하면 좋을 것이다.
비교유전체 분석을 통해서 특정 균주의 유전체에만 존재하는 염기서열을 추출할 필요가 종종 있다. 앞서 소개한 roary나 orthoMCL에서는 유전자 단위에서 모든 유전체에 공통으로 존재하는 core gene과 나머지에 해당하는 accessory gene을 추출해 주므로 후자를 활용하면 각 균주에만 존재하는 strain-specific gene을 찾을 수 있다. 그러나 이러한 특이적 영역을 반드시 유전자 단위로만 추출할 필요는 없다. 유전체 단위로 특이적인 영역을 찾아낸다면 qPCR 프라이머 등을 설계할 때 더욱 유용하게 활용할 수 있다. Panseq 프로그램을 novel region finder 모드로 실행하면 이러한 목적에 맞는 염기서열을 얻게 된다.
Panseq의 실행에 필요한 파라미터는 settings.txt(파일 이름은 중요하지 않음)에 미리 설정해 두어야 한다. Panseq는 query 파일을 일정 길이(fragmentationSize 파라미터로 변경 가능)로 잘라서 MUMmer 3 패키지의 nucmer를 사용하여 reference 염기서열과 비교한다.
유전체 서열 파일이 여러 contig로 구성된 경우에는 이것이 같은 균주에 속함을 명시적으로 나타내기 위하여 다음과 같이 서열 ID를 lcl|unique_identifier|sequence_id 포맷으로 바꾸어야 한다. 단일 서열로 구성된 경우는 수정할 필요는 없다.
>lcl|strain1|contig1 ATGCTTT… >lcl|strain1|contig2 ATGCTTT…
그러나 실제 상황에서는 서열 ID 뒤에 description 정보가 있으면 제대로 결과가 나오지 않는 것을 발견하였다. 따라서 서열 ID 라인을 일괄적으로 정리하는 작업이 필요하다. 염기서열 파일이 field1_field2_field3.fna일 때 field1을 새로 만들어질 서열 ID의 unique_identifier로 사용함과 동시에 description을 삭제하여 최종적으로 field1.fasta 파일에 저장하려면 다음과 같이 실행한다.
$ ls *fna | while read f
> do
> DATA=$(cut -d '_' -f 1,2 <<<$f)
> awk -v -s="lcl|${DATA}" '$1~/^>/{sub (/>/, ""); printf ">%s|%s\n", s, $1; next}{print}' $f > ${DATA}.fasta
> done
$ rm *fna
Panseq를 실행하려면 텍스트 파일(settings.txt)에 파라미터 값들을 먼저 저장해 두어야 한다. 다음은 novel region finding 모드로 실행할 때 필요한 최소한도의 settings.txt 파일 사례이다. 이를 참조하되 /usr/local/apps/Panseq/settings.txt 파일을 이용하여 적절히 변경하도록 한다. 특이적 영역을 추출할 유전체 염기서열 파일은 queryDirectory에, 비교를 할 나머지 유전체 파일은 referenceDirectory에 넣어둔다. MUMmer와 BLAST+ 및 Muscle은 PATH 환경변수에 저장되어 있어야 한다.
# 모든 디렉토리는 '/'로 끝나야 한다. queryDirectory <query sequence가 있는 디렉토리> referenceDirectory <reference sequence가 있는 디렉토리> baseDirectory <출력 파일을 저장할 디렉토리> numberOfCores 8 minimumNovelRegionSize 500 novelRegionFinderMode no_duplicates # or unique fragmentationSize 500 percentIdentityCutoff 95 runMode novel overwrite 1
Novel region finder는 “no_duplicates”와 “unique”의 두 가지 모드로 작동한다. no_duplicates 모드에서는 query로부터 reference에는 없는 영역을 찾아서 multi-fasta format으로 출력하는 반면, unique 모드에서는 각 query에는 고유하되 reference에는 없는 영역을 출력한다. 실행 방법은 다음과 같으며, 결과 파일은 refereneceDirectory/novelRegions.fasta이다.
$ /usr/local/apps/Panseq/lib/panseq settings.txt
Core/Accessary analysis 또는 loci selector 기능에 대해서는 매뉴얼을 참조하라.
기타 pan-genome 분석 도구
2018년에 논문으로 발표된 panX는 미생물 pan genome의 분석 및 시각화를 위한 도구이다. panX는 DIAMOND를 이용하여 단백질 서열의 빠른 정렬을 한 뒤 MCL로 클러스터링을 실시하며, 계통발생에 근거한 후처리를 실시한다. 더불어 core-genome의 SNP를 기반으로 트리를 작성하여 gene gain/loss도 추정한다. 가장 큰 특징은 시각화를 위한 별도의 프로그램(pan-genome-visualization)을 제공한다는 것이다. 사전에 계산된 pan genome은 panX 웹사이트에서 열람할 수 있다. 여기에서는 panX를 이용하여 직접 pan genome analysis를 실시하고 그 결과물을 시각화하는 방법을 알아본다.
입력파일은 RefSeq에서 입수하거나 사용자가 직접 만든 GenBank 파일로서 확장자는 gbk여야 한다. GenBank 파일을 data/TestSet 디렉토리에 모은 뒤 다음을 실행한다.
$ panX.py –fn data/TestSet –sl TestSet –t 24 > TestSet.log 2> TestSet.err
유전자의 클러스터링 결과물은 ./data/TestSet/allCluster_final.tsv을 비롯한 여러 파일로 저장되는데, 이를 직접 다루는 것보다는 pan-genome-visualization을 이용해서 웹브라우저에서 열람하고 그 환경 안에서 필요한 파일을 다운로드하는 것이 훨씬 바람직하다. 설치 및 데모용 결과 열람은 다음과 실행하면 된다. 설치 위치는 특별히 정해진 곳은 없지만 여기에서는 홈디렉토리에서 실시함을 가정하자. 아래의 명령어는 conda panX 환경과 무관하게 작동한다. 자바스크립트 런타임 환경인 Node.js와 소프트웨어 저장소 및 패키지 관리자인 npm은 관리자 권한으로 설치해 두어야 한다.
$ git clone https://github.com/neherlab/pan-genome-visualization $ cd pan-genome-visualization $ git submodule update --init $ npm install $ npm start
리눅스 명령 프롬프트에서 firefox를 실행한 뒤 http://localhost:8000/을 접속하거나. 혹은 8000번 포트로 tube server에 접속하면 된다(http://192.3.4.201:8000/). 사용자가 계산한 결과를 로드하려면 우선 nmp start를 실행했던 터미널에서 Ctrl+C를 입력하여 pan-genome-visualization을 중단한 뒤 panX.py를 실행했던 위치로 이동하여 다음을 실행한다(data 디렉토리가 있는 곳).
$ link-to-server.py –s TestSet –v ~/pan-genome-visualization $ cd ~/pan-genome-analysis $ npm start # http://localhost:8000/TestSet으로 접속
link-to-server.py를 실행하면 data/TestSet 디렉토리의 사본이 ~/pan-genome-visualization/public/dataset에 만들어지므로 npm을 종료한 뒤 재시작해도 그대로 접속 가능하다.
뒤에서 소개할 LS-BSR도 pan genome의 분석이 가능하다.
Scoary: pan-genome-wide association study
Genome-wide association study(GWAS)는 주로 인간의 SNP와 형질의 연관성을 알아보기 위한 대규모 집단 유전학 및 유전체학의 응용 분야로 알려져 있다. 박테리아의 GWAS, 즉 주로 임상적으로 중요한 형질(병원성이나 항생제 내성 등)과 SNP의 연관성을 탐구하는 연구는 매우 최근에 들어서 시작되었다.
Scoary는 앞서 설명한 Roary의 pan-genome analysis 결과물(“gene_presence_absence.csv” 파일)을 입력물로 하여 두 가지 카테고리로 나눌 수 있는 형질과 연관된 유전자를 찾아주는 프로그램이다. 앞서 소개한 Roary/Over-representation analysis(ORA)에서는 기능 증가와 연관된 유전자, 즉 그 유전자가 있음으로 인해서 enrich된 기능을 찾는 데에 초점을 맞추지만, Scoary에서는 비교 대상 게놈을 단지 binary category에 대해서 나누어 놓은 뒤 분석을 시작하므로 연구자가 관심을 갖는 특성은 예측된 유전자의 ‘없음’에서 비롯된 것일 수도 있다. 예를 들어서 항생제에 대한 세균의 내성은 항생제를 무력화하는 유전자가 있어서 생기는 일일 수도 있고, 항생제의 표적에 해당하는 유전자가 없어서 생기는 일일 수도 있다. Scoary라는 이름은 Roary에 대한 오마주이다.
Scoary의 입력물로는 gene_presence_absence.csv 파일과 더불어 trait table 파일(traits.csv)이 필요하다. Trait의 유무는 0과 1로 표현하되, 분석할 특성이 여러 개라면 콤마로 구분한다. 첫 컬럼은 균주명이고, 첫 줄은 제목에 해당한다. 다음은 trait table의 구조와 실행 사례이다.
읽어야 할 논문
Pandora: nucleotide-resolution bacterial pan-genomics with reference graphs Genome Biology (2021)