Mercurial > repos > mahtabm > ensembl

=head1 LICENSE

  Copyright (c) 1999-2011 The European Bioinformatics Institute and
  Genome Research Limited.  All rights reserved.

  This software is distributed under a modified Apache license.
  For license details, please see

    http://www.ensembl.org/info/about/code_licence.html

=head1 CONTACT

  Please email comments or questions to the public Ensembl
  developers list at <ensembl-dev@ebi.ac.uk>.

  Questions may also be sent to the Ensembl help desk at
  <helpdesk@ensembl.org>.

=head1 NAME

Bio::EnsEMBL::Funcgen::MAGE

=head1 SYNOPSIS

my $imp = Bio::EnsEMBL::Funcgen::Importer->new(%params);
$imp->register_experiment();


=head1 DESCRIPTION

B<This program> is a base main class for all MAGE type array importers(e.g. Nimblegen).

=cut

################################################################################

package Bio::EnsEMBL::Funcgen::Parsers::MAGE;

use Bio::EnsEMBL::Funcgen::Utils::EFGUtils qw(get_date open_file run_system_cmd);
use Bio::EnsEMBL::Utils::Exception qw( throw deprecate );
use Bio::EnsEMBL::Utils::Argument qw( rearrange );
use Bio::EnsEMBL::Funcgen::Utils::Helper;
use Bio::MAGE::XMLUtils;


use File::Path;
use strict;
use vars qw(@ISA);
@ISA = qw(Bio::EnsEMBL::Funcgen::Utils::Helper);


################################################################################

=head2 new

 Description : Constructor method

 Arg  [1]    : hash containing optional attributes:

 ReturnType  : Bio::EnsEMBL::Funcgen::MAGE
 Example     : my $Exp = Bio::EnsEMBL::Nimblegen->new(%params);
 Exceptions  : throws if mandatory params are not set or DB connect fails
 Caller      : General
 Status      : Medium - potential for %params names to change, remove %attrdata?

=cut

################################################################################

sub new{
  my ($caller) = shift;
  my $class = ref($caller) || $caller;
  my $self = $class->SUPER::new(@_);

  #This needs to be an Importer!
  throw("This is base class for the experiment Bio::EnsEMBL::Funcgen::Parsers, needs to inherit from Bio::EnsEMBL::Funcgen::Importer") if(! $self->isa("Bio::EnsEMBL::Funcgen::Importer"));


  #Are we not passing any Helper params?
  #Log file etc is set in the run script

  my ($write_mage, $no_mage, $vendor)
	= rearrange(['WRITE_MAGE', 'NO_MAGE', 'VENDOR'], @_);


  #$self->{'update_xml'} = $update_xml || 0;
  $self->{'write_mage'} = $write_mage || 0;
  $self->{'no_mage'} = $no_mage || 0;
  $self->{'vendor'} = $vendor;


  if ($self->vendor ne 'NIMBLEGEN'){
	$self->{'no_mage'} = 1;
	warn "Hardcoding no_mage for non-NIMBLEGEN imports";
  }


  if($self->{'no_mage'} && $self->{'write_mage'}){
	throw('-no_mage and -write_mage options are mutually exclusive, please select just one');
  }

  return $self;
}

=head2 process_experiment_config

  Example    : $self->init_experiment_import();
  Description: Initialises import by creating working directories
               and by storing the Experiemnt
  Returntype : none
  Exceptions : warns and throws depending on recover and Experiment status
  Caller     : general
  Status     : at risk - merge with register exeriment

=cut

#This is actually processing the tab2mage file & writing XML

sub process_experiment_config{
  my $self = shift;

   #Here, this is where we need to call the a Parser from Importer to do this for only MAGE experiments
  #validate_import?

  #This is only used for the first test below.
  my $exp_adaptor = $self->db->get_ExperimentAdaptor();
  my $xml = $exp_adaptor->fetch_mage_xml_by_experiment_name($self->name);# if $self->{'write_xml'};

  #DO NOT CHANGE THIS LOGIC!
  #write mage if we specify or we don't have a the final xml or the template
  #recovery is turned on to stop exiting when previously stored chips are found from the 'write_mage' run.
  #This does mean that if you import without running the write_mage step
  #you could potentially be overwriting someone elses experiment info
  #No way of getting around this, need to make warning obvious, add to end of log!!!
  #We always want to write and update xml and ResultSets if we're running the 2nd stage of the import
  #Why would we ever want to skip the validate process?
  #Leave for now as this is working as we want it
  #But propose to remove skip functionality

  if( ! $self->{'no_mage'}){

	if($self->{'write_mage'} || !( -f $self->get_config('tab2mage_file') || $xml)){
	  $self->{'write_mage'} = 1;
	  $self->backup_file($self->get_config('tab2mage_file'));
	}
	#elsif($xml && (! $self->{'update_xml'})){#Changed this so we always update
	#elsif(! $self->{'update_xml'}){


	  #Here, we need to always update_xml
	  #If we are doing the 2nd stage
	  #Currently this is skipping as we haven't explicitly set it
	  #To remove this...
	  #what we need to do is check that we don't test for update_xml,
	  # i.e. assuming that we're running the second stage of the import.
	  # Therefore we need a boolean to set whether it is the first stage..else update_xml implicit
	  # write mage is explicit flag
	  # Or if we have not tab2mage file?
	  # we can then override this explicitly with update_xml?
	  # WE're never likely edit the xml directly, so we always want to validate and update
	  # so write mage flag become update_experiment?  No this is no obvious behaviour
	  # We need to warn about removing the write_mage flag after we have updated it
	  # Otherwise we will never get to 2nd stage


	#No mage is still valid as we may want to jus import and experiment
	#Before receiving correct meta data
	#When we can then rerun the import with -write_mage to update the resultsets

	#  $self->{'recover'} = 1;
	#  $self->{'skip_validate'} = 1;
	#}
	elsif( -f $self->get_config('tab2mage_file')){#Run Tab2Mage

	  $self->backup_file($self->get_config('mage_xml_file'));
	  my $cmd = 'tab2mage.pl -e '.$self->get_config('tab2mage_file').' -k -t '.$self->get_dir('output').
		' -c -d '.$self->get_dir('results');

	  $self->log('Reading tab2mage file');
	  my $t2m_exit_code = run_system_cmd($cmd, 1);#no exit flag due to non-zero exit codes
	  warn "tab2mage exit code is  $t2m_exit_code";

	  if(! ($t2m_exit_code > -1) && ($t2m_exit_code <255)){
		throw("tab2mage failed.  Please check and correct:\t".$self->get_config('tab2mage_file')."\n...and try again");
	  }

	  $self->{'recover'} = 1;
	}
  }

  return;
}

=heead init_tab2mage_export

  Example    : $self->init_tab2mage_export;
  Description: Writes the standard experiment section of the tab2mage file
  Returntype : FileHandle
  Exceptions : ???
  Caller     : general
  Status     : at risk

=cut

sub init_tab2mage_export{
  my $self = shift;

  $self->backup_file($self->get_config('tab2mage_file')) if(-f $self->get_config('tab2mage_file'));

  my $t2m_file = open_file($self->get_config('tab2mage_file'), '>');

  #reformat this
  my $exp_section = "experiment section\ndomain\t".(split/@/, $self->contact())[1]."\naccession\t\n".
	"quality_control\tbiological_replicate\nexperiment_design_type\tbinding_site_identification\n".
	  "name\t".$self->name()."\nrelease_date\t\nsubmission_date\t\nsubmitter\t???\n".
		"submitter_email\t???\ninvestigator\t???\ninvestigator_email\t???\norganization\t???\naddress\t".
		  "???\npublication_title\t\nauthors\t\njournal\t\nvolume\t\nissue\t\npages\t\nyear\t\npubmed_id\t\n";

  my $protocol_section = "Protocol section\naccession\tname\ttext\tparameters\n";

  foreach my $protocol(sort (keys %{$self->get_config('protocols')})){
	$protocol_section .= $self->get_config('protocols')->{$protocol}->{'accession'}.
	  "\t".$self->get_config('protocols')->{$protocol}->{'name'}.
		"\t".$self->get_config('protocols')->{$protocol}->{'text'}."\t";

	$protocol_section .= (defined $self->get_config('protocols')->{$protocol}->{'parameters'}) ?
	  $self->get_config('protocols')->{$protocol}->{'parameters'}."\t\n" : "\t\n";
  }

  #File[raw]	Array[accession]	Array[serial]	Protocol[grow]	Protocol[treatment]	Protocol[extraction]	Protocol[labeling]	Protocol[hybridization]	Protocol[scanning]	BioSource	Sample	Extract	LabeledExtract	Immunoprecipitate	Hybridization	BioSourceMaterial	SampleMaterial	ExtractMaterial	LabeledExtractMaterial	Dye	BioMaterialCharacteristics[Organism]	BioMaterialCharacteristics[BioSourceType]	BioMaterialCharacteristics[StrainOrLine]	BioMaterialCharacteristics[CellType]	BioMaterialCharacteristics[Sex]	FactorValue[StrainOrLine]	FactorValue[Immunoprecipitate]


  #Need to do this bit better?
  #have array of fields.  We can then populate a hash in the read method based on field names, then use the array to print in order

  my $hyb_header = "\nHybridization section\n".join("\t", @{$self->hybridisation_fields()});

  print $t2m_file $exp_section."\n".$protocol_section."\n".$hyb_header."\n";

  return $t2m_file;
}


#Move to MAGE package?

sub hybridisation_fields{
  my $self = shift;

  return ['File[raw]', 'Array[accession]', 'Array[serial]',
		  (map 'Protocol['.$_.']', (sort (keys %{$self->get_config('protocols')}))),
		  'BioSource', 'Sample', 'Extract', 'LabeledExtract', 'Immunoprecipitate', 'Hybridization',
		  'BioSourceMaterial', 'SampleMaterial', 'ExtractMaterial', 'LabeledExtractMaterial',
		  'Dye', 'BioMaterialCharacteristics[Organism]', 'BioMaterialCharacteristics[BioSourceType]',
		  'BioMaterialCharacteristics[StrainOrLine]', 'BioMaterialCharacteristics[CellType]',
		  'BioMaterialCharacteristics[Sex]', 'FactorValue[StrainOrLine]', 'FactorValue[Immunoprecipitate]'];
}


#=head2 register_experiment
#
#  Example    : $imp->register_experiment()
#  Description: General control method, performs all data import and normalisations
#  Arg [1]    : optional - dnadb DBAdaptor
#  Returntype : none
#  Exceptions : throws if arg is not Bio::EnsEMBL::DBSQL::DBAdaptor
#  Caller     : general
#  Status     : Medium
#
#=cut

#write/validate_mage

sub write_validate_experiment_config{
  my $self = shift;


  if($self->{'write_mage'} || $self->{'no_mage'}){
	$self->read_data("array");

	if(! $self->{'no_mage'}){
	  $self->log("PLEASE CHECK AND EDIT AUTOGENERATED TAB2MAGE FILE:\t".$self->get_config('tab2mage_file'));
	  #we could make this print only if it was set by the user, not by the Importer
	  $self->log('REMEMBER TO REMOVE -write_mage FLAG BEFORE UPDATING');
	  exit;
	}
  }
  elsif(! $self->{'no_mage'}){#This should be a no_channel flag, set dependent on import mode(gff_chip, gff_chan)
	#Need to accomodate chip level imports in validate?

	if (! $self->{'skip_validate'}){

	  $self->log("Validating mage file:\t".$self->get_config('mage_xml_file'));


	  #Updating ResultSets:
	  #Given that we might want to add a chip to an experiment we will also need to update the tab2MAGE
	  #mage_xml and ResultSets accordingly.

	  #This should happen if we specify update_xml
	  #Should recovery also always force update?
	  #Considering the two run modes, write tab2mage & validate and import
	  #There is a subtle difference between recovery and update mage_xml
	  #Do we always run in recovery mode for the validate&import step?
	  #Yes we do, so can't guarantee the this means we want to update.
	  #So we need to change update_xml to update to reflect the changed functionality on ResultSets

	  #If we run an update without on then chips will be loaded but xml and ResultSets will not be altered :(
	  #If we're running the 2nd stage we should always be updating the xml anyway!!!!
	  #As there is no reason to rerun the validate & import step without it.(unless we're debugging maybe)
	  #So why should we ever run without it?

	  #To update ResultSets we validate as normal and then update where appropriate
	  #What has precedence? Replicate name?
	  #Update echip types as appropriate
	  #What if this invalidates original rsets?
	  #Then list sets not covered for removal by script?


	  my (%echips, @log);
	  my $rset_adaptor = $self->db->get_ResultSetAdaptor;
	  my $chan_anal = $self->db->get_AnalysisAdaptor->fetch_by_logic_name('RawValue');

	  #need to change this to default analysis
	  #There we issues with setting VSN_GLOG as an env var
	  #as this is tested for and the norm was skipped for some reason?
	  my $chip_anal = $self->db->get_AnalysisAdaptor->fetch_by_logic_name($self->norm_method());

	  #Try import sets like this first, so we know ther is new data
	  my $chan_rset = $self->get_import_ResultSet($chan_anal, 'channel');
	  my $rset =  $self->get_import_ResultSet($chip_anal, 'experimental_chip');


	  #Else get them anyway and log
	  if(! $rset){

		if($chan_rset){
		  $self->log('Identified partial Channel only import, updating MAGE-XML');
		}
		else{
		  ($chan_rset) = @{$rset_adaptor->fetch_all_by_name_Analysis($self->experiment->name.'_IMPORT', $chan_anal)};
		  #Don't need to test for >1 here as this has already been done in get_import_ResultSet
		  $self->log('All ExperimentalChips imported, updating MAGE-XML only');
		}

		($rset) = @{$rset_adaptor->fetch_all_by_name_Analysis($self->experiment->name.'_IMPORT', $chip_anal)};
	  }


  #This will never happen now due to the change tab2mage rules in init_experiment
  #Remove?
  if(! $rset){
	throw('Cannot find ResultSet, are you trying to import a new experiment which already has a tab2mage file present?  Try removing the file, or specifying the -write_mage flag to parse_and_import.pl');
  }

  if(! -l $self->get_dir('output').'/MAGE-ML.dtd'){
	system('ln -s '.$ENV{'EFG_DATA'}.'/MAGE-ML.dtd '.$self->get_dir('output').'/MAGE-ML.dtd') == 0 ||
	  throw('Failed to link MAGE-ML.dtd');
  }


  $self->log('VALIDATING MAGE XML');
  my $reader = Bio::MAGE::XML::Reader->new();
	  my $mage_xml ||= $self->get_config('mage_xml_file');
  $self->{'mage'} = $reader->read($mage_xml);

  #this should only ever return 1 for an import
  foreach my $mage_exp(@{$self->{'mage'}->getExperiment_package->getExperiment_list()}){

	if($mage_exp->getName() ne $self->name()){
	  $self->log('MAGE experiment name ('.$mage_exp->getName().') does not match import name ('.$self->name().')');
	}

	#add more experiment level validation here?

	foreach my $assay (@{$mage_exp->getBioAssays()}){

	  if($assay->isa('Bio::MAGE::BioAssay::PhysicalBioAssay')){#channel
		$self->log('Validating PhysicalBioAssay "'.$assay->getName()."'\n");#hyb name(this is the file name for measured assays

		my $bioassc = $assay->getBioAssayCreation();#This is a Hybridisation
		my $array = $bioassc->getArray();#this is an ArrayChip
		my $design_id = $array->getArrayDesign->getIdentifier();
		my $chip_uid = $array->getArrayIdentifier();


		foreach my $echip(@{$rset->get_ExperimentalChips()}){


		  if($echip->unique_id() eq $chip_uid){
			$self->log("Found ExperimentalChip:\t".$chip_uid);

			if(! exists $echips{$chip_uid}){
			  $echips{$chip_uid} = {(
									 total_biorep            => undef,
									 total_biotechrep        => undef,
									 experimental_biorep     => undef,
									 experimental_biotechrep => undef,
									 total_dye               => undef,
									 experimental_dye        => undef,
									 cell_type               => undef,
									 feature_type            => undef,
									)};
			}

			#Validate ArrayChip
			my ($achip) = @{$self->db->get_ArrayChipAdaptor->fetch_all_by_ExperimentalChips([$echip])};

			if($achip->design_id() ne $design_id){
			  push @log, "ArrayDesign Identifier (${design_id}) does not match ArrayChip design ID (".
				$achip->design_id().")\n\tSkipping channel and replicate validation";
			  #skip the channel/replicate validation here?
			}
			else {			#validate channels and replicate names

			  foreach my $src_biomat (@{$bioassc->getSourceBioMaterialMeasurements()}) { #Channel materials(X1)?
				my $biomat = $src_biomat->getBioMaterial();	#LabelledExtract (IP/Control)
				#we could sub this passing $echip and biomat?
				#messy to pass regexs and populate correct echip hash attrs
				#also messy to populate log
				#keeping nested loop also prevents further obfuscation
				#do we need to do all the defined checks, or maybe just the first one?
				#Then we can skip all following warning?

				foreach my $treat (@{$biomat->getTreatments()}) {
				  #As there is effectively one more level of material extraction for the IP channel
				  #this loop will returns materials an iteration out of sync for each channel

				  foreach my $ssrc_biomat (@{$treat->getSourceBioMaterialMeasurements()}) {	#Channel measurement(x1)
					my $sbiomat = $ssrc_biomat->getBioMaterial();
					#This will either be techrep name for control of IP name for experimental channel
					#SOM0035_BR1_TR2 IP  #Immunoprecicpitate
					#SOM0035_BR1_TR2     #Extract

					if ($sbiomat->getName() =~ /BR[0-9]+_TR[0-9]+$/o) { #Total

					  if (! defined $echips{$chip_uid}{'total_biotechrep'}) {
						$echips{$chip_uid}{'total_biotechrep'} = $sbiomat->getName();
					  }
					  else{
						push @log, "Found two TOTAL Channels on same chip with biotechreps:\t".$sbiomat->getName().
						  " and ".$echips{$chip_uid}{'total_biotechrep'};
					  }
					}else{#Experimental

					  #get feature type from assay
					  my $fv_ref = $assay->getBioAssayFactorValues();
					  if(! defined $fv_ref){
						throw('No FactorValues found, you must populate the "Immunoprecipitate" field. Maybe you forgot to specify -feature_type?');
					  }

					  my ($feature_type);

					  foreach my $fvalue(@{$fv_ref}){

						if($fvalue->getValue()->getCategory() eq 'Immunoprecipitate'){
						  $feature_type = $fvalue->getName();
						  $feature_type =~ s/anti\s*-\s*//;
						  $feature_type =~ s/\s*antibody\s*//;
						}
					  }
					  $echips{$chip_uid}{'feature_type'} = $feature_type;
					}

					foreach my $ttreat (@{$sbiomat->getTreatments()}) {

					  foreach my $tsrc_biomat (@{$ttreat->getSourceBioMaterialMeasurements()}) {
						my $tbiomat = $tsrc_biomat->getBioMaterial();
						#SOM0035_BR1_TR2     #Extract (exp)
						#SOM0035_BR1              #Sample (total)

						if ($tbiomat->getName() =~ /BR[0-9]+_TR[0-9]+$/o) { #experimental

						  if (! defined $echips{$chip_uid}{'experimental_biotechrep'}) {
							$echips{$chip_uid}{'experimental_biotechrep'} = $tbiomat->getName();
						  }
						  else{
							push @log, "Found two EXPERIMENTAL Channels on same chip with biotechreps:\t".$tbiomat->getName().
							  " and ".$echips{$chip_uid}{'experimental_biotechrep'};
						  }

						  my $dye = $biomat->getLabels()->[0]->getName();

						  foreach my $chan (@{$echip->get_Channels()}) {

							if ($chan->type() eq 'EXPERIMENTAL') {

							  if (uc($dye) ne uc($chan->dye())) {
								push @log, "EXPERIMENTAL channel dye mismatch:\tMAGE = ".uc($dye).' vs DB '.uc($chan->dye);
							  } else {
								$echips{$chip_uid}{'experimental_dye'} = uc($dye);
							  }
							}
						  }
						}
						else { #control

						  if (! defined $echips{$chip_uid}{'total_biorep'}) {
							$echips{$chip_uid}{'total_biorep'} = $tbiomat->getName();
						  }
						  else{
							push @log, "Found two TOTAL Channels on same chip with biotechreps:\t".$tbiomat->getName().
							  " and ".$echips{$chip_uid}{'total_biorep'};
						  }

						  my $dye = $biomat->getLabels()->[0]->getName();

						  foreach my $chan (@{$echip->get_Channels()}) {

							if ($chan->type() eq 'TOTAL') {

							  if (uc($dye) ne uc($chan->dye())) {
								push @log, "TOTAL channel dye mismatch:\tMAGE = ".uc($dye).' vs DB '.uc($chan->dye);
							  }
							  else {
								$echips{$chip_uid}{'total_dye'} = uc($dye);
							  }
							}
						  }
						}
						#could do one more iteration and get Source and FeatureType?
						#we should really extend this, and then update the EC cell_type and feature_types
						#these features might not be biotmats tho...need to check


						foreach my $ftreat (@{$tbiomat->getTreatments()}) {

						  foreach my $fsrc_biomat (@{$ftreat->getSourceBioMaterialMeasurements()}) {
							my $fbiomat = $fsrc_biomat->getBioMaterial();
							#EXPERIMENTAL - biorep
							#TOTAL        - source/cell type
							my $cell_type;

							if($fbiomat->getName() =~ /BR[0-9]+$/o){#EXPERIMETNAL

							  if(! defined $echips{$chip_uid}{'experimental_biorep'}){
								$echips{$chip_uid}{'experimental_biorep'} = $fbiomat->getName();
							  }
							  else{
								push @log, "Found two Experimental Channels on same chip with bioreps:\t".$fbiomat->getName().
								  " and ".$echips{$chip_uid}{'experimental_biorep'};
							  }


							  #last treatment/measurement/biomat level should go here
							  #as TOTAL channel does not have another level and will fail
							  foreach my $xtreat (@{$fbiomat->getTreatments()}) {

								foreach my $xsrc_biomat (@{$xtreat->getSourceBioMaterialMeasurements()}) {
								  my $xbiomat = $xsrc_biomat->getBioMaterial();

								  foreach my $char(@{$xbiomat->getCharacteristics()}){
									$cell_type = $char->getValue() if($char->getCategory() eq 'CellType');
								  }
								}
							  }

							}else{#this should be BioSource
							  #which should have CellType as characteristic
							  #we could change tab2mage and have this as a factor value,
							  #but don't want to start messing with "standard" format

							  foreach my $char(@{$fbiomat->getCharacteristics()}){
								$cell_type = $char->getValue() if($char->getCategory() eq 'CellType');
							  }
							}

							#can have cell_type validation here
							if(! defined $echips{$chip_uid}{'cell_type'}){
							  $echips{$chip_uid}{'cell_type'} = $cell_type;
							}
							elsif( $echips{$chip_uid}{'cell_type'} ne $cell_type){
							  push @log, "Found Channels on same chip (${chip_uid}) with different cell types:\t".
								$cell_type." and ".$echips{$chip_uid}{'cell_type'};
							}
						  }
						}
					  }
					}
				  }
				}
			  }
			}
		  }						#end of echip
		}						#end of foreach echip
	  }							#end of physbioassay
	}							#end of foreach assay
  }								#end of foreach exp


  #we should fail here with log before we update the result sets

   #we need to build rep names
  #we're currently using sample labels, in the tab2mage file
  #altho' previous sets have been using exp name
  #these have been manually patched afterwards

  #More desirable to have exp name as rset name, but no way of doing BR validation
  #based on sample label, if we don't have it in the tab2mage
  #if we change it in the DB then we need to update the tab2mage

  #no way to do this when generating tab2mage as the user hasn't yet defined the reps
  #we could just make reps based on sample labels
  #then we just assume that alterations made by the user are correct
  #as we can no longer validate using sample labels
  #can still validate using cell/feature type

  #no longer need vendor specific validation as this will be done in tab2mage generation


  #We need to validate reps here
  #then update ec records as appropriate and then create rsets

  my (%bio_reps, %tech_reps);
  my $ct_adaptor = $self->db->get_CellTypeAdaptor();
  my $ft_adaptor = $self->db->get_FeatureTypeAdaptor();

#select rs.*, ec.*, c.* from result_set rs, chip_channel cc, channel c, experimental_chip ec where rs.result_set_id=cc.result_set_id and cc.table_name='experimental_chip' and cc.table_id=ec.experimental_chip_id and cc.table_id=c.experimental_chip_id order by name;

  foreach my $echip (@{$rset->get_ExperimentalChips()}) {

	my ($biorep, $biotechrep);

	if (! exists $echips{$echip->unique_id()}) {
	  push @log, "No MAGE entry found for ExperimentalChip:\t".$echip->unique_id();
	}
	else {

	  foreach my $chan_type('total', 'experimental'){

		$biorep = $echips{$echip->unique_id()}{$chan_type.'_biorep'};
		$biotechrep = $echips{$echip->unique_id()}{$chan_type.'_biotechrep'};

		if (! defined $biotechrep) {
		  push @log, 'ExperimentalChip('.$echip->unique_id().') Extract field do not meet naming convention(SAMPLE_BRN_TRN)';
		}							#! defined biorep? will never occur at present
		elsif ($biotechrep !~ /\Q$biorep\E/) {
		  push @log, "Found Extract(techrep) vs Sample(biorep) naming mismatch\t${biotechrep}\tvs\t$biorep";
		}

		if ( ! $echips{$echip->unique_id()}{$chan_type.'_dye'}) {
		  push @log, "No ".uc($chan_type)." channel found for ExperimentalChip:\t".$echip->unique_id();
		}

	  }

	  #Is this is really implicit in the test above
	  if($echips{$echip->unique_id()}{'experimental_biorep'} ne $echips{$echip->unique_id()}{'total_biorep'}){
		push @log, "Found biorep mismatch between channels of ExperimentalChip ".$echip->unique_id().":\n".
		  "\tEXPERIMENTAL\t".$echips{$echip->unique_id()}{'experimental_biorep'}."\tTOTAL\t".
			$echips{$echip->unique_id()}{'total_biorep'};
	  }

	  #Is this is really implicit in the test above
	  if($echips{$echip->unique_id()}{'experimental_biotechrep'} ne $echips{$echip->unique_id()}{'total_biotechrep'}){
		push @log, "Found biotechrep mismatch between channels of ExperimentalChip ".$echip->unique_id().":\n".
		  "\tEXPERIMENTAL\t".$echips{$echip->unique_id()}{'experimental_biotechrep'}."\tTOTAL\t".
			$echips{$echip->unique_id()}{'total_biotechrep'};
	  }


	}


	#Now we need to validate ec has same feature/cell type as other ecs in this br
	#this does not handle import sets which ARE allowed to have same name but different types

	#warn "Processing ".$echip->unique_id()." $biorep $biotechrep";


	if(exists $bio_reps{$biorep}){


	  if(! defined $bio_reps{$biorep}{'cell_type'}){
		push @log, "Found undefined CellType for biorep $biorep";
	  }
	  elsif($bio_reps{$biorep}{'cell_type'}->name() ne  $echips{$echip->unique_id()}{'cell_type'}){
		push @log, "Found CellType mismatch between $biorep and ExperimentalChip ".$echip->unique_id();
	  }


	  if(! defined $bio_reps{$biorep}{'feature_type'}){
		push @log, "Found undefined FeatureType for biorep $biorep";
	  }
	  elsif($bio_reps{$biorep}{'feature_type'}->name() ne  $echips{$echip->unique_id()}{'feature_type'}){
		push @log, "Found FeatureType mismatch between $biorep and ExperimentalChip ".$echip->unique_id();
	  }

	  #warn "$biorep exists with\t".$bio_reps{$biorep}{'cell_type'}->name().' '.$bio_reps{$biorep}{'feature_type'}->name();

	  #We need to set the tech rep here too!
	  #Do we need to validate this also, as above.
	  #This would be overkill due to the inherant nature of the TR to BR relationship

	  if(! exists $tech_reps{$biotechrep}){
		 $tech_reps{$biotechrep}{'cell_type'} = $bio_reps{$biorep}{'cell_type'};
		 $tech_reps{$biotechrep}{'feature_type'} = $bio_reps{$biorep}{'feature_type'};
	  }


	}else{

	  #warn "Creating new BR $biorep and TR $biotechrep";

	  if(defined $echips{$echip->unique_id()}{'cell_type'}){

		my $cell_type = $ct_adaptor->fetch_by_name($echips{$echip->unique_id()}{'cell_type'});

		if(! defined $cell_type){
		  push @log, 'CellType '.$echips{$echip->unique_id()}{'cell_type'}.' does not exist in the database, please use the import_type.pl script';
		}else{
		  $bio_reps{$biorep}{'cell_type'} = $cell_type;
		  $tech_reps{$biotechrep}{'cell_type'} = $cell_type;
		#  warn "Setting ".$echip->unique_id()." $biorep $biotechrep ".$cell_type->name;
		}
	  }else{
		warn "No CellType specified for ExperimentalChip:\t".$echip->unique_id()."\n";
	  }


	  if(defined $echips{$echip->unique_id()}{'feature_type'}){
		my $feature_type = $ft_adaptor->fetch_by_name($echips{$echip->unique_id()}{'feature_type'});

		if(! defined $feature_type){
		  push @log, 'FeatureType '.$echips{$echip->unique_id()}{'feature_type'}.' does not exist in the database, please use the import_type.pl script';
		}
		else{
		  $bio_reps{$biorep}{'feature_type'} = $feature_type;
		  $tech_reps{$biotechrep}{'feature_type'} = $feature_type;

		  #warn "Setting ".$echip->unique_id()." $biorep $biotechrep ".$feature_type->name;
		}
	  }else{
		warn "No FeatureType specified for ExperimentalChip:\t".$echip->unique_id()."\n";
	  }
	}

	push @{$tech_reps{$biotechrep}{'echips'}}, $echip->unique_id();
	push @{$bio_reps{$biorep}{'echips'}}, $echip->unique_id();
  }


  if (@log) {
	$self->log("MAGE VALIDATION REPORT\n::\t".join("\n::\t", @log));
	throw("MAGE VALIDATION FAILED\nPlease correct tab2mage file and try again:\t".$self->get_config('tab2mage_file'));
  } else {
	$self->log('MAGE VALDIATION SUCCEEDED');
  }


  #we also need to build the tech rep results sets(not displayable)
  #do we need to have result sets for each biorep too?
  #update ExperimentalChip replicate info
  my (%rsets);
  my %types = (
			   feature => {},
			   cell    => {},
			  );


  #This needs to update and split the import/top level sets so they are of same types
  #update ec type here as we have ec context
  #careful not to update multiple times, just once for each ec

  my $eca = $self->db->get_ExperimentalChipAdaptor();


  foreach my $echip (@{$rset->get_ExperimentalChips()}) {
	my ($cell_type, $feature_type);

	#Set biorep info and rset
	foreach my $biorep (keys %bio_reps){

	  foreach my $chip_uid(@{$bio_reps{$biorep}{'echips'}}){

		if($chip_uid eq $echip->unique_id()){
		  $echip->biological_replicate($biorep);
		  $cell_type = $bio_reps{$biorep}{'cell_type'};
		  $feature_type = $bio_reps{$biorep}{'feature_type'};

		  if(! defined $rsets{$biorep}){

			$rsets{$biorep} = Bio::EnsEMBL::Funcgen::ResultSet->new
			  (
			   -NAME         => $biorep,#this may not be unique, prepend with exp name? Force method to use Experiment_and_name?
			   -ANALYSIS     => $rset->analysis(),
			   -TABLE_NAME   => 'experimental_chip',
			   -FEATURE_TYPE => $feature_type,
			   -CELL_TYPE    => $cell_type,
			  );

			#record cell and feature types
			$types{'feature'}{$feature_type->name()} = $feature_type;
			$types{'cell'}{$cell_type->name()} = $cell_type;
			$self->log("Created BioRep ResultSet:\t".$rsets{$biorep}->log_label);
		  }

		  $rsets{$biorep}->add_table_id($echip->dbID(), $rset->get_chip_channel_id($echip->dbID()));
		}
	  }
	}

	#reset echip types
	$echip->feature_type($feature_type);
	$echip->cell_type($cell_type);


	#set tech rep info and rset
	foreach my $techrep(keys %tech_reps){

	  foreach my $chip_uid(@{$tech_reps{$techrep}{'echips'}}){

		if($chip_uid eq $echip->unique_id()){
		  $echip->technical_replicate($techrep);

		  if(! defined $rsets{$techrep}){
			$rsets{$techrep} = Bio::EnsEMBL::Funcgen::ResultSet->new
			  (
			   -NAME       => $techrep,#this may not be unique, prepend with exp name? Force method to use Experiment_and_name?
			   -ANALYSIS   => $rset->analysis(),
			   -TABLE_NAME => 'experimental_chip',
			   -FEATURE_TYPE => $tech_reps{$techrep}{'feature_type'},
			   -CELL_TYPE    => $tech_reps{$techrep}{'cell_type'},
			  );

			$self->log("Created TechRep ResultSet:\t".$rsets{$techrep}->log_label);
		  }
		  $rsets{$techrep}->add_table_id($echip->dbID(), $rset->get_chip_channel_id($echip->dbID()));
		}
	  }
	}

	$echip->adaptor->update_replicate_types($echip);#store rep info
  }


  ### Reset/Update/Clean import sets type fields
  my $sql;

  if(scalar keys %{$types{'feature'}} >1){
	$self->log('Resetting IMPORT FeatureType to NULL for multi-FeatureType Experiment');
	$sql = "UPDATE result_set set feature_type_id='NULL' where result_set_id in (".$rset->dbID().', '.$chan_rset->dbID().')';

  }else{
	my ($ftype) = values %{$types{'feature'}};

	if(! defined $rset->feature_type()){
	  $self->log('Updating IMPORT FeatureType to '.$ftype->name());
	  $sql = "UPDATE result_set set feature_type_id=".$ftype->dbID()." where result_set_id in (".$rset->dbID().', '.$chan_rset->dbID().')';
	}
	elsif($rset->feature_type->dbID ne $ftype->dbID()){
	  $self->log('WARNING: FeatureType mismatch.  Updating IMPORT FeatureType('.$rset->feature_type->name().') to match meta('.$ftype->name.')');
	  $sql = "UPDATE result_set set feature_type_id=".$ftype->dbID()." where result_set_id in (".$rset->dbID().', '.$chan_rset->dbID().')';

	}
  }

  $self->db->dbc->do($sql) if $sql;

  undef $sql;

  if(scalar keys %{$types{'cell'}} >1){
	$self->log('Resetting IMPORT CellType to NULL for multi-CellType Experiment');
	my $sql = "UPDATE result_set set cell_type_id='NULL' where result_set_id in (".$rset->dbID().', '.$chan_rset->dbID().')';
  }else{
	my ($ctype) = values %{$types{'cell'}};

	if(! defined $rset->cell_type()){
	  $self->log('Updating IMPORT CellType to '.$ctype->name());
	  $sql = "UPDATE result_set set cell_type_id=".$ctype->dbID()." where result_set_id in (".$rset->dbID().', '.$chan_rset->dbID().')';
	}
	elsif($rset->cell_type->dbID ne $ctype->dbID()){
	  $self->log('WARNING: CellType mismatch.  Updating IMPORT CellType('.$rset->cell_type->name().') to match meta('.$ctype->name.')');
	  $sql = "UPDATE result_set set cell_type_id=".$ctype->dbID()." where result_set_id in (".$rset->dbID().', '.$chan_rset->dbID().')';
	}
  }

  $self->db->dbc->do($sql) if $sql;

  ### Generate new top level sets here based on br type combos
  #we risk duplicating sets here if import set is set to one cell/featuretype
  #duplicate anyway, as import is really just for easy tracking of all chips during import

  my %toplevel_sets;
  my $toplevel_cnt = 1;
  #could tidy up toplevel_sets implmentation

  foreach my $new_rset(values %rsets){

	my $ftype_name = (defined $new_rset->{'feature_type'}) ? $new_rset->{'feature_type'}->name() : undef;
	my $ctype_name = (defined $new_rset->{'cell_type'}) ? $new_rset->{'cell_type'}->name() : undef;

	if(! exists $toplevel_sets{$ftype_name}){
	  $toplevel_sets{$ftype_name} = {};
	  $toplevel_sets{$ftype_name}{'feature_type'} = $new_rset->{'feature_type'};
	}


	if(! exists $toplevel_sets{$ftype_name}{$ctype_name}){
	  $toplevel_sets{$ftype_name}{$ctype_name}{'cell_type'} = $new_rset->{'cell_type'};
	  $toplevel_sets{$ftype_name}{$ctype_name}{'rsets'} = [$new_rset];
	}else{
	  push @{$toplevel_sets{$ftype_name}{$ctype_name}{'rsets'}}, $new_rset;
	}
  }


  #build toplevel sets for each feature/cell type combo using constituent rsets
  foreach my $ftype_name(keys %toplevel_sets){

	foreach my $ctype_name(keys %{$toplevel_sets{$ftype_name}}){

	  next if $ctype_name eq 'feature_type';#skip feature type

	  #we need to give these a different key so we're not overwriting in the rset hash
	  $rsets{$self->experiment->name().'_'.$toplevel_cnt} = Bio::EnsEMBL::Funcgen::ResultSet->new
		(
		 -NAME       => $self->experiment->name(),
		 -ANALYSIS   => $rset->analysis(),
		 -TABLE_NAME => 'experimental_chip',
		 -FEATURE_TYPE => $toplevel_sets{$ftype_name}{'feature_type'},
		 -CELL_TYPE    => $toplevel_sets{$ftype_name}{$ctype_name}{'cell_type'},
		);

	  $self->log("Created toplevel ResultSet for:\t". $rsets{$self->experiment->name().'_'.$toplevel_cnt}->log_label);

	  #add consituent table ids
	  foreach my $new_rset(@{$toplevel_sets{$ftype_name}{$ctype_name}{'rsets'}}){

		foreach my $ec_id(@{$new_rset->table_ids()}){

		  #Only add it if it has not already been added
		  if(!  $rsets{$self->experiment->name().'_'.$toplevel_cnt}->get_chip_channel_id($ec_id)){
			$rsets{$self->experiment->name().'_'.$toplevel_cnt}->add_table_id($ec_id, $new_rset->get_chip_channel_id($ec_id));
		  }
		}
	  }
	  $toplevel_cnt++;
	}
  }

  #ResultSet update strategy
  #To avoid messyness in resolving result_set differences
  #Simply delete all that are not used as supporting sets
  #and load new ones, log old supporting rsets for manual
  #reassignment and rollback.
  #If we have clash between an old set and a new set, rename old
  #set and log
  #We might not always have the previous data files.
  #But we might want to maintain all the previous rsets and just add a new one
  #At present this would require acquiring the previous Tab2Mage file
  #and adding the new data to it.
  #We could do with a way to merge data already in the DB with new meta data to form a new Tab2Mage file
  #and validate that


  my @previous_rep_sets;
  my @supporting_rset_dsets;


  #Get non-import Sets
  map {push @previous_rep_sets, $_ if $_->name !~ /_IMPORT$/}
	@{$rset_adaptor->fetch_all_by_Experiment_Analysis($self->experiment, $chip_anal)};


  #rollback_ResultSet if possible?
	  #This is just checking if they are supporting, not actually rolling them back
  if(@previous_rep_sets){
	$self->log('Found previously stored ResultSets');

	foreach my $prev_rset(@previous_rep_sets){
	  #This should not rollback anything, just return skipped sets
	  #i.e. sets which have a product feature set
	  #It also used to delete the supporting set records which maybe important for redefining the DataSet below
	  my $rset_dset = $self->rollback_ResultSet($prev_rset);
	  push @supporting_rset_dsets, $rset_dset if @$rset_dset;
	}
  }

  #Note: If we remove chips from an experiment, they are only removed from the non-import sets
  #To fully remove them, you need to use the rollback_experiment.pl script with -chip_ids
  #can we log this in get_import_ResultSet?

  $self->log('Storing ResultSets');

  #Store new tech, biol and toplevel type rsets
  foreach my $new_rset(values %rsets){
	my $replace_txt;

	#Rename old set if we have a name/anal/type clash
	foreach my $prs(@supporting_rset_dsets){

	  my ($pset, $dset) = @$prs;

	  if($pset->log_label eq $new_rset->log_label){
		my $new_name = "OLD_".$rset->log_label;
		$self->log("Found update supporting ResultSet clash, renaming to:\t${new_name}");
		$self->unlink_ResultSet_DataSet($rset, $dset, $new_name);

		#This pset dbID has already been removed
		#Will get updated with new rset dbID when updating DataSet
		$replace_txt = 'Proposed ResultSet(dbID) replacement for DataSet('.$dset->name."):\t".$pset->dbID.' > ';
	  }
	}


	$new_rset->add_status('DAS_DISPLAYABLE');
	my ($new_rset) = @{$rset_adaptor->store($new_rset)};

	if(defined $replace_txt){
	  $self->log($replace_txt.$new_rset->dbID);
	}
  }

  my $xml_file = open_file($self->get_config('mage_xml_file'));

  #slurp in changing separator to null so we get it all in one string.
  $self->experiment->mage_xml(do{ local ($/); <$xml_file>});
  close($xml_file);

  $self->experiment($self->db->get_ExperimentAdaptor->update_mage_xml_by_Experiment($self->experiment()));
  	}
  }

  return;
}


1;
author	mahtabm
date	Thu, 11 Apr 2013 02:01:53 -0400
parents
children