iTOL 그림 새로 그리기 - 최적의 genus delineation 방법은?

녹슬고 무뎌진 칼날 갈기! 숫돌이 어디로 갔나...

2019년에 Bergey's Manual of Systematic Bacteriology에 공동 저자로서 새로운 속(genus)에 관한 챕터를 투고한 일이 있었다. 리뷰 과정이 도대체 어떻게 돌아가는지 답이 없다가 무려 4년이 지난 2023년이 되어서 리뷰 의견이 왔다. 이 원고를 들고 있던 편집자가 은퇴를 한 것이 지연의 사유였다. 이건 해도 너무 하는 것 아닌가! 그 사이에 이 속에 해당되는 새로운 종이 보고된 것도 있으니 이를 추가하여 다시 phylogenomic analysis에 대한 그림을 만들어야 했다. 작년까지 iTOL을 유료 회원으로 매우 편리하게 사용해 왔었지만 외부로 파견 근무를 나오면서 구독 갱신을 하지 못하여 free access로 전환되었다. Tree 업로드는 되지만 annotation 정보는 저장이 되지 않는 불편함을 감수해야 한다. 그래서 대부분의 설정 사항을 파일로 잘 저장해서 필요한 때마다 drag & drop으로 밀어 넣어 해결하기로 했다. 

Genus level의 구별에 대해서는 별로 심각하게 생각해 본 일이 없었다. 그저 POCP(percentage of conserved proteins) 분석 정도를 통해서 대략적인 genus 경계를 제시했었다. 그러나 이런 목적이라면 AAI(amino acid identity)가 더 낫다고 하니 분석을 새로 실시해야 한다. Kostas lab의 AAI calculator는 여러 샘플을 한꺼번에 다루기가 매우 불편하다. 더 나은 방법이 없을까? 

AAI는 68% 정도를 genus 구별의 threshold로 여기면 된다고 한다. POCP는 대략 50%가 기준인 것으로 알려져 있으나 그 폭이 너무 크다. 자료 출처: MyTaxa: an advanced taxonomic classifier for genomic and metagenomic sequences (2014).

지금은 CJ바이오사인스의 대표가 된 천종식 교수께서 서울대 재직 당시 개발했던 EzAAI pipeline을 쓰기로 했다. Gene prediction, identity calculation, clustering(newick format file 생성)까지 한 번에 해 주니 이보다 편할 수는 없다. 일단 글 쓰다가 천 대표께 문자 한번 보내고...

 EzAAI와 ezTree를 써서 트리 파일을 각각 만든 다음, 이를 iTOL에 먼저 업로드한 뒤 annotation file을 순차적으로 올려서 그림을 그려 보았다. iTOL 데이터 파일을 어떻게 만들더라? 단편적으로 기록해 둔 메모를 찾아서 적용하되 오랜만에 Perl 코드 조각을 새로 쓰기도 하였다. 두 가지 방법으로 만든 tree의 토폴로지 상에는 크게 다른 점이 없다.

해상도는 의도적으로 낮게...

Genome으로부터 계산되는 여러 수치를 다양한 계층의 택소노미 구별에 쓸 수 있는데, 미생물 그룹에 따라서 보편적인 값을 적용하기 어렵다는 것이 큰 문제이다. Species 수준에서는 그나마 ANI와 dDDH가 "Gold Standard"로 인식되면서 꽤 좁은 threshold를 적용할 수 있지만, 연구자에게 중요한 genus가 특히 어렵다. Family보다 더 높은 단계의 taxnomy 레벨은 연구자의 가슴에 별로 와 닿지 않는다. 이런 문제에 천착한 논문으로는 Next-generation systematics: an innovative approach to resolve the structure of complex prokaryotic taxa(2016) 같은 것이 있다. 최근 논문을 찾다가 해양과기원 권개경 박사팀에서 작년에 Frontiers in Microbiology에 발표한 논문을 발견하였다. 요즘도 멋진 헤어스타일 유지하고 계신지 모르겠다.

Establishment of genomie base criteria for classification of the family Desulfovibrionaceae and proposal of two novel genera, Alkalidesulfovibrio gen. nov. and Salidesulfovibrio gen. no. (2022)

중요한 문장 몇 개를 발췌해 보았다.

• POCP와 AAI는 genus delineation에 유용하다.
• Genus delineation을 위한 AAI cutoff는 60-80%으로서 비교적 폭이 넓지만 POCP(대략 50%)보다는 낫다고 여겨진다.
• POCP = 50%를 단일 threshold로 하면 monophyly rule for taxon delineation을 위배할 수 있다(참조: The species concept for prokaryotes, FEMS Microbiology Reviews 2001. 이 논문은 genome sequencing이 보편화되기 전에 나온 것이다. 즉, species demarcation을 위한 ANI 개념이 나오기 전의 논문이다.).
• 클러스터링 작업을 자동으로 반복하고 그 결과를 평가하여 60~80% AAI 사이의 최적 임계값을 결정하는 파이썬 스크립트를 제작하였다.

Supplementary Figure 2(링크의 data sheet 2).

필요한 부분만 빨리 발췌하여 읽은 것이기에 중요한 메시지를 혹시 놓치지는 않았나 모르겠다. 다음은 methanotroph와 관련한 다른 논문(2018년도, 링크)에서 인용한 것.

Another tool proposed for delimiting taxonomic ranks at the genus level is AAI (Konstantinidis and Tiedje, 2005). However, prokaryotic taxa exhibit a range of AAI values, making distinct boundaries difficult to define (Konstantinidis and Tiedje, 2005; Luo et al., 2014). AAI comparisons by Luo et al. (2014) of related but different genera typically ranged from 60 to 80%.

동일 family에 속하는 서로 다른 genus에 대하여 서로 다른 AAI threshold를 적용해야 한다는 뜻이렸다! 비교적 최근에 나온 참고할 논문을 하나 더 소개하는 것으로 오늘 글을 마치고자 한다.

A genus definitio for Bacteria and Archaea based on a standard genome relatedness index (2020).

[하루에 한 R] pair 형태의 데이터를 matrix로 전환할 때 주의할 점

오늘 쓰려는 주제는 R 활용 기법과 직접적인 관련은 없음을 미리 밝힌다.

한 쌍의 유전체 사이에 어떤 계산을 하여 수치를 얻었다고 가정하자. 기본적으로 얻어지는 결과물은 다음과 같이 pair 형태일 것이다.

친절한 프로그램이라면 matrix 형태로 전환한 파일을 제공해 줄 것이고, 이를 R에서 데이터 프레임으로 읽어들여서 여러 가지의 후속 분석을 실시할 수 있다. 만약 pairwise data만 위 그림과 같이 제공한다면 R에서 reshape 패키지를 이용하여 matrix 형태로 전환하면 된다.

그런데 여기에서 주의할 점이 하나 있다. 바로 자기 자신에 대한 데이터를 포함하고 있느냐를 고려해야 한다. 위에서 보인 그림 중 첫번째 것은 dRep의 결과물이다. FASTA 파일 여러개를 제공하면 프로그램이 알아서 모든 쌍에 대한 ANI를 계산하는데, 이때 자기 자신에 대한 것도 포함한다. 계산값을 보여주는 첫 줄이 그러한 것으로서 alignment coverage와 ani 전부 1,0으로 나타난다.

이것이 전부가 아니다. A->B, B->A 각 경우에 대한 값이 전부 존재하는가? 그렇다면 동일한가? 결과가 없는 것은 어떻게 할 것인가? 의외로 생각할 것이 많다. 그리고 각 경우에 대해서 처리하는 R 코드가 조금씩 달라진다.

위에 보인 그림에서 아래의 것은 POCP.sh를 사용하여 percentage of conserved protein을 계산한 것이다. 이 스크립트에서는 유전체 서열 한 쌍과 스레드 수만 인수로 제공해야 한다. 따라서 query와 subject 서열의 쌍을 별도로 계산하여 넣어 주어야 한다. 나는 Perl의 Algorithm::Combinatorics 모듈을 사용했으므로 당연히 자기 자신에 대한 계산을 하지 않게 만들었다(샘플 코드는 여기를 참조). 계산을 하나마나 POCP = 100%가 나올 것이기 때문이다. 하지만 이를 고려하지 않고 R의 reshape 패키지를 그냥 돌리면 매트릭스 형태가 좀 이상해진다. 전부 1(=100%)로 채워지는 대각선에 대한 배려가 없기 때문이다. 10개의 유전체 서열을 전부 비교하여도 R로 읽어들어 전환한 다음에는 10x10의 매트릭스에서 하나가 부족해진다. query = subject에 해당하는 대각선도 제대로 존재하지 않는다.

이 문제를 해결하려면 pairwise 결과 파일의 뒷부분에 자기 자신에 대한 분석값을 넣는 것이 좋다. A, B, C,...J의 10개 유전체에 대해서 POCP를 계산한다면 다음의 10줄을 넣기만 하면 된다. 자기 자신에 대한 POCP 값은 당연히 100%이므로 이를 세번째 컬럼에 입력해 넣는다.

A A 100.0
B B 100.0
C C 100.0
..
J J 100.0

이렇게 하여 reshape의 cast() 함수를 돌려 처리를 하면 10x10의 매트릭스, 정확히 말하자면 데이터 프레임(df라 하자)이 생길 것이다. A->B, B->A 두 경우에 대해서 한쪽으로만 계산을 하였으니 데이터 프레임으로 전환된 뒤의 구조를 확인해 보면 대각선은 전부 100이고 나머지 절만만 채워진 상태일 것이다.

NA를 전부 0으로 치환한 다음 df와 t(df)를 더하면 대각선 절반에 대한 값이 채워진다. 그러나 대각선은 전부 200이 될 것이다. 이는 또 적절히 처리하면 된다.

오늘의 교훈은 무엇인가? 데이터 프레임을 한번쯤은 엑셀 등에서 열어보고 어떤 상태인지를 눈으로 확인하는 것이 좋다.

2019년 3월 14일 업데이트

Algorithm::Combinatorics의 combination_with_repetition() 서브루틴을 사용하면 자기 자신에 대한 검색 결과도 만들 수 있다. 이것을 사용하여 POPC.sh를 실행하는 Perl 코드는 다음과 같다.

#!/usr/bin/perl

use strict; use warnings;
use Algorithm::Combinatorics qw(combinations_with_repetition);

open LIST, 'list';
my @temp = ();
while (<LIST>) {
    chomp;
    push @temp, $_;
}
my $faa_ref = \@temp;

my $iter = combinations_with_repetition($faa_ref, 2);

while (my $c = $iter->next) {
    my $out = 'pocp' . '_' . $c->[0] . '_vs_' . $c->[1];
    print "Running for $out...\n";
    system("POCP.sh $c->[0] $c->[1] 24 > $out");
}