Genome Informatics Laboratory at KRIBB

User Tools

Site Tools

Trace:
bioinfo:average_nucleotide_identity_ani_의_계산

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
bioinfo:average_nucleotide_identity_ani_의_계산 [2023/06/26 12:55] – [dRep을 이용한 빠른 ANI 계산] hyjeongbioinfo:average_nucleotide_identity_ani_의_계산 [2023/08/11 16:48] (current) – [쌍 형태의 ANI 자료를 매트릭스로 전환하기] hyjeong
Line 21: Line 21:
 [[https://github.com/MrOlm/drep|dRep]]은 많은 수의 세균 유전체 서열을 빠르게 비교하는 도구이다. dRep 'compare' 명령의 첫 단계에서는 유전체 거리를 산출하는 신속한 방법([[https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0997-x|Mash]]: MinHash distance를 이용)을 이용하여 primary cluster를 만들고, 두 번째 단계에서는 ANI를 계산하는 느리지만 정확한 방법([[https://pubmed.ncbi.nlm.nih.gov/26150420/|gANI]] or nucmer 등, [[https://drep.readthedocs.io/en/latest/choosing_parameters.html?highlight=ANIn#overview-of-genome-comparison-algorithms|선택할 수 있는 전체 알고리즘]])을 사용하여 대략적으로 species 단위에 해당하는 secondary clustering을 실시한다. Primary clustering에만 의존할 수 없는 이유는 genome completeness가 낮은 상황에서는 Mash가 잘 대처하지 못하기 때문이다. 원래 dRep은 연관성이 있는 개별 메타게놈 데이터를 독립적으로 조립한 다음 이를 서로 비교하여 거의 동일한 게놈끼리 것끼리 그룹을 지어서 대표적인 것만을 추려내는(de-replication; dRep dereplicate) 용도로 개발된 것이다. 그러나 dRep compare 명령을 사용하여 종 수준의 클러스터링 결과를 PDF 보고서로 얻을 수 있다. 뿐만 아니라 pair-wise ANI 값도 결과로 제공된다. 이를 매트릭스로 전환하려면 R에서 추가 작업을 해야 한다.  [[https://github.com/MrOlm/drep|dRep]]은 많은 수의 세균 유전체 서열을 빠르게 비교하는 도구이다. dRep 'compare' 명령의 첫 단계에서는 유전체 거리를 산출하는 신속한 방법([[https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0997-x|Mash]]: MinHash distance를 이용)을 이용하여 primary cluster를 만들고, 두 번째 단계에서는 ANI를 계산하는 느리지만 정확한 방법([[https://pubmed.ncbi.nlm.nih.gov/26150420/|gANI]] or nucmer 등, [[https://drep.readthedocs.io/en/latest/choosing_parameters.html?highlight=ANIn#overview-of-genome-comparison-algorithms|선택할 수 있는 전체 알고리즘]])을 사용하여 대략적으로 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를 계산한다. 따라서 dRep의 pairwise ANI를 매트릭스로 전환하면 서로 다른 클러스터 사이에서는 ANI가 정의되지 않는다.+다음의 그림은 서로 다른 //Clostridium// 종의 genome을 한꺼번에 클러스터링한 결과를 비교한 것이다. Pyani(왼쪽)는 모든 유전체쌍에 대하여 ANI를 계산하지만 dRep(오른쪽)은 동일 primary cluster에 속하는 것 안에서만 ANI(secondary clustering)를 계산한다. 따라서 dRep의 pairwise ANI를 매트릭스로 전환하면 서로 다른 클러스터 사이에서는 ANI가 정의되지 않는다.
  
 {{ :bioinfo:pyani_vs_drep.png?400 |}} {{ :bioinfo:pyani_vs_drep.png?400 |}}
Line 53: Line 53:
   $ awk -F"," -v OFS="\t" 'NR!=1{print $1, $2, $4}' Ndb.csv > pairwise.txt   $ awk -F"," -v OFS="\t" 'NR!=1{print $1, $2, $4}' Ndb.csv > pairwise.txt
      
-실습을 통하여 dRep 계산과 MASH matrix를 만드는 방법을 알아보자. {{:bioinfo:ani_r_exercise.zip|}}을 다운로드하여 압축을 풀고 accessions.txt 파일에 수록된 정보를 이용하여 51개의 유전체 파일을 다운로드하고 압축해제 및 파일명 단순화를 한 뒤 dRep v3.4.3을 실행하고 이에 따른 후속 작업을 하는 명령어를 다음에 보였다.+실습을 통하여 dRep 계산과 MASH matrix를 만드는 방법을 알아보자. {{:bioinfo:ani_r_exercise.zip|}}을 다운로드하여 압축을 푼다. accessions.txt 파일에 수록된 정보를 이용하여 //Paenibacillus polymyxa// 및 관련 species에 속하는 51개 균주의 유전체 파일을 다운로드하고 압축해제 및 파일명 단순화를 한 뒤 dRep v3.4.3을 실행한 다음, 이에 따른 후속 작업을 하는 명령어를 다음에 보였다. 특별히 지정하지 않으면 secondary clustering algorithm으로는 [[https://github.com/ParBLiSS/FastANI|fastANI]]가 쓰인다.
  
   # accessions.txt의 각 줄을 분리하여 acc1,acc2,acc3,,,accN 형태의 문자열로 만드는 트릭을 눈여겨보라.   # accessions.txt의 각 줄을 분리하여 acc1,acc2,acc3,,,accN 형태의 문자열로 만드는 트릭을 눈여겨보라.
Line 91: Line 91:
      
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
-   
      
       ..:: dRep compare finished ::..       ..:: dRep compare finished ::..
Line 98: Line 97:
   Figures.............................. /media/sf_Shared_Folder/ani_r_exercise/drep_outdir/figures/   Figures.............................. /media/sf_Shared_Folder/ani_r_exercise/drep_outdir/figures/
   Warnings............................. /media/sf_Shared_Folder/ani_r_exercise/drep_outdir/log/warnings.txt   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 ===== ===== Genome-to-Genome Distance Calculator =====
bioinfo/average_nucleotide_identity_ani_의_계산.1687751749.txt.gz · Last modified: by hyjeong