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/05 16:03] – [사전 준비 사항] hyjeongcog_assignment [2021/08/10 12:51] (current) – [Functional classification] hyjeong
Line 1: Line 1:
 ====== NCBI COG software를 이용한 query protein의 COG assignment 방법 ====== ====== 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.]]에 게재되었다.+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를 만들고자 한다면 Readme 파일의 앞부분도 읽어야 한다. 참고로 이 실행 예제는 COG2014 data 디렉토리에서 실행하는 것을 전제로 한다.+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/   * [NCBI의 COG 공식 웹사이트] https://www.ncbi.nlm.nih.gov/COG/
Line 9: Line 9:
   * [2003년 COG의 2014년도 업데이트, data] ftp://ftp.ncbi.nih.gov/pub/COG/COG2014/data   * [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)]]   * [소프트웨어 다운로드] [[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 정보가 부여되었으며, 명칭 자체도 많이 다듬어졌다고 한다. 개정판에는 총 711개의 유전체에서 예측한 단백질 서열을 대상으로 작성한 4631개의 COG가 수록되었다. 2003년판의 4873 COG 중 242개가 제거되었는데(주로 yeast), 새로 추가된 것은 없다. Functional classification 체계가 더욱 세분화되고, 기능을 잘 모르던 COG에 대해서 좀 더 구체적인 annotation 정보가 부여되었으며, 명칭 자체도 많이 다듬어졌다고 한다.
Line 25: Line 27:
   ...   ...
      
-세번째 필드가 바로 현재의 GenBank genome data에서 protein_id로 표시되는 값이다.약간 성가시지만 이 값을 참조하면 된다. +세번째 필드가 바로 현재의 GenBank genome data에서 protein_id로 표시되는 값이다. 약간 성가시지만 이 값을 참조하면 된다. 
  
 {{ :protein_id.png?direct&400 |}} {{ :protein_id.png?direct&400 |}}
  
 RefSeq에서는 이 키워드로 찾지 못한다. 왜냐하면 RefSeq의 protein_id는 WP_042447374.1 형식이기 때문이다. 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 파일을 참조하여 작성한 것임). COG 2014년 개정판 데이터와 COG software 2012년판을 다운로드하여 압축을 푼다. 소프트웨어 압축파일을 풀면 다음의 네 가지 디렉토리가 만들어진다(본 과정은 패키지에 수록된 Readme.2012.04.txt 파일을 참조하여 작성한 것임).
Line 37: Line 48:
 각 디렉토리로 들어가서 make를 실행하면 디렉토리 이름에 해당하는 실행파일이 생긴다. gcc 4.4.7로는 잘 빌드되지만 gcc 5.3.0(Linuxbrew)로는 오류가 발생하였다. 만들어진 실행파일을 $PATH에 위치시킨다. 각 디렉토리로 들어가서 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 계산 작업에 직접적으로 필요하지는 않다. +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**+  * **prot2003-2014.fa** BLASTDB로 전환된 다음에는 필요하지 않다.
   * cog2003-2014.csv => **COGs.csv**라는 이름의 심볼릭 링크를 만들든지 복사를 해라.   * cog2003-2014.csv => **COGs.csv**라는 이름의 심볼릭 링크를 만들든지 복사를 해라.
   * **COG.p2o.csv**: 다음과 같이 <protein id>,<genome id>로 만들어진 파일이다. cog2003-2014.csv 파일의 첫번째와 두번째 컬럼을 awk 명령으로 떼어내면 된다. 아래에 예를 든 awk 명령어는 별로 마음에 들지 않는다.   * **COG.p2o.csv**: 다음과 같이 <protein id>,<genome id>로 만들어진 파일이다. cog2003-2014.csv 파일의 첫번째와 두번째 컬럼을 awk 명령으로 떼어내면 된다. 아래에 예를 든 awk 명령어는 별로 마음에 들지 않는다.
Line 61: Line 72:
   $ grep '>' Gen2.fa | sed -e 's/^>//' | awk 'OFS=","{print $1,"Gen2"}' > file2   $ grep '>' Gen2.fa | sed -e 's/^>//' | awk 'OFS=","{print $1,"Gen2"}' > file2
   $ cat file1 file2 > GenQuery.p2o.csv   $ 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는 일정한 곳에 보관하여 계속 재사용하는 것이 좋을 것이다. Query protein set와 COG protein에 대해서 blast DB를 만든다. 사용한 fasta file의 이름과 DB name이 다른 것에 유의한다. 혼동을 피하기 위해서 이렇게 한 것이니 각자 자율적으로 결정해도 된다. prot2003-2014.fa에서 만들어진 DB는 일정한 곳에 보관하여 계속 재사용하는 것이 좋을 것이다.
Line 75: Line 120:
 ===== COGnitor process ===== ===== COGnitor process =====
 ==== BLAST 검색 세 차례 ==== ==== BLAST 검색 세 차례 ====
-query의 자체 검색(filter OFF), 그리고 query -> COGs에 대한 검색(filter OFF and ON)을 실시한다. blastp가 아니라 psiblast를 이용한다는 점이 의외였다. Multiple processor를 쓴다면 -num_threads <NUM>을 해도 좋다. 세번째 psiblast 검색에 가장 많은 시간이 걸린다. 이를 단축할만한 획기적인 방법이 있으면 정말 좋을 것이다.+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 Query -show_gis -outfmt 7 -num_descriptions 10 -num_alignments 10 -dbsize 100000000 -comp_based_stats F -seg no -out BLASTss/QuerySelf.tab
Line 95: Line 140:
  
 ==== COGNITOR 실행 ==== ==== 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   $ 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
 +  #
      
-모든 것이 끝났다. **GenQuery.COG.csv**에 최종적으로 할당된 COG 정보가 수록된다. 주의할 점은 하나의 Query protein에 대하여 복수의 COG가 부여될 수 있. 이러한 경우에는 5번째 필드인 cognitor-score를 가지고 1등을 선별해야 할 것이다.+  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   $ cat GenQuery.COG.csv
Line 108: Line 167:
   ...   ...
  
 +[[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를 의미하는 수치이다. 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 파일 아니겠는가?  마찬가지로 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를 실시하는 스크립트를 만들어 보자.
cog_assignment.1515135814.txt.gz · Last modified: 2021/03/17 13:09 (external edit)