[Bio] / Babel / bin / build_md52id_for_ncbi.pl Repository:
ViewVC logotype

Annotation of /Babel/bin/build_md52id_for_ncbi.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (view) (download) (as text)

1 : wilke 1.1 use strict;
2 :     use warnings;
3 :     use vars qw($opt_f $opt_d $opt_t $opt_s);
4 :     use Getopt::Std;
5 :     use Digest::MD5;
6 :     #use Bio::SearchIO;
7 :    
8 :    
9 :    
10 :     getopts('f:d:t:s:');
11 :    
12 :     # opt_t target directory
13 :     # opt_d destination directory
14 :    
15 :     my @files;
16 :     my $destination_dir = $opt_d || "/tmp/";
17 :     my $source = $opt_s || '';
18 :    
19 :     # parameter check
20 :    
21 :     unless (-d $destination_dir){
22 :     print STDERR "No destination directory $destination_dir\n";
23 :     exit;
24 :     }
25 :    
26 :     if ($opt_t and -d $opt_t){
27 :    
28 :     my @list =`find $opt_d -name fasta`;
29 :     foreach my $var (@list){
30 :     chomp $var;
31 :     push @files, $var;
32 :     }
33 :     print STDERR join "\n" , @files , "\n";
34 :    
35 :     }
36 :     elsif($opt_f and -f $opt_f){
37 :    
38 :     push @files , $opt_f;
39 :     }
40 :    
41 :    
42 :    
43 :     # read files
44 :     my $id_hash = {};
45 :     my $source_hash = {};
46 :    
47 :     $id_hash = read_fasta( $id_hash , @files , $source_hash);
48 :     open(SOURCES , ">$destination_dir/sources") or die "Can't open file sources\n";
49 :     foreach my $s (keys %$source_hash){
50 :     print SOURCES "$s\t" . $source_hash->{$s} . "\n";
51 :     }
52 :     close(SOURCES);
53 :    
54 :     # my $md5 = Digest::MD5::md5_hex( uc $sequence );
55 :    
56 :    
57 :     sub read_fasta{
58 :     my ($hash , $files , $source_hash) = @_;
59 :    
60 :     my $default = $/;
61 :    
62 :    
63 :     open(NR , ">$destination_dir/md5.fasta") or die "Can't open nr\n";
64 :     open(FILE , ">$destination_dir/md52id2func") or die "Can't open md2id2func";
65 :    
66 :    
67 :     open(FUNC , ">$destination_dir/assigned_functions") or die "Can't open file assigned_functions\n";
68 :     open(ORG , ">$destination_dir/org.table") or die "Can't open file org.table\n";
69 :     open(FAS , ">$destination_dir/fasta") or die "Can't open fasta\n";
70 :     open(ALL , ">$destination_dir/$source.faa") or die "Can't open $source.faa\n";
71 :     open(GI , ">$destination_dir/gi2id") or die "Can't open gi2id\n";
72 :    
73 :     foreach my $file (@files){
74 :    
75 :     my ($path,$t_source) = $file =~ /([\w\/]+)\/(\w+)\/[\w\.]+$/;
76 :     $source = $t_source unless ($source);
77 :    
78 :     print "Reading $file\n";
79 :     print "$source\n$path\n";
80 :    
81 :     my $ids = {};
82 :    
83 :    
84 :     open(FASTA , $file ) or die "Can't open file $file \n!";
85 :     # set line end
86 :     $/="\n>";
87 :    
88 :     my $count = 0;
89 :     while(my $line = <FASTA>){
90 :    
91 :     my @lines = split "\n" , $line;
92 :     my $end = pop @lines;
93 :     my $id_line = shift @lines;
94 :     my $fasta = join "" , @lines;
95 :    
96 :     next unless ($fasta);
97 :     my $md5 = Digest::MD5::md5_hex( uc $fasta );
98 :     my @entries;
99 :     if ($id_line =~/>/){
100 :     @entries = split ">gi" , $id_line ;
101 :     }
102 :     else{
103 :     push @entries , $id_line ;
104 :     }
105 :    
106 :     unless ($hash->{ $md5 }) { $hash->{ $md5 } = 0} ;
107 :    
108 :     # print STDERR $id_line , "\n\n";
109 :    
110 :     foreach my $entry (@entries){
111 :    
112 :     next unless ($entry);
113 :    
114 :     my ($gi , $source , $id) = $entry =~/^g{0,1}i{0,1}\|(\d+)\|(\w+)\|([^\s]+)/;
115 :    
116 :     my ($func , $org) = $entry =~/[^\s]+\s([^\[]+)\[?([^\]\[]*)\]?/;
117 :    
118 :     print GI "$gi\t$source|$id\n";
119 :    
120 :     my @ids = split /\|/ , $id ;
121 :     push @ids , $gi ;
122 :    
123 :     #print STDERR "$gi \t $source \t @ids \t $func \t $org\n";
124 :    
125 :     $source_hash->{$source}++;
126 :    
127 :    
128 :     unless($org){
129 :     # print STDERR "Error , no organism name\n";
130 :     # print STDERR $id_line , "\n";
131 :     $org = "organism not parsed/found in ncbi nr"
132 :    
133 :     }
134 :     unless ($gi){
135 :     print STDERR "No GI:\t$entry\n";
136 :     next;
137 :     }
138 :     unless ($id){
139 :     print STDERR "No ID:\t$entry\n";
140 :     next;
141 :     }
142 :    
143 :     # next if ($id =~/^sp/);
144 :    
145 :    
146 :    
147 :     foreach my $id (@ids){
148 :    
149 :     print FILE "$md5\t$id\t" . $func ."\t". $org ."\t$source\n";
150 :    
151 :     # for build nr
152 :     my $prefix = $source ;
153 :     $prefix = "gi" if ($id =~/^(\d+)$/);
154 :     print ALL ">$prefix|$id\t$func\t$org\n$fasta\n";
155 :     print FAS ">$prefix|$id\n$fasta\n";
156 :     print FUNC "$prefix|$id\t$func\n";
157 :     print ORG "$prefix|$id\t$org\n";
158 :     }
159 :    
160 :     unless ( $hash->{ $md5 } ) {
161 :     print NR ">lcl|$md5\n$fasta\n";
162 :     $hash->{ $md5 }++;
163 :     }
164 :     $count++;
165 :     #exit if ($count > 3);
166 :     }
167 :     }
168 :    
169 :     close(FASTA);
170 :    
171 :     # set line end back to default
172 :     $/=$default;
173 :     }
174 :    
175 :     close(FILE);
176 :     close(NR);
177 :    
178 :     close(ORG);
179 :     close(FUNC);
180 :     close(FAS);
181 :    
182 :     $/=$default;
183 :     return $hash;
184 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3