annotatingTranscriptomes/annotationsPileup.pl at master · z0on/annotatingTranscriptomes · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
#!/usr/bin/perl

my $usage="
Annotations pileup: takes fasta file with contigs as main identifiers
(the first word after \'>\' symbol),
adds to the headers whatever it finds in the supplied space- or tab-delimited
lookup tables:
contig2isogroup (c2i)
isogroup2genename (i2gn)
isogroup2GO (i2go)
isogroup2KEGG (i2k)

arguments:

required:
fasta=[file name]
c2i=[file name]

optional:
i2gn=[file name]
i2go=[file name]
i2k=[file name]

";

my $fasta;
if ("@ARGV"=~/fasta=(\S+)/) { $fasta=$1;}
else { die $usage;}
my $c2i;
if ("@ARGV"=~/c2i=(\S+)/) { $c2i=$1;}
else { die $usage;}
my $i2gn;
if ("@ARGV"=~/i2gn=(\S+)/) { $i2gn=$1;}
my $i2go;
if ("@ARGV"=~/i2go=(\S+)/) { $i2go=$1;}
my $i2k;
if ("@ARGV"=~/i2k=(\S+)/) { $i2k=$1;}

my $e1;
my $e2;
my @ee;

my %c2iso={};
open TT, $c2i or die "cannot open contigs to isogroups table $c2i\n";
while (<TT>) {
	chomp;
	($e1,$e2)=split(/\s/,$_);
	$c2iso{$e1}=$e2;
#print "$e1|$c2iso{$e2}|c2iso\n";
}
close TT;
my %i2gene={};
if ($i2gn){
	open TT, $i2gn or die "cannot open isogroups to genenames table $i2gn\n";
	while (<TT>) {
		chomp;
		($e1,@ee)=split(/\s/,$_);
#		$e2=join("_",@ee);
		$i2gene{$e1}="@ee";
#print "$e1|$i2gene{$e1}|i2gene\n";
	}
}
close TT;
my %i2GO={};
if ($i2go){
	open TT, $i2go or die "cannot open isogroups to GO table $i2go\n";
	while (<TT>) {
		chomp;
		($e1,$e2)=split(/\s/,$_);
		$i2GO{$e1}=$e2;
#print "$e1|$i2GO{$e1}|i2GO\n";
	}
}
close TT;
my %i2kegg={};
if ($i2k){
	open TT, $i2k or die "cannot open isogroups to KEGG table $i2k\n";
	while (<TT>) {
		chomp;
		($e1,$e2)=split(/\s/,$_);
		$i2kegg{$e1}=$e2;
#print "$e1|$i2kegg{$e1}|i2kegg\n";
	}
}
close TT;

open FF, $fasta or die "cannot open fasta file $fasta\n";
my $cname="";
my $seq="";
my @noiso=();
my $newname;

while (<FF>) {
	chomp;
	if ($_=~/^>(\S+)/) {
		$newname=$1;
		if ($seq && $cname) {
			my $iso="";
			if ($c2iso{$cname}) {
				$iso=$c2iso{$cname};
				print ">$cname gene=$iso";
			}
			else {
				print ">$cname gene=NA ";
				push @noiso, $cname;
			}
			if ($i2gn && $i2gene{$iso}) {
				print " Gene=$i2gene{$iso}" unless ($i2gene{$iso}=~/annot missing/);
			}
			if ($i2go && $i2GO{$iso}) {
				print " GOterms=$i2GO{$iso}";
			}
			if ($i2k && $i2kegg{$iso}) {
				print " KEGGterms=$i2kegg{$iso}";
			}
			print "\n$seq\n";
		}
		$cname=$newname;
		$seq="";
	}
	else {$seq.=$_;}
}
if ($seq && $cname) {
	my $iso="";
	if ($c2iso{$cname}) {
		$iso=$c2iso{$cname};
		print ">$cname gene=$iso";
	}
	else {
		print ">$cname gene=NA ";
		$noiso++;
	}
	if ($i2gn && $iso && $i2gene{$iso}) {
		print " Gene=$i2gene{$iso}";
	}
	if ($i2go && $iso && $i2GO{$iso}) {
		print " GOterms=$i2GO{$iso}";
	}
	if ($i2k && $iso && $i2kegg{$iso}) {
		print " KEGGterms=$i2kegg{$iso}";
	}
	print "\n$seq\n";
}

if ($noiso[0]) {
	print STDERR "\n",$#noiso+1," contigs without isogroup designation:
	@noiso
	check whether your $c2i table is up to date
	";
}