写在前面
午间，收到朋友求助，问及如何进行序列的 KOG 注释...Emmm，这个可就太久远了。事实上，无论是 COG 还是 KOG 数据库，基本都没怎么更新过。其中 KOG 某种程度属于 COG 的拓展，后者主要用于原核生物，前者用于真核生物。尽管COG在2014年（我刚接触生信的时候）更新了一般，但KOG就没更新过，停留在 20年前。这个数据库主要是基于 7 个主要真核生物构建，其中植物应该就拟南芥... 收到朋友的求助，我还是非常惊讶，怎么这个年代还有人做KOG注释。然而，无论如何，还是看看。
下述为操作记录，比较简单
下载数据
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
结果示例
对应的后续可以做一下结果可视化，比如大家都喜欢进行归类信息，做个柱形图云云，具体参考「fun.txt」的分类整理，然后绘图就可以了（PS: 注意，有一些Group同时归属于两个大字母，比如 T 或者 U）
写在最后
啊，怀念小时候。2014年我刚接触生物信息的时候，这些东西看起来都是高大上，只要会整，那就.... 现在来看...Emmm
网友评论