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

Annotation of /FigKernelPackages/ExpressionDir.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.1 package ExpressionDir;
2 :    
3 :     use Data::Dumper;
4 :     use strict;
5 :     use SeedAware;
6 :     use File::Copy;
7 :     use File::Spec::Functions;
8 :     use base 'Class::Accessor';
9 :     use Carp;
10 :     use Fcntl ':seek';
11 :    
12 :     __PACKAGE__->mk_accessors(qw(genome_dir expr_dir error_dir));
13 :    
14 :     =head3 new
15 :    
16 :     my $edir = ExpressionDir->new($genome_dir);
17 :    
18 :     Create a new ExpressionDir object to be associated with the given genome directory.
19 :     If a subdirectory ExpressionData does not yet exist, one is created.
20 :    
21 :     =cut
22 :    
23 :     sub new
24 :     {
25 :     my($class, $genome_dir) = @_;
26 :    
27 :     my $edir = catfile($genome_dir, "ExpressionData");
28 :     if (! -d $edir)
29 :     {
30 :     mkdir($edir);
31 :     }
32 :     my $errdir = catfile($genome_dir, 'errors');
33 :     if (! -d $errdir)
34 :     {
35 :     mkdir($errdir);
36 :     }
37 :    
38 :    
39 :     my $self = {
40 :     genome_dir => $genome_dir,
41 :     expr_dir => $edir,
42 :     error_dir => $errdir,
43 :     };
44 :     return bless $self, $class;
45 :     }
46 :    
47 :     sub parse_probe_format_1
48 :     {
49 :     my($self, $in_file, $out_file) = @_;
50 :    
51 :     my($fh);
52 :    
53 :     open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
54 :     my $l = <$fh>;
55 :     chomp $l;
56 :     $l =~ s/\r//;
57 :     my @hdrs = split(/\t/, $l);
58 :     my %hdrs;
59 :     $hdrs{$hdrs[$_]} = $_ for 0..$#hdrs;
60 :    
61 :     my $x_col = $hdrs{"Probe X"};
62 :     my $y_col = $hdrs{"Probe Y"};
63 :     my $seq_col = $hdrs{"Probe Sequence"};
64 :     if (!(defined($x_col) && defined($y_col) && defined($seq_col)))
65 :     {
66 :     close($fh);
67 :     return undef;
68 :     }
69 :    
70 :     my $out;
71 :     open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
72 :    
73 :     while (<$fh>)
74 :     {
75 :     chomp;
76 :     s/\r//g;
77 :     my @flds = split(/\t/,$_);
78 :     my($x,$y,$seq);
79 :     $x = $flds[$x_col];
80 :     $y = $flds[$y_col];
81 :     $seq = $flds[$seq_col];
82 :     my $id = "$x\_$y";
83 :     print $out "$id\t$seq\n";
84 :     }
85 :     close($fh);
86 :     close($out);
87 :     return 1;
88 :     }
89 :    
90 :     sub parse_probe_format_2
91 :     {
92 :     my($self, $in_file, $out_file) = @_;
93 :    
94 :     my($fh);
95 :    
96 :     local $/ = "\n>";
97 :    
98 :     open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
99 :     my $l = <$fh>;
100 :     chomp $l;
101 :     $l =~ s/\r//;
102 :    
103 :     if ($l !~ /^>?\S+:(\d+):(\d+);\s+Interrogation_Position=\d+;\s+Antisense;\n([ACGT]+)/s)
104 :     {
105 :     close($fh);
106 :     return undef;
107 :     }
108 :     seek($fh, 0, SEEK_SET);
109 :    
110 :     my $out;
111 :     open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
112 :    
113 :     while (<$fh>)
114 :     {
115 :     chomp;
116 :    
117 :     if ($_ =~ /^>?\S+:(\d+):(\d+);\s+Interrogation_Position=\d+;\s+Antisense;\n([ACGT]+)/s)
118 :     {
119 :     if (length($3) < 15)
120 :     {
121 :     close($fh);
122 :     confess "Bad length at line $. of $in_file";
123 :     }
124 :     print $out "$1\_$2\t$3\n";
125 :     }
126 :     else
127 :     {
128 :     confess "Bad input at line $. of $in_file";
129 :     }
130 :     }
131 :     close($out);
132 :     close($fh);
133 :     return 1;
134 :     }
135 :    
136 :     #
137 :     # Our "native" format, used for passing through pre-parsed data.
138 :     #
139 :     sub parse_probe_format_3
140 :     {
141 :     my($self, $in_file, $out_file) = @_;
142 :    
143 :     my($fh);
144 :    
145 :     open($fh, "<", $in_file) or confess "Cannot open $in_file for reading: $!";
146 :     my $l = <$fh>;
147 :     chomp $l;
148 :     $l =~ s/\r//;
149 :    
150 :     if ($l !~ /^\d+_\d+\t[ACGT]+$/)
151 :     {
152 :     close($fh);
153 :     return undef;
154 :     }
155 :     seek($fh, 0, SEEK_SET);
156 :    
157 :     my $out;
158 :     open($out, ">", $out_file) or confess "Cannot open $out for writing: $!";
159 :    
160 :     while (<$fh>)
161 :     {
162 :     if ($_ =~ /^\d+_\d+\t[ACGT]+$/)
163 :     {
164 :     print $out $_;
165 :     }
166 :     else
167 :     {
168 :     confess "Bad input at line $. of $in_file";
169 :     }
170 :     }
171 :     close($out);
172 :     close($fh);
173 :     return 1;
174 :     }
175 :    
176 :     sub compute_probe_to_peg
177 :     {
178 :     my($self, $probes) = @_;
179 :    
180 :     my $my_probes = catfile($self->expr_dir, "probes.in");
181 :    
182 :     copy($probes, $my_probes) or confess "Cannot copy $probes to $my_probes: $!";
183 :    
184 :     my $probes_fasta = catfile($self->expr_dir, "probes.fasta");
185 :    
186 :     #
187 :     # Attempt to translate probe file.
188 :     #
189 :     my $success;
190 :     for my $meth (qw(parse_probe_format_1 parse_probe_format_2 parse_probe_format_3))
191 :     {
192 :     if ($self->$meth($my_probes, $probes_fasta))
193 :     {
194 :     print STDERR "Translated $probes to $probes_fasta using $meth\n";
195 :     $success = 1;
196 :     last;
197 :     }
198 :     else
199 :     {
200 :     print STDERR "Failed to translate $probes to $probes_fasta using $meth\n";
201 :     }
202 :     }
203 :     if (!$success)
204 :     {
205 :     confess "Could not translate $probes\n";
206 :     }
207 :    
208 :     my $peg_probe_table = catfile($self->expr_dir, 'peg.probe.table');
209 :     my $probe_occ_table = catfile($self->expr_dir, 'probe.occ.table');
210 :    
211 :     system_with_redirect([executable_for("make_probes_to_genes"),
212 :     $probes_fasta,
213 :     catfile($self->genome_dir, 'contigs'),
214 :     catfile($self->genome_dir, 'Features', 'peg', 'tbl'),
215 :     $peg_probe_table,
216 :     $probe_occ_table],
217 :     { stderr => catfile($self->expr_dir, 'problems') });
218 :    
219 :    
220 :     system_with_redirect([executable_for("remove_multiple_occurring_probes"),
221 :     $peg_probe_table,
222 :     $probe_occ_table],
223 :     { stdout => catfile($self->expr_dir, 'peg.probe.table.no.multiple') } );
224 :    
225 :     $self->make_missing_probes($peg_probe_table, $probes_fasta,
226 :     catfile($self->expr_dir, 'probe.no.match'));
227 :     }
228 :     sub make_missing_probes
229 :     {
230 :     my($self, $probe_table, $probes, $output) = @_;
231 :     open(MATCH,"<", $probe_table) or die "Cannot open $probe_table: $!";
232 :     open(PROBES,"<", $probes) or die "Cannot open $probes: $!";
233 :     open(OUTPUT, ">", $output) or die "Cannot open $output: $!";
234 :     my %locations;
235 :     while(<MATCH>)
236 :     {
237 :     chomp;
238 :     my($peg,$loc)=split "\t";
239 :     $locations{$loc} = $peg;
240 :     }
241 :    
242 :     while(<PROBES>)
243 :     {
244 :     chomp;
245 :     my($loc,$seq) = split "\t";
246 :     print OUTPUT $loc, "\n" if ! exists $locations{$loc};
247 :     }
248 :     close(MATCH);
249 :     close(PROBES);
250 :     close(OUTPUT);
251 :     }
252 :    
253 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3