[Bio] / FigKernelScripts / svr_genbank_to_seed.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/svr_genbank_to_seed.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : gdpusch 1.1 # -*- perl -*-
2 :     # This is a SAS Component.
3 :     ########################################################################
4 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
5 :     # for Interpretations of Genomes. All Rights Reserved.
6 :     #
7 :     # This file is part of the SEED Toolkit.
8 :     #
9 :     # The SEED Toolkit is free software. You can redistribute
10 :     # it and/or modify it under the terms of the SEED Toolkit
11 :     # Public License.
12 :     #
13 :     # You should have received a copy of the SEED Toolkit Public License
14 :     # along with this program; if not write to the University of Chicago
15 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
16 :     # Genomes at veronika@thefig.info or download a copy from
17 :     # http://www.theseed.org/LICENSE.TXT.
18 :     ########################################################################
19 :    
20 :     # usage: svr_genbank_to_seed OrgDir genbank.file
21 : gdpusch 1.4 my $usage = "svr_genbank_to_seed OrgDir genbank.file";
22 : gdpusch 1.1
23 :     use strict;
24 :     use warnings;
25 :    
26 :     use FIGV;
27 :     use gjogenbank;
28 :     use gjoseqlib;
29 :    
30 :     use File::Path;
31 :     use Data::Dumper;
32 :    
33 : gdpusch 1.4 my $orgID;
34 :     my $org_dir;
35 :     my $genbank_file;
36 :     my $extension = 0;
37 :     use Getopt::Long;
38 :     my $rc = GetOptions("orgID=s" => \$orgID,
39 :     "orgDir=s" => \$org_dir,
40 :     "gbk" => \$genbank_file,
41 :     "extension=i" => \$extension,
42 :     );
43 :     if (! $rc) { print STDERR "$usage\n"; exit }
44 :    
45 :    
46 :     if (!defined($org_dir) && !defined($genbank_file) && (@ARGV == 2)) {
47 :     ($org_dir, $genbank_file) = @ARGV;
48 :     }
49 :    
50 :     if (!defined($orgID) && ($org_dir =~ m/(\d+\.\d+)$/o)) {
51 :     $orgID = $1;
52 :     }
53 :    
54 :    
55 :     if (-d $org_dir) {
56 :     warn "WARNING: orgDir='$org_dir' already exists\n";
57 :     }
58 :     else {
59 :     ($rc = mkpath($org_dir)) || die qq(Could not create orgDir='$org_dir', rc=$rc);
60 :     }
61 :    
62 :     my $accession = gjogenbank::parse_next_genbank($genbank_file);
63 :     my $db_xref = $accession->{FEATURES}->{source}->[0]->[1]->{db_xref}->[0];
64 :     if (!defined($orgID)) {
65 :     if (defined($db_xref)) {
66 :     if ($db_xref =~ m/^taxon:(\d+)$/) {
67 :     $orgID = "$1.$extension";
68 :     warn "WARNING: Using orgID='$orgID' based on GenBank db_xref='$db_xref'";
69 :     }
70 :     else {
71 :     die "ERROR: No orgID and could not parse db_xref='$db_xref' --- aborting";
72 :     }
73 :     }
74 :     else {
75 :     die "ERROR: No orgID and no db_xref --- aborting";
76 :     }
77 :     }
78 :    
79 :     if (defined($orgID)) {
80 :     open(GENOME_ID, q(>), "$org_dir/GENOME_ID")
81 :     || die "Could not write-open file '$org_dir/GENOME_ID'";
82 :     print GENOME_ID ($orgID, "\n");
83 :     close(GENOME_ID);
84 :     }
85 :     else {
86 :     die "ERROR: no orgID found on command-line, at end of orgDir path, or within GenBank file --- aborting";
87 :     }
88 :    
89 : gdpusch 1.1
90 :     my $figV = FIGV->new($org_dir);
91 :    
92 :     my $contigs_file = $org_dir.q(/contigs);
93 :     open(my $contigs_fh, q(>), $contigs_file)
94 :     || die qq(Could not write-open contigs file \'$contigs_file\');
95 :    
96 : gdpusch 1.4 do {
97 : gdpusch 1.2 my $contig_id = $accession->{ACCESSION}->[0];
98 : gdpusch 1.1 my $contig_dna = $accession->{SEQUENCE};
99 : gdpusch 1.4
100 :     my $this_xref = $accession->{FEATURES}->{source}->[0]->[1]->{db_xref}->[0];
101 :     if (defined($this_xref) && ($this_xref ne $db_xref)) {
102 :     warn "WARNING: db_xref mismatch for contig='$contig_id': '$this_xref' differs from '$db_xref'\n";
103 :     }
104 :    
105 : gdpusch 1.1 $figV->display_id_and_seq( $contig_id, \$contig_dna, $contigs_fh);
106 :    
107 :     foreach my $cds (@ { $accession->{FEATURES}->{CDS} }) {
108 : gdpusch 1.3 # die Dumper($cds);
109 : gdpusch 1.1 my $gb_loc = gjogenbank::location( $cds, $accession );
110 :     my $locus = gjogenbank::genbank_loc_2_seed($contig_id, $gb_loc);
111 :     my $func = gjogenbank::product( $cds ) || q();
112 :     my $translation = gjogenbank::CDS_translation($cds);
113 : gdpusch 1.3 my $pseudo = defined( $cds->[1]->{pseudo}->[0] );
114 : gdpusch 1.1
115 :     my $gene_name = defined($cds->[1]->{gene}->[0]) ? $cds->[1]->{gene}->[0] : q();
116 :     my $locus_tag = defined($cds->[1]->{locus_tag}->[0]) ? $cds->[1]->{locus_tag}->[0] : q();
117 :     my $protein_id = defined($cds->[1]->{protein_id}->[0]) ? $cds->[1]->{protein_id}->[0] : q();
118 :    
119 :     my @db_xrefs = defined($cds->[1]->{db_xref}->[0]) ? @ { $cds->[1]->{db_xref} } : ();
120 :    
121 :     my @gi_nums = map { m/GI\:(\d+)/o ? (q(gi|).$1) : () } @db_xrefs;
122 :     my @gene_nums = map { m/GeneID\:(\d+)/o ? (q(GeneID|).$1) : () } @db_xrefs;
123 :    
124 :     my @aliases = grep { $_ } ($gene_name, $locus_tag, $protein_id, @gi_nums, @db_xrefs, @gene_nums);
125 : gdpusch 1.2 my $aliases = join(q(,), @aliases);
126 : gdpusch 1.1
127 : gdpusch 1.3 if ($locus) {
128 :     my $fid;
129 :     if ($translation) {
130 : gdpusch 1.4 if ($fid = $figV->add_feature(q(Initial Import), $orgID, q(peg), $locus, $aliases, $translation)) {
131 : gdpusch 1.3 if ($func) {
132 :     $figV->assign_function($fid, q(master:Initial Import), $func);
133 :     }
134 :     }
135 :     }
136 :     elsif ($pseudo) {
137 :     my $sequence = gjogenbank::ftr_seq( $cds, $contig_dna);
138 : gdpusch 1.4 if ($fid = $figV->add_feature(q(Initial Import), $orgID, q(pseudo), $locus, $aliases, $sequence)) {
139 : gdpusch 1.3 if ($func) {
140 :     if ($func !~ m/pseudogene/i) {
141 : gdpusch 1.4 $func .= " \# pseudogene";
142 : gdpusch 1.3 }
143 :    
144 :     $figV->assign_function($fid, q(master:Initial Import), $func);
145 :     }
146 :     }
147 :     }
148 :     else {
149 :     die (qq(Could not add feature\n), Dumper($cds));
150 :     }
151 :     }
152 :     else {
153 : gdpusch 1.4 warn (qq(Could not parse CDS feature in accession '$contig_id':\n),
154 :     Dumper($cds),
155 :     qq(\n));
156 : gdpusch 1.3 }
157 :     }
158 :    
159 :     foreach my $rna (map { $_ ? @$_ : () }
160 :     map { my $x = $accession->{FEATURES}->{$_}
161 : gdpusch 1.4 } qw( rRNA tRNA misc_RNA ncRNA )
162 :     ) {
163 : gdpusch 1.3 # die Dumper($rna);
164 :    
165 :     my $gb_loc = gjogenbank::location( $rna, $accession );
166 :     my $locus = gjogenbank::genbank_loc_2_seed($contig_id, $gb_loc);
167 :     my $func = gjogenbank::product( $rna ) || q();
168 :     my $sequence = gjogenbank::ftr_seq( $rna, $contig_dna);
169 :    
170 :     my $gene_name = defined($rna->[1]->{gene}->[0]) ? $rna->[1]->{gene}->[0] : q();
171 :     my $locus_tag = defined($rna->[1]->{locus_tag}->[0]) ? $rna->[1]->{locus_tag}->[0] : q();
172 :    
173 :     my @db_xrefs = defined($rna->[1]->{db_xref}->[0]) ? @ { $rna->[1]->{db_xref} } : ();
174 :    
175 :     my @gi_nums = map { m/GI\:(\d+)/o ? (q(gi|).$1) : () } @db_xrefs;
176 :     my @gene_nums = map { m/GeneID\:(\d+)/o ? (q(GeneID|).$1) : () } @db_xrefs;
177 :    
178 :     my @aliases = grep { $_ } ($gene_name, $locus_tag, @gi_nums, @db_xrefs, @gene_nums);
179 :     my $aliases = join(q(,), @aliases);
180 :    
181 :     if ($locus && defined($func) && $sequence) {
182 : gdpusch 1.4 if (my $fid = $figV->add_feature(q(Initial Import), $orgID, q(rna), $locus, $aliases, $sequence)) {
183 : gdpusch 1.1 if ($func) {
184 : gdpusch 1.3 if ($func =~ m/23S\s+(ribosomal)?\s+RNA/i) {
185 : gdpusch 1.4 $func = "LSU rRNA \#\# 23S rRNA, large subunit ribosomal RNA";
186 : gdpusch 1.3 }
187 :     elsif ($func =~ m/16S\s+(ribosomal)?\s+RNA/i) {
188 : gdpusch 1.4 $func = "SSU rRNA \#\# 16S rRNA, small subunit ribosomal RNA";
189 : gdpusch 1.3 }
190 :     elsif ($func =~ m/5S\s+(ribosomal)?\s+RNA/i) {
191 : gdpusch 1.4 $func = "5S rRNA \#\# 5S ribosomal RNA";
192 : gdpusch 1.3 }
193 :    
194 : gdpusch 1.1 $figV->assign_function($fid, q(master:Initial Import), $func);
195 :     }
196 :     }
197 :     else {
198 : gdpusch 1.3 die (qq(Could not add feature\n), Dumper($rna));
199 : gdpusch 1.1 }
200 :     }
201 :     else {
202 : gdpusch 1.4 warn (qq(Could not parse RNA feature in accession '$contig_id': locus='$locus', func='$func', sequence='$sequence'\n),
203 :     Dumper($rna),
204 :     qq(\n));
205 : gdpusch 1.1 }
206 :     }
207 : gdpusch 1.4 } until (not defined($accession = gjogenbank::parse_next_genbank($genbank_file)));
208 : gdpusch 1.1 exit(0);

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3