Toolbar

Summary

Bio::EnsEMBL::Compara::Production::EPOanchors::ExonerateAnchors

Package variables

Included modules

Inherit

Bio::EnsEMBL::Hive::Process

Synopsis

$exonate_anchors->fetch_input();
$exonate_anchors->run();
$exonate_anchors->write_output(); writes to database

Description

Given a database with anchor sequences and a target genome. This modules exonerates
the anchors against the target genome. The required information (anchor batch size,
target genome file, exonerate parameters are provided by the analysis, analysis_job
and analysis_data tables

Methods

Methods description

None available.

Methods code

sub analysis_data_id {

	my $self = shift;
	if (@_) {
		$self->{_analysis_data_id} = shift;
	}
	return $self->{_analysis_data_id};

}

sub analysis_id {

	my $self = shift;
	if (@_) {
		$self->{_analysis_id} = shift;
	}
	return $self->{_analysis_id};

}

sub anchor_sequences_mlssid {

	my $self = shift;
	if (@_){
		$self->{_anchor_sequences_mlssid} = shift;
	}
	$self->{_anchor_sequences_mlssid};

}

sub ancs_from_to {

	my $self = shift;
	if (@_) {
		$self->{_ancs_from_to} = shift;
	}
	return $self->{_ancs_from_to};

}

sub configure_defaults {

 	my $self = shift;
	$self->analysis->program("/usr/local/ensembl/bin/exonerate-1.0.0") unless $self->analysis->program;
	$self->get_parameters( "{ bestn=>11, gappedextension=>\"no\", softmasktarget=>\"no\", percent=>75, showalignment=>\"no\", model=>\"affine:local\", }")
		unless $self->exonerate_options;
  	return 1;

}

sub exo_file {

	my $self = shift;
	if (@_) {
		$self->{_exo_file} = shift;
	}
	return $self->{_exo_file};

}

sub exonerate_mlssid {

	my $self = shift;
	if (@_){
		$self->{_exonerate_mlssid} = shift;
	}
	$self->{_exonerate_mlssid};

}

sub exonerate_options {

	my $self = shift;
	if (@_){
		$self->{_exonerate_options} = shift;
	}
	$self->{_exonerate_options};

}

sub fetch_input {

	my ($self) = @_;
	$self->configure_defaults();
	$self->get_parameters($self->parameters);
	#create a Compara::DBAdaptor which shares the same DBI handle with $self->db (Hive DBAdaptor)
	$self->{'comparaDBA'} = Bio::EnsEMBL::Compara::Production::DBSQL::DBAdaptor->new(-DBCONN=>$self->db->dbc);
	$self->{'comparaDBA'}->dbc->disconnect_if_idle();
	$self->{'hiveDBA'} = Bio::EnsEMBL::Hive::DBSQL::DBAdaptor->new(-DBCONN => $self->{'comparaDBA'}->dbc);
	$self->{'hiveDBA'}->dbc->disconnect_if_idle();
	$self->get_input_id($self->input_id);
	my $anchor_seq_adaptor = $self->{comparaDBA}->get_AnchorSeqAdaptor();
	my $analysis_data_adaptor = $self->{hiveDBA}->get_AnalysisDataAdaptor();
	my $target_genome_files = eval $analysis_data_adaptor->fetch_by_dbID($self->analysis_data_id);
	$self->target_file( $target_genome_files->{target_genomes}->{ $self->target_genome_db_id } );
	my $anchors = $anchor_seq_adaptor->get_anchor_sequences($self->ancs_from_to, $self->anchor_sequences_mlssid);
	my $query_file = $self->worker_temp_directory  . "anchors." . join ("-", @{$self->ancs_from_to});
	open F, ">$query_file" || throw("Couldn't open $query_file");
	foreach my $anchor_seq( @{ $anchors } ) {
		print F ">", join(":", @{$anchor_seq}[0..5]), "\n", $anchor_seq->[-1], "\n";
	}
	$self->query_file($query_file);
	return 1;

}

sub get_input_id {

	my $self = shift;
	my $input_id_string = shift;

	return unless($input_id_string);
	print("parsing input_id string : ",$input_id_string,"\n");
	
	my $params = eval($input_id_string);
	return unless($params);
	
	if(defined($params->{'ancs_from_to'})) {
		$self->ancs_from_to($params->{'ancs_from_to'});
	}
	if(defined($params->{'target_genome'})) {
		$self->target_genome_db_id($params->{'target_genome'});
	}
	if(defined($params->{'analysis_data_id'})) {
		$self->analysis_data_id($params->{'analysis_data_id'});
	}
	if(defined($params->{'analysis_id'})) {
		$self->analysis_id($params->{'analysis_id'});
	}
	if(defined($params->{'exonerate_mlssid'})) {
		$self->exonerate_mlssid($params->{'exonerate_mlssid'});
	}
	if(defined($params->{'anchor_sequences_mlssid'})) {
		$self->anchor_sequences_mlssid($params->{'anchor_sequences_mlssid'});
	}
	return 1;
}

1;

}

sub get_parameters {

	my $self = shift;
	my $param_string = shift;
	
	return unless($param_string);
	my $params = eval($param_string);
	$self->exonerate_options($params);

}

sub merge_overlapping_target_regions {

 #merge overlapping target regions hit by different seqs in the same anchor	my $self = shift;
	my $mapped_anchors = shift;
	my $HIT_NUMS;
	foreach my $anchor(sort {$a <=> $b} keys %{$mapped_anchors}) {
	        foreach my $targ_info(sort keys %{$mapped_anchors->{$anchor}}) {
	                @{$mapped_anchors->{$anchor}{$targ_info}} = sort {$a->[0] <=> $b->[0]} @{$mapped_anchors->{$anchor}{$targ_info}};
	                for(my$i=0;$i<@{$mapped_anchors->{$anchor}{$targ_info}};$i++) {
	                        my $anc_look_up_name = join(":", $anchor, $targ_info, $mapped_anchors->{$anchor}{$targ_info}->[$i]->[0]);
				if($i < @{$mapped_anchors->{$anchor}{$targ_info}} - 1) {
		                        if($mapped_anchors->{$anchor}{$targ_info}->[$i]->[1] >= $mapped_anchors->{$anchor}{$targ_info}->[$i+1]->[0]) {  
		                                unless($mapped_anchors->{$anchor}{$targ_info}->[$i]->[2] eq 
							$mapped_anchors->{$anchor}{$targ_info}->[$i+1]->[2]) {       
		                                        print STDERR "possible palindromic sequences: $anchor ", 
								"$mapped_anchors->{$anchor}{$targ_info}->[$i]->[2] ", 
								$mapped_anchors->{$anchor}{$targ_info}->[$i+1]->[2], "\n";
		                                        $mapped_anchors->{$anchor}{$targ_info}->[$i]->[2] = 0;
		                                }       
		                                if($mapped_anchors->{$anchor}{$targ_info}->[$i]->[1] < 
							$mapped_anchors->{$anchor}{$targ_info}->[$i+1]->[1]) {
		                                        $mapped_anchors->{$anchor}{$targ_info}->[$i]->[1] = 
								$mapped_anchors->{$anchor}{$targ_info}->[$i+1]->[1];
		                                }       
		                                $mapped_anchors->{$anchor}{$targ_info}->[$i]->[3] += $mapped_anchors->{$anchor}{$targ_info}->[$i+1]->[3];
		                                $mapped_anchors->{$anchor}{$targ_info}->[$i]->[3] /= 2; # simplistic scoring
						#count the organisms from which the anchor seqs were derived 
		                                $HIT_NUMS->{$anc_look_up_name}{anc_orgs}{$mapped_anchors->{$anchor}{$targ_info}->[$i+1]->[4]}++;
						#count number of anchor seqs that map
						$HIT_NUMS->{$anc_look_up_name}{seq_nums}++;
		                                splice(@{$mapped_anchors->{$anchor}{$targ_info}}, $i+1, 1);
		                                $i--;   
						next;
		                        }       
				}
				$HIT_NUMS->{$anc_look_up_name}{anc_orgs}{$mapped_anchors->{$anchor}{$targ_info}->[$i]->[4]}++;
				$HIT_NUMS->{$anc_look_up_name}{seq_nums}++;
	                }       
	        }       
	}
	return $HIT_NUMS;

}

sub process_exonerate_hits {

	my $self = shift;
	my($hits, $target2dnafrag, $hit_numbers) = @_;
	my($records_to_load);
	foreach my $anchor_id (sort keys %{$hits}) {
		foreach my $targ_dnafrag_info (sort keys %{$hits->{$anchor_id}}) {
			my $dnafrag_id = $target2dnafrag->{$targ_dnafrag_info};
			foreach my $hit_position (@{$hits->{$anchor_id}->{$targ_dnafrag_info}}) {
				my $index = join(":", $anchor_id, $targ_dnafrag_info, $hit_position->[0]);
				my $number_of_org_hits = keys %{$hit_numbers->{$index}->{anc_orgs}};
				my $number_of_seq_hits = $hit_numbers->{$index}->{seq_nums};
				push(@{$records_to_load}, join(":", $self->exonerate_mlssid, $anchor_id, $dnafrag_id, 
							@{$hit_position}[0..3], $number_of_org_hits, $number_of_seq_hits));
			}
		}
	}
	return $records_to_load;

}

sub query_file {

	my $self = shift;
	if (@_) {
		$self->{_query_file} = shift;
	}
	return $self->{_query_file};

}

sub run {

	my ($self) = @_;
	my $program = $self->analysis->program;
	my $query_file = $self->query_file;
	my $target_file = $self->target_file;
	my $option_st;
	foreach my $option(sort keys %{ $self->exonerate_options }) {
		$option_st .= " --" . $option . " " . $self->exonerate_options->{$option}; 
	}
	my $command = join(" ", $program, $option_st, $query_file, $target_file); 
	print $command, "\n";
	my $exo_fh;
	open( $exo_fh, "$command |" ) or throw("Error opening exonerate command: $? $!"); #run exonerate	
	$self->exo_file($exo_fh);
	return 1;

}

sub target_file {

	my $self = shift;
	if (@_){
		$self->{_target_file} = shift;
	}
	return $self->{_target_file};

}

sub target_genome_db_id {

	my $self = shift;
	if (@_){
		$self->{_target_genome_db_id} = shift;
	}
	return $self->{_target_genome_db_id};

}

sub write_output {

	my ($self) = @_;
	my $anchor_align_adaptor = $self->{'comparaDBA'}->get_AnchorAlignAdaptor();
	my $exo_fh = $self->exo_file;
	my ($hits, $target2dnafrag);
	while(my $mapping = <$exo_fh>){ 
		next unless $mapping =~/^vulgar:/;
		my($anchor_info, $targ_strand, $targ_info, $targ_from, $targ_to, $score) = (split(" ",$mapping))[1,8,5,6,7,9];
		($targ_from, $targ_to) = ($targ_to, $targ_from) if ($targ_from > $targ_to); #exonerate can switch these around
		$targ_strand = $targ_strand eq "+" ? "1" : "-1";
		$targ_from++; #modify the exonerate start position
		my($anchor_name, $anc_org) = split(":", $anchor_info);
		push(@{$hits->{$anchor_name}{$targ_info}}, [ $targ_from, $targ_to, $targ_strand, $score, $anc_org ]);
		$target2dnafrag->{$targ_info}++;
	}
	foreach my $target_info (sort keys %{$target2dnafrag}) {
		my($coord_sys, $dnafrag_name) = (split(":", $target_info))[0,2];
		$target2dnafrag->{$target_info} = $anchor_align_adaptor->fetch_dnafrag_id(
							$coord_sys, $dnafrag_name, $self->target_genome_db_id);
	}
	my $hit_numbers = $self->merge_overlapping_target_regions($hits);
	my $records = $self->process_exonerate_hits($hits, $target2dnafrag, $hit_numbers);	
	$anchor_align_adaptor->store_exonerate_hits($records);
	return 1;

}

General documentation

This modules is part of the Ensembl project
Email compara@ebi.ac.uk

This modules is part of the EnsEMBL project ()
Questions can be posted to the ensembl-dev mailing list: ensembl-dev@ebi.ac.uk

The rest of the documentation details each of the object methods.
Internal methods are usually preceded with a _