File Coverage

blib/lib/GenOO/GeneCollection/Factory/GTF.pm
Criterion Covered Total %
statement 77 79 97.4
branch 24 32 75.0
condition 16 30 53.3
subroutine 10 10 100.0
pod 0 1 0.0
total 127 152 83.5


line stmt bran cond sub pod time code
1             # POD documentation - main docs before the code
2              
3             =head1 NAME
4              
5             GenOO::GeneCollection::Factory::GTF - Factory to create GeneCollection from a GTF file
6              
7             =head1 SYNOPSIS
8              
9             Creates GenOO::GeneCollection containing genes from a GTF file
10             Preferably use it through the generic GenOO::GeneCollection::Factory
11              
12             my $factory = GenOO::GeneCollection::Factory->new('GTF',{
13             file => 'sample.gtf'
14             });
15              
16             =head1 DESCRIPTION
17              
18             # An instance of this class is a concrete factory for the creation of a
19             # L<GenOO::GeneCollection> containing genes from a GTF file. It offers the
20             # method "read_collection" (as the consumed role requires) which returns
21             # the actual L<GenOO::GeneCollection> object in the form of
22             # L<GenOO::RegionCollection::Type::DoubleHashArray>. The latter is the
23             # implementation
24             # of the L<GenOO::RegionCollection> class based on the complex data
25             # structure L<GenOO::Data::Structure::DoubleHashArray>.
26              
27             =head1 EXAMPLES
28              
29             # Create a concrete factory
30             my $factory_implementation = GenOO::GeneCollection::Factory->new('GTF',{
31             file => 'sample.gtf'
32             });
33              
34             # Return the actual GenOO::GeneCollection object
35             my $collection = $factory_implementation->read_collection;
36             print ref($collection) # GenOO::GeneCollection::Type::DoubleHashArray
37              
38             =cut
39              
40             # Let the code begin...
41              
42             package GenOO::GeneCollection::Factory::GTF;
43             $GenOO::GeneCollection::Factory::GTF::VERSION = '1.5.2';
44              
45             #######################################################################
46             ####################### Load External modules #####################
47             #######################################################################
48 1     1   4345 use Modern::Perl;
  1         8  
  1         23  
49 1     1   186 use autodie;
  1         3  
  1         8  
50 1     1   5601 use Moose;
  1         3  
  1         7  
51 1     1   7640 use namespace::autoclean;
  1         2  
  1         18  
52              
53              
54             #######################################################################
55             ####################### Load GenOO modules #####################
56             #######################################################################
57 1     1   129 use GenOO::RegionCollection::Factory;
  1         3  
  1         75  
58 1     1   12 use GenOO::Transcript;
  1         4  
  1         42  
59 1     1   7 use GenOO::Gene;
  1         7  
  1         48  
60 1     1   7 use GenOO::Data::File::GFF;
  1         2  
  1         1178  
61              
62              
63             #######################################################################
64             ####################### Interface attributes ######################
65             #######################################################################
66             has 'file' => (
67             isa => 'Maybe[Str]',
68             is => 'ro'
69             );
70              
71              
72             #######################################################################
73             ########################## Consumed Roles #########################
74             #######################################################################
75             with 'GenOO::RegionCollection::Factory::Requires';
76              
77              
78             #######################################################################
79             ######################## Interface Methods ########################
80             #######################################################################
81             sub read_collection {
82 6     6 0 36455 my ($self) = @_;
83              
84 6         235 my $genes = $self->_read_gtf($self->file);
85              
86 6         362 return GenOO::RegionCollection::Factory->create('RegionArray', {
87             array => $genes
88             })->read_collection;
89             }
90              
91             #######################################################################
92             ######################### Private methods ##########################
93             #######################################################################
94             sub _read_gtf {
95 6     6   32 my ($self, $file)=@_;
96              
97 6         37 my %transcripts;
98             my %transcript_splice_starts;
99 6         0 my %transcript_splice_stops;
100 6         0 my %genes;
101              
102 6         240 my $gff = GenOO::Data::File::GFF->new(file => $file);
103              
104 6         119 while (my $record = $gff->next_record){
105 5808 100 100     161051 next unless (($record->feature eq 'exon') or ($record->feature eq 'start_codon') or ($record->feature eq 'stop_codon'));
      100        
106 3522 50       123735 my $tid = $record->attribute('transcript_id')
107             or die "transcript_id attribute must be defined\n";
108 3522 50       120520 my $gid = $record->attribute('gene_id')
109             or die "gene_id attribute must be defined\n";
110              
111 3522 50       95311 if ($record->strand == 0){
112 0         0 next;
113             }
114              
115             # Get transcript with id or create a new one. Update coordinates if required
116 3522         7920 my $transcript = $transcripts{$tid};
117 3522 100       8492 if (not defined $transcript) {
118 348         9577 $transcript = GenOO::Transcript->new(
119             id => $tid,
120             chromosome => $record->rname,
121             strand => $record->strand,
122             start => $record->start,
123             stop => $record->stop,
124             splice_starts => [$record->start], # will be re-written later
125             splice_stops => [$record->stop], # will be re-written later
126             );
127 348         3042 $transcripts{$tid} = $transcript;
128 348         1237 $transcript_splice_starts{$tid} = [];
129 348         1001 $transcript_splice_stops{$tid} = [];
130 348 100       1308 if (!exists $genes{$gid}) {
131 156         867 $genes{$gid} = [];
132             }
133 348         707 push @{$genes{$gid}}, $transcript;
  348         1055  
134             }
135             else {
136 3174 100       84800 $transcript->start($record->start) if ($record->start < $transcript->start);
137 3174 50       85664 $transcript->stop($record->stop) if ($record->stop > $transcript->stop);
138             }
139              
140 3522 100       97232 if ($record->feature eq 'exon') {
    100          
    50          
141 2982         4834 push @{$transcript_splice_starts{$tid}}, $record->start;
  2982         81003  
142 2982         4801 push @{$transcript_splice_stops{$tid}}, $record->stop;
  2982         79614  
143             }
144             elsif ($record->feature eq 'start_codon') {
145 270 100 33     7354 if ($record->strand == 1 and
    50 66        
      33        
      33        
146             (!defined $transcript->coding_start or
147             $record->start < $transcript->coding_start)) {
148              
149 168         4623 $transcript->coding_start($record->start);
150             }
151             elsif ($record->strand == -1 and
152             (!defined $transcript->coding_stop or
153             $record->stop > $transcript->coding_stop)) {
154              
155 102         2838 $transcript->coding_stop($record->stop);
156             }
157             }
158             elsif ($record->feature eq 'stop_codon') {
159 270 100 33     7196 if ($record->strand == 1 and
    50 66        
      33        
      33        
160             (!defined $transcript->coding_stop or
161             $record->stop > $transcript->coding_stop)) {
162              
163 168         4584 $transcript->coding_stop($record->stop);
164             }
165             elsif ($record->strand == -1 and
166             (!defined $transcript->coding_start or
167             $record->start < $transcript->coding_start)) {
168              
169 102         2803 $transcript->coding_start($record->start);
170             }
171             }
172             }
173              
174 6         204 foreach my $tid (keys %transcripts) {
175             $transcripts{$tid}->set_splice_starts_and_stops(
176 348         1450 $transcript_splice_starts{$tid}, $transcript_splice_stops{$tid});
177             }
178              
179 6         62 my @genes;
180 6         126 GENE: foreach my $gid (keys %genes) {
181 156         4040 my $gene = GenOO::Gene->new(name => $gid);
182 156         261 my @gene_transcripts = sort {$a->start <=> $b->start} @{$genes{$gid}};
  300         7862  
  156         702  
183 156         346 my $tr = shift @gene_transcripts;
184 156         522 $gene->add_transcript($tr);
185 156         3750 $tr->gene($gene);
186 156         331 foreach $tr (@gene_transcripts) {
187             # if transcrpt doesn't overlap previous ones then skip the gene.
188 192 50       746 if (!$gene->overlaps($tr)) {
189 0         0 next GENE;
190             }
191 192         637 $gene->add_transcript($tr);
192 192         4485 $tr->gene($gene);
193             }
194 156         406 push @genes, $gene;
195             }
196              
197 6         333 return \@genes;
198             }
199              
200             1;