Table of Contents

Average Nucleotide Identity(ANI)의 계산

온라인에서 ANI 계산을 할 수 있는 곳으로 잘 알려진 것은 Kostas lab의 ANI calculator와 EzBioCloud(구 천랩)의 ANI Calculator가 있다. 후자의 프로그램은 BLAST가 아닌 USEARCH 기반의 OrthoANIu 알고리즘을 사용하여 상동 유전체 절편 사이의 평균 identity를 구하는 것이 특징이다.

명령행 환경에서 두 개의 유전체 서열(FASTA 파일) 사이의 ANI를 간단히 계산하려면 ANI.pl Perl 스크립트를 사용한다. 이를 사용하려면 legacy BLAST에 포함된 blastall과 formatdb의 전체 경로를 옵션으로 지정해야 한다. 만약 두 실행파일이 PATH 환경변수에 포함되어 있다면 프로그램의 이름만 제공하면 된다. 

$ mkdir result # 사전에 결과 디렉토리를 만들어야 함
$ ANI.pl –bl blastall –fd formatdb –qr genome_1.fa –sb genome_2.fa –od result 
/path/to/genome_1.fa VS /path/to/genome_2.fa
ANI: 65.6929759299781

다수의 유전체 염기서열 파일에 대한 ANI 매트릭스 계산은 JSpecies가 그 원조라 할 수 있다. EzBioCloud에서 제공하는 독립형 JAVA 프로그램인 OrthoANIu tool도 유용할 것이다.

pyani

Pyani는 복수의 유전체 서열 파일에 대하여 ANI 및 tetranucleotide frequency 계산과 더불어 heatmap 이미지를 생성해 준다. 사용 가능한 쓰레드를 전부 동원하여 계산을 하므로 실행 시간은 매우 빠르다. Stable version(0.2.x)는 여전히 legacy blast와 MUMmer를 alignment engine으로 사용한다. 

$ average_nucleotide_identity.py –o OUTDIR –m ANIb –g –i FASTA_DIR

결과 디렉토리(OUTDIR)에 만들어지는 ANIb_percent_identity.tab 파일은 R에서 읽어 들여서 hierarchical clustering 등 추가적인 분석이나 조작을 할 수 있다. Row와 column은 각각 독립적으로 정렬된 상태라서 대각선에 대하여 대칭을 이루고 있지 않으니 주의해야 한다. Percent identity, alignment length 및 alignment coverage를 이용한 heatmap 그래픽 파일을 제공해 주는데, 문서에 따르면 흔히 쓰이는 average/UPGMA(unweight pair group method with arithmetic mean) 방법이 아닌 single-linkage clustering에 의한 결과를 보여주는 것으로 여겨진다.

dRep을 이용한 빠른 ANI 계산

dRep은 많은 수의 세균 유전체 서열을 빠르게 비교하는 도구이다. dRep 'compare' 명령의 첫 단계에서는 유전체 거리를 산출하는 신속한 방법(Mash: MinHash distance를 이용)을 이용하여 primary cluster를 만들고, 두 번째 단계에서는 ANI를 계산하는 느리지만 정확한 방법(gANI or nucmer 등, 선택할 수 있는 전체 알고리즘)을 사용하여 대략적으로 species 단위에 해당하는 secondary clustering을 실시한다. Primary clustering에만 의존할 수 없는 이유는 genome completeness가 낮은 상황에서는 Mash가 잘 대처하지 못하기 때문이다. 원래 dRep은 연관성이 있는 개별 메타게놈 데이터를 독립적으로 조립한 다음 이를 서로 비교하여 거의 동일한 게놈끼리 것끼리 그룹을 지어서 대표적인 것만을 추려내는(de-replication; dRep dereplicate) 용도로 개발된 것이다. 그러나 dRep compare 명령을 사용하여 종 수준의 클러스터링 결과를 PDF 보고서로 얻을 수 있다. 뿐만 아니라 pair-wise ANI 값도 결과로 제공된다. 이를 매트릭스로 전환하려면 R에서 추가 작업을 해야 한다.

다음의 그림은 서로 다른 Clostridium 종의 genome을 한꺼번에 클러스터링한 결과를 비교한 것이다. Pyani(왼쪽)는 모든 유전체쌍에 대하여 ANI를 계산하지만 dRep(오른쪽)은 동일 primary cluster에 속하는 것 안에서만 ANI(secondary clustering)를 계산한다. 따라서 dRep의 pairwise ANI를 매트릭스로 전환하면 서로 다른 클러스터 사이에서는 ANI가 정의되지 않는다.

gANI(웹 버전 및 다운로드)는 별도로 설치해야 한다. 실행파일 이름은 ANIcalculator이다. 이 프로그램은 오픈소스가 아니라서 dRep의 저자는 goANI라는 이름으로 이를 따로 구현했다고 한다.

dRep의 결과물을 저장할 OUTDIR은 compare 또는 dereplicate 명령어 바로 뒤에 지정해야 한다. 특별히 설정하지 않으면 6개의 프로세서를 동시에 사용한다. 다음의 실행 사례는 두 번째 클러스터링(–S_algorithm)에서는 gANI를 사용하되 0.965를 기준으로 종 클러스터를 나누겠다는 의미이다. 두 번째 클러스터링 알고리즘의 기본 설정은 ANImf(align whole genomes with nucmer; filter alignment; compare aligned region)이며, –S_ani의 기본 수치는 0.99이다. gANI(whole-genome based ANI)는 유전체 서열이 아니라 상동 유전자 사이의 ANI 수치(최소 96.5%)와 AF(aligned fraction; 0.6 이상)을 감안하여 종을 구분하는 방법에서 제안된 것이다. 

$ dRep compare –f # for help
$ dRep compare OUTDIR –p 16 --S_ani 0.965 --S_algorithm gANI –g path/to/genomes/*fasta

쌍 형태의 ANI 자료를 매트릭스로 전환하기

dRep 결과에 포함되는 data_tables 디렉토리에는 추가 분석에 쓸 수 있는 텍스트 파일이 들어있다. 

./data_tables/Bdb.csv  # Sequence locations and filenames
./data_tables/Cdb.csv  # Genomes and cluster designations
./data_tables/Chdb.csv # CheckM results for Bdb
./data_tables/Mdb.csv  # Raw results of MASH comparisons
./data_tables/Ndb.csv  # Raw results of ANIn comparisons
./data_tables/Sdb.csv  # Scoring information
./data_tables/Wdb.csv  # Winning genomes
./data_tables/Widb.csv # Winning genomes' checkM information

예를 들어 Ndv.csv에는 다음의 정보가 수록된다. Secondary clustering algorithm에 따라서 ANI 정보를 담고 있는 컬럼의 위치가 달라지므로 유의해야 한다. 

query,reference,ani,alignment_coverage,primary_cluster # gANI 
query,reference,alignment_length,similarity_errors,ref_coverage,querry_coverage,ani,reference_length,querry_length,alignment_coverage,primary_cluster # ANImf

이를 R 환경(R-project)에서 적절히 조작하여 매트릭스 형태로 전환하면 계층적 클러스터링(hierarchical clustering), heatmap 작성, iTOL에서 쓰이는 tree annotation용 데이터셋 작성 등에 유용하게 쓸 수 있다. 먼저 Ndb.csv의 에서 첫 줄은 제외하고 나머지로부터 query, reference 및 ANI에 해당하는 값을 추출하여 pairwise.txt 파일에 저장한다. 만약 primary cluster 형성에 사용한 pairwise mash value를 사용하고 싶다면 Mdb.csv에서 컬럼 1(genome1), 컬럼 2(genome2) 및 컬럼 4(similarity)를 추출하면 된다. 

$ awk -F"," -v OFS="\t" 'NR!=1{print $1, $2, $4}' Ndb.csv > pairwise.txt

실습을 통하여 dRep 계산과 MASH matrix를 만드는 방법을 알아보자. ani_r_exercise.zip을 다운로드하여 압축을 푼다. accessions.txt 파일에 수록된 정보를 이용하여 Paenibacillus polymyxa 및 관련 species에 속하는 51개 균주의 유전체 파일을 다운로드하고 압축해제 및 파일명 단순화를 한 뒤 dRep v3.4.3을 실행한 다음, 이에 따른 후속 작업을 하는 명령어를 다음에 보였다. 특별히 지정하지 않으면 secondary clustering algorithm으로는 fastANI가 쓰인다. 

# accessions.txt의 각 줄을 분리하여 acc1,acc2,acc3,,,accN 형태의 문자열로 만드는 트릭을 눈여겨보라.
$ ncbi-genome-download -F fasta -A $(paste -sd, accessions.txt) --flat-output -o genomes bacteria 
$ cd genomes; gzip -d *gz
$ $ for i in $(ls *fna)
> do
> mv $i $(cut -d_ -f1,2 <<<$i).fna
> done
$ cd ..
$ dRep compare drep_outdir -g genomes/*fna
***************************************************
    ..:: dRep compare Step 1. Cluster ::..
***************************************************
    
Running primary clustering
Running pair-wise MASH clustering
3 primary clusters made
Running secondary clustering
Running 1883 fastANI comparisons- should take ~ 13.0 min
Step 4. Return output
***************************************************
    ..:: dRep compare Step 2. Evaluate ::..
***************************************************
    
will provide warnings about clusters
2 warnings generated: saved to /media/sf_Shared_Folder/ani_r_exercise/drep_outdir/log/warnings.txt
***************************************************
    ..:: dRep compare Step 3. Analyze ::..
***************************************************
    
making plots 1, 2, 3, 4
Plotting primary dendrogram
Plotting secondary dendrograms
Plotting MDS plot
Plotting scatterplots

$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

    ..:: dRep compare finished ::..

Genome comparison data............... /media/sf_Shared_Folder/ani_r_exercise/drep_outdir/data_tables/
Figures.............................. /media/sf_Shared_Folder/ani_r_exercise/drep_outdir/figures/
Warnings............................. /media/sf_Shared_Folder/ani_r_exercise/drep_outdir/log/warnings.txt

$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

$ awk -F"," -v OFS="\t" 'NR!=1{print $1, $2, $4}' drep_outdir/data_tables/Mdb.csv > pairwise.txt

# 데이터 파일(pairwise.txt)이 만들어졌으므로 R console에서 후속 작업을 진행한다.
> library(reshape)
> d = read.table(file="pairwise.txt",sep="\t")
> head(d) # 데이터 확인
> d2 = cast(d,V1~V2)
> row.names(d2) = d2[,1]
> d2 = d2[,-1]
# NA value(결측치)가 있는지 확인하라. Secondary clustering 결과를 사용한다면 결측치가 있을 것이다.
> sum(is.na(d2))
# NA를 0으로 대체한다. ANI 값이라면 이렇게 하는 것이 타당할 것이다.
> d2[is.na(d2)] = 0
# 대칭 여부를 확인한다.
> isSymmetric(as.matrix(d2))
[1] TRUE
# 대칭이 아니라면 다음과 같이 조작하는 것이 나중의 분석을 위해 바람직할 것이다.
> d3 = (d2 + t(d2))/2
> dim(d3)
> View(d3)
> write.table(d3,"dRep_ANI_matrix.txt",sep="\t")

drep_outdir/data_tables/Cdb.csv 파일로부터 각 secondary cluster에 몇 개의 genome이 속하는지 알아보자. 

$ awk -F, 'NR!=1{print $2}' drep_outdir/data_tables/Cdb.csv | sort | uniq -c 
     21 1_1
      2 1_2
      9 1_3
     11 1_4
      1 2_1
      1 2_2
      1 2_3
      3 3_1
      1 3_2
      1 3_3

가장 멤버의 수가 많은 cluster 1_1에는 21개의 genome이 속한다.

모든 유전체에 대한 pairwise ANI 수치를 빠짐없이 얻고 싶다면 pyani를 이용하는 것이 나을 것이다. Pyani의 결과물에는 ANI 매트릭스가 포함되어 있어서 추가적인 계산은 필요하지 않다. 단, 계산 시간이 더 많이 걸리는 것은 감안해야 한다. dRep의 결과 파일을 조작하여 매트릭스를 만드는 이유는 iTOL 서버 등에서 재사용할 수 있는 tree 파일을 얻기 위함이다. 클러스터링 정보는 이미 PDF로 작성된 리포트 파일로 제공하고 있으니 그것을 활용하면 된다. 반면 pyani는 ANI 매트릭스와 heatmap을 제공하지만 정작 어느 유전체가 어느 클러스터에 속하는지를 나타낸 별도의 파일을 만들어내지는 않는다.

Genome-to-Genome Distance Calculator

DSMZ의 Genome-to-Genome Distance Calculator(GGDC)는 유전체 서열을 이용하여 in silico DNA-DNA hybridization(DDH) 수치를 추정해 준다. ANI는 원래 DDH 실험을 대신하기 위해 고안된 것임을 상기하라. Query 및 reference의 genome FASTA file을 로컬 머신에서 업로드하는 것 외에도 accession을 입력할 수 있다는 것이 특징이다. 올바른 accession number의 입력 방법은 GGDC FAQ의 How should GenBank accession numbers be specified?와 'How can I check whether or not my accession number is usable GenBank accession number?를 읽어보기 바란다.

GGDC는 standalone version이 존재하지 않고, ANI 계산보다 시간이 오래 걸린다.