User Tools

Site Tools


cog_assignment

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
cog_assignment [2018/01/06 09:29] – [개선 아이디어] hyjeongcog_assignment [2021/08/10 12:51] (current) – [Functional classification] hyjeong
Line 1: Line 1:
 +====== NCBI COG software를 이용한 query protein의 COG assignment 방법 ======
  
 +COG(clusters of orthologous groups)의 개요는 블로그에 포스팅한 [[http://blog.genoglobe.com/2016/11/cogclusters-of-orthologous-groups-2014.html|COG 2014년 개정판 음미하기]]를 참조하라. 개정판에 대한 논문은 2015년 초에 [[https://www.ncbi.nlm.nih.gov/pubmed/25428365|Nucl. Acids Res.]]에 게재되었다. 2020년 업데이트에 대한 글은 [[https://blog.genoglobe.com/2021/08/cog-2020.html|반가운 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
 +  * [소프트웨어 다운로드] [[ftp://ftp.ncbi.nih.gov/pub/wolf/COGs/COGsoft|NCBI]] [[https://sourceforge.net/projects/cogtriangles/files/|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 정보가 부여되었으며, 명칭 자체도 많이 다듬어졌다고 한다.
 +
 +==== 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로 표시되는 값이다. 약간 성가시지만 이 값을 참조하면 된다. 
 +
 +{{ :protein_id.png?direct&400 |}}
 +
 +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 형식은 [[https://www.ncbi.nlm.nih.gov/protein/AAM01497.1/|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
 +  ...
 +  
 +==== 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
 +
 +파일의 구조를 간단히 그림으로 설명하면 다음과 같다. 그림 아래에서 **GenQueryquery.p2o.csv**라는 파일명은 **GenQuery.p2o.csv**로 고쳐야 한다. 그림 파일의 원본이 없어서 볼썽사납지만 그대로 둔다.
 +{{ :cog.png?600 |}}
 +
 +이상의 작업을 하는 것이 너무 성가셔서 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
 +  
 +
 +===== COGnitor process =====
 +==== BLAST 검색 세 차례 ====
 +query의 자체 검색(filter OFF), 그리고 query -> COGs에 대한 검색(filter OFF and ON)을 실시한다. blastp가 아니라 psiblast를 이용한다는 점이 의외였다. Multiple processor를 쓴다면 -num_threads <NUM>을 해도 좋다. 두번째 및 세번째 psiblast 검색에 가장 많은 시간이 걸린다(-seg no로 필터링을 하지 않는 것이 좀 더 빠름. 당연하지 않은가? [[http://www.biology.wustl.edu/gcg/seg.html|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
 +==== 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 실행 ====
 +이 과정에서는 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|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 파일 아니겠는가? 
 +
 +==== Functional classification ====
 +COG의 각 functional class에 몇 개씩의 유전자가 분포하는지를 집계하고 싶다면 [[cogclass2018.pl|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 파일의 끝부분에 다음과 같은 집계표가 나온다.
 +{{ :cog_class.png?600 |}}
 +===== 개선 아이디어 =====
 +Gen1.fa, Gen2.fa... 파일을 인수로 공급하면 일괄적으로 COG assignment를 실시하는 스크립트를 만들어 보자.