寫在前面
午間，收到朋友求助，問及如何進行序列的 KOG 注釋...Emmm，這個可就太久遠了。事實上，無論是 COG 還是 KOG 數(shù)據(jù)庫，基本都沒怎么更新過。其中 KOG 某種程度屬于 COG 的拓展，后者主要用于原核生物，前者用于真核生物。盡管COG在2014年（我剛接觸生信的時候）更新了一般，但KOG就沒更新過，停留在 20年前。這個數(shù)據(jù)庫主要是基于 7 個主要真核生物構建，其中植物應該就擬南芥... 收到朋友的求助，我還是非常驚訝，怎么這個年代還有人做KOG注釋。然而，無論如何，還是看看。
下述為操作記錄，比較簡單
下載數(shù)據(jù)
wget https://ftp.ncbi.nih.gov/pub/COG/KOG/fun.txt
wget https://ftp.ncbi.nih.gov/pub/COG/KOG/kog
wget https://ftp.ncbi.nih.gov/pub/COG/KOG/kyva
# wget https://ftp.ncbi.nih.gov/pub/COG/KOG/kyva=gb
# wget https://ftp.ncbi.nih.gov/pub/COG/KOG/lse
# wget https://ftp.ncbi.nih.gov/pub/COG/KOG/pa
# wget https://ftp.ncbi.nih.gov/pub/COG/KOG/readme
wget https://ftp.ncbi.nih.gov/pub/COG/KOG/twog
序列比對
建庫
# makeblastdb -in kyva -input_type fasta -dbtype prot
比對（邏輯上可以用 ghostz 或者 diamond 替代）
# blastx -db kyva -query transcripts.fa -out transcripts.2.kog.blast.tab -evalue 1e-5 -qcov_hsp_perc 33 -outfmt 6 -max_target_seqs 20 -num_threads 40
實在太慢，我用 diamond
diamond makedb --in kyva -d kyva.diamond
diamond blastx --db kyva.diamond --query transcripts.fa --out transcripts.2.kog.diamond.tab --threads 20
過濾
perl -F'\t' -lane 'next if $seen{$F[0]}++;next if $F[-2]>1e-5;print' transcripts.2.kog.diamond.tab > transcripts.2.kog.diamond.tab.filtered
KOG 關系映射
perl -lane 'next unless $_;if(/^\[([A-Z]+)\]\s*(.*?)$/){$anno=$2;@flag=split "",$1;}else{print join qq{\t},$F[1],$1,$2}' kog > kog.parsed.tab
表格合并
perl -F'\t' -lane 'if(!$flag){$anno{$F[0]}=q{[}.$F[1].q{]}.$F[2]}else{print join qq{\t},$F[0],qq{$F[1]/$F[-2]/$anno{$F[1]}}}$flag=1 if eof ARGV' kog.parsed.tab transcripts.2.kog.diamond.tab.filtered > kog.anno.tab
結果示例
image
對應的后續(xù)可以做一下結果可視化，比如大家都喜歡進行歸類信息，做個柱形圖云云，具體參考「fun.txt」的分類整理，然后繪圖就可以了（PS: 注意，有一些Group同時歸屬于兩個大字母，比如 T 或者 U）
image
image
