User Tools

Site Tools


Sidebar

This is the sidebar. Without it, the main text is too wide!

cog_assignment

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

COG(clusters of orthologous groups)의 개요는 블로그에 포스팅한 COG 2014년 개정판 음미하기를 참조하라. 개정판에 대한 논문은 2015년 초에 Nucl. Acids Res.에 게재되었다.

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

자료 다운로드

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

COG 데이터 사용 시 주의할 점

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 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
  • 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
...

Query sequence의 서열 ID 정비

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

Query sequence file의 준비

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

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

COGnitor process

BLAST 검색 세 차례

query의 자체 검색(filter OFF), 그리고 query → COGs에 대한 검색(filter OFF and ON)을 실시한다. blastp가 아니라 psiblast를 이용한다는 점이 의외였다. Multiple processor를 쓴다면 -num_threads <NUM>을 해도 좋다. 세번째 psiblast 검색에 가장 많은 시간이 걸린다. 이를 단축할만한 획기적인 방법이 있으면 정말 좋을 것이다.

$ 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

Preparation of "sequence universe"

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

BLAST 결과물 처리

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

COGNITOR 실행

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

모든 것이 끝났다. 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
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 파일 아니겠는가?

개선 아이디어

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

cog_assignment.txt · Last modified: 2018/02/19 08:40 by hyjeong