Mercurial > repos > mahtabm > ensembl

comparison variant_effect_predictor/Bio/EnsEMBL/Funcgen/Importer.pm @ 0:1f6dce3d34e0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author mahtabm
date Thu, 11 Apr 2013 02:01:53 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:1f6dce3d34e0
1
2 =head1 LICENSE
3
4 Copyright (c) 1999-2012 The European Bioinformatics Institute and
5 Genome Research Limited. All rights reserved.
6
7 This software is distributed under a modified Apache license.
8 For license details, please see
9
10 http://www.ensembl.org/info/about/code_licence.html
11
12 =head1 CONTACT
13
14 Please email comments or questions to the public Ensembl
15 developers list at <ensembl-dev@ebi.ac.uk>.
16
17 Questions may also be sent to the Ensembl help desk at
18 <helpdesk@ensembl.org>.
19
20 =head1 NAME
21
22 Bio::EnsEMBL::Funcgen::Importer
23
24 =head1 SYNOPSIS
25
26 my $imp = Bio::EnsEMBL::Funcgen::Importer->new(%params);
27 $imp->register_experiment();
28
29
30 =head1 DESCRIPTION
31
32 B<This program> is the main class coordinating import of tiling array design and experimental data.
33 It utilises several underlying parser classes specific to array vendor or import file type.
34
35 =cut
36
37 ################################################################################
38
39 package Bio::EnsEMBL::Funcgen::Importer;
40
41 use Bio::EnsEMBL::Funcgen::Utils::EFGUtils qw(get_date open_file run_system_cmd);
42 use Bio::EnsEMBL::Utils::Exception qw( throw deprecate );
43 use Bio::EnsEMBL::Utils::Argument qw( rearrange );
44 use Bio::EnsEMBL::Funcgen::Experiment;
45 use Bio::EnsEMBL::Funcgen::DBSQL::DBAdaptor;
46 use Bio::EnsEMBL::DBSQL::DBAdaptor;
47
48 use File::Path;
49 use strict;
50 use vars qw(@ISA);
51
52
53 ################################################################################
54
55 =head2 new
56
57 Description : Constructor method
58
59 Arg [1] : hash containing optional attributes:
60 -name Name of Experiment(dir)
61 -format of array e.g. Tiled(default)
62 -vendor name of array vendor
63 -description of the experiment
64 -pass DB password
65 -host DB host
66 -user DB user
67 -port DB port
68 -registry_host Host to load registry from
69 -registry_port Port for registry host
70 -registry_user User for registry host
71 -registry_pass Password for registry user
72 -ssh Flag to set connection over ssh via forwarded port to localhost (default = 0); remove?
73 -group name of experimental/research group
74 ` -location of experimental/research group
75 -contact e/mail address of primary contact for experimental group
76 -species
77 -assembly Genome assembly version i.e. 36 for NCBI36
78 -recover Recovery flag (default = 0)
79 -data_dir Root data directory (default = $ENV{'EFG_DATA'})
80 -output_dir review these dirs ???????
81 -input_dir ?????????
82 -import_dir ???????
83 -norm_dir ??????
84 -fasta dump FASTA flag (default =0)
85 -array_set Flag to treat all chip designs as part of same array (default = 0)
86 -array_name Name for array set
87 -array_file Path of array file to import for sanger ENCODE array
88 -result_set_name Name to give the raw and normalised result sets (default uses experiment and analysis name)
89 -norm_method Normalisation method (Nimblegen default = VSN_GLOG);
90 -dbname Override for autogeneration of funcgen dbaname
91 -reg_config path to local registry config file (default = ~/ensembl.init || undef)
92 -design_type MGED term (default = binding_site_identification) get from meta/MAGE?
93
94 -farm Flag to submit jobs to farm e.g. normalisation jobs
95 -batch_job Flag to signify that this Importer is running as a prepared batch/farm job
96 -prepared Flag to signify result files have been previously imported in prepare mode
97 and file names will differ to those record in InputSubset
98
99
100 #-use_defaults This changes some mandatory parameters to optional, instead using either DEFAULT or the input file name for the following options -name, -input_set, -feature_type, -cell_type etc ???
101
102 -verbose
103 ReturnType : Bio::EnsEMBL::Funcgen::Importer
104 Example : my $Exp = Bio::EnsEMBL::Importer->new(%params);
105 Exceptions : throws if mandatory params are not set or DB connect fails
106 Caller : General
107 Status : Medium - potential for %params names to change, remove %attrdata?
108
109 =cut
110
111 ################################################################################
112
113 sub new{
114 my ($caller) = shift;
115
116 #my $reg = "Bio::EnsEMBL::Registry";
117 my $class = ref($caller) || $caller;
118
119 #$user, $host, $port, $pass, $dbname, $db, $assm_version, $reg_config, $reg_db, $ucsc_coords, $verbose, $release, $reg_host, $reg_port, $reg_user, $reg_pass
120 #$ftype_name, $ctype_name,
121
122 my ($name, $format, $vendor, $group, $location, $contact,
123 $array_name, $array_set, $array_file, $data_dir, $result_files,
124 $exp_date, $desc, $design_type, $output_dir, $input_dir,
125 $batch_job, $farm, $prepared,
126 $norm_method, $old_dvd_format, $parser_type)
127 = rearrange(['NAME', 'FORMAT', 'VENDOR', 'GROUP', 'LOCATION', 'CONTACT',
128 'ARRAY_NAME', 'ARRAY_SET', 'ARRAY_FILE', 'DATA_DIR', 'RESULT_FILES',
129 'EXPERIMENT_DATE', 'DESCRIPTION',
130 'DESIGN_TYPE', 'OUTPUT_DIR', 'INPUT_DIR', #to allow override of defaults
131 'BATCH_JOB', 'FARM', 'PREPARED', 'NORM_METHOD',
132 'OLD_DVD_FORMAT', 'PARSER'], @_);
133
134
135
136 #### Define parent parser class based on vendor
137 throw("Mandatory argument -vendor not defined") if ! defined $vendor;
138
139 #This will override the default Vendor Parser type
140 #Evals simply protect from messy errors if parser type not found
141 my $parser_error;
142 my $vendor_parser = ucfirst(lc($vendor));
143
144
145 #WARNING evaling these parsers to enable pluggability hides errors in parser
146 #use a perl -MBio::EnsEMBL::Funcgen::Parsers:ParserType to debug
147 #get rid of all this case guessing and force correct parser name usage?
148
149
150 #WARNING
151 #Dynamic setting of ISA in this way reports the resultant object as Importer, when
152 #some throws/methods are actually in other base/custom Parsers
153 #This can seem a little counterintuitive, but allows plugability
154 #With out the need for separate control scripts
155
156
157 #Change this to be set and required/detected in the parse_and_import.pl script
158 #Then we can have Importer.pm as the base class and get rid of this.
159 #as well as set_config methods?
160
161
162 eval {require "Bio/EnsEMBL/Funcgen/Parsers/${vendor_parser}.pm";};
163
164 if ($@) {
165 #Don't warn/throw yet as we might have a standard parser format
166
167 $parser_error .= "There is no valid parser for the vendor your have specified:\t".$vendor.
168 "\nMaybe this is a typo or maybe you want to specify a default import format using the -parser option\n".$@;
169 }
170
171
172
173 if (defined $parser_type) {
174
175 #try normal case first
176 eval {require "Bio/EnsEMBL/Funcgen/Parsers/${parser_type}.pm";};
177
178 if ($@) {
179 $parser_type = ucfirst(lc($parser_type));
180
181 #Now eval the new parser
182 eval {require "Bio/EnsEMBL/Funcgen/Parsers/${parser_type}.pm";};
183
184 if ($@) {
185
186 #Might be no default
187 my $txt = "There is no valid parser for the -parser format your have specified:\t".$parser_type."\n";
188
189 if (! $parser_error) {
190 $txt .= "Maybe this is a typo or maybe you want run with the default $vendor_parser parser\n";
191 }
192
193 throw($txt.$@);
194 }
195
196 #warn about over riding vendor parser here
197 if (! $parser_error) {
198 #Can't log this as we haven't blessed the Helper yet
199 warn("WARNING\t::\tYou are over-riding the default ".$vendor." parser with -parser ".$parser_type);
200 }
201 }
202 } else {
203 throw($parser_error) if $parser_error;
204 $parser_type = $vendor_parser;
205 }
206
207
208 #we should now really set parser_type as an attrtibute?
209 unshift @ISA, 'Bio::EnsEMBL::Funcgen::Parsers::'.$parser_type;
210 #change this to be called explicitly from the load script?
211
212 #### Create object from parent class
213
214 my $self = $class->SUPER::new(@_);
215
216 #### Set vars and test minimum mandatory params for any import type
217
218 $self->{'name'} = $name || throw('Mandatory param -name not met'); #This is not mandatory for array design import
219 ## $self->{'user'} = $user || $ENV{'EFG_WRITE_USER'};
220 $self->vendor(uc($vendor)); #already tested
221 $self->{'format'} = uc($format) || 'TILED'; #remove default?
222 $self->group($group) if $group;
223 $self->location($location) if $location;
224 $self->contact($contact) if $contact;
225 ## $species || throw('Mandatory param -species not met');
226 $self->array_name($array_name) if $array_name;
227 $self->array_set($array_set) if $array_set;
228 $self->array_file($array_file) if $array_file;
229 $self->{'data_dir'} = $data_dir || $ENV{'EFG_DATA'};
230 $self->result_files($result_files)if $result_files;
231 $self->experiment_date($exp_date) if $exp_date;
232 $self->description($desc) if $desc; #experiment
233 ## $self->feature_set_description($fset_desc) if $fset_desc;
234
235 #$assm_version || throw('Mandatory param -assembly not met');
236 #Only required if setting DB by params e.g db not passed or generated from reg
237 #i.e. most of the time
238 #Why was this made mandatory?
239 #Default to dnadb||efgdb assm from the dbname
240
241 $self->{'design_type'} = $design_type || 'binding_site_identification'; #remove default?
242 $self->{'output_dir'} = $output_dir if $output_dir; #config default override
243 $self->{'input_dir'} = $input_dir if $input_dir; #config default override
244 $self->farm($farm) if $farm;
245 $self->batch_job($batch_job);
246 $self->prepared($prepared);
247 ## $self->{'ssh'} = $ssh || 0;
248 ## $self->{'_dump_fasta'} = $fasta || 0;
249 #$self->{'recover'} = $recover || 0; Now in BaseImporter
250 #check for ~/.ensembl_init to mirror general EnsEMBL behaviour
251 ## $self->{'reg_config'} = $reg_config || ((-f "$ENV{'HOME'}/.ensembl_init") ? "$ENV{'HOME'}/.ensembl_init" : undef);
252 $self->{'old_dvd_format'} = $old_dvd_format || 0;
253 ## $self->{'ucsc_coords'} = $ucsc_coords || 0;
254 ## $self->{'verbose'} = $verbose || 0;
255 ## $self->{'release'} = $release;
256
257
258
259 ## if($reg_host && $self->{'reg_config'}){
260 ## warn "You have specified registry parameters and a config file:\t".$self->{'reg_config'}.
261 ## "\nOver-riding config file with specified paramters:\t${reg_user}@${reg_host}:$reg_port";
262 ## }
263
264
265 #Will a general norm method be applicable for all imports?
266 #Already casued problems with Bed imports... remove?
267 #Could set NORM_METHOD in Parser!!
268 #warn "Need to fully implement norm_method is validate_mage, remove ENV NORM_METHOD?";
269 $self->{'norm_method'} = $norm_method; # || $ENV{'NORM_METHOD'};
270
271 #if ($self->vendor ne 'NIMBLEGEN'){
272 # $self->{'no_mage'} = 1;
273 # warn "Hardcoding no_mage for non-NIMBLEGEN imports";
274 # }
275
276
277 # if($self->{'no_mage'} && $self->{'write_mage'}){
278 # throw('-no_mage and -write_mage options are mutually exclusive, please select just one');
279 # }
280
281 ## #### Set up DBs and load and reconfig registry
282 ##
283 ## ### Load Registry
284 ## #Can we load the registry using the assembly version, then just redefine the efg DB?
285 ## #We have problems here if we try and load on a dev version, where no dev DBs are available on ensembldb
286 ## #Get the latest API version for the assembly we want to use
287 ## #Then load the registry from that version
288 ## #Then we can remove some of the dnadb setting code below?
289 ## #This may cause problems with API schema mismatches
290 ## #Can we just test whether the current default dnadb contains the assembly?
291 ## #Problem with this is that this will not have any other data e.g. genes etc
292 ## #which may be required for some parsers
293 ##
294 ## #How does the registry pick up the schema version??
295 ##
296 ## #We should really load the registry first given the dnadb assembly version
297 ## #Then reset the eFG DB as appropriate
298 ##
299 ##
300 ## if ($reg_host || ! defined $self->{'_reg_config'}) {
301 ## #defaults to current ensembl DBs
302 ## $reg_host ||= 'ensembldb.ensembl.org';
303 ## $reg_user ||= 'anonymous';
304 ##
305 ## #Default to the most recent port for ensdb
306 ## if( (! $reg_port) &&
307 ## ($reg_host eq 'ensdb-archive') ){
308 ## $reg_port = 5304;
309 ## }
310 ##
311 ## #This will try and load the dev DBs if we are using v49 schema or API?
312 ## #Need to be mindful about this when developing
313 ## #we need to tip all this on it's head and load the reg from the dnadb version!!!!!!!
314 ##
315 ## my $version_text= ($self->{'release'}) ? 'version '.$self->{'release'} : 'current version';
316 ## $self->log("Loading $version_text registry from $reg_user".'@'.$reg_host);
317 ##
318 ## #Note this defaults API version, hence running with head code
319 ## #And not specifying a release version will find not head version
320 ## #DBs on ensembldb, resulting in an exception from reset_DBAdaptor below
321 ## $reg->load_registry_from_db(
322 ## -host => $reg_host,
323 ## -user => $reg_user,
324 ## -port => $reg_port,
325 ## -pass => $reg_pass,
326 ## #-host => "ens-staging",
327 ## #-user => 'ensro',
328 ## -db_version => $self->{'release'},#51
329 ## -verbose => $self->verbose,
330 ## );
331 ##
332 ## throw('Not sensible to set the import DB as the default eFG DB from '.$reg_host.', please define db params') if ((! $dbname) && (! $db));
333 ## }
334 ## else{
335 ## $self->log("Loading registry from:\t".$self->{'_reg_config'});
336 ## $reg->load_all($self->{'_reg_config'}, 1);
337 ## }
338 ##
339 ##
340 ## #Need to test the DBs here, as we may not have loaded any!
341 ## #get_alias wil fail otherwise?
342 ## #This is a cyclical dependancy as we need alias to get species which we use to grab the DB
343 ## #alias is dependant on core DB being loaded with relevant meta entries.
344 ## #revise this when we split the Importer
345 ##
346 ##
347 ## #Validate species
348 ## my $alias = $reg->get_alias($species) || throw("Could not find valid species alias for $species\nYou might want to clean up:\t".$self->get_dir('output'));
349 ## $self->species($alias);
350 ## $self->{'param_species'} = $species;#Only used for dir generation
351 ##
352 ##
353 ## #SET UP DBs
354 ## if($db){
355 ## #db will have been defined before reg loaded, so will be present in reg
356 ##
357 ## if(! (ref($db) && $db->isa('Bio::EnsEMBL::Funcgen::DBSQL::DBAdaptor'))){
358 ## $self->throw('-db must be a valid Bio::EnsEMBL::Funcgen::DBSQL::DBAdaptor');
359 ## }
360 ## }
361 ## else{ #define eFG DB from params or registry
362 ##
363 ## if($reg_db){#load eFG DB from reg
364 ##
365 ## if($dbname){
366 ## throw("You cannot specify DB params($dbname) and load from the registry at the same time.");
367 ## }
368 ##
369 ## $self->log('WARNING: Loading eFG DB from Registry');
370 ## $db = $reg->get_DBAdaptor($self->species(), 'funcgen');
371 ## throw("Unable to retrieve ".$self->species." funcgen DB from the registry") if ! $db;
372 ## }
373 ## else{#resets the eFG DB in the custom or generic registry
374 ##
375 ## $dbname || throw('Must provide a -dbname if not using default custom registry config');
376 ## #$user || throw('Must provide a -user parameter');#make this default to EFG_WRITE_USER?
377 ## $pass || throw('Must provide a -pass parameter');
378 ##
379 ## #remove this and throw?
380 ## if(! defined $host){
381 ## $self->log('WARNING: Defaulting to localhost');
382 ## $host = 'localhost';
383 ## }
384 ##
385 ## $port ||= 3306;
386 ## my $host_ip = '127.0.0.1';#is this valid for all localhosts?
387 ##
388 ## if ($self->{'ssh'}) {
389 ## $host = `host localhost`; #mac specific? nslookup localhost wont work on server/non-PC
390 ## #will this always be the same?
391 ##
392 ## if (! (exists $ENV{'EFG_HOST_IP'})) {
393 ## warn "Environment variable EFG_HOST_IP not set for ssh mode, defaulting to $host_ip for $host";
394 ## } else {
395 ## $host_ip = $ENV{'EFG_HOST_IP'};
396 ## }
397 ##
398 ## if ($self->host() ne 'localhost') {
399 ## warn "Overriding host ".$self->host()." for ssh connection via localhost($host_ip)";
400 ## }
401 ## }
402 ##
403 ##
404 ## #data version is only used if we don't want to define the dbname
405 ## #This should never be guessed so don't need data_version here
406 ## #$dbname ||= $self->species()."_funcgen_".$self->data_version();
407 ##
408 ##
409 ## #Remove block below when we can
410 ## my $dbhost = ($self->{'ssh'}) ? $host_ip : $host;
411 ##
412 ## #This isn't set yet!?
413 ## #When we try to load, say 49, when we only have 48 on ensembldb
414 ## #This fails because there is not DB set for v49, as it is not on ensembl DB
415 ## #In this case we need to load from the previous version
416 ## #Trap this and suggest using the -schema_version/release option
417 ## #Can we autodetect this and reload the registry?
418 ## #We want to reload the registry anyway with the right version corresponding to the dnadb
419 ## #We could either test for the db in the registry or just pass the class.
420 ##
421 ## $db = $reg->reset_DBAdaptor($self->species(), 'funcgen', $dbname, $dbhost, $port, $self->user, $pass,
422 ## {
423 ## -dnadb_host => $reg_host,
424 ## -dnadb_port => $reg_port,
425 ## -dnadb_assembly => $assm_version,
426 ## -dnadb_user => $reg_user,
427 ## -dnadb_pass => $reg_pass,
428 ## });
429 ##
430 ##
431 ## #ConfigRegistry will try and set this
432 ## #This will fail if there is already one in the registry as it will try
433 ## #and defined a new unique species so as not to overwrite the original
434 ## #e.g. homo_sapiens1
435 ##
436 ## #This is why it was orignally written backwards as we can't easily dynamically redefine
437 ## #an adaptor in the registry without ConfigRegistry trying to change the name
438 ## #the very act of creating a new db to redefine the registry with causes ConfigRegistry
439 ## #to try and register it with a unique species name
440 ##
441 ## #Delete the old funcgen DB from the registry first
442 ## #$reg->remove_DBAdaptor($self->species, 'funcgen');
443 ##
444 ## #ConfigRegistry will automatically configure this new db
445 ##
446 ## #$db = Bio::EnsEMBL::Funcgen::DBSQL::DBAdaptor->new(
447 ## # -user => $user,
448 ## # -host => ($self->{'ssh'}) ? $host_ip : $host,
449 ## # -port => $port,
450 ## # -pass => $pass,
451 ## # #we need to pass dbname else we can use non-standard dbs
452 ## # -dbname => $dbname,
453 ## # -species => $self->species(),
454 ## # -group => 'funcgen',
455 ## # );
456 ##
457 ##
458 ## #if we get a species like homo_sapiens1 here
459 ## #This is because ConfigRegistry is try to make the dbname different between the
460 ## #one already present and the one you're trying to add
461 ## }
462 ## }
463 ##
464 ##
465 ## ### VALIDATE DNADB
466 ## #This is now done in DBAdaptor
467 ##
468 ## #We can change this to just use the assembly version
469 ## #we could even have the wordy assmelby version from the meta table
470 ## #do the standard ensembl subs
471 ## #s/[A-Za-z]//g;
472 ## #s/\.//g;
473 ## #And then validate?
474 ## #Just stick to number version for now.
475 ##
476 ##
477 ## #Now we need to set the dnadb_host params to avoid ensembldb defaults
478 ## #This should check the registry first
479 ## #Then load from the registry db?
480 ## #If we have a multi host registry config file this isn't going to work!
481 ##
482 ## #Is this required anymore as the DBAdaptor handles this?
483 ## #Not if we pass a db with an incorrect dnadb attached.
484 ##
485 ## #if($db->_get_schema_build($db->dnadb()) !~ /_[0-9]+_${assm_version}[a-z]*$/){
486 ### my $warning = "WARNING: dnadb does not match assembly_version $assm_version. Using ensembldb.enembl.org to define the dnadb";
487 ### $warning .= ' rather than the reg_config' if (defined $self->{'_reg_config'});
488 ##
489 ## #We need to account for reg_config DBs which may have custom info in
490 ## #So try reg_config host first, then try ensembldb with warning
491 ## #Could have a reg_config only flag for core dbs
492 ## #Need to implement more params in set_dnadb_by_assembly_version
493 ### $self->log($warning);
494 ##
495 ### $db->set_dnadb_by_assembly_version($assm_version);
496 ### }
497 ##
498 ##
499 ##
500 ##
501 ## #Test connections
502 ## $self->db($db);
503 ## $db->dbc->db_handle;
504 ## $db->dnadb->dbc->db_handle;
505 ## #Set re/disconnect options
506 ##
507 ## #These really need setting dependant on the import parser
508 ## $db->dbc->disconnect_when_inactive(1);
509 ## $db->dnadb->dbc->disconnect_when_inactive(1);
510
511
512
513 ## ### Check analyses/feature_type/cell_type
514 ## if($feature_analysis){
515 ## my $fanal = $self->db->get_AnalysisAdaptor->fetch_by_logic_name($feature_analysis);
516 ## throw("The Feature Analysis $feature_analysis does not exist in the database") if(!$fanal);
517 ## $self->feature_analysis($fanal);
518 ##
519 ## #This currently fails before the config gets loaded!
520 ## #Need to load config before this validation!
521 ## }
522 ##
523 ## if($ctype_name){
524 ## my $ctype = $self->db->get_CellTypeAdaptor->fetch_by_name($ctype_name);
525 ## throw("The CellType $ctype_name does not exist in the database") if(!$ctype);
526 ## $self->cell_type($ctype);
527 ## }
528 ##
529 ## if ($ftype_name) {
530 ## my $ftype = $self->db->get_FeatureTypeAdaptor->fetch_by_name($ftype_name);
531 ## throw("The FeatureType $ftype_name does not exist in the database") if(!$ftype);
532 ## $self->feature_type($ftype);
533 ## }
534
535
536 #Set config here instead?
537 #So we can check all mandatory params
538 #Set vendor specific attr dependent vars
539
540 #Generic input dir
541 $self->{'input_dir'} ||= $self->get_dir("data").'/input/'.$self->{'param_species'}.'/'.$self->vendor().'/'.$self->name();
542
543 if (! -d $self->get_dir('input')) {
544
545 if (@{$self->result_files}) {
546 #This is really InputSet specific
547 #Could go in init_experiment_import
548 $self->log("Processing files:\n\t\t".join("\n\t\t",@{$self->result_files}));
549 } else {
550 throw('input_dir is not defined or does not exist ('.$self->get_dir('input').')');
551 }
552 }
553
554 #Parser specific config
555 $self->set_config();
556
557
558 $self->debug(2, "Importer class instance created.");
559 $self->debug_hash(3, \$self);
560
561 return ($self);
562 }
563
564 =head2 registry_host
565
566 Example : my $reg_host = $imp->registry_host;
567 Description: Accessor for registry host attribute
568 Returntype : string e.g. ensembldb.ensembl.org
569 Exceptions : None
570 Caller : general
571 Status : at risk
572
573 =cut
574
575 sub registry_host{
576 return $_[0]->{'reg_host'};
577 }
578
579 =head2 registry_user
580
581 Example : my $reg_user = $imp->registry_user;
582 Description: Accessor for registry user attribute
583 Returntype : string e.g. ensembldb.ensembl.org
584 Exceptions : None
585 Caller : general
586 Status : at risk
587
588 =cut
589
590 sub registry_user{
591 return $_[0]->{'reg_user'};
592 }
593
594 =head2 registry_port
595
596 Example : my $reg_port = $imp->registry_port;
597 Description: Accessor for registry port attribute
598 Returntype : string e.g. ensembldb.ensembl.org
599 Exceptions : None
600 Caller : general
601 Status : at risk
602
603 =cut
604
605 sub registry_port{
606 return $_[0]->{'reg_port'};
607 }
608
609 =head2 registry_pass
610
611 Example : my $reg_pass = $imp->registry_pass;
612 Description: Accessor for registry pass attribute
613 Returntype : string e.g. ensembldb.ensembl.org
614 Exceptions : None
615 Caller : general
616 Status : at risk
617
618 =cut
619
620 sub registry_pass{
621 return $_[0]->{'reg_pass'};
622 }
623
624
625 #init method kept separate from new due to differing madatory check and set up
626
627 =head2 init_array_import
628
629 Example : $self->init_import();
630 Description: Initialises import by creating working directories
631 and by storing the Experiemnt
632 Returntype : none
633 Exceptions : warns and throws depending on recover and Experiment status
634 Caller : general
635 Status : at risk - merge with register_array_design
636
637 =cut
638
639 sub init_array_import{
640
641 my ($self) = shift;
642
643 # we need to define which paramters we'll be storing
644 #use the logic names of the analyses as the field headers
645
646 #need to test for vendor here
647
648 #Sanger, NIMBLEGEN(no design_id issue, could get from the ndf, but we want it in the DesignNotes.txt)
649 #Then we can change the Array/Chip generation to soley use the DesignNotes.txt rather than SampleKey
650 #which is experiment specific
651 #or eFG format.
652
653 $self->create_output_dirs('caches', 'fastas');
654
655
656 }
657
658
659 =head2 init_experiment_import
660
661 Example : $self->init_import();
662 Description: Initialises import by creating working directories
663 and by storing the Experiemnt
664 Returntype : none
665 Exceptions : warns and throws depending on recover and Experiment status
666 Caller : general
667 Status : at risk - merge with register exeriment
668
669 =cut
670
671 sub init_experiment_import{
672 my ($self) = shift;
673
674 #Change this to take config mandatory params?
675 #No specific to exp import
676 #Name is used in set_config anyway
677 #Currently we only have array and experiment import, both of which should have names
678 #Make mandatory?
679
680 foreach my $tmp ("group", "data_dir") { #name now generically mandatory
681 throw("Mandatory arg $tmp not been defined") if (! defined $self->{$tmp});
682 }
683 #Should we separate path on group here too, so we can have a dev/test group?
684
685 #Create output dirs
686 #This should be moved to the Parser to avoid generating directories which are needed for different imports
687 $self->create_output_dirs('raw', 'norm', 'caches', 'fastas');
688 throw("No result_files defined.") if (! defined $self->result_files());
689
690 #Log input files
691 #if (@{$self->result_files()}) {
692 # $self->log("Found result files arguments:\n\t".join("\n\t", @{$self->result_files()}));
693 #}
694 #This is done in new
695
696 #check for cell||feature and warn if no met file supplied?
697
698 if ($self->norm_method) {
699 my $norm_anal = $self->db->get_AnalysisAdaptor->fetch_by_logic_name($self->norm_method);
700
701 #should we list the valid analyses?
702 throw($self->norm_method.' is not a valid analysis') if ! $norm_anal;
703 $self->norm_analysis($norm_anal);
704 } else {
705 $self->log('WARNING: No normalisation analysis specified');
706 }
707
708 #warn "Need to check env vars here or in Parser or just after set_config?";
709 #Need generic method for checking ENV vars in Helper
710 #check for ENV vars?
711 #R_LIBS
712 #R_PATH if ! farm
713 #R_FARM_PATH
714
715 $self->validate_group(); #import experimental_group
716
717 #Get experiment
718 my $exp_adaptor = $self->db->get_ExperimentAdaptor();
719 my $exp = $exp_adaptor->fetch_by_name($self->name()); #, $self->group());
720 $self->process_experiment_config if $self->can('process_experiment_config'); #Parsers::MAGE::process_experiment_config
721
722 #Moved MAGE support form here to MAGE.pm
723
724 #Recovery now set so deal with experiment
725 if ($self->recovery() && ($exp)) {
726 $self->log("Using previously stored Experiment:\t".$exp->name);
727 } elsif ((! $self->recovery()) && $exp) {
728 throw("Your experiment name is already registered in the database, please choose a different \"name\", this will require renaming you input directory, or specify -recover if you are working with a failed/partial import.");
729 #can we skip this and store, and then check in register experiment if it is already stored then throw if not recovery
730 } else { # (recover && exp) || (recover && ! exp)
731
732
733 $exp = Bio::EnsEMBL::Funcgen::Experiment->new(
734 -EXPERIMENTAL_GROUP => $self->{egroup},
735 -NAME => $self->name(),
736 -DATE => $self->experiment_date(),
737 -PRIMARY_DESIGN_TYPE => $self->design_type(),
738 -DESCRIPTION => $self->description(),
739 -ADAPTOR => $self->db->get_ExperimentAdaptor(),
740 );
741
742 ($exp) = @{$exp_adaptor->store($exp)};
743 }
744
745
746 $self->experiment($exp);
747
748 #remove and add specific report, this is catchig some Root stuff
749 #$self->log("Initiated efg import with following parameters:\n".Data::Dumper::Dumper(\$self));
750
751 return;
752 }
753
754
755 #Move this to new or init_experiment
756
757 =head2 validate_group
758
759 Example : $self->validate_group();
760 Description: Validates groups details
761 Returntype : none
762 Exceptions : throws if insufficient info defined to store new Group and is not already present
763 Caller : general
764 Status : Medium - check location and contact i.e. group name clash?
765
766 =cut
767
768 sub validate_group{
769 my ($self) = shift;
770
771 my $egroup = $self->db->get_ExperimentalGroupAdaptor->fetch_by_name($self->group);
772
773 if (! defined $egroup) {
774
775 throw("validate_group does not yet fully support ExperimentalGroups, please add manually");
776
777 #if ($self->location() && $self->contact()) {
778 # $self->db->import_group($self->group(), $self->location, $self->contact());
779 #} else {
780 # throw("Group ".$self->group()." does not exist, please specify a location and contact to register the group");
781 #}
782 }
783
784 $self->{egroup} = $egroup;
785
786 return;
787 }
788
789 =head2 create_output_dirs
790
791 Example : $self->create_output_dirs();
792 Description: Does what it says on the tin, creates dirs in
793 the root output dir foreach @dirnames, also set paths in self
794 Arg [1] : mandatory - list of dir names
795 Returntype : none
796 Exceptions : none
797 Caller : general
798 Status : Medium - add throw?
799
800 =cut
801
802 sub create_output_dirs{
803 my ($self, @dirnames) = @_;
804
805 #output dir created in control script
806 #avoids errors when logs generated first
807
808
809 foreach my $name (@dirnames) {
810
811 if ($name eq 'caches') {
812 $self->{"${name}_dir"} = $ENV{'EFG_DATA'}."/${name}/".$self->db->dbc->dbname() if(! defined $self->{"${name}_dir"});
813 } elsif ($name eq 'fastas') {
814 $self->{"${name}_dir"} = $ENV{'EFG_DATA'}."/${name}/" if(! defined $self->{"${name}_dir"});
815 } else {
816 $self->{"${name}_dir"} = $self->get_dir('output')."/${name}" if(! defined $self->{"${name}_dir"});
817 }
818
819 if (! (-d $self->get_dir($name) || (-l $self->get_dir($name)))) {
820 $self->log("Creating directory:\t".$self->get_dir($name));
821 #This did not throw with mkdir!!
822 mkpath $self->get_dir($name) || throw('Failed to create directory: '. $self->get_dir($name));
823 chmod 0744, $self->get_dir($name);
824 }
825 }
826
827 return;
828 }
829
830 =head2 vendor
831
832 Example : $imp->vendor("NimbleGen");
833 Description: Getter/Setter for array vendor
834 Arg [1] : optional - vendor name
835 Returntype : string
836 Exceptions : none
837 Caller : general
838 Status : Stable
839
840 =cut
841
842 sub vendor{
843 my ($self) = shift;
844
845 if (@_) {
846 $self->{'vendor'} = shift;
847 $self->{'vendor'} = uc($self->{'vendor'});
848 }
849
850 return $self->{'vendor'};
851 }
852
853
854 =head2 feature_type
855
856 Example : $imp->feature_type($ftype);
857 Description: Getter/Setter for Experiment FeatureType
858 Arg [1] : optional - Bio::EnsEMBL::Funcgen::FeatureType
859 Returntype : Bio::EnsEMBL::Funcgen::FeatureType
860 Exceptions : Throws if arg is not valid or stored
861 Caller : general
862 Status : at risk
863
864 =cut
865
866 sub feature_type{
867 my ($self) = shift;
868
869 if (@_) {
870 my $ftype = shift;
871
872 #do we need this as we're checking in new?
873 if (! ($ftype->isa('Bio::EnsEMBL::Funcgen::FeatureType') && $ftype->dbID())) {
874 throw("Must pass a valid stored Bio::EnsEMBL::Funcgen::FeatureType");
875 }
876
877 $self->{'feature_type'} = $ftype;
878 }
879
880 return $self->{'feature_type'};
881 }
882
883 =head2 feature_analysis
884
885 Example : $imp->feature_analysis($fanal);
886 Description: Getter/Setter for Analysis used for creating the imported Features
887 Arg [1] : optional - Bio::EnsEMBL::Analysis
888 Returntype : Bio::EnsEMBL::Analysis
889 Exceptions : Throws if arg is not valid or stored
890 Caller : general
891 Status : at risk
892
893 =cut
894
895 sub feature_analysis{
896 my ($self) = shift;
897
898 if (@_) {
899 my $fanal = shift;
900
901 #do we need this as we're checking in new?
902 if (! (ref ($fanal) && $fanal->isa('Bio::EnsEMBL::Analysis') && $fanal->dbID())) {
903 throw("Must pass a valid stored Bio::EnsEMBL::Analysis");
904 }
905
906 $self->{'feature_analysis'} = $fanal;
907 }
908
909 return $self->{'feature_analysis'};
910 }
911
912 =head2 norm_analysis
913
914 Example : $imp->norm_analysis($anal);
915 Description: Getter/Setter for the normalisation analysis
916 Arg [1] : optional - Bio::EnsEMBL::Analysis
917 Returntype : Bio::EnsEMBL::Analysis
918 Exceptions : Throws if arg is not valid or stored
919 Caller : general
920 Status : at risk
921
922 =cut
923
924 sub norm_analysis{
925 my ($self) = shift;
926
927 if (@_) {
928 my $anal = shift;
929
930 #do we need this as we're checking in new?
931 if (! (ref($anal) && $anal->isa('Bio::EnsEMBL::Analysis') && $anal->dbID())) {
932 throw("Must pass a valid stored Bio::EnsEMBL::Analysis");
933 }
934
935 $self->{'norm_analysis'} = $anal;
936 }
937
938 return $self->{'norm_analysis'};
939 }
940
941
942
943 =head2 cell_type
944
945 Example : $imp->cell_type($ctype);
946 Description: Getter/Setter for Experiment CellType
947 Arg [1] : optional - Bio::EnsEMBL::Funcgen::CellType
948 Returntype : Bio::EnsEMBL::Funcgen::CellType
949 Exceptions : Throws if arg is not valid or stored
950 Caller : general
951 Status : at risk
952
953 =cut
954
955 sub cell_type{
956 my ($self) = shift;
957
958 if (@_) {
959 my $ctype = shift;
960
961 #do we need this as we're checking in new?
962 if (! ($ctype->isa('Bio::EnsEMBL::Funcgen::CellType') && $ctype->dbID())) {
963 throw("Must pass a valid stored Bio::EnsEMBL::Funcgen::CellType");
964 }
965
966 $self->{'cell_type'} = $ctype;
967 }
968
969 return $self->{'cell_type'};
970 }
971
972
973 ##=head2 ucsc_coords
974 ##
975 ## Example : $start += 1 if $self->ucsc_coords;
976 ## Description: Getter for UCSC coordinate usage flag
977 ## Returntype : boolean
978 ## Exceptions : none
979 ## Caller : general
980 ## Status : at risk
981 ##
982 ##=cut
983 ##
984 ##sub ucsc_coords{
985 ## my $self = shift;
986 ## return $self->{'ucsc_coords'};
987 ##}
988
989
990
991 =head2 array_file
992
993 Example : my $array_file = $imp->array_file();
994 Description: Getter/Setter for sanger/design array file
995 Arg [1] : optional - path to adf or gff array definition/mapping file
996 Returntype : string
997 Exceptions : none
998 Caller : general
999 Status : at risk
1000
1001 =cut
1002
1003 sub array_file{
1004 my ($self) = shift;
1005 $self->{'array_file'} = shift if(@_);
1006 return $self->{'array_file'};
1007 }
1008
1009 =head2 array_name
1010
1011 Example : my $array_name = $imp->array_name();
1012 Description: Getter/Setter for array name
1013 Arg [1] : optional string - name of array
1014 Returntype : string
1015 Exceptions : none
1016 Caller : general
1017 Status : at risk
1018
1019 =cut
1020
1021 sub array_name{
1022 my ($self) = shift;
1023 $self->{'array_name'} = shift if(@_);
1024 return $self->{'array_name'};
1025 }
1026
1027
1028 =head2 array_set
1029
1030 Example : $imp->array_set(1);
1031 Description: Getter/Setter for array set flag
1032 Arg [1] : optional boolean - treat all array chips as the same array
1033 Returntype : boolean
1034 Exceptions : none
1035 Caller : general
1036 Status : at risk
1037
1038 =cut
1039
1040 sub array_set{
1041 my ($self) = shift;
1042 $self->{'array_set'} = shift if(@_);
1043 return $self->{'array_set'};
1044 }
1045
1046
1047 =head2 add_Array
1048
1049 Arg [1] : Bio::EnsEMBL::Funcgen::Array
1050 Example : $self->add_Array($array);
1051 Description: Setter for array elements
1052 Returntype : none
1053 Exceptions : throws if passed non Array or if more than one Array set
1054 Caller : Importer
1055 Status : At risk - Implement multiple arrays? Move to Experiment?
1056
1057 =cut
1058
1059 sub add_Array{
1060 my $self = shift;
1061
1062 #do we need to check if stored?
1063 if (! $_[0]->isa('Bio::EnsEMBL::Funcgen::Array')) {
1064 throw("Must supply a Bio::EnsEMBL::Funcgen::Array");
1065 } elsif (@_) {
1066 push @{$self->{'arrays'}}, @_;
1067 }
1068
1069 throw("Does not yet support multiple array imports") if(scalar (@{$self->{'arrays'}}) > 1);
1070 #need to alter read_probe data at the very least
1071
1072 return;
1073 }
1074
1075
1076
1077 =head2 arrays
1078
1079 Example : foreach my $array(@{$imp->arrays}){ ...do an array of things ...};
1080 Description: Getter for the arrays attribute
1081 Returntype : ARRAYREF
1082 Exceptions : none
1083 Caller : general
1084 Status : at risk
1085
1086 =cut
1087
1088 sub arrays{
1089 my $self = shift;
1090
1091 if (! defined $self->{'arrays'}) {
1092 $self->{'arrays'} = $self->db->get_ArrayAdaptor->fetch_all_by_Experiment($self->experiment());
1093 }
1094
1095 return $self->{'arrays'};
1096 }
1097
1098
1099
1100 =head2 location
1101
1102 Example : $imp->vendor("Hinxton");
1103 Description: Getter/Setter for group location
1104 Arg [1] : optional - location
1105 Returntype : string
1106 Exceptions : none
1107 Caller : general
1108 Status : Stable
1109
1110 =cut
1111
1112
1113 sub location{
1114 my ($self) = shift;
1115 $self->{'location'} = shift if(@_);
1116 return $self->{'location'};
1117 }
1118
1119
1120 =head2 contact
1121
1122 Example : my $contact = $imp->contact();
1123 Description: Getter/Setter for the group contact
1124 Arg [1] : optional - contact name/email/address
1125 Returntype : string
1126 Exceptions : none
1127 Caller : general
1128 Status : Stable
1129
1130 =cut
1131
1132 sub contact{
1133 my ($self) = shift;
1134 $self->{'contact'} = shift if(@_);
1135 return $self->{'contact'};
1136 }
1137
1138 =head2 name
1139
1140 Example : $imp->name('Experiment1');
1141 Description: Getter/Setter for the experiment name
1142 Arg [1] : optional - experiment name
1143 Returntype : string
1144 Exceptions : none
1145 Caller : general
1146 Status : Stable
1147
1148 =cut
1149
1150
1151 sub name{
1152 my ($self) = shift;
1153 $self->{'name'} = shift if(@_);
1154 return $self->{'name'};
1155 }
1156
1157 =head2 result_files
1158
1159 Example : $imp->result_files(\@files);
1160 Description: Getter/Setter for the result file paths
1161 Arg [1] : Listref of file paths
1162 Returntype : Listref
1163 Exceptions : none
1164 Caller : general
1165 Status : At risk
1166
1167 =cut
1168
1169
1170 sub result_files{
1171 my ($self) = shift;
1172 $self->{'result_files'} = shift if(@_);
1173 return $self->{'result_files'};
1174 }
1175
1176
1177
1178
1179
1180 =head2 experiment_date
1181
1182 Example : $imp->experiment_date('2006-11-02');
1183 Description: Getter/Setter for the experiment date
1184 Arg [1] : optional - date string in yyyy-mm-dd
1185 Returntype : string
1186 Exceptions : none
1187 Caller : general
1188 Status : At risk
1189
1190 =cut
1191
1192
1193
1194 sub experiment_date{
1195 my ($self) = shift;
1196
1197 if (@_) {
1198 my $date = shift;
1199
1200 if ($date !~ /[0-9]{4}-[0-9]{2}[0-9]{2}/o) {
1201 throw('Parameter -experiment_date needs to fe in the format: YYYY-MM-DD');
1202 }
1203
1204 $self->{'experiment_date'} = $date;
1205 } elsif ($self->vendor() eq "nimblegen" && ! defined $self->{'experiment_date'}) {
1206 $self->{'experiment_date'} = &get_date("date", $self->get_config("chip_file")),
1207 }
1208
1209 return $self->{'experiment_date'};
1210 }
1211
1212
1213
1214 =head2 group
1215
1216 Example : my $exp_group = $imp->group();
1217 Description: Getter/Setter for the group name
1218 Arg [1] : optional - group name
1219 Returntype : string
1220 Exceptions : none
1221 Caller : general
1222 Status : Stable
1223
1224 =cut
1225
1226 sub group{
1227 my ($self) = shift;
1228 $self->{'group'} = shift if(@_);
1229 return $self->{'group'};
1230 }
1231
1232
1233 =head2 description
1234
1235 Example : $imp->description("Human chrX H3 Lys 9 methlyation");
1236 Description: Getter/Setter for the experiment element
1237 Arg [1] : optional - experiment description
1238 Returntype : string
1239 Exceptions : none
1240 Caller : general
1241 Status : Stable
1242
1243 =cut
1244
1245 sub description{
1246 my $self = shift;
1247
1248 if (@_) {
1249 $self->{'description'} = shift;
1250 }
1251
1252 return $self->{'description'};
1253 }
1254
1255 ##=head2 feature_set_description
1256 ##
1257 ## Example : $imp->description("ExperimentalSet description");
1258 ## Description: Getter/Setter for the FeatureSet description for an
1259 ## InputSet import e.g. preprocessed GFF/Bed data
1260 ## Arg [1] : optional - string feature set description
1261 ## Returntype : string
1262 ## Exceptions : none
1263 ## Caller : general
1264 ## Status : At risk
1265 ##
1266 ##=cut
1267 ##
1268 ##sub feature_set_description{
1269 ## my $self = shift;
1270 ##
1271 ## $self->{'feature_set_description'} = shift if @_;
1272 ##
1273 ## return $self->{'feature_set_description'};
1274 ##}
1275
1276 =head2 format
1277
1278 Example : $imp->format("Tiled");
1279 Description: Getter/Setter for the array format
1280 Arg [1] : optional - array format
1281 Returntype : string
1282 Exceptions : none
1283 Caller : general
1284 Status : Stable
1285
1286 =cut
1287
1288 sub format{
1289 my ($self) = shift;
1290 $self->{'format'} = shift if(@_);
1291 return $self->{'format'};
1292 }
1293
1294 =head2 experiment
1295
1296 Example : my $exp = $imp->experiment();
1297 Description: Getter/Setter for the Experiment element
1298 Arg [1] : optional - Bio::EnsEMBL::Funcgen::Experiment
1299 Returntype : Bio::EnsEMBL::Funcgen::Experiment
1300 Exceptions : throws if arg is not an Experiment
1301 Caller : general
1302 Status : Stable
1303
1304 =cut
1305
1306 sub experiment{
1307 my ($self) = shift;
1308
1309 if (@_) {
1310
1311 if (! $_[0]->isa('Bio::EnsEMBL::Funcgen::Experiment')) {
1312 throw("Must pass a Bio::EnsEMBL::Funcgen::Experiment object");
1313 }
1314
1315 $self->{'experiment'} = shift;
1316 }
1317
1318 return $self->{'experiment'};
1319 }
1320
1321 ##=head2 db
1322 ##
1323 ## Example : $imp->db($funcgen_db);
1324 ## Description: Getter/Setter for the db element
1325 ## Arg [1] : optional - Bio::EnsEMBL::Funcgen::DBSQL::DBAdaptor
1326 ## Returntype : Bio::EnsEMBL::Funcgen::DBSQL::DBAdaptor
1327 ## Exceptions : throws if arg is not an DBAdaptor
1328 ## Caller : general
1329 ## Status : Stable
1330 ##
1331 ##=cut
1332 ##
1333 ##sub db{
1334 ## my $self = shift;
1335 ##
1336 ## if (defined $_[0] && $_[0]->isa("Bio::EnsEMBL::Funcgen::DBSQL::DBAdaptor")) {
1337 ## $self->{'db'} = shift;
1338 ## } elsif (defined $_[0]) {
1339 ## throw("Need to pass a valid Bio::EnsEMBL::Funcgen::DBSQL::DBAdaptor");
1340 ## }
1341 ##
1342 ## return $self->{'db'};
1343 ##}
1344 ##
1345 ##=head2 pass
1346 ##
1347 ## Example : $imp->pass("password");
1348 ## Description: Getter/Setter for the db password
1349 ## Arg [1] : optional - db password
1350 ## Returntype : string
1351 ## Exceptions : none
1352 ## Caller : general
1353 ## Status : Stable
1354 ##
1355 ##=cut
1356 ##
1357 ##
1358 ##sub pass{
1359 ## my $self = shift;
1360 ## $self->{'pass'} = shift if(@_);
1361 ## return $self->{'pass'};
1362 ##}
1363 ##
1364 ##=head2 pass
1365 ##
1366 ## Example : $imp->host("hoastname");
1367 ## Description: Getter/Setter for the db hostname
1368 ## Arg [1] : optional - db hostname
1369 ## Returntype : string
1370 ## Exceptions : none
1371 ## Caller : general
1372 ## Status : Stable
1373 ##
1374 ##=cut
1375 ##
1376 ##sub host{
1377 ## my $self = shift;
1378 ## $self->{'host'} = shift if(@_);
1379 ## return $self->{'host'};
1380 ##}
1381 ##
1382 ##=head2 port
1383 ##
1384 ## Example : $imp->port(3306);
1385 ## Description: Getter/Setter for the db port number
1386 ## Arg [1] : optional - db port number
1387 ## Returntype : int
1388 ## Exceptions : none
1389 ## Caller : general
1390 ## Status : Stable
1391 ##
1392 ##=cut
1393 ##
1394 ##sub port{
1395 ## my $self = shift;
1396 ## $self->{'port'} = shift if(@_);
1397 ## return $self->{'port'};
1398 ##}
1399 ##
1400 ##=head2 user
1401 ##
1402 ## Example : $imp->user("user_name");
1403 ## Description: Getter/Setter for the db user name
1404 ## Arg [1] : optional - db user name
1405 ## Returntype : string
1406 ## Exceptions : none
1407 ## Caller : general
1408 ## Status : Stable
1409 ##
1410 ##=cut
1411 ##
1412 ##sub user{
1413 ## my $self = shift;
1414 ## $self->{'user'} = shift if(@_);
1415 ## return $self->{'user'};
1416 ##}
1417
1418 ##=head2 dump_fasta
1419 ##
1420 ## Example : if($self->dump_fasta()){...do fasta dump...}
1421 ## Description: Getter/Setter for the dump_fasta flag
1422 ## Arg [1] : optional - 0 or 1
1423 ## Returntype : boolean
1424 ## Exceptions : none
1425 ## Caller : self
1426 ## Status : Stable
1427 ##
1428 ##=cut
1429 ##
1430 ##
1431 ##sub dump_fasta{
1432 ## my $self = shift;
1433 ## $self->{'_dump_fasta'} = shift if @_;
1434 ## return $self->{'_dump_fasta'};
1435 ##}
1436 ##
1437
1438
1439 ##=head2 species
1440 ##
1441 ## Example : $imp->species("homo_sapiens");
1442 ## Description: Getter/Setter for species
1443 ## Arg [1] : optional - species name(alias?)
1444 ## Returntype : string
1445 ## Exceptions : none ? throw if no alias found?
1446 ## Caller : general
1447 ## Status : Medium - may move reg alias look up to this method
1448 ##
1449 ##=cut
1450 ##
1451 ##sub species{
1452 ## my $self = shift;
1453 ##
1454 ## #should we do reg alias look up here?
1455 ## #Will we ever want to redefine species?
1456 ## #Change to just getter?
1457 ##
1458 ## $self->{'species'} = shift if(@_);
1459 ##
1460 ## return $self->{'species'};
1461 ##}
1462
1463 =head2 get_dir
1464
1465 Example : $imp->get_dir("import");
1466 Description: Retrieves full path for given directory
1467 Arg [1] : mandatory - dir name
1468 Returntype : string
1469 Exceptions : none
1470 Caller : general
1471 Status : at risk - move to Helper?
1472
1473 =cut
1474
1475 sub get_dir{
1476 my ($self, $dirname) = @_;
1477 return $self->get_data("${dirname}_dir");
1478 }
1479
1480 =head2 norm_method
1481
1482 Example : my $norm_method = $imp->norm_method()
1483 Description: Getter/Setter for normalisation method
1484 Arg [1] : mandatory - method name
1485 Returntype : string
1486 Exceptions : none ? throw if no analysis with logic name
1487 Caller : general
1488 Status : At risk - restrict to logic_name and validate against DB, allow multiple
1489
1490 =cut
1491
1492 #Move to Nimblegen?
1493 #Do we ever want to normalise other data?
1494
1495 sub norm_method{
1496 my $self = shift;
1497
1498 if (@_) {
1499 $self->{'norm_method'} = shift;
1500 } elsif (! defined $self->{'norm_method'}) {
1501 $self->{'norm_method'}= $self->get_config('norm_method');
1502 }
1503
1504 return $self->{'norm_method'};
1505 }
1506
1507
1508 =head2 get_config
1509
1510 Arg [1] : mandatory - name of the data element to retrieve from the config hash
1511 Example : %dye_freqs = %{$imp->get_config('dye_freqs')};
1512 Description: returns data from the definitions hash
1513 Returntype : various
1514 Exceptions : none
1515 Caller : Importer
1516 Status : at risk - replace with direct calls in the inherited Defs class?
1517
1518 =cut
1519
1520
1521 sub get_config{
1522 my ($self, $data_name) = @_;
1523 return $self->get_data('config', $data_name); #will this cause undefs?
1524 }
1525
1526
1527
1528
1529
1530 =head2 register_experiment
1531
1532 Example : $imp->register_experiment()
1533 Description: General control method, performs all data import and normalisations
1534 Arg [1] : optional - dnadb DBAdaptor
1535 Returntype : none
1536 Exceptions : throws if arg is not Bio::EnsEMBL::DBSQL::DBAdaptor
1537 Caller : general
1538 Status : Medium
1539
1540 =cut
1541
1542
1543 #Need to split this method?
1544 #Pre farm process
1545 #Define and store all sets,
1546 #pre process input file once if required
1547 #How are we going to be able to tell wether this has been done successfully?
1548 #runner will catch error, therefore safe to assume that file is complete if farm job
1549 #unless we are not running with farm
1550
1551 #farm specific processes
1552 #actually parse and store data
1553
1554 #Can we separate these for different import types?
1555
1556 sub register_experiment{
1557 my ($self) = shift;
1558
1559 #Need to check for dnadb passed with adaptor to contructor
1560 if (@_) {
1561
1562 if ( ! $_[0]->isa("Bio::EnsEMBL::DBSQL::DBAdaptor")) {
1563 throw("You need to pass a valid dnadb adaptor to register the experiment");
1564 }
1565 $self->db->dnadb($_[0]);
1566 } elsif ( ! $self->db) {
1567 throw("You need to pass/set a DBAdaptor with a DNADB attached of the relevant data version");
1568 }
1569
1570 #This could still be the default core db for the current version
1571 #warn here if not passed DB?
1572 #These should be vendor independent, only the read methods should need specific order?
1573
1574 $self->init_experiment_import();
1575 #can we just have init here instead?
1576
1577
1578 #This could do with a rewrite to move some things to the parsers
1579 #$self::SUPER->register_experiment
1580
1581 $self->write_validate_experiment_config if $self->can('write_validate_experiment_config');
1582
1583 #This is too array specific!
1584 #Can we have an array of import_methods in the config?
1585 #foreach my $method(@{$self->get_config{'import_methods'}}){
1586 #$self->method;
1587 #}
1588 #We're already doing this in read_data for each of these data_types
1589
1590 #Need to be able to run this separately, so we can normalise previously imported sets with different methods
1591 #should be able t do this without raw data files e.g. retrieve info from DB
1592 #Is this implemented?
1593
1594 $self->read_data("probe");
1595 $self->read_data("results");
1596
1597
1598 my $norm_method = $self->norm_method();
1599
1600 if (defined $norm_method) {
1601 warn "norm method is $norm_method";
1602
1603 $self->R_norm($norm_method);
1604 #change this to $self->$norm_method
1605 #so we can have non-R_norm normalisation
1606 }
1607
1608
1609 return;
1610 }
1611
1612 #Move array specific ones to Nimblegen.pm?
1613 #Also used by ArrayDesign and Sanger.pm
1614 #So need to create Array.pm baseclass, which is a Helper.
1615
1616 =head2 store_set_probes_features
1617
1618 Arg [1] : mandatory - array chip id
1619 Arg [2] : optional - Bio::EnsEMBL::Funcgen::ProbeSet
1620 Arg [3] : mandatory - hashref of keys probe id, values are
1621 hash of probe/features with values
1622 Bio::EnsEMBL::Funcgen::Probe/Features for a given
1623 probe set if defined.
1624 Example : $self->store_set_probes_features($ac->dbID(), $ops, \%pfs);
1625 Description: Stores probe set, probes and probe features
1626 Returntype : none
1627 Exceptions : none
1628 Caller : self
1629 Status : Medium
1630
1631 =cut
1632
1633 sub store_set_probes_features{
1634 my ($self, $ac_id, $pf_hash, $ops) = @_;
1635
1636 ### Deal with ProbeSets
1637 if ($ops) {
1638 $ops->size(scalar(keys %$pf_hash));
1639 ($ops) = $self->db->get_ProbeSetAdaptor->store($ops);
1640 }
1641
1642 #If we're going to validate fully, we need to check for probes in this probeset on this array chip
1643 #Update size if we have any new probes
1644 #Overkill? Only do on recover? Do not read if array chip is IMPORTED
1645 #This does not make any attempt to validate probes/set vs previously stored data
1646
1647 for my $probe_id (keys %$pf_hash) {
1648
1649 #set probeset in probe and store
1650 #the process corresponding feature
1651 my $probe = $pf_hash->{$probe_id}->{'probe'};
1652 $probe->probeset($ops) if $ops;
1653 ($probe) = @{$self->db->get_ProbeAdaptor->store($probe)};
1654
1655 #Can't use get_all_Arrays here as we can't guarantee this will only ever be the array we've generated
1656 #Might dynamically load array if non-present
1657 #This is allowing multiple dbIDs per probe??? Is this wrong?
1658 #$self->cache_probe_info($probe->get_probename(), $probe->dbID());###########Remove as we're now importing all then resolving
1659
1660
1661 foreach my $feature (@{$pf_hash->{$probe_id}->{'features'}}) {
1662 $feature->probe($probe);
1663 ($feature) = @{$self->db->get_ProbeFeatureAdaptor->store($feature)};
1664 }
1665 }
1666
1667 undef $ops; #Will this persist in the caller?
1668 undef %{$pf_hash};
1669
1670 return;
1671 }
1672
1673
1674 =head2 cache_slice
1675
1676 Arg [0] : string - region_name e.g. X
1677 Arg [1] : optional - coordinate system name e.g. supercontig, defaults to chromosome
1678 Example : my $slice = $self->cache_slice(12);
1679 Description: Caches or retrieves from cache a given slice
1680 Returntype : Bio::EnsEMBL::Slice
1681 Exceptions : throws f no region name specified
1682 Caller : self
1683 Status : At risk
1684
1685 =cut
1686
1687 sub cache_slice{
1688 my ($self, $region_name, $cs_name, $total_count) = @_;
1689
1690 throw("Need to define a region_name to cache a slice from") if ! $region_name;
1691 $self->{'slice_cache'} ||= {};
1692 $region_name =~ s/chr//;
1693 $region_name = "MT" if $region_name eq "M";
1694
1695 if (! exists ${$self->{'seen_slice_cache'}}{$region_name}) {
1696 my $slice = $self->slice_adaptor->fetch_by_region($cs_name, $region_name);
1697
1698 #Set seen cache so we don't try this again
1699 $self->{seen_slice_cache}{$region_name} = $slice;
1700
1701 if (! $slice) {
1702 warn("-- Could not generate a slice for ${cs_name}:$region_name\n");
1703 } else {
1704
1705 my $sr_name = $slice->seq_region_name; #In case we passed a slice name
1706
1707 if (@{$self->{seq_region_names}}) {
1708 return if ! grep(/^${sr_name}$/, @{$self->{seq_region_names}}); #not on required slice
1709 }
1710 }
1711
1712 $self->{'slice_cache'}->{$region_name} = $slice;
1713 }
1714
1715 if ($total_count && exists ${$self->{'seen_slice_cache'}}{$region_name}) {
1716 #This is an InputSet specific method
1717 $self->count('total_features') if $self->can('count');
1718 }
1719
1720 #Only return if exists to avoid creating hash key
1721 return (exists $self->{'slice_cache'}->{$region_name}) ? $self->{'slice_cache'}->{$region_name} : undef;
1722 }
1723
1724 =head2 slice_cache
1725
1726 Example : my @seen_slices = values(%{$self->slice_cache});;
1727 Description: Returns the slice cache i.e. all the Slices seen in the data filtered
1728 by the defined slices. This method can be used to run only the appropriate
1729 slice jobs after a prepare stage.
1730 Returntype : Hashref of seq_region name Bio::EnsEMBL::Slice pairs
1731 Exceptions : None
1732 Caller : self
1733 Status : At risk
1734
1735 =cut
1736
1737
1738 sub slice_cache{
1739 my $self = shift;
1740
1741 return $self->{'slice_cache'};
1742 }
1743
1744
1745
1746
1747 =head2 cache_probe_info
1748
1749 Arg [0] : mandatory - probe name
1750 Arg [1] : mandatory - probe dbID
1751 Arg [2] : optioanl int - x coord of probe on array
1752 Arg [3] : optional int - y coord of probe on array
1753 Example : $self->cache_probe_info("Probe1", $probe->dbID());
1754 Or for result files which do not have X & Y, we need to cache
1755 X & Y from the design files: $self->cache_probe_info('Probe2', $probe->dbID(), $x, $y);
1756 Description: Setter for probe cache values
1757 Returntype : none
1758 Exceptions : throws is cache conflict encountered
1759 Caller : self
1760 Status : At risk - merge with following?
1761
1762 =cut
1763
1764 sub cache_probe_info{
1765 my ($self, $pname, $pid, $x, $y) = @_;
1766
1767 throw('Deprecated, too memory expensive, now resolving DB duplicates and using Tied File cache');
1768 throw("Must provide a probe name and id") if (! defined $pname || ! defined $pid);
1769
1770
1771 #do we need to loop through the file here?
1772 #if we make sure we're testing for a previous dbID before storing probes then we don't need to do this
1773 #we can catch the error when we get the probe id as we can check for >1 id for the same probe name
1774 #if (defined $self->{'_probe_cache'}->{$pname} && ($self->{'_probe_cache'}->{$pname}->[0] != $pid)) {
1775 # throw("Found two differing dbIDs for $pname, need to sort out redundant probe entries");
1776 #}
1777
1778 $self->{'_probe_cache'}->{$pname} = (defined $x && defined $y) ? [$pid, $x, $y] : [$pid];
1779
1780 return;
1781 }
1782
1783
1784 =head2 get_probe_id_by_name_Array
1785
1786 Arg [1] : mandatory - probe name
1787 Example : $pid = $self->get_probe_id_by_name($pname);
1788 Description: Getter for probe cache values
1789 Returntype : int
1790 Exceptions : none
1791 Caller : self
1792 Status : At risk - merge with previous, move to importer?
1793
1794 =cut
1795
1796 sub get_probe_id_by_name_Array{
1797 my ($self, $name, $array) = @_;
1798
1799 #this is only ever called for fully imported ArrayChips, as will be deleted if recovering
1800 $self->resolve_probe_data() if(! exists $self->{'_probe_cache'}{$array->name()});
1801
1802 #we want to cycle through the given cache starting from the last position or 0.
1803 #we don't want to have to test for the size of the cache each time as this will be quite expensive
1804 #so we should store sizes separately along with current position
1805
1806
1807 my ($pid, $line);
1808
1809 #check current line
1810 if ($line = $self->{'_probe_cache'}{$array->name()}{'current_line'}) {
1811 if ($line =~ /^\Q${name}\E\t/) {
1812 $pid = (split/\t/o, $line)[1];
1813 }
1814 }
1815
1816
1817 if (! $pid) {
1818 while ($line = $self->{'_probe_cache'}{$array->name()}{'handle'}->getline()) {
1819
1820 if ($line =~ /^\Q${name}\E\t/) {
1821 $pid = (split/\t/o, $line)[1];
1822 $self->{'_probe_cache'}{$array->name()}{'current_line'} = $line;
1823 last;
1824 }
1825 }
1826 }
1827
1828 #do not remove this
1829 if (! $pid) {
1830 throw("Did not find probe name ($name) in cache, cache may need rebuilding, results may need sorting, or do you have an anomolaous probe?")
1831 } else {
1832 chomp $pid;
1833 }
1834
1835 return $pid;
1836 }
1837
1838 =head2 get_probe_cache_by_Array
1839
1840 Arg[1] : Bio::EnsEMBL::Funcgen::Array
1841 Arg[2] : boolean - from db flag, only to be used by Importer->resolve_probe_data !
1842 Example : $self->get_probe_cache_by_Array();
1843 Description: Gets the probe info cache which is an array tied to a file
1844 Returntype : Boolean - True if cache has been generated and set successfully
1845 Exceptions : none
1846 Caller : general
1847 Status : At risk
1848
1849 =cut
1850
1851 #from db flag should only be used by importer
1852 #this is because there is no guarantee that it will be resolved unless
1853 #called by resolve_probe_data
1854 #which then renames the file and resets the handle
1855 #can we clean this up and protect/hide this functionality?
1856 #can we check the cache file name in the get methods and throw if it contains unresolved?
1857 #make this private?
1858
1859 sub get_probe_cache_by_Array{
1860 my ($self, $array, $from_db) = @_;
1861
1862 my $msg = "Getting probe cache for ".$array->name();
1863 $msg .= " from DB" if $from_db;
1864 $self->log($msg); #, 1);
1865
1866 if (! ($array && $array->isa('Bio::EnsEMBL::Funcgen::Array') && $array->dbID())) {
1867 throw('Must provide a valid stored Bio::EnsEMBL::Funcgen::Array object');
1868 }
1869
1870 my $set = 0;
1871 my $cache_file = $self->get_dir('caches').'/'.$array->name().'.probe_cache';
1872
1873 ### Generate and resolve fresh cache from DB
1874 if ($from_db) {
1875
1876 $cache_file .= '.unresolved'; #This will be renamed by the caller if it is resolved
1877
1878 if (exists $self->{'_probe_cache'}{$array->name()}) {
1879 $self->log('Rebuilding probe_cache from DB for '.$array->name(), 1);
1880
1881
1882 #untie @{$self->{'_probe_cache'}{$array->name()}{'entries'}};
1883 #close($self->{'_probe_cache'}{$array->name()}{'handle'});#do we need to do this?
1884 delete $self->{'_probe_cache'}{$array->name()}; #implicitly closes
1885 $self->log('Deleted old cache', 1);
1886 } else {
1887 $self->log('Building probe_cache from DB for '.$array->name(), 1);
1888 }
1889
1890 #Move this to ProbeAdaptor?
1891 #This is where we'd set the unique key for a vendor and resolves duplicates based on the key
1892 my $cmd = 'SELECT name, probe_id from probe WHERE array_chip_id IN ('.join(',', @{$array->get_array_chip_ids()}).') ORDER by name, probe_id';
1893 $cmd = 'mysql '.$self->db->connect_string()." -e \"$cmd\" >".$cache_file;
1894 run_system_cmd($cmd);
1895
1896 }
1897
1898 ### Set cache
1899 if (-f $cache_file) {
1900 $self->log('MD5 check here?',1);
1901 $self->{'_probe_cache'}{$array->name()}{'current_line'} = undef;
1902 $self->{'_probe_cache'}{$array->name()}{'handle'} = open_file($cache_file);
1903
1904 #can we do a select count instead? and do this instead of the MD5?
1905 #$cmd = "wc -l $cache_file";
1906 #my $size = `$cmd`;
1907
1908 $set = 1;
1909 } else {
1910 warn 'Failed to get probe cache for array:'.$array->name();
1911 }
1912
1913 return $set;
1914 }
1915
1916
1917 #should reorganise these emthods to split reading the array data, and the actual data
1918 #currently:
1919 #meta reads array and chip data
1920 #probe reads probe_set, probes, which should definitely be in array, probe_feature? and results
1921 #native data format may not map to these methods directly, so may need to call previous method if required data not defined
1922
1923
1924 =head2 read_data
1925
1926 Example : $self->read_data("probe")
1927 Description: Calls each method in data_type array from config hash
1928 Arg [1] : mandatory - data type
1929 Returntype : none
1930 Exceptions : none
1931 Caller : self
1932 Status : At risk
1933
1934 =cut
1935
1936 sub read_data{
1937 my($self, $data_type) = @_;
1938
1939 map {my $method = "read_${_}_data"; $self->$method()} @{$self->get_config("${data_type}_data")};
1940 return;
1941 }
1942
1943
1944 =head2 design_type
1945
1946 Example : $self->design_type("binding_site_identification")
1947 Description: Getter/Setter for experimental design type
1948 Arg [1] : optional - design type
1949 Returntype : string
1950 Exceptions : none
1951 Caller : general
1952 Status : At risk
1953
1954 =cut
1955
1956 sub design_type{
1957 my $self = shift;
1958 return $self->{'design_type'};
1959 }
1960
1961
1962 =head2 get_chr_seq_region_id
1963
1964 Example : $seq_region_id = $self->get_seq_region_id('X');
1965 Description: Calls each method in data_type array from config hash
1966 Arg [1] : mandatory - chromosome name
1967 Arg [2] : optional - start value
1968 Arg [3] : optional - end value
1969 Returntype : int
1970 Exceptions : none
1971 Caller : self
1972 Status : At risk
1973
1974 =cut
1975
1976 #convinience wrapper method
1977 #could we use the seq region cache instead?
1978 #this seems like a lot of overhead for getting the id
1979 sub get_chr_seq_region_id{
1980 my ($self, $chr, $start, $end) = @_;
1981 #what about strand info?
1982
1983 #do we need the start and stop?
1984
1985 #use start and stop to prevent problems with scaffodl assemblies, i.e. >1 seq_region_id
1986 #my $slice = $self->slice_adaptor->fetch_by_region("chromosome", $chr, $start, $end);
1987 #we could pass the slice back to the slice adaptor for this, to avoid dbid problems betwen DBs
1988
1989 ###would need to implement other cs's here
1990
1991 return $self->slice_adaptor->fetch_by_region("chromosome", $chr, $start, $end)->get_seq_region_id();
1992 }
1993
1994
1995
1996 =head2 vsn_norm
1997
1998 Example : $self->vsn_norm();
1999 Description: Convinience/Wrapper method for vsn R normalisation
2000 Returntype : none
2001 Exceptions : none
2002 Caller : general
2003 Status : At risk
2004
2005 =cut
2006
2007 #Have Norm class or contain methods in importer?
2008 #Need to have analysis set up script for all standard analyses.
2009
2010 sub vsn_norm{
2011 my $self = shift;
2012 return $self->R_norm("VSN_GLOG");
2013 }
2014
2015
2016 =head2 farm
2017
2018 Arg [1] : Boolean
2019 Example : $importer->farm(1);
2020 Description: Flag to turn farm submission on
2021 Returntype : Boolean
2022 Exceptions : Throws is argument not a boolean
2023 Caller : general
2024 Status : At risk
2025
2026 =cut
2027
2028
2029 sub farm{
2030 my ($self, $farm) = @_;
2031
2032 $self->{'farm'} ||= undef; #define farm
2033
2034 if (defined $farm) {
2035 throw("Argument to farm must be a boolean 1 or 0") if(! ($farm == 1 || $farm == 0));
2036 $self->{'farm'} = $farm;
2037 }
2038
2039 return $self->{'farm'};
2040
2041 }
2042
2043 =head2 batch_job
2044
2045 Arg [1] : Boolean
2046 Example : $importer->batch_job(1);
2047 Description: Flag to turn on batch_job status
2048 Returntype : Boolean
2049 Exceptions : Throws is argument not a boolean
2050 Caller : general
2051 Status : At risk
2052
2053 =cut
2054
2055
2056 sub batch_job{
2057 my ($self, $batch_job) = @_;
2058
2059 #$self->{'batch_job'} ||= undef;
2060
2061 if (defined $batch_job) {
2062 throw("Argument to batch_job must be a boolean 1 or 0") if(! ($batch_job == 1 || $batch_job == 0));
2063 $self->{'batch_job'} = $batch_job;
2064 }
2065
2066 return $self->{'batch_job'};
2067
2068 }
2069
2070 =head2 prepared
2071
2072 Arg [1] : Boolean
2073 Example : $importer->prepared(1);
2074 Description: Flag to turn on prepared file status
2075 This signifies that the files have been previously imported
2076 using prepare mode and may not match the InputSubset names
2077 Returntype : Boolean
2078 Exceptions : None
2079 Caller : general
2080 Status : At risk
2081
2082 =cut
2083
2084
2085 sub prepared{
2086 my ($self, $prepared) = @_;
2087 $self->{'prepared'} = $prepared if (defined $prepared);
2088 return $self->{'prepared'};
2089 }
2090
2091
2092 =head2 R_norm
2093
2094 Example : $self->R_norm(@logic_names);
2095 Description: Performs R normalisations for given logic names
2096 Returntype : none
2097 Exceptions : Throws if R exits with error code or if data not not valid for analysis
2098 Caller : general
2099 Status : At risk
2100
2101 =cut
2102
2103 sub R_norm{
2104 my ($self, @logic_names) = @_;
2105 #This currently normalises a single two colour chip at a time
2106 #rather than normalising across a set of chip
2107 #also does in sets of analyses
2108 #good for keeping data separate, but not efficient in terms of querying
2109 #convert to use one script which only queries for each echip once, then does each anal
2110
2111
2112 my $aa = $self->db->get_AnalysisAdaptor();
2113 my $rset_adaptor = $self->db->get_ResultSetAdaptor();
2114 my $ra_id = $aa->fetch_by_logic_name("RawValue")->dbID();
2115 my %r_config = (
2116 "VSN_GLOG" => {( libs => ['vsn'],
2117 #norm_method => 'vsn',
2118 )},
2119
2120 "T.Biweight" => {(
2121 libs => ['affy'],
2122 #norm_method => 'tukey.biweight',
2123 )},
2124 );
2125
2126
2127 foreach my $logic_name (@logic_names) {
2128 #throw("Not yet implemented TukeyBiweight") if $logic_name eq "Tukey_Biweight";
2129
2130 #this has already been chcecked and set as the norm_analysis
2131 #need to resolve this multi method approach
2132 my $norm_anal = $aa->fetch_by_logic_name($logic_name);
2133
2134
2135
2136 #This only creates an RSet for the IMPORT set
2137 #So if we re-run with a different analysis
2138 #tab2mage will have already been validated
2139 #So RSet generation will be skipped
2140 #We need to recreate the each non-import RSet for this norm analysis
2141
2142 #This also means the RSets are being created before the data has been imported
2143 #This avoids having to parse tab2mage each time but means we have an uncertain status of these Rsets
2144
2145 my $rset = $self->get_import_ResultSet($norm_anal, 'experimental_chip');
2146
2147 my @chips = ();
2148
2149 if (! $rset) {
2150 $self->log("All ExperimentalChips already have status:\t${logic_name}");
2151 } else { #Got data to normalise and import
2152 my @dbids;
2153 my $R_file = $self->get_dir("norm")."/${logic_name}.R";
2154 my $job_name = $self->experiment->name()."_${logic_name}";
2155 my $resultfile = $self->get_dir("norm")."/result.${logic_name}.txt";
2156 my $outfile = $self->get_dir("norm")."/${logic_name}.out";
2157
2158 #How do we get farm job output i.e. run time memusage
2159 #from interactive job?
2160 #This assumes R_PATH
2161 my $errfile = $self->get_dir("norm")."/${logic_name}.err";
2162
2163 #Let's build this better so we capture the farm output aswell as the job output.
2164 my $cmdline = "$ENV{'R_PATH'} --no-save < $R_file"; # >$errfile 2>&1";
2165 #-K option waits for job to complete
2166 my $bsub = "bsub -K -J $job_name ".$ENV{'R_BSUB_OPTIONS'}.
2167 " -e $errfile -o $outfile $ENV{'R_FARM_PATH'} CMD BATCH $R_file";
2168
2169 #Can we separate the out and err for commandline?
2170 my $r_cmd = (! $self->farm()) ? "$cmdline >$outfile 2>&1" : $bsub;
2171
2172 $self->backup_file($resultfile); #Need to do this as we're appending in the loop
2173
2174 #setup qurey
2175 #warn "Need to add host and port here";
2176 #Set up DB, defaults and libs for each logic name
2177 my $query = "options(scipen=20);library(RMySQL);library(Ringo);";
2178 #scipen is to prevent probe_ids being converted to exponents
2179 #Ringo is for default QC
2180
2181 #foreach my $ln(@logic_names){
2182
2183 foreach my $lib (@{$r_config{$logic_name}{'libs'}}) {
2184 $query .= "library($lib);";
2185 }
2186 #}
2187
2188 $query .= "con<-dbConnect(dbDriver(\"MySQL\"), host=\"".$self->db->dbc->host()."\", port=".$self->db->dbc->port().", dbname=\"".$self->db->dbc->dbname()."\", user=\"".$self->db->dbc->username()."\"";
2189
2190 #should use read only pass here as we are printing this to file
2191 $query .= ", pass=\"".$self->db->dbc->password."\")\n";
2192
2193
2194 #Build queries for each chip
2195 foreach my $echip (@{$self->experiment->get_ExperimentalChips()}) {
2196
2197
2198 #should implement logic name here?
2199 #can't as we need seperate ResultSet for each
2200
2201
2202 if ($echip->has_status($logic_name)) {
2203 $self->log("ExperimentalChip ".$echip->unique_id()." already has status:\t$logic_name");
2204 } else {
2205
2206 #warn "Need to roll back here if recovery, as norm import may fail halfway through";
2207
2208 push @chips, $echip;
2209 my $cc_id = $rset->get_chip_channel_id($echip->dbID());
2210
2211 #if ($self->recovery()){
2212 # $self->log('Rolling back results for ExperimentalChip('.$echip->dbID().") $logic_name");
2213 # $self->db->rollback_results($cc_id) if $self->recovery();
2214 # }
2215
2216 $self->log("Building $logic_name R cmd for ".$echip->unique_id());
2217 @dbids = ();
2218
2219 foreach my $chan (@{$echip->get_Channels()}) {
2220
2221 if ($chan->type() eq "EXPERIMENTAL") {
2222 push @dbids, $chan->dbID();
2223 } else {
2224 unshift @dbids, $chan->dbID();
2225 }
2226 }
2227
2228 throw("vsn does not accomodate more than 2 channels") if (scalar(@dbids > 2) && $logic_name eq "VSN_GLOG");
2229
2230 #should do some of this with maps?
2231 #HARDCODED metric ID for raw data as one
2232
2233 #Need to get total and experimental here and set db_id accordingly
2234 #can probably do this directly into one df
2235
2236 $query .= "c1<-dbGetQuery(con, 'select r.probe_id as PROBE_ID, r.score as CONTROL_score, r.X, r.Y from result r, chip_channel c, result_set rs where c.table_name=\"channel\" and c.table_id=${dbids[0]} and c.result_set_id=rs.result_set_id and rs.analysis_id=${ra_id} and c.chip_channel_id=r.chip_channel_id')\n";
2237
2238 $query .= "c2<-dbGetQuery(con, 'select r.probe_id as PROBE_ID, r.score as EXPERIMENTAL_score, r.X, r.Y from result r, chip_channel c, result_set rs where c.table_name=\"channel\" and c.table_id=${dbids[1]} and c.result_set_id=rs.result_set_id and rs.analysis_id=${ra_id} and c.chip_channel_id=r.chip_channel_id')\n";
2239
2240 #Can we define some of the basic structures here and reuse in the QC and each norm method?
2241
2242
2243 #Is this going to eat up memory?
2244 #can we strip out and separate the data from c1 and c2 into RGList and
2245 #individual vector for probe_ids, then rm c1 and c2 to free up memory
2246
2247 #create RGList object
2248 $query .= "R<-as.matrix(c1['CONTROL_score'])\nG<-as.matrix(c2['EXPERIMENTAL_score'])\n";
2249 $query .= "genes<-cbind(c1['PROBE_ID'], c1['X'], c1['Y'])\n";
2250 $query .= "testRG<-new('RGList', list(R=R, G=G, genes=genes))\n";
2251
2252
2253 #QC plots here before doing norm
2254
2255 #open pdf device
2256 $query .= "pdf('".$self->get_dir('norm').'/'.$echip->unique_id."_QC.pdf', paper='a4', height = 15, width = 9)\n";
2257 #set format
2258 $query .= "par(mfrow = c(2,2), font.lab = 2)\n";
2259
2260 #Channel densisties
2261 #These need limma or Ringo
2262 $query .= "plotDensities(testRG)\n";
2263
2264 #MvA Plot
2265
2266 $query .= 'meanLogA<-((log(testRG$R, base=exp(2)) + log(testRG$G, base=exp(2)))/2)'."\n";
2267 $query .= 'logIntRatioM<-(log(testRG$R, base=exp(2)) - log(testRG$G, base=exp(2)))'."\n";
2268 $query .= "yMin<-min(logIntRatioM)\n";
2269 $query .= "yMax<-max(logIntRatioM)\n";
2270
2271
2272
2273 #Need to validate yMax here
2274 #If is is Inf then we need to sort the vector and track back until we find the high real number
2275 #count number of Infs and note on MvA plot
2276 $query .= "infCount<-0\n";
2277 $query .= "if( yMax == Inf){; sortedM<-sort(logIntRatioM); lengthM<-length(logIntRatioM); indexM<-lengthM\n"
2278 ."while (yMax == Inf){; indexM<-(indexM-1); yMax<-sortedM[indexM];}; infCount<-(lengthM-indexM);}\n";
2279
2280 #
2281 $query .= "if(infCount == 0){\n";
2282 $query .= 'plot(meanLogA, logIntRatioM, xlab="A - Average Log Ratio",ylab="M - Log Ratio",pch=".",ylim=c(yMin,yMax), main="'.$echip->unique_id.'")'."\n";
2283 $query .= "} else {\n";
2284 $query .= 'plot(meanLogA, logIntRatioM, xlab="A - Average Log Ratio",ylab="M - Log Ratio",pch=".",ylim=c(yMin,yMax), main="'.$echip->unique_id.'", sub=paste(infCount, " Inf values not plotted"));'."}\n";
2285
2286
2287 #$query .= 'plot(log(testRG$R*testRG$G, base=exp(2))/2, log(testRG$R/testRG$G, base=exp(2)),xlab="A",ylab="M",pch=".",ylim=c(-3,3), main="'.$echip->unique_id."\")\n";
2288
2289 #Plate plots
2290 $query .= 'image(testRG, 1, channel = "green", mycols = c("black", "green4", "springgreen"))'."\n";
2291 $query .= 'image(testRG, 1, channel = "red", mycols = c("black", "green4", "springgreen"))'."\n";
2292
2293 $query .= "dev.off()\n";
2294 #Finished QC pdf printing
2295
2296
2297
2298 #The simple preprocess step of Ringo is actually vsn, so we can nest these in line
2299
2300
2301 ### Build Analyses cmds ###
2302
2303 if ($logic_name eq 'T.Biweight') {
2304
2305 #log2 ratios
2306 $query .= 'lr_df<-cbind((log(c2["EXPERIMENTAL_score"], base=exp(2)) - log(c1["CONTROL_score"], base=exp(2))))'."\n";
2307
2308 #Adjust using tukey.biweight weighted average
2309 #inherits first col name
2310 $query .= 'norm_df<-(lr_df["EXPERIMENTAL_score"]-tukey.biweight(as.matrix(lr_df)))'."\n";
2311 $query .= 'formatted_df<-cbind(rep.int(0, length(c1["PROBE_ID"])), c1["PROBE_ID"], sprintf("%.3f", norm_df[,1]), rep.int('.$cc_id.', length(c1["PROBE_ID"])), c1["X"], c1["Y"])'."\n";
2312
2313 } elsif ($logic_name eq 'VSN_GLOG') {
2314 #could do this directly
2315 $query .= "raw_df<-cbind(c1[\"CONTROL_score\"], c2[\"EXPERIMENTAL_score\"])\n";
2316 #variance stabilise
2317 $query .= "norm_df<-vsn(raw_df)\n";
2318
2319
2320 #do some more calcs here and print report?
2321 #fold change exponentiate? See VSN docs
2322 #should do someplot's of raw and glog and save here?
2323 #set log func and params
2324 #$query .= "par(mfrow = c(1, 2)); log.na = function(x) log(ifelse(x > 0, x, NA));";
2325 #plot
2326 #$query .= "plot(exprs(glog_df), main = \"vsn\", pch = \".\");".
2327 # "plot(log.na(exprs(raw_df)), main = \"raw\", pch = \".\");";
2328 #FAILS ON RAW PLOT!!
2329 #par(mfrow = c(1, 2))
2330 #> meanSdPlot(nkid, ranks = TRUE)
2331 #> meanSdPlot(nkid, ranks = FALSE)
2332
2333
2334 #Now create table structure with glog values(diffs)
2335 #3 sig dec places on scores(doesn't work?!)
2336 $query .= 'formatted_df<-cbind(rep.int(0, length(c1["PROBE_ID"])), c1["PROBE_ID"], sprintf("%.3f", (exprs(norm_df[,2]) - exprs(norm_df[,1]))), rep.int('.$cc_id.', length(c1["PROBE_ID"])), c1["X"], c1["Y"])'."\n";
2337
2338 }
2339 #load back into DB
2340 #c3results<-cbind(rep("", length(c3["probe_id"])), c3["probe_id"], c3["c3_score"], rep(1, length(c3["probe_id"])), rep(1, length(c3["probe_id"])))
2341 #may want to use safe.write here
2342 #dbWriteTable(con, "result", c3results, append=TRUE)
2343 #dbWriteTable returns true but does not load any data into table!!!
2344
2345 $query .= "write.table(formatted_df, file=\"${resultfile}\", sep=\"\\t\", col.names=FALSE, row.names=FALSE, quote=FALSE, append=TRUE)\n";
2346
2347 #tidy up here??
2348 }
2349 }
2350
2351 $query .= "q();";
2352
2353 open(RFILE, ">$R_file") || throw("Cannot open $R_file for writing");
2354 print RFILE $query;
2355 close(RFILE);
2356
2357 my $submit_text = "Submitting $logic_name job";
2358 $submit_text .= ' to farm' if $self->farm;
2359 $self->log("${submit_text}:\t".localtime());
2360 run_system_cmd($r_cmd);
2361 $self->log("Finished $logic_name job:\t".localtime());
2362 $self->log('See '.$self->get_dir('norm').' for ExperimentalChip QC files');
2363
2364 #Now load file and update status
2365 #Import directly here to avoid having to reparse all results if we crash!!!!
2366 $self->log("Importing:\t$resultfile");
2367 $self->db->load_table_data("result", $resultfile);
2368 $self->log("Finishing importing:\t$resultfile");
2369
2370
2371 foreach my $echip (@chips) {
2372 $echip->adaptor->store_status($logic_name, $echip);
2373 }
2374
2375 #Recreate all non-import RSets for analysis if not already present
2376 #
2377
2378 my $rset_a = $self->db->get_ResultSetAdaptor();
2379 my %seen_rsets;
2380
2381 foreach my $anal_rset (@{$rset_a->fetch_all_by_Experiment($self->experiment)}) {
2382 next if($anal_rset->name =~ /_IMPORT$/o);
2383 next if(exists $seen_rsets{$anal_rset->name});
2384 next if $anal_rset->analysis->logic_name eq $norm_anal->logic_name;
2385 $seen_rsets{$rset->name} = 1;
2386 $anal_rset->analysis($norm_anal);
2387 $anal_rset->{'dbID'} = undef;
2388 $anal_rset->{'adaptor'} = undef;
2389
2390 #add the chip_channel_ids from the new anal IMPORT set
2391 foreach my $table_id (@{$anal_rset->table_ids}) {
2392 $anal_rset->{'table_id_hash'}{$table_id} = $rset->get_chip_channel_id($table_id);
2393 }
2394
2395 $self->log('Adding new ResultSet '.$anal_rset->name.' with analysis '.$norm_anal->logic_name);
2396 $rset_a->store($anal_rset);
2397 }
2398
2399
2400
2401 }
2402 }
2403
2404 return;
2405 }
2406
2407 #can we sub this? args: table_name, logic_name
2408 #also use result_set_name
2409 #would also clean all data for result set if recovery
2410 #return would be result_set
2411 #Can we extend this to incorporate InputSet parser define_sets?
2412
2413 sub get_import_ResultSet{
2414 my ($self, $anal, $table_name) = @_;
2415
2416 if (!($anal && $anal->isa("Bio::EnsEMBL::Analysis") && $anal->dbID())) {
2417 throw("Must provide a valid stored Bio::EnsEMBL::Analysis");
2418 }
2419
2420 $self->log("Getting import $table_name ResultSet for analysis:\t".$anal->logic_name());
2421
2422 my ($rset, @new_chip_channels);
2423 my $result_adaptor = $self->db->get_ResultSetAdaptor();
2424 my $logic_name = $anal->logic_name;
2425 my $status = ($logic_name eq "RawValue") ? "IMPORTED" : $logic_name;
2426
2427 if (($logic_name) eq 'RawValue' && ($table_name eq 'experimental_chip')) {
2428 throw("Cannot have an ExperimentalChip ResultSet with a RawValue analysis, either specify 'channel' or another analysis");
2429 }
2430
2431 #Build IMPORT Set for $table_name
2432 foreach my $echip (@{$self->experiment->get_ExperimentalChips()}) {
2433
2434 #clean chip import and generate rset
2435
2436 if ($table_name eq 'experimental_chip') {
2437
2438 if ($echip->has_status($status)) { #this translates to each channel have the IMPORTED_RawValue status
2439 $self->log("ExperimentalChip(".$echip->unique_id().") already has status:\t".$status);
2440 } else {
2441 $self->log("Found ExperimentalChip(".$echip->unique_id().") without status $status");
2442
2443 push @new_chip_channels, $echip;
2444 }
2445
2446 } else { #channel
2447
2448 foreach my $chan (@{$echip->get_Channels()}) {
2449
2450 if ($chan->has_status($status)) { #this translates to each channel have the IMPORTED_RawValue status
2451 $self->log("Channel(".$echip->unique_id()."_".$self->get_config('dye_freqs')->{$chan->dye()}.") already has status:\t".$status);
2452 } else {
2453 $self->log("Found Channel(".$echip->unique_id()."_".$self->get_config('dye_freqs')->{$chan->dye()}.") without status $status");
2454 push @new_chip_channels, $chan;
2455 }
2456 }
2457 }
2458
2459 if (( ! $rset) && @new_chip_channels) {
2460 my(@tmp) = @{$result_adaptor->fetch_all_by_name_Analysis($self->name()."_IMPORT", $anal)};
2461
2462 if (scalar(@tmp) > 1) {
2463 throw('Found more than one IMPORT ResultSet for '.$self->name().'_IMPORT with analysis '.$logic_name);
2464 }
2465
2466 $rset = shift @tmp;
2467
2468
2469 #do we need to throw here if not recovery?
2470 #what if we want the import result set elsewhere during the first import?
2471
2472 #if ($self->recovery()) {
2473
2474 #fetch by anal and experiment_id
2475 #Need to change this to result_set.name!
2476 # warn("add chip set handling here");
2477
2478 #my @tmp = @{$result_adaptor->fetch_all_by_Experiment_Analysis($self->experiment(), $anal)};
2479 #throw("Found more than one ResultSet for Experiment:\t".$self->experiment->name()."\tAnalysis:\t".$anal->logic_name().')' if (scalar(@tmp) >1);
2480 #$rset = $tmp[0];
2481
2482 #warn "fetching rset with ".$self->name()."_IMPORT ". $anal->logic_name;
2483
2484 #$rset = $result_adaptor->fetch_by_name_Analysis($self->name()."_IMPORT", $anal);
2485 warn("Warning: Could not find recovery ResultSet for analysis ".$logic_name) if ! $rset;
2486 #}
2487
2488 if (! $rset) {
2489 $self->log("Generating new ResultSet for analysis ".$logic_name);
2490
2491 $rset = Bio::EnsEMBL::Funcgen::ResultSet->new
2492 (
2493 -analysis => $anal,
2494 -table_name => $table_name,
2495 -name => $self->name()."_IMPORT",
2496 -feature_type => $self->feature_type(),
2497 -cell_type => $self->cell_type(),
2498 );
2499
2500 #These types should be set to NULL during the MAGE-XML validation if we have more than one type in an experiment
2501
2502 ($rset) = @{$result_adaptor->store($rset)};
2503 }
2504 }
2505 }
2506
2507 #do we need this here as we're rolling back in the read methods?
2508 #we only want to roll back those chips/channels which have not been registered
2509
2510 if ($self->recovery()) {
2511
2512 my $ec_adaptor = $self->db->get_ExperimentalChipAdaptor();
2513
2514 foreach my $cc (@new_chip_channels) {
2515
2516 #only roll back if already part of import set
2517 #Not previously registered if not
2518 if ($rset->contains($cc) && $rset->get_chip_channel_id($cc->dbID())) {
2519
2520 if ($table_name eq 'channel') {
2521 my $chan_name = $ec_adaptor->fetch_by_dbID($cc->experimental_chip_id())->unique_id()."_".
2522 $self->get_config('dye_freqs')->{$cc->dye()};
2523 $self->log("Rolling back results for $table_name:\t".$chan_name);
2524
2525 } else {
2526 $self->log("Rolling back results for $table_name:\t".$cc->unique_id);
2527 }
2528
2529 $self->rollback_results([$rset->get_chip_channel_id($cc->dbID())]);
2530 }
2531 }
2532 }
2533
2534
2535 #check whether it is present in the ResultSet and add if not
2536 if ($rset) {
2537 #ids will already be present if not rset i.e. already imported
2538
2539 foreach my $cc (@new_chip_channels) {
2540 $rset->add_table_id($cc->dbID()) if(! $rset->contains($cc));
2541 }
2542 }
2543
2544
2545
2546 if ($rset) {
2547 $result_adaptor->store_chip_channels($rset);
2548 } else {
2549 $self->log("All ExperimentalChips have status:\t$status");
2550 }
2551
2552 #this only returns a result set if there is new data to import
2553 return $rset;
2554 }
2555
2556
2557
2558 =head2 resolve_probe_data
2559
2560 Example : $self->resolve_probe_data();
2561 Description: Resolves DB probe duplicates and builds local probe cache
2562 Returntype : none
2563 Exceptions : ????
2564 Caller : general
2565 Status : At risk
2566
2567 =cut
2568
2569 sub resolve_probe_data{
2570 my $self = shift;
2571
2572 $self->log("Resolving probe data", 1);
2573
2574 warn "Probe cache resolution needs to accomodate probesets too!";
2575
2576 foreach my $array (@{$self->arrays()}) {
2577 my $resolve = 0;
2578
2579 if ($self->get_probe_cache_by_Array($array)) { #cache already generated
2580
2581 #check if we have any new unresolved array chips to add to the cache
2582 foreach my $achip (@{$array->get_ArrayChips()}) {
2583
2584 if ($achip->has_status('RESOLVED')) {
2585 $self->log("ArrayChip has RESOLVED status:\t".$achip->design_id()); #, 1);
2586 next;
2587 } else {
2588 $self->log("Found un-RESOLVED ArrayChip:\t".$achip->design_id());
2589 $resolve = 1;
2590 last;
2591 }
2592 }
2593 } else { #no cache file
2594 $resolve = 1;
2595 $self->log('No probe cache found for array '.$array->name());
2596 }
2597
2598 if ($resolve) {
2599 $self->log('Resolving array duplicates('.$array->name().') and rebuilding probe cache.', 1);
2600 $self->get_probe_cache_by_Array($array, 1); #get from DB
2601
2602 #we need ot make sure we mark cache as unresolved, so we don't use it by mistake.
2603
2604
2605
2606 my ($line, $name, $pid, @pids);
2607 #my $index = 0;
2608 my $tmp_name = '';
2609 my $tmp_id = '';
2610
2611 #miss the header
2612
2613 while ($line = $self->{'_probe_cache'}{$array->name}{'handle'}->getline()) {
2614 ($name, $pid) = split/\t/o, $line;
2615
2616 if ($name eq $tmp_name) {
2617
2618 if ($pid != $tmp_id) {
2619 push @pids, $pid;
2620 #should reset to pid here if we have x y data else undef
2621 #ignore this and force result to have x y
2622 }
2623
2624 #can't do this naymore unless we figure out how to move the line pointer
2625 #would still need to sed the file anyway, better to regen from DB?
2626 #undef $self->{'_probe_cache'}{$array->name}{'entries'}->[$i];#delete true or to be resolved duplicate
2627 } elsif ($name ne $tmp_name) { #new probe
2628 $self->tidy_duplicates(\@pids) if(scalar(@pids) > 1);
2629 $tmp_name = $name;
2630 $tmp_id = $pid;
2631 @pids = ($pid);
2632 #$index = $i + 1;
2633 }
2634 }
2635
2636 $self->tidy_duplicates(\@pids) if(scalar(@pids) > 1);
2637
2638 #rename resovled cache and reset cache handle
2639 my $cmd = 'mv '.$self->get_dir('caches').'/'.$array->name().'.probe_cache.unresolved '.
2640 $self->get_dir('caches').'/'.$array->name().'.probe_cache';
2641
2642 run_system_cmd($cmd);
2643 $self->get_probe_cache_by_Array($array); #This sets the caches
2644
2645
2646 #warn "Only generate MD5 here, as this is guranteed to be correct";
2647
2648 foreach my $achip (@{$array->get_ArrayChips()}) {
2649
2650 if (! $achip->has_status('RESOLVED')) {
2651 $self->log("Updating ArrayChip to RESOLVED status:\t".$achip->design_id());
2652 $achip->adaptor->store_status('RESOLVED', $achip);
2653 }
2654 }
2655
2656 $self->log('Finished building probe cache for '.$array->name(), 1);
2657 }
2658 }
2659
2660 $self->log('Finished resolving probe data', 1);
2661
2662 return;
2663 }
2664
2665
2666 sub tidy_duplicates{
2667 my ($self, $pids) = @_;
2668
2669 my $pfa = $self->db->get_ProbeFeatureAdaptor();
2670 my ($feature, %features);
2671
2672 foreach my $dup_id (@$pids) {
2673
2674 foreach $feature(@{$pfa->fetch_all_by_Probe_id($dup_id)}) {
2675 #can we safely assume end will be same too?
2676 push @{$features{$feature->seq_region_name().':'.$feature->start()}}, $feature;
2677 }
2678 }
2679
2680 my (@reassign_ids, @delete_ids);
2681
2682 foreach my $seq_start_key (keys %features) {
2683 my $reassign_features = 1;
2684
2685 foreach $feature(@{$features{$seq_start_key}}) {
2686
2687 if ($feature->probe_id() == $pids->[0]) {
2688 $reassign_features = 0;
2689 } else {
2690 push @delete_ids, $feature->dbID();
2691 }
2692 }
2693
2694 #This assumes that we actually have at least one element to every seq_start_key array
2695 if ($reassign_features) {
2696 my $new_fid = pop @delete_ids;
2697 push @reassign_ids, $new_fid;
2698 }
2699 }
2700
2701 #resolve features first so we don't get any orphaned features if we crash.
2702 $pfa->reassign_features_to_probe(\@reassign_ids, $pids->[0]) if @reassign_ids;
2703 $pfa->delete_features(\@delete_ids) if @delete_ids;
2704
2705 return;
2706 }
2707
2708 1;