annotate ClassAnnotate.pl @ 0:07745c0958dd draft

Uploaded

1 #!/usr/bin/perl -w
2 #Filename:
3 #Author: Chen Tingting
4 #Email: [e-mail]
5 #Date: 2014/5/13
6 #Modified:
7 #Description: cluster annotate
8 my $version=1.00;
9 
10 use strict;
11 use Getopt::Long;
12 
13 my %opts;
14 GetOptions(\%opts,"i=s","g=s","n=s","r=s","p=s","o=s","t=s","l=s","h");
15 if (!(defined $opts{i} and defined $opts{g} and defined $opts{n} and defined $opts{r} and defined $opts{p} and defined $opts{o} and defined $opts{t}  and defined $opts{l}) || defined $opts{h}) { #necessary arguments
16 &usage;
17 }
18 
19 #my %gene;
20 my %gene1;
21 open IN,"<$opts{g}";
22 open OUT ,">$opts{l}";
23 print OUT "#ID\tchr\tstart\tend\tstrand\n";
24 my $n=1;
25 while (my $aline=<IN>) {
26 	chomp $aline;
27 	next if($aline=~/^\#/);
28 	my @tmp=split/\t/,$aline;
29 	my $ID;
30 	if ($tmp[2] eq "gene") {
31 		$tmp[0]=~s/Chr\./Chr/;
32 		#$tmp[0]=~s/Chr/chr/;
33 		my @infor=split/;/,$tmp[8];
34 		for (my $i=0;$i<@infor ;$i++) {
35 			if ($infor[$i]=~/Alias\=(\S+)$/) {
36 				$ID=$1;
37 				last;
38 			}
39 			else {
40 				$ID="unknown$n";
41 				$n++;
42 			}
43 		}
44 		#$gene{$tmp[0]}{$ID}=[$tmp[3],$tmp[4],$tmp[6]];#$gene{chr}{geneID}=[start,end,strand]
45 		push @{$gene1{$ID}},[$tmp[3],$tmp[4],$tmp[0]];
46 		print OUT "$ID\t$tmp[0]\t$tmp[3]\t$tmp[4]\t$tmp[6]\n";
47 	}
48 }
49 #while (my $aline=<IN>) {
50 #	chomp $aline;
51 #	next if($aline=~/^\#/);
52 #	my @tmp=split/\t/,$aline;#ID chr start end strand
53 #	push @{$gene1{$tmp[0]}},[$tmp[2],$tmp[3],$tmp[1]];
54 #	#$gene{$tmp[1]}{$tmp[0]}=[$tmp[2],$tmp[3],$tmp[1]];
55 #}
56 close IN;
57 close OUT;
58 
59 my %nat;
60 open TMP,">$opts{t}";
61 print TMP "#NAT_ID\tGene\tStrand\tChr\tGene_start\tGene_end\tAntiGene\tStrand\tChr\tAntiGene_start\tAntiGene_end\tType1\tType2\tNATS1_start\tNATS1_end\tNATS2_start\tNATS2_end\n";
62 
63 open IN,"<$opts{n}";
64 my $nat_n=1;
65 while (my $aline=<IN>) {
66 	next if($aline=~/^\#/);#osj LOC_Os05g02659 - LOC_Os01g24200 + trans 1559 1802 660 905 246 100nt -
67 	chomp $aline;
68 	my @arr=split /\t/,$aline;
69 	my ($ns,$ne,$ns2,$ne2)=(0,0,0,0);
70 	if ($arr[11]=~/Nearby/) {
71 		($ns,$ne)=&nearby($gene1{$arr[1]}[0][0],$gene1{$arr[1]}[0][1],$gene1{$arr[3]}[0][0],$gene1{$arr[3]}[0][1]);
72 		push @{$nat{$gene1{$arr[1]}[0][2]}},[$ns,$ne,$arr[5],$arr[11],"NATs".$nat_n];
73 		print TMP "NATs$nat_n\t$arr[1]\t$arr[2]\t$gene1{$arr[1]}[0][2]\t$gene1{$arr[1]}[0][0]\t$gene1{$arr[1]}[0][1]\t$arr[3]\t$arr[4]\t$gene1{$arr[3]}[0][2]\t$gene1{$arr[3]}[0][0]\t$gene1{$arr[3]}[0][1]\t$arr[5]\t$arr[11]\t$ns\t$ne\t$ns\t$ne\n";
74 		$nat_n++;
75 	}else{
76 		$ns=$gene1{$arr[1]}[0][0]+$arr[6]-1;
77 		$ne=$gene1{$arr[1]}[0][0]+$arr[7]-1;
78 		$ns2=$gene1{$arr[3]}[0][1]-$arr[9]+1;
79 		$ne2=$gene1{$arr[3]}[0][1]-$arr[8]+1;
80 		push @{$nat{$gene1{$arr[1]}[0][2]}},[$ns,$ne,$arr[5],$arr[11],"NATs$nat_n"."_1"];#start,end,class1,class2
81 		push @{$nat{$gene1{$arr[3]}[0][2]}},[$ns2,$ne2,$arr[5],$arr[11],"NATs$nat_n"."_2"];
82 		print TMP "NATs$nat_n\t$arr[1]\t$arr[2]\t$gene1{$arr[1]}[0][2]\t$gene1{$arr[1]}[0][0]\t$gene1{$arr[1]}[0][1]\t$arr[3]\t$arr[4]\t$gene1{$arr[3]}[0][2]\t$gene1{$arr[3]}[0][0]\t$gene1{$arr[3]}[0][1]\t$arr[5]\t$arr[11]\t$ns\t$ne\t$ns2\t$ne2\n";
83 		$nat_n++;
84 	}
85 }
86 close IN;
87 close TMP;
88 
89 my %repeat;
90 open IN,"<$opts{'r'}";
91 my $first=<IN>;
92 $first=<IN>;
93 $first=<IN>;
94 while (my $aline=<IN>) {
95 	chomp $aline;
96 	$aline=~s/^\s+//;
97 	my @tmp=split/\s+/,$aline;
98 	$tmp[4]=~s/chr0/Chr/;
99 	$tmp[4]=~s/chr/Chr/;
00 	push @{$repeat{$tmp[4]}},[$tmp[5],$tmp[6],$tmp[10]];#start,end,class
01 	#print "$tmp[4]\t$tmp[5]\t$tmp[6]\t$tmp[10]\n";
02 }
03 close IN;
04 
05 my %phase;
06 open IN,"<$opts{'p'}";
07 while (my $aline=<IN>) {
08 	chomp $aline;
09 	next if($aline=~/^\#/);
10 	my @tmp=split/\t/,$aline;
11 	if ($tmp[5]>=25) {
12 		$phase{$tmp[0]}=$tmp[5];
13 	}
14 }
15 close IN;
16 
17 my $filein=$opts{'i'};
18 my $fileout=$opts{'o'};
19 open IN,"<$filein"; #input file  
20 open OUT,">$fileout"; #output file  
21 while (my $aline=<IN>) {
22 	chomp $aline;
23 	if($aline=~/^\#/){
24 		print OUT "$aline\tPhase\tLong\tRepeatClass\tNatClass1\tNatClass2\tNatID\n";
25 		next;
26 	}
27 	my @tmp=split/\t/,$aline;
28 	my @inf=split/\:/,$tmp[0];
29 	my @pos=split/\-/,$inf[1];
30 	my $chr=$inf[0];
31 	$chr=~s/Chr0/Chr/;
32 	my $start=$pos[0];
33 	my $end=$pos[1];
34 	#=========Repeat============
35 	my @repeat;
36 	if (defined(@{$repeat{$chr}})) {
37 		my @r_array=sort {$a->[0] <=> $b->[0]} @{$repeat{$chr}};
38 		for (my $i=0;$i<@r_array ;$i++) {
39 			if ($start<$r_array[$i][0] && $end>$r_array[$i][0]) {
40 				push @repeat,$r_array[$i][2];
41 			}
42 			elsif($start>$r_array[$i][0] && $start<$r_array[$i][1]){
43 				push @repeat,$r_array[$i][2];
44 
45 			}
46 			elsif($end<$r_array[$i][0]){
47 				last;
48 			}
49 			else{next;}
50 		}
51 	}
52 	my $repeat;
53 	if (@repeat==0) {
54 		$repeat="\/";
55 	}
56 	else{
57 		$repeat=join ";",@repeat;
58 	}
59 	#=========nat===============
60 	my @nat1;#class 1
61 	my @nat2;#class 2
62 	my @nat;#nat ID
63 	#foreach my $chr (keys %nat) {
64 		my @n_array=sort {$a->[0] <=> $b->[0] } @{$nat{$chr}};
65 		for (my $i=0;$i<@n_array ;$i++) {
66 			if ($start<$n_array[$i][0] && $end>$n_array[$i][0]) {
67 				push @nat1,$n_array[$i][2];
68 				push @nat2,$n_array[$i][3];
69 				push @nat,$n_array[$i][4];
70 			}
71 			elsif($start>$n_array[$i][0] && $start<$n_array[$i][1]){
72 				push @nat1,$n_array[$i][2];
73 				push @nat2,$n_array[$i][3];
74 				push @nat,$n_array[$i][4];
75 			}
76 			elsif($end<$n_array[$i][0]){
77 				last;
78 			}
79 			else{next;}
80 		}
81 	#}
82 
83 	my $nat1;
84 	my $nat2;
85 	my $nat;
86 	if (@nat1==0) {
87 		$nat1="\/";
88 	}
89 	else{
90 		$nat1=join ";",@nat1;
91 	}
92 	if (@nat2==0) {
93 		$nat2="\/";
94 	}
95 	else{
96 		$nat2=join ";",@nat2;
97 	}
98 	if (@nat==0) {
99 		$nat="\/";
00 	}
01 	else{
02 		$nat=join ";",@nat;
03 	}
04 	#========phase==============
05 	my $phase="\/";
06 	if (defined($phase{$tmp[0]})) {
07 		$phase="phase";
08 	}
09 	#=========long===============
10 	my $long="\/";
11 	if ($tmp[1] eq "\>30nt" and $tmp[2]>=0.5) {
12 		$long="long";
13 	}
14 	print OUT "$aline\t$phase\t$long\t$repeat\t$nat1\t$nat2\t$nat\n";
15 }
16 
17 close IN;
18 close OUT;
19 
20 sub nearby{
21 	my @p=@_;
22 	my ($s,$e)=(0,0);
23 	if ($p[1]<$p[2]) {
24 		$s=$p[1];
25 		$e=$p[2];
26 	}else{
27 		$s=$p[3];
28 		$e=$p[0];
29 		}
30 	return ($s,$e);
31 }
32 
33 sub usage{
34 print <<"USAGE";
35 Version $version
36 Usage:
37 $0 -i -o -g -n -r -p -t -l
38 options:
39 -i input file
40 	-g gff file
41 	-n NATs file
42 	-r repeat file
43 	-p phase file
44 -o output file
45 -t nat output file
46 -l genelist output file
47 -h help
48 USAGE
49 exit(1);
50 }
51