COG database annotation

Step 1
wget -bc -r -np -nH -nd  ftp://ftp.ncbi.nih.gov/pub/COG/COG2020/data/fasta

wget -bc ftp://ftp.ncbi.nih.gov/pub/COG/COG2020/data/cog-20.cog.csv

wget -bc ftp://ftp.ncbi.nih.gov/pub/COG/COG2020/data/cog-20.def.tab

---

Step 2

cat fasta/*fa.gz > COG_2020.fa.gz

seqkit rmdup COG_2020.fa.gz > COG_2020_rmdup.fa.gz

diamond makedb --in COG_2020_rmdup.fa.gz -d COG_2020_rmdup.dmnd

---

Step3

diamond blastx -d COG_2020_rmdup.dmnd -b 6 -q test_file.fa -f 5 -o test_file.xml  #for nucleotide sequence

diamond blastp -d COG_2020_rmdup.dmnd -b 6 -q test_file.fa -f 5 -o test_file.xml  #for amino acid sequence

---

Step4

perl gene_cog_annotaion.pl -x test_file.xml -c cog-20.cog.csv -t cog-20.def.tab

result file: test_file_cog_func_anno.tsv, test_file.tsv

use strict;

use warnings;

use Getopt::Long;

use Bio::SearchIO;

my $ver = 0.77;

my %opts;

GetOptions(\%opts, "x=s", "t=s", "c=s", "h" );

if(!defined($opts{x}) || !defined($opts{c}) || !defined($opts{t}) || defined($opts{h})){

        print <<"      Usage End.";

#-------------------------------------------------------------------------------------------------------------------------------------------

        Description: COG database annotation from blast xml format out

      example:perl $0 -x test_file.xml -c cog-20.cog.csv -t cog-20.def.tab

      Version: $ver 

        General options:

                -x          blast/diamond output xml file                                          <infile>      xml format[must]

                -c          Comma-delimited plain text file assigning proteins to COGs(cog-20.cog.csv)      <infile>      csv format[must]

                -t          Tab-delimited plain text file with COG descriptions(cog-20.def.tab)            <infile>      tsv format[must]

                -h          Help document

#-------------------------------------------------------------------------------------------------------------------------------------------

      Usage End.

        exit;

}

my $blastxml = $opts{x};

my $csv = $opts{c};

my $tsv = $opts{t};

my ($name) = $blastxml=~/(\S+?)\.xml/;

open TSV,">$name.tsv";

print TSV "qseqid\tsseqid\thit_desc\thit_rank\tpercent_identity\tquery_start\tquery_end\tsubject_start\tsubject_end\tevalue\tbits\n";

my $searchio = Bio::SearchIO->new( -format => 'blastxml', -file => $blastxml );

my %info;

while ( my $result = $searchio->next_result() ) { #perldoc Bio::Search::Result::ResultI

my $id = $result->query_name();

my $query_desc = $result->query_description();

my $dbname = $result->database_name();

my $size = $result->database_letters();

my $num_entries = $result->database_entries();

#print "#dbname\t$dbname\nnum_entries\t$num_entries\nsize\t$size";

while( my $hit = $result->next_hit ) {  #process the Bio::Search::Hit::HitI object

my $hit_name = $hit->name();

my $hit_desc = $hit->description();

        my $len = $hit->length();

my $rank = $hit->rank(); #the Nth hit for a specific query

print TSV "$query_desc\t";

print TSV "$hit_name\t$hit_desc\t$rank\t";

while( my $hsp = $hit->next_hsp ) {    # process the Bio::Search::HSP::HSPI object

my $query_start = $hsp->start('query');

my $query_end = $hsp->end('query');

my $subject_start = $hsp->start('subject');

my $subject_end = $hsp->end('subject');

my $evalue = $hsp->evalue;

#my $num_identical = $obj->num_identical(); #returns the number of identical residues in the alignment

my $percent_identity = sprintf "%.2f",$hsp->percent_identity(); #Returns the calculated percent identity for an HSP(0,100)

#my $num_conserved = $obj->num_conserved($newval)  #returns the number of conserved residues in the alignment

my $query_len = $hsp->length('query');          #length of query seq (without gaps)

my $hit_len = $hsp->length('hit');              #length of hit seq (without gaps)

my $total_len = $hsp->length('total');          #length of alignment (with gaps)

my $query_gaps = $hsp->gaps('query');            #num conserved / length of query seq (without gaps)

my $hit_gaps = $hsp->gaps('hit');                #num conserved / length of hit seq (without gaps)

my $total_gaps = $hsp->gaps('tatla');            #num conserved / length of alignment (with gaps)

my $hit = $hsp->hit;

my $score = $hsp->score();

my $bits = $hsp->bits();

print TSV "$percent_identity\t$query_start\t$query_end\t$subject_start\t$subject_end\t$evalue\t$bits\n";

$info{$query_desc}="$hit_name\t$hit_desc\t$evalue" if $rank==1;

}

}

}

#-----------------------------------------------------------------------------------------------------------------------------------------------

my %cog_cog=&cog_csv;

my %cog_def=&cog_tab;

open F, ">$name\_cog_func_anno.tsv";

print F "#query_id\tprot_id\tcog_id\tgene_id\tcog_membership_class\tcog_functional_category\tcog_name\tcog_functional_pathway\n";

for my $query_id(sort keys %info){

my ($prot_id, $prot_disc, $evalue) = split/\t/, $info{$query_id};

print F "$query_id\t$prot_id\t";

my ($cog_id, $gene_id,$cog_membership_class) = split/\t/, $cog_cog{$prot_id};

print F "$cog_id\t$gene_id\t$cog_membership_class\t";

my ($cog_functional_category,$cog_name,$cog_functional_pathway)=split/\t/,$cog_def{$cog_id};

print F "$cog_functional_category\t$cog_name\t$cog_functional_pathway\n";

}

#-----------------------------------------------------------------------------------------------------------------------------------------------

sub cog_csv{

open IN, $csv;

my %cog_cog;

while(<IN>){

chomp;

my ($gene_id, $protein_id, $cog_id, $cog_membership_class) = (split/\,/,$_)[0, 2, 6, 8];

$cog_cog{$protein_id} = "$cog_id\t$gene_id\t$cog_membership_class";

}

return %cog_cog;

close;

}

sub cog_tab{

open IN, $tsv;

my %cog_def;

while(<IN>){

chomp;

my ($cog_id, $cog_functional_category,$cog_name,$cog_functional_pathway) = (split/\t/,$_)[0, 1, 2, 4];

$cog_functional_pathway = $cog_functional_pathway?$cog_functional_pathway:"NA";

$cog_def{$cog_id} = "$cog_functional_category\t$cog_name\t$cog_functional_pathway";

}

return %cog_def;

close IN;

}