Raw content of BioMart::Dataset::GenomicAlign
<![CDATA[# $Id: GenomicAlign.pm,v 1.2 2006/11/25 18:11:31 arek Exp $
#
# BioMart module for BioMart::Dataset::GenomicSequence
# You may distribute this module under the same terms as perl itself
# POD documentation - main docs before the code
=head1 NAME
BioMart::Dataset::GenomicAlign
=head1 SYNOPSIS
A hidden Dataset containing sequence attributes that can be imported to other
visible Datasets which are compatible with its required data input, based
on the presence of one or more importable-exportable relationships.
=head1 DESCRIPTION
Dataset providing Align Sequence attributes, which can be imported into
other Datasets. AlignSequence is itself not a visible Dataset.
=head1 AUTHOR - Arek Kasprzyk, Darin London
=head1 CONTACT
This module is part of the BioMart project http://www.biomart.org
Questions can be posted to the mart-dev mailing list:
mart-dev@ebi.ac.uk
=head1 METHODS
=head1 Developer Notes
The peptide translation algorithm is taken directly
from the CodonTable module that is part of the
BioPerl project. For more information about the
BioPerl project, visit:
http://www.bioperl.org
=cut
package BioMart::Dataset::GenomicAlign;
# implements Dataset interface
use strict;
use warnings;
use base qw(BioMart::DatasetI);
use BioMart::Dataset::GenomicSequence::DNAAdaptor;
use BioMart::Configuration::ConfigurationTree;
use BioMart::Configuration::AttributeTree;
use BioMart::Configuration::FilterTree;
use BioMart::Configuration::AttributeGroup;
use BioMart::Configuration::FilterGroup;
use BioMart::Configuration::AttributeCollection;
use BioMart::Configuration::FilterCollection;
use BioMart::Configuration::Attribute;
use BioMart::Configuration::AttributeList;
use BioMart::Configuration::BooleanFilter;
use BioMart::Configuration::ValueFilter;
use BioMart::Configuration::FilterList;
#change this to change the size of each batch
use constant BATCHSIZE => 100;
my $i;
sub _new {
my ($self, @param) = @_;
$self->SUPER::_new(@param);
my $i=1; ##
$self->attr('dna', undef);
$self->attr('dnaparams', undef);
######### modified
$self->attr('recipe', 'raw_sequence'); ## this will hold a subRef
$self->attr('ignore', undef);
$self->attr('ignore_row', undef);
$self->attr('seq_edits', undef); #will hold seq_edits to be applied after sequence has been produced
#$self->attr('codon_table_id', 1); #codon table defaults to 1
$self->attr('seq_name', undef); #this is linked to the Attribute->name, determines which sequence recipe to run
$self->attr('translate', 0); # set to true for peptide
#$self->attr('downstream_flank', 0);
#$self->attr('upstream_flank', 0);
$self->attr('importable', undef);
$self->attr('lastPkey', undef);
$self->attr('importable_names', undef); # initialized when first row processed in first batch for a given query
$self->attr('importable_indices', undef); # initialized when first row processed in first batch for a given query
$self->attr('returnRow_indices', undef); # initialized when first row processed in first batch for a given query
$self->attr('returnRow', undef);
$self->attr('batchIndex', 0); #increment each time a new pkey is encountered.
$self->attr('seq_species', undef);
#process rows until this equals batchSize
$self->attr('locations', {}); # not used by all sequences
$self->attr('outRow', undef); # not used by all sequences
#attributes calculated over sequence locations
$self->attr('calc_location', undef);
$self->attr('sequence', undef);
}
#private methods
sub _ignoreRow {
my ($self, $curRow) = @_;
my $ignore = $self->get('ignore');
return 0 unless ($ignore);
my $ignore_row = $self->get('ignore_row');
my $test = $self->_getLocationFrom($curRow, $ignore_row);
#if the actual value is false, return false, else, return ignore for the value
return $test->{ $ignore_row } && $ignore->{ $test->{ $ignore_row } };
}
sub _editSequence {
my ($self, $seqref) = @_;
my $seq_edits = $self->get('seq_edits');
if ($$seqref && $seq_edits) {
foreach my $seq_edit (split /\;/, $seq_edits) {
my ($start, $end, $alt_seq) = split /\,/, $seq_edit;
my $len = $end - $start + 1;
substr($$seqref, $start - 1, $len) = $alt_seq;
}
}
}
sub _initializeDNAAdaptor {
my $self = shift;
my $dna_params = $self->getConfigurationTree()->optionalParameters;
unless ($dna_params) {
BioMart::Exception::Configuration->throw("GenomicSequence Dataset requires optional_parameters to be set in the DatasetConfig\n");
}
my $dna = {};
## from the has_mmu_align_sequence dataset
## need to parse the attributes - to get different DNAadaptor for each species
## Attribute : optional parameters
## hsapiens_sequence,hsapiens_genomic_sequence__dna_chunks__main,chr_name,chr_start,sequence,100000;mmus_oriented_raw_sequence,mmusculus_genomic_sequence__dna_chunks__main,chr_name,chr_start,sequence,100000
##
foreach my $dna_params4specie ( split /\;/, $dna_params ){
#print STDERR ("##GenomicAlign _initializeDNAAdaptor ".localtime(time)."\n");
# print "\n+In subroutine __dnaAdaptor dna_params4specie $dna_params4specie\n";
## hsa_oriented_raw_sequence,hsapiens_genomic_sequence__dna_chunks__main,chr_name,chr_start,sequence,100000
my ($attribute_name, $dnatablename, $chunk_name_fieldname, $chunk_start_fieldname, $seqfieldname,$chunk_size) = split /\,/, $dna_params4specie ;
# print "\n\n\nstoring dna for $attribute_name $dnatablename $chunk_name_fieldname $chunk_start_fieldname $seqfieldname $chunk_size\n\n\n";
# print "chunk_size: $chunk_size\n";
$dna->{$attribute_name} = BioMart::Dataset::GenomicSequence::DNAAdaptor->new('seq_name' => $attribute_name,
'dna_tablename' => $dnatablename, ##hsapiens_genomic_sequence__dna_chunks__main
'seq_fieldname' => $seqfieldname, ##sequence
'chunk_name_fieldname' => $chunk_name_fieldname, ##chr_name
'chunk_start_fieldname' => $chunk_start_fieldname, ##chr_start
'chunk_size' => $chunk_size, ##100000
'configurator' => $self->getParam('configurator')
);
unless ($dna->{$attribute_name}) {
BioMart::Exception::Configuration->throw("Couldnt connect to DNAAdaptor for $attribute_name\n\n");
}
}
# print "dna: ",$dna,"\n";
$self->set('dna', $dna);
# print "\n+Out subroutine: __dnaAdaptor GenomicAlign\n";
}
sub __processNewQuery {
my ($self, $query) = @_;
#print STDERR ("##GenomicAlign 1_processNewQuery ".localtime(time)."\n");
# print "\nname:",$self->name," GenomicAlign\n";
# print "\nquery: ",$query," GenomicAlign\n";
# print "\ntableau: ",$query->getAllAttributes($self->name)," GenomicAlign\n";
my $attribute = $query->getAllAttributes($self->name)->[0];
my $seq_name = $attribute->name;
$self->set('seq_name', $seq_name);
##$self->set('translate', ($seq_name =~ m/peptide$/));
#print STDERR ("##GenomicAlign 2_processNewQuery ".localtime(time)."\n");
my $ignore = 'IGNORE';
if ($seq_name =~ m/raw/){
$self->set('recipe', '_rawSequences');
} else {
BioMart::Exception::Configuration->throw("Unsupported sequence name $seq_name recieved by GenomicAlign\n");
}
#print STDERR ("##GenomicAlign 3_processNewQuery ".localtime(time)."\n");
#$self->set('downstream_flank', 0);
#$self->set('upstream_flank', 0);
$self->set('importable', undef);
$self->set('lastPkey', undef);
$self->set('importable_indices', undef);
$self->set('returnRow_indices', undef);
$self->set('locations', {});
$self->set('outRow', undef);
$self->set('calc_location', undef);
$self->set('sequence', undef);
#determine which BaseSequenceA object to create
my $filters = $query->getAllFilters($self->name);
#print STDERR ("##GenomicAlign 4_processNewQuery ".localtime(time)."\n");
foreach my $filt (@{$filters}) {
if ($filt->isa("BioMart::Configuration::FilterList")) {
if ($filt->linkName) {
if ($self->get('importable') ) {
BioMart::Exception::Configuration->throw("Recieved two importables, can only work with one\n");
} else {
$self->set('importable', $filt);
}
} else {
BioMart::Exception::Configuration->throw("Recieved invalid linkName ".$filt->linkName."\n");
}
} else {
#must be a downstream or upstream valueFilter
unless ($filt->isa("BioMart::Configuration::ValueFilter")) {
BioMart::Exception::Configuration->throw("Recieved unknown filter ".$filt->name." in GenomicSequence Dataset!\n");
}
if ($self->get($filt->name)) {
BioMart::Exception::Configuration->throw("Recieved two ".$filt->name." flanking filters in GenomicSequence Dataset\n");
}
#could still be some strange ValueFilter that is not upstream or downstream, but not likely
#will throw an exception if this is the case
my $table = $filt->getTable;
my $row = $table->nextRow;
my $value = $row->[0];
if ($value) {
$self->set($filt->name, $value);
}
}
}
#print STDERR ("##GenomicAlign 5_processNewQuery ".localtime(time)."\n");
unless ($self->get('importable')) {
BioMart::Exception::Configuration->throw("No Importable Recieved in GenomicAlign\n");
}
}
sub _continueWithBatch {
my ($self, $batchSize, $rtable) = @_;
#always true if underlying table is an AttributeTable and it has rows
my $continue = ($rtable->isa("BioMart::ResultTable"))
? $rtable->inCurrentBatch()
: $rtable->hasMoreRows;
if ($continue && $batchSize) {
my $batchIndex = $self->get('batchIndex');
$continue = ($batchIndex < $batchSize);
}
return $continue;
}
sub _incrementBatch {
my $self = shift;
my $batchIndex = $self->get('batchIndex');
$batchIndex++;
$self->set('batchIndex', $batchIndex);
}
sub _initializeIndices {
my ($self, $numFields) = @_;
my $importable_names = [];
my $returnRow_indices = {};
my $importable_indices = {};
my $filts = $self->get('importable')->getAllFilters;
#define where the importable fields are in rtable
my $index = 0;
foreach my $filt (@{$filts}) {
push @{$importable_names}, $filt->name;
$importable_indices->{$filt->name} = $index;
$index++;
}
#define where fields needing to be merged into final returnRow are in rtable
my $resultIndex = 0;
while ($index < $numFields) {
$returnRow_indices->{$index} = $resultIndex;
$index++;
$resultIndex++;
}
$self->set('importable_indices', $importable_indices);
$self->set('returnRow_indices', $returnRow_indices);
$self->set('importable_names', $importable_names);
}
sub _initializeReturnRow {
my ($self, $curRow) = @_;
#This does __NOT__ concatenate fields that are many<->one with the pkey
my $returnRow = [];
foreach my $val (@{$curRow}) {
push @{$returnRow}, $val;
}
return $returnRow;
}
sub _processRow {
my ($self, $atable, $curRow) = @_;
# if this is the very first row for a new query, initialize the indices using its length as numFields
unless ($self->get('importable_indices')) {
if ($self->get('exhausted')) {
$atable->addRow(["No Sequence Returned"]);
} else {
my $numFields = @{$curRow};
$self->_initializeIndices($numFields);
}
}
my $method = $self->get('recipe');
$self->$method($atable, $curRow);
}
sub _calcSeqOverLocations {
my ($self, $this_location) = @_;
$this_location->{start} || return; # Sanity check
$this_location->{end} || return; # Sanitt check
my $calc_location = $self->get('calc_location');
if ($calc_location) {
$calc_location->{"start"} = $this_location->{"start"} if ($this_location->{"start"} < $calc_location->{"start"});
$calc_location->{"end"} = $this_location->{"end"} if ($this_location->{"end"} > $calc_location->{"end"});
} else {
$calc_location = {};
foreach my $key (keys %{$this_location}) {
$calc_location->{$key} = $this_location->{$key};
}
}
$self->set('calc_location', $calc_location);
}
sub _getLocationFrom {
my ($self, $curRow, @expectedFields) = @_;
my $importable_indices = $self->get('importable_indices');
my $location = {};
foreach my $expectedField (@expectedFields) {
$location->{$expectedField} = ( exists( $importable_indices->{$expectedField} ) ) ? $curRow->[ $importable_indices->{$expectedField} ] : undef;
}
return $location;
}
######################################################
sub _processSequence { ############################
my ($self, $location, $attribute_name, $count) = @_;
#@species_attribute_name contains hsa et mmu
#print STDERR ("##GenomicAlign _processSequence starting ".localtime(time)."\n");
#print "attribute_name $attribute_name\n";
#print "count $count\n";
# my $i=1;
my $seq = '';
my $dna = $self->get('dna')->{$attribute_name};
my $chr = $location->{'chr'.$count};
my $start = $location->{'start'.$count};
my $end = $location->{'end'.$count};
my $strand = $location->{'strand'.$count};
#------------------------
#my $seq = '';
#my $dna = $self->get('dna')->{$attribute_name};
#my $chr = $location->{'chr1'};
#$chr = $location->{'chr2'} unless (defined $chr);
#my $start = $location->{'start1'};
#$start = $location->{'start2'} unless (defined $start);
#my $end = $location->{'end1'};
#$end = $location->{'end2'} unless (defined $end);
#my $strand = $location->{'strand1'};
#$strand = $location->{'strand2'} unless (defined $strand);
#$strand = 1 unless (defined $strand);
#print "coucou $chr $start $end $strand\n";
if ($strand < 0) {
#print STDERR ("##GenomicAlign NO rc \n".$dna->getSequence( $chr, $start, $end )."\n");
#print STDERR ("#GenomicAlign start rc \n"._rc( $dna->getSequence( $chr, $start, $end ))."\n");
$seq .= $self->_rc( $dna->getSequence( $chr, $start, $end ) );
#print STDERR ("##GenomicAlign start end ".localtime(time)."\n");
#print STDERR ("##GenomicAlign RC \n".$seq."\n");
} else {
#print STDERR ("##GenomicAlign start getSequence ".localtime(time)."\n");
$seq .= $dna->getSequence( $chr, $start, $end );
#print STDERR ("##GenomicAlign start getSequence ".localtime(time)."\n");
}
$i++;
#print STDERR ("##GenomicAlign 3_processSequence $attribute_name ".localtime(time)."\n");
if (length($seq)) {
return $seq;
}
# print STDERR ("##GenomicAlign 4_processSequence $attribute_name ".localtime(time)."\n");
return undef;
}
sub _addRow {
my ($self, $atable, $outRow) = @_;
$atable->addRow($outRow);
$self->_incrementBatch;
}
sub _rc{
my ($self, $seq) = @_;
#print STDERR "GenomicAlign reverse start ".localtime(time)."\n";
$seq = reverse($seq);
#print STDERR "GenomicAlign reverse end tr start ".localtime(time)."\n";
$seq =~ tr/YABCDGHKMRSTUVyabcdghkmrstuv/RTVGHCDMKYSAABrtvghcdmkysaab/;
#print STDERR "GenomicAlign tr end _rc end ".localtime(time)."\n";
return $seq;
}
#interface methods
sub _getConfigurationTree {
my $self = shift;
return $self->getParam('configurator')->getConfigurationTree($self->virtualSchema, $self->name);
}
sub _getExportables {
my ($self, $linkName) = @_;
my $exportables = $self->get('exportables');
if ($linkName) {
return [ $exportables->{$linkName} ];
}
my $ref = [];
push @{$ref}, values %{$exportables};
return $ref;
}
sub _getImportables {
my ($self, $linkName) = @_;
my $importables = $self->get('importables');
if ($linkName) {
#return [ $importables->{$linkName} ];
return $importables->{$linkName};
}
my $ref = [];
push @{$ref}, values %{$importables};
return $ref;
}
sub _getResultTable {
my ($self, @param) = @_;
$self->set('batchIndex', 0);
local($^W) = 0; # prevent "odd number of elements" warning with -w.
my(%param) = @param;
my $query = $param{'query'};
my $atable = $param{'table'};
my $batch_size = $param{'batch_size'};
if ($self->serverType eq "web"){
my $batch_start = $param{'batch_start'} || 0;
my $location = $self->getParam('configurator')->get('location');
my $xml = $query->toXML($batch_start,$batch_size,0);
foreach my $el($location->getResultSet("","POST",$xml)){
if ($el =~ /No Sequence Returned/) {
$self->_setExhausted(1);
last;
}
my @clean=split(/\t/,$el);
$atable->addRow([@clean]);
}
return $atable;
} else {
$self->_initializeDNAAdaptor($query->
getInterfaceForDataset($self->name));
}
# print STDERR ("##BATCH SIZE FOR ".$self->name." IS ".$batch_size."\n");
my $importable = $self->get('importable');
my $rtable = $importable->getTable();
my $has_rows = $rtable->hasMoreRows;
while ($has_rows && $self->_continueWithBatch($batch_size, $rtable)) {
$self->_processRow( $atable, $rtable->nextRow );
}
# the last and final call to GenomicSequence, after the call which exhausts the importable,
# will result in the last sequence being processed and added to the resultTable.
# the next call after this returns undef.
unless ($has_rows) {
$self->_setExhausted(1);
# $self->_processRow($atable);
}
$importable->setTable($rtable);
$self->set('importable', $importable);
my $dna = $self->get('dna');
foreach my $attribute_name (keys %$dna) {
$dna->{$attribute_name}->close;
}
return $atable;
}
### sequence __recipes__
sub _rawSequencesOriginal { ## SHOULD BE REMOVED
my ($self, $atable, $curRow) = @_; ## AS IT WAS FROM GENOMICSEQUENCE.PM
my $rank = 1;
if ($curRow) {
my $importable_indices = $self->get('importable_indices');
my $locations = {};
my $location = $self->_getLocationFrom($curRow, "chr", "start", "end");
$location->{"strand"} = ( exists( $importable_indices->{"strand"} ) ) ? $curRow->[ $importable_indices->{"strand"} ] : 1;
$locations->{$rank} = $location if ($location->{"start"});
my $sequence = $self->_processSequence($locations);
$self->_editSequence(\$sequence);
if ($sequence) {
$self->_addRow($atable, $self->_initializeReturnRow($curRow), $sequence);
}
}
#else there is no last entry
}
sub _rawSequences {
my ($self, $atable, $curRow) = @_;
my $rank = 1;
my $overall_count = 0;
my $local_count = 0;
my $species_numbers = 0;
my $count = 1;
my $n = 0;
my @importable_names = @{$self->get('importable_names')};
my $dna_params = $self->getConfigurationTree()->optionalParameters;
my @species_dna_params = split(/\;/, $dna_params);
my @species_attribute_name;
# print STDERR ("##GenomicAlign 1_rawSequences ".localtime(time)."\n");
foreach my $sdp (@species_dna_params) {
my ($attribute_name) = split(/\,/,$sdp);
push @species_attribute_name, $attribute_name;
#print "## $attribute_name\n";
}
# print STDERR ("##GenomicAlign 2_rawSequences ".localtime(time)."\n");
my $initRow = $self->_initializeReturnRow($curRow);
# print STDERR ("##GenomicAlign 3_rawSequences ".localtime(time)."\n");
while (my $attribute_name = shift @species_attribute_name){
my $importable_indices = $self->get('importable_indices');
# print STDERR ("##GenomicAlign 4_rawSequences ".localtime(time)."\n");
my ($name, $start, $end, $strand);
foreach my $var (\$name, \$start, \$end, \$strand) {
$$var = shift @importable_names;
last if (defined $strand);
#print "strand $strand";
#shift @importable_names;
}
# print STDERR ("##GenomicAlign 5_rawSequences $attribute_name".localtime(time)."\n");
my $location = $self->_getLocationFrom($curRow, ($name, $start, $end, $strand));
# print STDERR ("##GenomicAlign 6_rawSequences $attribute_name ".localtime(time)."\n");
#print "$attribute_name $name, $start, $end, $strand\n";
my $sequence = $self->_processSequence($location, $attribute_name, $count);
# print STDERR ("##GenomicAlign 7_rawSequences $attribute_name ".localtime(time)."\n");
if ($sequence) {
push @{$initRow}, $sequence;
}
## IMPORTANT remove the length as the coordinate aer as folow
## name, start, end, strand, length
shift @importable_names ;
# print STDERR ("##GenomicAlign 8_rawSequences $attribute_name ".localtime(time)."\n");
$count++;
}
$self->_addRow($atable, $initRow);
# print STDERR ("##GenomicAlign 9_rawSequences ".localtime(time)."\n");
}
1;
]]>