午间,收到朋友求助,问及如何进行序列的 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
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
联系客服