NCBI COG software를 이용한 query protein의 COG assignment 방법

COG(clusters of orthologous groups)의 개요는 블로그에 포스팅한 COG 2014년 개정판 음미하기를 참조하라. 개정판에 대한 논문은 2015년 초에 Nucl. Acids Res.에 게재되었다. 2020년 업데이트에 대한 글은 반가운 COG 2020에 작성하였다.

NCBI가 공개한 COG software는 직접 COG를 만드는 것 외에도 기존의 COG 체계에 대해 query protein을 검색하여 알맞은 COG 정보를 부여하는 기능을 제공한다. 그런데 설명 파일이 다소 난해하여 한동안 쓸 생각을 하지 않고 있다가 약간의 시행착오를 거쳐서 성공에 이르렀기에 설명 자료를 남긴다. 본 페이지에서는 COG software에 포함된 Readme.2012.04.txt 파일의 3.2 COGnitor 섹션만을 실제 실행에 옮기면서 현실에 맞게 작성한 것이다. 나만의 COG를 만들고자 한다면 '3.1. Making COGs' 섹션도 읽어야 한다. 참고로 이 실행 예제는 COG2014 data 디렉토리에서 실행하는 것을 전제로 한다.

• [NCBI의 COG 공식 웹사이트] https://www.ncbi.nlm.nih.gov/COG/
• [2003년 COG의 2014년도 업데이트, html] ftp://ftp.ncbi.nih.gov/pub/COG/COG2014/static/lists/homeCOGs.html
• [2003년 COG의 2014년도 업데이트, data] ftp://ftp.ncbi.nih.gov/pub/COG/COG2014/data
• [소프트웨어 다운로드] NCBI SourceForge (COGsoft)
• [데이터 파일 위치(로컬)] /data/Utilities/DB/COG/COG2014
• [COG 2020] https://ftp.ncbi.nih.gov/pub/COG/COG2020/data/

개정판에는 총 711개의 유전체에서 예측한 단백질 서열을 대상으로 작성한 4631개의 COG가 수록되었다. 2003년판의 4873 COG 중 242개가 제거되었는데(주로 yeast), 새로 추가된 것은 없다. Functional classification 체계가 더욱 세분화되고, 기능을 잘 모르던 COG에 대해서 좀 더 구체적인 annotation 정보가 부여되었으며, 명칭 자체도 많이 다듬어졌다고 한다.

prot2003-2014.fa 파일의 모습은 다음과 같다.

>gi|103485499|ref|YP_615060.1| chromosomal replication initiation protein [Sphingopyxis alaskensis RB2256]
MSGDAAALWPRVAEGLRRDLGARTFDHWLKPVRFADYCALSGVVTLETASRFSANWINERFGDRLELAWRQQLPAVRSVS
...

문제는 여기에 표시된 서열 ID(gi number and protein accession)이 이제는 유효하지 않다는 것이다. NCBI 웹사이트의 검색창에 이를 넣으면 더 이상 유효하지 않다는 메시지와 함께 현재의 레코드로 연결이 되지만 일괄적인 작업은 어렵다. 현재의 GenBank 정보와의 연결성을 알고 싶다면 prot2003-2014.gi2gbk.tab 파일을 열어보라.

103485499       YP_615060       ABF51727
103485500       YP_615061       ABF51728
...

세번째 필드가 바로 현재의 GenBank genome data에서 protein_id로 표시되는 값이다. 약간 성가시지만 이 값을 참조하면 된다.

RefSeq에서는 이 키워드로 찾지 못한다. 왜냐하면 RefSeq의 protein_id는 WP_042447374.1 형식이기 때문이다.

COG 2020의 cog-20.fa 파일은 다음과 같은 모습이다.

>AAM01497_1 Glutamate-1-semialdehyde aminotransferase [Methanopyrus kandleri AV19]
MGYEDEFPESLELFKRAERVMPGGVSSPVRRFDPYPFYVERAEGSRLYTVDGHVLIDYCLAFGPLILGHAHPEVVEAVVER
...

AAM01497_1라는 protein ID 형식은 AAM01497.1로 바꾸어야 할 것 같다. 왜냐하면 이것이 현재 통용되는 accession이기도 하고, 아래에 일부를 보인 cog-20.cog.csv 파일의 세 번째 컬럼에서도 그렇게 쓰이기 때문이다. Perl script를 사용하여 원본 cog-20.fa의 서열 ID 라인을 수정하여 blast DB를 만들었다.

MK0280,GCA_000007185.1,AAM01497.1,430,1-430,430,COG0001,COG0001,0,570.0,1.0e-200,432,3-431
MA_0581,GCA_000007345.1,AAM04025.1,424,1-424,424,COG0001,COG0001,0,574.0,1.0e-200,432,3-426
...

COG 2014년 개정판 데이터와 COG software 2012년판을 다운로드하여 압축을 푼다. 소프트웨어 압축파일을 풀면 다음의 네 가지 디렉토리가 만들어진다(본 과정은 패키지에 수록된 Readme.2012.04.txt 파일을 참조하여 작성한 것임).

COGcognitor  COGlse  COGmakehash  COGreadblast  COGtriangles

각 디렉토리로 들어가서 make를 실행하면 디렉토리 이름에 해당하는 실행파일이 생긴다. gcc 4.4.7로는 잘 빌드되지만 gcc 5.3.0(Linuxbrew)로는 오류가 발생하였다. 만들어진 실행파일을 $PATH에 위치시킨다.

COG data 파일 중 당장 필요한 것은 다음과 같다. /data/Utilities/DB/COG/COG2014에 있는 것을 현재의 작업 디렉토리로 가지고 온다. 아니면 이 디렉토리로 query sequence file을 복사하여 거기에서 작업을 해라! 이렇게 하는 것이 더 편리하다. 그 대신 이전에 작업한 흔적은 잘 지워야 한다. fun2003-2014.tab(one-letter functional category)과 cognames2003-2014.tab(COG 번호, function code, COG 명칭) 파일은 COG 계산 작업에 직접적으로 필요하지는 않다.

• prot2003-2014.fa BLASTDB로 전환된 다음에는 필요하지 않다.
• cog2003-2014.csv ⇒ COGs.csv라는 이름의 심볼릭 링크를 만들든지 복사를 해라.
• COG.p2o.csv: 다음과 같이 <protein id>,<genome id>로 만들어진 파일이다. cog2003-2014.csv 파일의 첫번째와 두번째 컬럼을 awk 명령으로 떼어내면 된다. 아래에 예를 든 awk 명령어는 별로 마음에 들지 않는다.

$ awk -F "," OFS='","{print $1,$2}' cog2003-2014.csv > COG.p2o.csv
$ cat COG.p2o.csv
158333741,Acaryochloris_marina_MBIC11017_uid58167
158339504,Acaryochloris_marina_MBIC11017_uid58167
...

COG 분석을 할 genome에 대해서 단백질 서열을 수록한 multi-fasta 파일이 각각 필요하다(Gen1.fa, Gen2.fa, Gen3.fa…). 서열 ID는 “lcl|<text-id>” 또는 “gi|<num-id>“의 형식을 갖추어야 한다. 점(.)이나 대쉬(-)는 포함할 수 있으나 권장하지 않으며, 쉼표(,), 세미콜론(;), 콜론(:) 등의 non-alphanumeric character는 있어서는 안된다. 대소문자를 구별하지 않는 것에도 유의하라.

Gen1.fa, Gen2.fa라는 두 개의 multi-fasta file이 있다고 가정하자. 즉 Gen1이라는 genome의 protein set와 Gen2라는 genome의 protein set가 있다는 뜻이다. 두 개의 파일은 하나로 합쳐 놓는다.

$ cat Gen1.fa Gen2.fa > query.fa

query에게도 <protein id>,<genome id> 파일이 있어야 한다. 역시 마음에 들지 않는 awk 코맨드이지만 다음과 같이 하면 될 것이다. 단, GenQuery.p2o.csv 파일에 수록된 서열 ID 시작 부분에서 “lcl|” 또는 “gi|”를 제거해야 한다.

$ grep '>' Gen1.fa | sed -e 's/^>//' | awk 'OFS=","{print $1,"Gen1"}' > file1
$ grep '>' Gen2.fa | sed -e 's/^>//' | awk 'OFS=","{print $1,"Gen2"}' > file2
$ cat file1 file2 > GenQuery.p2o.csv

파일의 구조를 간단히 그림으로 설명하면 다음과 같다. 그림 아래에서 GenQueryquery.p2o.csv라는 파일명은 GenQuery.p2o.csv로 고쳐야 한다. 그림 파일의 원본이 없어서 볼썽사납지만 그대로 둔다.

이상의 작업을 하는 것이 너무 성가셔서 Gen1,fa, Gen2.fa.. 파일을 인수로 주면 query.fa 및 GenQuery.p2o.csv 파일을 만드는 스크립트 process_COG_query.pl를 하나 작성해 보았다.

#!/bin/sh

# usage : process_COG_query.sh query1.fa query2.fa ...

for query in "$@"
do
    cat $query >> tmp.query.fa
    prefix=${query%.*}
    grep '>' $query | sed -e 's/^>//' | awk -v gen="$prefix" 'OFS=","{print $1,gen}' >> tmp.file
done

sed -e 's/^lcl|//' tmp.file > GenQuery.p2o.csv

awk '/^>/ {
          split($0, field, "|")
	  if (field[1] == ">lcl")
	      {
	      print
	      }
	  else
	      {
	      gsub(">",">lcl|", $0);
	      print 
	      }
	  ;next}
	 {print}' tmp.query.fa > query.fa

rm tmp.query.fa tmp.file

Query protein set와 COG protein에 대해서 blast DB를 만든다. 사용한 fasta file의 이름과 DB name이 다른 것에 유의한다. 혼동을 피하기 위해서 이렇게 한 것이니 각자 자율적으로 결정해도 된다. prot2003-2014.fa에서 만들어진 DB는 일정한 곳에 보관하여 계속 재사용하는 것이 좋을 것이다.

$ makeblastdb -in query.fa -dbtype prot -out Query
$ makeblastdb -in  prot2003-2014.fa -dbtype prot -out COGs

마지막으로 BLASTcogn, BLASTconv, BLASTff, BLASTno, BLASTss라는 디렉토리를 만든다.

$ mkdir BLASTcogn BLASTconv BLASTff BLASTno BLASTss

query의 자체 검색(filter OFF), 그리고 query → COGs에 대한 검색(filter OFF and ON)을 실시한다. blastp가 아니라 psiblast를 이용한다는 점이 의외였다. Multiple processor를 쓴다면 -num_threads <NUM>을 해도 좋다. 두번째 및 세번째 psiblast 검색에 가장 많은 시간이 걸린다(-seg no로 필터링을 하지 않는 것이 좀 더 빠름. 당연하지 않은가? SEG 프로그램을 실행하지 않으므로!). 이를 단축할만한 획기적인 방법이 있으면 정말 좋을 것이다. 아래 사례에서는 -num_threads int_value 옵션이 빠져 있으니 그대로 실행하면 시간이 대단히 많이 걸림에 유의하라.

$ psiblast -query query.fa -db Query -show_gis -outfmt 7 -num_descriptions 10 -num_alignments 10 -dbsize 100000000 -comp_based_stats F -seg no -out BLASTss/QuerySelf.tab
$ psiblast -query query.fa -db COGs -show_gis -outfmt 7 -num_descriptions 1000 -num_alignments 1000 -dbsize 100000000 -comp_based_stats F -seg no -out BLASTno/QueryCOGs.tab
$ psiblast -query query.fa -db COGs -show_gis -outfmt 7 -num_descriptions 1000 -num_alignments 1000 -dbsize 100000000 -comp_based_stats T -seg yes -out BLASTff/QueryCOGs.tab

• psiblast warning: The parameter -num_descriptions is ignored for output formats > 4 . Use -max_target_seqs to control output

COG software의 Readme 파일에는 거창하게 나왔지만, 사실은 별 것이 아니다. 모든 서열에게 숫자로 이루어진 ID를 부여하는 것이다. ./BLASTcogn/hash.csv file이 만들어진다.

$ cat GenQuery.p2o.csv COG.p2o.csv > tmp.p2o.csv
$ COGmakehash -i=tmp.p2o.csv -o=./BLASTcogn -s="," -n=1

$ COGreadblast -d=./BLASTcogn -u=./BLASTno -f=./BLASTff -s=./BLASTss -e=0.1 -q=2 -t=2

이 과정에서는 COGs.csv라는 파일이 필요하다. 2003-2014 버전과 2020 버전의 컬럼 수가 다르므로 COG software를 쓰려면 2020 버전 파일을 고쳐야 한다. 여기서 약간의 시행착오를 겪었다. 두 버전의 컬럼 설명이 완전히 1:1 대응하듯이 일치하지를 않아서 그야말로 '대충' 전환을 하였다.

$ COGcognitor -i=./BLASTcogn -t=COGs.csv -q=GenQuery.p2o.csv -o=GenQuery.COG.csv

다음은 전환용 스크립트. 작동의 완전성을 보장하지 못한다. 특히 마지막 컬럼은 2003-2014 버전의 설명에 의하면 'The membership-class field indicates the nature of the match between the sequence and the COG consensus'라 하여 0부터 3까지의 값을 갖는데, 나는 무조건 0(the domain matches the COG consensus)로 해 버렸다.

#!/usr/bin/perl
#

open COG2020, 'cog-20.cog.csv';
open OUT, '>COGs.csv';
while (<COG2020>) {
    chomp;
    my @temp = split /,/, $_;
    if ($temp[12] =~ /^(\d+)\-(\d+)$/) {
       ($start, $end) = ($1, $2);
    }
    print OUT join ',', $temp[2], $temp[1], $temp[2], $temp[3], $start, $end, $temp[6], 0 . ",\n";
}

모든 것이 끝났다. GenQuery.COG.csv에 최종적으로 할당된 COG 정보가 수록된다. 주의할 점은 하나의 Query protein에 대하여 복수의 COG가 부여될 수 있다. 이러한 경우에는 5번째 필드인 cognitor-score를 가지고 1등을 선별해야 할 것이다.

$ cat GenQuery.COG.csv
6666666.214180.peg.4,650,1,650,3763.34,COG3505
6666666.214180.peg.8,788,1,788,6665.55,COG3451
6666666.214180.peg.9,829,1,480,0,-1
6666666.214180.peg.9,829,481,829,3475.96,COG0741
6666666.214180.peg.24,271,1,271,1221.63,COG1192
...

findBestFromCOGs.pl이라는 스크립트를 사용하여 GenQuery.COG.csv를 처리하면 각 protein query에 대한 best COG를 출력한다.

$ findBestFromCOGs.pl GenQuery.COG.csv > GenQuery.COG.csv.bestHit
$ cat GenQuery.COG.csv.bestHit
AND37645.1,448,1,448,32900.1,COG0593,sinlge
AND37646.1,378,1,378,17170.9,COG0592,sinlge
AND37647.1,69,1,69,2273.97,COG2501,sinlge
AND37648.1,372,1,372,15085.1,COG1195,sinlge
AND37649.1,82,1,82,0,-1,none
...(생략)...

COGNITOR data(COGcognitor 프로그램의 아웃풋을 COG software Readme 파일 섹션 2.11에서는 이렇게 부름)의 형식은 ”<prot-id>,<prot-length>,<prot-start>,<prot-end>,<cognitor-score>,<cluster-id>,…“이다, <cognitor-score>란 주어진 query protein fragment를 해당 COG cluster에 할당한 것에 대한 상대적인 confidence를 의미하는 수치이다.

마찬가지로 Clustering data(Readme 파일 섹션 2.10)란, COG를 만드는 프로그램인 COGtriangle의 결과물이다. 바로 COG data package에 포함된 cog2003-2014.csv 파일 아니겠는가?

COG의 각 functional class에 몇 개씩의 유전자가 분포하는지를 집계하고 싶다면 COGclass2018.pl 스크립트를 이용하라. 하나의 COG에 대하여 복수의 functional class가 있는 경우는 무작위로 하나만 선택한다. 이게 최선이 아닌 줄은 알지만 어쩔 수가 없어서…

$ COGclass2018.pl GenQuery.COG.csv.bestHit > COGclass
...
COG5503 KV (random selection) ===> K # STDERR 출력
COG5524 CT (random selection) ===> C # STDERR 출력

COGclass 파일의 끝부분에 다음과 같은 집계표가 나온다.

Gen1.fa, Gen2.fa… 파일을 인수로 공급하면 일괄적으로 COG assignment를 실시하는 스크립트를 만들어 보자.