[Bio] / FigKernelPackages / gjogff.pm Repository:
ViewVC logotype

Diff of /FigKernelPackages/gjogff.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1, Sun Jan 2 19:02:50 2011 UTC revision 1.2, Mon Jan 24 00:05:10 2011 UTC
# Line 7  Line 7 
7  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
8  #  Parse the descriptions of the features into key-value pairs  #  Parse the descriptions of the features into key-value pairs
9  #  #
10  #     ( \%features_by_contig, \@dna ) = read_gff( \*FH );  #     ( \%features_by_contig, \@dna ) = read_gff( \*FH, \%options );
11    #
12    #  Options:
13    #
14    #     contig_ids => \%ids
15    #     contig_ids => \@ids
16    #     dna        => \@dna
17  #  #
18  #  Features are:  #  Features are:
19  #  #
20  #     [ $type, $loc, $id, $name, \%keys_values ]  #     [ $type, $loc, $id, $name, \%keys_values ]
21  #  #
22  #     Location is Sapling style contig_beg±length  #     Location is Sapling style contig_begĀ±length
23  #  #
24  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
25  sub read_gff  sub read_gff
26  {  {
27        my $opts = @_ && $_[ 0] && ref $_[ 0] eq 'HASH' ? shift :
28                   @_ && $_[-1] && ref $_[-1] eq 'HASH' ? pop   :
29                                                          {};
30    
31        my $contig_ids = $opts->{ contig_ids } && ref $opts->{ contig_ids } eq 'HASH'  ? $opts->{ contig_ids }                                     :
32                         $opts->{ contig_ids } && ref $opts->{ contig_ids } eq 'ARRAY' ? { map { $_ => 1 } @{ $opts->{ contig_ids } } }            :
33                         $opts->{ contig_len } && ref $opts->{ contig_len } eq 'HASH'  ? $opts->{ contig_len }                                     :
34                         $opts->{ contig_len } && ref $opts->{ contig_len } eq 'ARRAY' ? { map { $_ => 1 } @{ $opts->{ contig_len } } }            :
35                         $opts->{ dna }        && ref $opts->{ dna }        eq 'ARRAY' ? { map { $_->[0] => length $_->[2] } @{ $opts->{ dna } } } :
36                                                                                         undef;
37    
38      my ( $fh, $close ) = input_filehandle( $_[0] );      my ( $fh, $close ) = input_filehandle( $_[0] );
39    
40      my $dna;      my $dna;
41      my %ftrs;      my %ftrs;
42      local $_;      local $_;
43    
44        my %no_contig;
45      while ( defined( $_ = <$fh> ) )      while ( defined( $_ = <$fh> ) )
46      {      {
47          if ( m/^##FASTA/i ) { $dna = 1; last }  #  GFF DNA introduction          if ( m/^##FASTA/i ) { $dna = 1; last }  #  GFF DNA introduction
# Line 39  Line 57 
57    
58      my @dna;      my @dna;
59      @dna = gjoseqlib::read_fasta( $fh ) if $dna;      @dna = gjoseqlib::read_fasta( $fh ) if $dna;
60        $contig_ids ||= { map { $_->[0] => 1 } @dna }  if @dna;
61    
62      close( $fh ) if $close;      close( $fh ) if $close;
63    
# Line 47  Line 66 
66      #     [ $type, $loc, $id, $name, $frame, \%keys_values ]      #     [ $type, $loc, $id, $name, $frame, \%keys_values ]
67      #      #
68    
69      my $features_by_contig = process_gff_features( \%ftrs );      my $features_by_contig = process_gff_features( \%ftrs, $contig_ids );
70    
71      wantarray ? ( $features_by_contig, \@dna ) : $features_by_contig;      wantarray ? ( $features_by_contig, \@dna ) : $features_by_contig;
72  }  }
# Line 102  Line 121 
121  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
122  #  Integrate and sort feature information:  #  Integrate and sort feature information:
123  #  #
124  #     \%features_by_contig = process_gff_features( \%features );  #     \%features_by_contig = process_gff_features( \%features, \%contig_ids );
125  #  #
126  #  Features are:  #  Features are:
127  #  #
128  #     [ $type, $loc, $id, $name, \%keys_values ]  #     [ $type, $loc, $id, $name, \%keys_values ]
129  #  #
130  #     Location is Sapling style contig_beg±length  #     Location is Sapling style contig_begĀ±length
131  #  #
132  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
133  sub process_gff_features  sub process_gff_features
134  {  {
135      my $ftrs = shift || {};      my $ftrs = shift || {};
136      my @ftrs;      my $contig_ids = $_[0] && ref $_[0] eq 'HASH' ? shift : {};
137        my $have_ids   = keys %$contig_ids;
138        my $have_len   = ( values %$contig_ids )[0] > 1;
139    
140        my @ftrs;
141        my %missing_contig;              #  Not currently reporting
142      foreach my $id ( keys %$ftrs )      foreach my $id ( keys %$ftrs )
143      {      {
144          my $types = $ftrs->{ $id };          my $types = $ftrs->{ $id };
# Line 127  Line 150 
150              #              0     1       2       3       4     5       6         7          8              #              0     1       2       3       4     5       6         7          8
151              #  $datum = [ $id, $name, $contig, $left, $right, $dir, $frame, \%keys_values, $mid ]              #  $datum = [ $id, $name, $contig, $left, $right, $dir, $frame, \%keys_values, $mid ]
152    
153              my %names = map { $_->[1] => 1 } @$data;              #
154              if ( keys %names > 1 )              #  Contig anlysis:
155                #
156                #  Number of contigs  contig_data   ids_matching      response
157                #  -----------------------------------------------------------------
158                #          1              yes            0         contig missing
159                #          1              yes            1              okay
160                #          1              no                            okay
161                #         >1              yes            0         contig missing
162                #         >1              yes            1           keep the 1
163                #         >1              yes           >1         multiple contigs
164                #         >1              no                       multiple contigs
165                #  -----------------------------------------------------------------
166                #
167    
168                my %contigs     = map { $_->[2] => 1 } @$data;
169                my @contigs     = keys %contigs;
170                my @have_contig = grep { $contig_ids->{ $_ } } @contigs;
171                if ( @contigs == 1 )
172              {              {
173                  print STDERR "Feature '$id' of type '$type' has multiple names: ",                  $missing_contig{ $contigs[0] } = 1 if $have_ids && ! @have_contig;
174                                join( ', ', sort keys %names ), "\n";              }
175                elsif ( @have_contig == 1 )
176                {
177                    #  Filter raw data to the contig we have:
178                    @$data = grep { $_->[2] eq $have_contig[0] } @$data;
179                }
180                elsif ( $have_ids && ! @have_contig )
181                {
182                    $missing_contig{ $contigs[0] } = 1;
183                }
184                else
185                {
186                    print STDERR "Feature '$id' of type '$type' has multiple contigs: ",
187                                  join( ', ', sort @contigs ), "\n";
188                  $okay = 0;                  $okay = 0;
189              }              }
190    
191              my %contigs = map { $_->[2] => 1 } @$data;              my %names = map { $_->[1] => 1 } @$data;
192              if ( keys %contigs > 1 )              if ( keys %names > 1 )
193              {              {
194                  print STDERR "Feature '$id' of type '$type' has multiple contigs: ",                  print STDERR "Feature '$id' of type '$type' has multiple names: ",
195                                join( ', ', sort keys %contigs ), "\n";                                join( ', ', sort keys %names ), "\n";
196                  $okay = 0;                  $okay = 0;
197              }              }
198    
# Line 189  Line 242 
242                  }                  }
243              }              }
244    
245              #  The JGI files record 0 frame in all forward genes, so kill this.              if ( $have_len && $contig_ids->{ $cont } )
246                {
247                    if ( max( map { $_->[4] } @$data ) > $contig_ids->{ $cont } )
248                    {
249                        print STDERR "Feature $type: $id $name extends beyond end of the contig DNA.\n";
250                        next;
251                    }
252                }
253    
254                #  The JGI files list frame 0 in all forward genes, so kill this.
255              #  $keys_values{ codon_start } = [ $frame =~ /\d/ ? $frame+1 : 1 ];              #  $keys_values{ codon_start } = [ $frame =~ /\d/ ? $frame+1 : 1 ];
256    
257              push @ftrs, [ [ $type, $loc, $id, $name, \%keys_values ], $cont, $mid ];              push @ftrs, [ [ $type, $loc, $id, $name, \%keys_values ], $cont, $mid ];
# Line 202  Line 264 
264          push @{ $features_by_contig{ $_->[1] } }, $_->[0];          push @{ $features_by_contig{ $_->[1] } }, $_->[0];
265      }      }
266    
267        if ( $have_ids && keys %missing_contig )
268        {
269            print STDERR "No data for contigs:\n";
270            print STDERR map { "    $_\n" } sort keys %missing_contig;
271        }
272    
273      \%features_by_contig;      \%features_by_contig;
274  }  }
275    
276    
277    sub max { my $max = shift; foreach ( @_ ) { $max = $_ if $_ > $max } $max }
278    
279    
280  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
281  #  Deal with %.. escaped text  #  Deal with %.. escaped text
282  #  #

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3