Raw content of BioMart::Dataset::GenomicSequence #$Id: GenomicSequence.pm,v 1.21.2.1 2008/08/13 12:59:45 syed 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::GenomicSequence =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 Genomic Sequence attributes, which can be imported into other Datasets. GenomicSequence is itself not a visible Dataset. =head1 AUTHOR - Arek Kasprzyk, Syed Haider, Darin London, Damian Smedley =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::GenomicSequence; # implements Dataset interface use strict; use warnings; use base qw(BioMart::DatasetI); use Log::Log4perl; use Data::Dumper; 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; #change this if want a different default Codon Table ID use constant DEFAULTCODONTABLEID => 1; #exon will build their ignore keys if necessary use constant IGNORE => { q(gene_exon) => {} }; use vars qw(@NAMES @TABLES $CODONS $TRCOL %IUPAC_DNA); BEGIN { @NAMES = #id ( 'Standard', #1 'Vertebrate Mitochondrial',#2 'Yeast Mitochondrial',# 3 'Mold, Protozoan, and CoelenterateMitochondrial and Mycoplasma/Spiroplasma',#4 'Invertebrate Mitochondrial',#5 'Ciliate, Dasycladacean and Hexamita Nuclear',# 6 '', '', 'Echinoderm Mitochondrial',#9 'Euplotid Nuclear',#10 '"Bacterial"',# 11 'Alternative Yeast Nuclear',# 12 'Ascidian Mitochondrial',# 13 'Flatworm Mitochondrial',# 14 'Blepharisma Nuclear',# 15 'Chlorophycean Mitochondrial',# 16 '', '', '', '', 'Trematode Mitochondrial',# 21 'Scenedesmus obliquus Mitochondrial', #22 'Thraustochytrium Mitochondrial' #23 ); @TABLES = qw( FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNKKSS**VVVVAAAADDEEGGGG FFLLSSSSYY**CCWWTTTTPPPPHHQQRRRRIIMMTTTTNNKKSSRRVVVVAAAADDEEGGGG FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNKKSSSSVVVVAAAADDEEGGGG FFLLSSSSYYQQCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG '' '' FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNNKSSSSVVVVAAAADDEEGGGG FFLLSSSSYY**CCCWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG FFLLSSSSYY**CC*WLLLSPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNKKSSGGVVVVAAAADDEEGGGG FFLLSSSSYYY*CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNNKSSSSVVVVAAAADDEEGGGG FFLLSSSSYY*QCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG FFLLSSSSYY*LCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG '' '' '' '' FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNNKSSSSVVVVAAAADDEEGGGG FFLLSS*SYY*LCC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG FF*LSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG ); my @nucs = qw(t c a g); my $x = 0; ($CODONS, $TRCOL) = ({}, {}); for my $i (@nucs) { for my $j (@nucs) { for my $k (@nucs) { my $codon = "$i$j$k"; $CODONS->{$codon} = $x; $TRCOL->{$x} = $codon; $x++; } } } %IUPAC_DNA = ( A => [qw(A)], C => [qw(C)], G => [qw(G)], T => [qw(T)], U => [qw(U)], M => [qw(A C)], R => [qw(A G)], W => [qw(A T)], S => [qw(C G)], Y => [qw(C T)], K => [qw(G T)], V => [qw(A C G)], H => [qw(A C T)], D => [qw(A G T)], B => [qw(C G T)], X => [qw(G A T C)], N => [qw(G A T C)] ); } sub _new { my ($self, @param) = @_; $self->SUPER::_new(@param); $self->attr('dna', undef); $self->attr('dnaparams', undef); $self->attr('recipe', undef); #this will hold a subRef $self->attr('ignore', undef); $self->attr('ignore_row', undef); $self->attr('seq_edits', undef); $self->attr('codon_table_id', undef); #codon table defaults to undef $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_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 seen $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); $self->attr('attribute_merge_required', 0); # last set of structure rows in a batch are not processed as they may be incomplete # set of say all exons, so store them somewhere for next time around and also store the # corresponding results for merging $self->attr('rowsFromLastBatch', undef); # these help us solve the problem when gene coding and flank sequnces are requested, and # for more than one transcripts coming in different batches, it repeats the same sequence # with different header if transcript atts are requested. # $self->attr('rowsFromLastBatch', undef), works for this case too. $self->attr('seqStorageHash', undef); $self->attr('onHoldSeqsKeys', undef); #-- new GS development --- $self->attr('exon_idHash', undef); #------------------------- } #private methods sub _rc{ my ($self, $seq) = @_; $seq = reverse($seq); $seq =~ tr/YABCDGHKMRSTUVyabcdghkmrstuv/RTVGHCDMKYSAABrtvghcdmkysaab/; return $seq; } 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 _translate { my ($self, $seq) = @_; BioMart::Exception::Configuration->throw("Calling translate without a seq argument!") unless defined $seq; return '' unless $seq; my $id = $self->get('codon_table_id') || DEFAULTCODONTABLEID; my ($partial) = 0; $partial = 2 if length($seq) % 3 == 2; $seq = lc $seq; $seq =~ tr/u/t/; my $protein = ""; if ($seq =~ /[^actg]/ ) { #ambiguous chars for (my $i = 0; $i < (length($seq) - 2 ); $i+=3) { my $triplet = substr($seq, $i, 3); if (exists $CODONS->{$triplet}) { $protein .= substr($TABLES[$id-1], $CODONS->{$triplet},1); } else { $protein .= $self->_translate_ambiguous_codon($triplet); } } } else { # simple, strict translation for (my $i = 0; $i < (length($seq) - 2 ); $i+=3) { my $triplet = substr($seq, $i, 3); if (exists $CODONS->{$triplet}) { $protein .= substr($TABLES[$id-1], $CODONS->{$triplet}, 1); } else { $protein .= 'X'; } } } if ($partial == 2) { # 2 overhanging nucleotides my $triplet = substr($seq, ($partial -4)). "n"; if (exists $CODONS->{$triplet}) { my $aa = substr($TABLES[$id-1], $CODONS->{$triplet},1); $protein .= $aa; } else { $protein .= $self->_translate_ambiguous_codon($triplet, $partial); } } return $protein; } sub _translate_ambiguous_codon { my ($self, $triplet, $partial) = @_; $partial ||= 0; my $id = $self->get('codon_table_id') || DEFAULTCODONTABLEID; my $aa; my @codons = _unambiquous_codons($triplet); my %aas =(); foreach my $codon (@codons) { $aas{substr($TABLES[$id-1],$CODONS->{$codon},1)} = 1; } my $count = scalar keys %aas; if ( $count == 1 ) { $aa = (keys %aas)[0]; } elsif ( $count == 2 ) { if ($aas{'D'} and $aas{'N'}) { $aa = 'B'; } elsif ($aas{'E'} and $aas{'Q'}) { $aa = 'Z'; } else { $partial ? ($aa = '') : ($aa = 'X'); } } else { $partial ? ($aa = '') : ($aa = 'X'); } return $aa; } sub _unambiquous_codons{ my ($value) = @_; my @nts = (); my @codons = (); my ($i, $j, $k); @nts = map { $IUPAC_DNA{uc $_} } split(//, $value); for my $i (@{$nts[0]}) { for my $j (@{$nts[1]}) { for my $k (@{$nts[2]}) { push @codons, lc "$i$j$k"; } } } return @codons; } 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; } } # important to clear if two seq_edit sequences come one after another $self->set('seq_edits', ""); } sub _initializeDNAAdaptor { my ($self,$interface) = @_; my $dna_params = $self->getConfigurationTree($interface)-> optionalParameters; unless ($dna_params) { BioMart::Exception::Configuration->throw("GenomicSequence Dataset requires optional_parameters to be set in the DatasetConfig\n"); } my ($dnatablename, $chunk_name_fieldname, $chunk_start_fieldname, $seqfieldname,$chunk_size) = split /\,/, $dna_params; my $dna = BioMart::Dataset::GenomicSequence::DNAAdaptor->new( 'seq_name' => $self->name, 'dna_tablename' => $dnatablename, 'seq_fieldname' => $seqfieldname, 'chunk_name_fieldname' => $chunk_name_fieldname, 'chunk_start_fieldname' => $chunk_start_fieldname, 'chunk_size' => $chunk_size, 'configurator' => $self->getParam('configurator'), ); unless ($dna) { BioMart::Exception::Configuration->throw("Couldnt connect to DNAAdaptor\n"); } $self->set('dna', $dna); } sub __processNewQuery { my ($self, $query) = @_; my $attribute = $query->getAllAttributes($self->name)->[0]; my $seq_name = $attribute->name; # hack to keep webservices working for 0_5 originating query XML if ($seq_name eq 'pkey'){ $attribute = $query->getAllAttributes($self->name)->[-1]; $seq_name = $attribute->name; } $self->set('seq_name', $seq_name); $self->set('translate', ($seq_name =~ m/peptide$/)); my $ignore = IGNORE; if ($seq_name =~ m/(coding|cdna|peptide)$/) { # this is not actually going to ignore anything, but is # simply used to determine translation table for each gene # without creating a second instance variable $self->set('ignore_row', "type"); $self->set('recipe', '_codingCdnaPeptideSequences'); } elsif ($seq_name =~ m/(transcript_exon_intron|transcript_flank|coding_transcript_flank)$/) { $self->set('recipe', '_transcript_exonIntronFlankSequences'); } elsif ($seq_name =~ m/(gene_exon_intron|gene_flank|coding_gene_flank)$/) { $self->set('recipe', '_gene_exonIntronFlankSequences'); } elsif ($seq_name =~ m/raw/){ $self->set('recipe', '_rawSequences'); } elsif ($seq_name =~ m/(gene_exon|transcript_exon|transcript_intron)$/) { # set the system to ignore rows with duplicate pkeys # for gene_exon $self->set('ignore', $ignore->{$1}); #undef for transcript_exon $self->set('ignore_row', "pkey"); $self->set('recipe', '_exonSequences'); } elsif ($seq_name =~ m/utr$/ or $seq_name =~ m/intergenic$/) { $self->set('recipe', '_utrSequences'); } elsif ($seq_name =~ m/snp$/) { $self->set('recipe', '_snpSequences'); } else { BioMart::Exception::Configuration->throw("Unsupported sequence name $seq_name recieved by GenomicSequence\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); $self->set('rowsFromLastBatch', undef); $self->set('seqStorageHash', undef); $self->set('onHoldSeqsKeys', undef); #-- new GS development --- $self->set('exon_idHash', undef); #------------------------- #determine which BaseSequenceA object to create my $filters = $query->getAllFilters($self->name); 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); } } } unless ($self->get('importable')) { BioMart::Exception::Configuration->throw("No Importable Recieved in GenomicSequence\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 $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}) { $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); } sub _initializeReturnRow { my ($self, $curRow) = @_; # this method used to handle the structure attributes as FASTA headers # no longer necessary with new attribute merging code in DatasetI.pm # if hashed attributes exist from previous structure dataset then just # return $curRow, otherwise return [] # my $importable = $self->get('importable'); # my $rtable = $importable->getTable(); return $self->get('attribute_merge_required') ? $curRow : []; } 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; } # (New-1) hack for MartBuilder built marts. different set of exportables importables # need to calculate start/end w.r.t old system # here we tamper start and end values of $location so they correspond to # the old system whereby they meant coding_start/end or 3/5utr_start/end # this is executed only for coding_gene_flank, coding_transcript_flank # 5utr, 3utr, coding and peptide sequence types. my $newCalculationsRequired = 0; foreach my $key (keys %$importable_indices) { $newCalculationsRequired = 1 if ($key eq 'coding_start_offset'); } if($newCalculationsRequired) { my $pkey = $curRow->[$importable_indices->{'pkey'}]; my $rank = $curRow->[$importable_indices->{'exon_id'}]; my $start_rank = $curRow->[$importable_indices->{'start_exon_id'}]; my $end_rank = $curRow->[$importable_indices->{'end_exon_id'}]; my $exon_start = $curRow->[$importable_indices->{'start'}]; my $exon_end = $curRow->[$importable_indices->{'end'}]; my $strand = $curRow->[$importable_indices->{'strand'}]; my $coding_start_offset = $curRow->[$importable_indices->{'coding_start_offset'}]; my $coding_end_offset = $curRow->[$importable_indices->{'coding_end_offset'}]; my ($new_start, $new_end); #foreach (%$location) { # print "<br/>KEY=$_ : VALUE=", $location->{$_}; #} #exit; my $exon_idHash = $self->get('exon_idHash'); if ($start_rank && $end_rank) { if ($self->get('seq_name') eq '5utr') { if ($rank == $start_rank) { if ($strand == 1) { $new_start = $exon_start; $new_end = $exon_start + $coding_start_offset - 2; } if ($strand == -1) { $new_start = $exon_end - $coding_start_offset + 2; $new_end = $exon_end; } # encountered coding START exon. storing as transcriptId.start_exonId to make it unique $self->_setCodingExonFlag($pkey.$start_rank, 1); } else { if ($exon_idHash->{$pkey.$start_rank} && $exon_idHash->{$pkey.$start_rank} == 1) { $new_start = undef; $new_end = undef; } else { $new_start = $exon_start; $new_end = $exon_end; } } } elsif ($self->get('seq_name') eq '3utr') { if ($rank == $end_rank) { if ($strand == 1) { $new_start = $exon_start + $coding_end_offset; $new_end = $exon_end; } if ($strand == -1) { $new_start = $exon_start; $new_end = $exon_end - $coding_end_offset; } # encountered coding END exon. storing as transcriptId.end_exonId to make it unique $self->_setCodingExonFlag($pkey.$end_rank, 1); } else { if ($exon_idHash->{$pkey.$end_rank} && $exon_idHash->{$pkey.$end_rank} == 1) { $new_start = $exon_start; $new_end = $exon_end; } else { $new_start = undef; $new_end = undef; } } } if ($self->get('seq_name') =~ m/^coding.*?|peptide/) { # this includes sequence types : coding_gene_flank, coding_transcript_flank, coding, peptide if ($rank == $start_rank && $rank == $end_rank) { # coding start and finish on the same exon if ($strand == 1) { $new_start = $exon_start + $coding_start_offset - 1; $new_end = $exon_start + $coding_end_offset - 1; } if ($strand == -1) { $new_start = $exon_end - $coding_end_offset + 1; $new_end = $exon_end - $coding_start_offset + 1; } } elsif ($rank == $start_rank) { if ($strand == 1) { $new_start = $exon_start + $coding_start_offset - 1; $new_end = $exon_end; } if ($strand == -1) { $new_start = $exon_start; $new_end = $exon_end - $coding_start_offset + 1; } # encountered coding START exon. storing as transcriptId.start_exonId to make it unique $self->_setCodingExonFlag($pkey.$start_rank, 1); } elsif ($rank == $end_rank) { if ($strand == 1) { $new_start = $exon_start; $new_end = $exon_start + $coding_end_offset - 1; } if ($strand == -1) { $new_start = $exon_end - $coding_end_offset + 1; $new_end = $exon_end; } # encountered coding END exon. clear ranscriptId.start_exonId stored earlier $self->_setCodingExonFlag($pkey.$start_rank, 0); } else { if ($exon_idHash->{$pkey.$start_rank} && $exon_idHash->{$pkey.$start_rank} == 1) { $new_start = $exon_start; $new_end = $exon_end; } else { $new_start = undef; $new_end = undef; } } } # sanity check if ($new_start && $new_end && ($new_start > $new_end)) { $new_start = undef; $new_end = undef; } # update the coordinates to be returned by this function - now similar to old coordinates $location->{'start'} = $new_start; $location->{'end'} = $new_end; } else { # no coding start/end set the vals to NULL $location->{'start'} = undef; $location->{'end'} = undef; } } return $location; } sub _setCodingExonFlag { my ($self, $key, $val) = @_; # the key value is pkey.start_exonId OR pkey.end_exonId # doesnt matter if its start or end exon id as the corresponding # value is always the same for each exon. my $exon_idHash = $self->get('exon_idHash'); $exon_idHash->{$key} = $val; $self->set('exon_idHash', $exon_idHash); } sub _modFlanks { my ($self, $location, $shift, $flank) = @_; $location->{start} || return $location; # Sanity check $location->{end} || return $location; # Sanity check if ($shift == 1) { # shift for flanks only - if user accidentally chooses both flanks, # assume upstream as the original martview if ($self->get('upstream_flank') && $self->get('downstream_flank')){ BioMart::Exception::Usage->throw("For this sequence option choose upstream OR downstream gene flanking sequence, NOT both, as makes no sense to simply concatenate them together.\n"); } if ($self->get('upstream_flank')) { if ($location->{"strand"} < 0) { $location->{"start"} = $location->{"end"} + 1; $location->{"end"} += $self->get('upstream_flank'); } else { $location->{"end"} = $location->{"start"} - 1; $location->{"start"} -= $self->get('upstream_flank'); } } elsif ($self->get('downstream_flank')) { if ($location->{"strand"} < 0) { $location->{"end"} = $location->{"start"} - 1; $location->{"start"} -= $self->get('downstream_flank'); } else { $location->{"start"} = $location->{"end"} + 1; $location->{"end"} += $self->get('downstream_flank'); } } else { BioMart::Exception::Usage->throw("Requests for flank sequence must be accompanied by an upstream_flank or downstream_flank request\n"); } } elsif ($shift == 2) #-- this is specific for: flanks + coding|cdna [ben] { if ($self->get('upstream_flank') && ($flank eq "up")) { #-- the flank variable is used here, in case the user ask for both upstream #-- and downstream flank. With $flank you are sure, you don't enter in #-- the upstream loop when you should get in the downstream loop. #-- or you don't enter twice the upstream loop when you ask for both up+down #-- or you get a single exon gene/transcript if ($location->{"strand"} < 0) { $location->{"end"} += $self->get('upstream_flank'); } else { $location->{"start"} -= $self->get('upstream_flank'); } } elsif ($self->get('downstream_flank') && ($flank eq "down")) { if ($location->{"strand"} < 0) { $location->{"start"} -= $self->get('downstream_flank'); } else { $location->{"end"} += $self->get('downstream_flank'); } } } else #-- if shift = 0 { if ($location->{"strand"} < 0) { $location->{"start"} -= $self->get('downstream_flank'); $location->{"end"} += $self->get('upstream_flank'); } else { $location->{"start"} -= $self->get('upstream_flank'); $location->{"end"} += $self->get('downstream_flank'); } } #sometimes users request more flanking sequence than is avaiable $location->{"start"} = 1 if ($location->{"start"} < 1); return $location; } sub _processSequence { my ($self, $locations) = @_; my $seq = ''; my $temp_Seq = ''; my $first_coding_exon_flag = 0; my $dna = $self->get('dna'); foreach my $rank (sort { $a <=> $b } keys %{$locations}) { my $location = $locations->{$rank}; my $chr = $location->{'chr'}; my $start = $location->{'start'}; my $end = $location->{'end'}; my $strand = exists( $location->{'strand'}) ? $location->{'strand'} : 1; my $phase = $location->{'phase'} || 0; if ($first_coding_exon_flag == 0) { if ($strand < 0) { $temp_Seq = $self->_rc( $dna-> getSequence( $chr, $start, $end ) ); } else { $temp_Seq = $dna->getSequence( $chr, $start, $end ); } if($temp_Seq) { # incase its not the first coding exon, # undef is returned by DNAAdapter if ($phase > 0) { # copying Ns in beginning, exactly as the # value of phase of first coding exon. $seq = 'N'x$phase; } $seq .= $temp_Seq; $first_coding_exon_flag = 1; } } else { if ($strand < 0) { $seq .= $self->_rc( $dna->getSequence( $chr, $start, $end ) ); } else { $seq .= $dna->getSequence( $chr, $start, $end ); } } } if (length($seq)) { return $seq; } return undef } sub _addRow { my ($self, $atable, $outRow, $sequence) = @_; push @{$outRow}, $sequence; $atable->addRow($outRow); $self->_incrementBatch; } #interface methods sub _getConfigurationTree { my ($self,$interface,$dsCounter)=@_;; return $self->getParam('configurator')->getConfigurationTree( $self->virtualSchema, $self->name, $interface, $dsCounter); } 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); my $logger=Log::Log4perl->get_logger(__PACKAGE__); $logger->info("QUERY XML: $xml"); 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)); } my $importable = $self->get('importable'); my $rtable = $importable->getTable(); my $attribute_count = @{$query->getAllAttributes}; if ($rtable->hashedResults || $attribute_count > 1){ $self->set('attribute_merge_required','1'); } my $has_rows = $rtable->hasMoreRows; my %avoidDuplication = (); NEXTROW: while ($has_rows && $self->_continueWithBatch($batch_size, $rtable)) { my $curRow = $rtable->nextRow; my $rowAsString = $self->toString($curRow); # avoid using "@A" eq "@B" type comparison, it floods error log next NEXTROW if (!$rowAsString || exists $avoidDuplication{$rowAsString} ) ; # no point retrieving sequence for same coordinates $avoidDuplication{$rowAsString} = ''; $self->_processRow( $atable, $curRow); #print "<BR> BAtch: $batch_size [ ", $self->_continueWithBatch($batch_size, $rtable), " ] [ ", $self->get('batchIndex'), " ]"; } # 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); ## additional logic to send all ORPHAN sequences at the end where the total number ## of transcripts encountered is not equal transcript_count. that may happen when you request a ## GENE type sequence with a filter restricting some transcripts. consequently you donot see all the ## transcripts and the sequence should not be kept by GS, so lets release them all at the end my $storageHash = $self->get('seqStorageHash'); if ($storageHash) { foreach my $fkey (%$storageHash) { my $sequence = $storageHash->{$fkey}->{'seq'}; my $outRow = $storageHash->{$fkey}->{'outRow'}; $self->_addRow($atable, $outRow, $sequence); } } } $importable->setTable($rtable); $self->set('importable', $importable); $self->get('dna')->close; return $atable; } ### sequence __recipes__ sub _codingCdnaPeptideSequences { my ($self, $atable, $curRow) = @_; # Determine this and last primary keys my $importable_indices = $self->get('importable_indices'); # Get the primary sequence ID from this row. Use DUMMY if missing my $pkey = $curRow ? ($curRow->[$importable_indices->{"pkey"}] || 'DUMMY') : undef; my $lastPkey = $self->get('lastPkey') || $pkey; my $locations = $self->get('locations'); my $outRow = $self->get('outRow'); if( ( ! defined $pkey ) or ( $pkey ne $lastPkey ) ){ # Start of new row, or end of results; Dump the current sequence #-- NOTE: #-- upstream and downstream sequences are added only once, and #-- not for every exons, as it used to be before. So _modFlank is #-- ran only at the very end, and once (once if upsream, once if #-- downstream) to modify $locations. [ben] #-- This is very important especially if you ask for both up AND downstream flanks #-- so you are sure that you enter only once in each _modFlank loop, and not twice #-- in the upstream as it was before. my ($up, $down); if ($self->get('upstream_flank')){$up = "up"} if ($self->get('downstream_flank')){$down = "down"} #-- get exons rank info my @ranks = sort { $a <=> $b } keys %{$locations}; my $array_size = scalar @ranks; my ($firstExon,$lastExon); if ($array_size == 1){ #-- if the gene/transcript has 1 exon only #-- so you can get a downstream flank even if rank = 1 $firstExon = shift @ranks; # firstExon and lastExon are both set to 1 $lastExon = $firstExon ; }else { #-- if the gene/transcript has more than 1 exon $firstExon = shift @ranks; #get first exon $lastExon = pop @ranks; #get last exon } #-- then run _modflank for the first and/or last exons if ($self->get('upstream_flank')) { my $location = $locations->{$firstExon}; $location = $self->_modFlanks($location, 2, $up); $locations->{$firstExon} = $location if ($location->{"start"}); } if ($self->get('downstream_flank')) { my $location = $locations->{$lastExon}; $location = $self->_modFlanks($location, 2, $down); $locations->{$lastExon} = $location if ($location->{"start"}); } my $sequence; if( grep{ $locations->{$_}->{"start"} } keys %$locations ) { $sequence = $self->_processSequence($locations); $sequence = $self->_translate($sequence) if ($self->get('translate')); $self->_editSequence(\$sequence); } if ($sequence) { $self->_addRow($atable, $outRow, $sequence); } else { $self->_addRow($atable, $outRow, "Sequence unavailable"); } $locations = {}; $outRow = undef; } # End sequence dumping if ($curRow) { # Update the location corresponding to this row my $rank = $curRow->[ $importable_indices->{"rank"} ]; # Requesting for phase info as well, to fix the bug of additional # Ns in the beginning - [syed] my $location = $self->_getLocationFrom($curRow, "chr", "start", "end", "strand", "phase", "codon_table_id", "seq_edits"); $self->set('codon_table_id',$location->{"codon_table_id"}); # $self->set('seq_edits',$location->{"seq_edits"}); # (New-2) hack for MartBuilder built marts. different set of exportables importables # need to store seq_edits as 12,12,U;14,14:U which is old formart # new exportable passes it in ensembl core format 12 12 U per row $self->set_seq_edits($location->{"seq_edits"}); # $location = $self->_modFlanks($location, 0); #-- This was soo wrong - it added a flank (upstream and/or downstream) #-- for every exons, instead of once, at the 5' or 3' end - [ben] $locations->{$rank} = $location if ($location->{"start"} && $rank); } $outRow ||= $self->_initializeReturnRow($curRow); $self->set('locations', $locations); $self->set('lastPkey', $pkey); $self->set('outRow', $outRow); } sub set_seq_edits { my ($self, $seq_edits) = @_; if ($seq_edits) { my $temp_str; # old format 12,12,U;14,14,M if ($seq_edits =~ m/\,/) { $temp_str = $seq_edits; } else { # new format 12 12 U $temp_str = $self->get('seq_edits'); $seq_edits =~ s/\s/\,/g; if ($temp_str) { $temp_str = $temp_str.';'.$seq_edits if($temp_str !~ m/$seq_edits/); } else { $temp_str = $seq_edits; } } $self->set('seq_edits', $temp_str); } else { $self->set('seq_edits', ""); } } sub _transcript_exonIntronFlankSequences { my ($self, $atable, $curRow) = @_; my $rank = 1; # Determine this and last primary keys my $importable_indices = $self->get('importable_indices'); # Get the primary sequence ID from this row. Use DUMMY if missing my $pkey = $curRow ? ($curRow->[$importable_indices->{"pkey"}] || 'DUMMY') : undef; my $lastPkey = $self->get('lastPkey') || $pkey; my $outRow = $self->get('outRow'); if( ( ! defined $pkey ) or ( $pkey ne $lastPkey ) ){ # Start of new row, or end of results; Dump the current sequence my $shift = ($self->get('seq_name') =~ m/flank/); my $location = $self->_modFlanks( $self->get('calc_location'), $shift ); $self->set('calc_location', undef); # Reset location cache my $sequence; if ($location->{"start"}) { my $locations = { $rank => $location }; $sequence = $self->_processSequence($locations); $self->_editSequence(\$sequence); } if ($sequence) { $self->_addRow($atable, $outRow, $sequence); } else { $self->_addRow($atable, $outRow, "Sequence unavailable"); } $outRow = undef; } # End sequence dumping if ($curRow) { # Update the location corresponding to this row my $location = $self->_getLocationFrom($curRow, "chr", "start", "end", "strand"); $self->_calcSeqOverLocations( $location ); } $outRow ||= $self->_initializeReturnRow($curRow); $self->set('lastPkey', $pkey); $self->set('outRow', $outRow); } sub _gene_exonIntronFlankSequences { my ($self, $atable, $curRow) = @_; my $rank = 1; # Determine this and last primary keys my $importable_indices = $self->get('importable_indices'); # Get the primary sequence ID from this row. Use DUMMY if missing my $pkey = $curRow ? ($curRow->[$importable_indices->{"pkey"}] || 'DUMMY') : undef; my $lastPkey = $self->get('lastPkey') || $pkey; my $outRow = $self->get('outRow'); if( ( ! defined $pkey ) or ( $pkey ne $lastPkey ) ) { # Start of new row, or end of results; Dump the current sequence my $shift = ($self->get('seq_name') =~ m/flank/); my $location = $self->_modFlanks( $self->get('calc_location'), $shift ); $self->set('calc_location', undef); # Reset location cache my $sequence; if ($location->{"start"}) { my $locations = { $rank => $location }; $sequence = $self->_processSequence($locations); $self->_editSequence(\$sequence); } my $storageHash = $self->get('seqStorageHash'); # This IF statement is just there to keep those in sink who have Ensembl 42 or earlier # and this new code. They should not expect transcript_pkey and transcript_count from XML if (exists $storageHash->{$lastPkey}) { # this is well tricky. lastDS is 2 when its a two DS query adn its the second DS's GS. # we do not wait for all transcripts to arrive because in a two DS query HUMAN/MSD, there is a chance # few transcripts are absent for not having any PDBids. So we should not stop such sequences. # The sequences are meant to stop till all the transcripts are seen when its a single DS query if (($storageHash->{$lastPkey}->{'transcriptCount'} >= $storageHash->{$lastPkey}->{'totalTranscripts'}) || $self->lastDS() == 2) { # use the old sequence from previous batch, the new one in some case (f_gene,c_f_gene,u_gene is different) $sequence = $storageHash->{$lastPkey}->{'seq'} if ($storageHash->{$lastPkey}->{'seq'}); if ($sequence) { $self->_addRow($atable, $outRow, $sequence); } else { $self->_addRow($atable, $outRow, "Sequence unavailable"); } delete $storageHash->{$lastPkey}; # over with this Gene. } else { # sequences gets an update, as and when a transcript of a gene # arrives in a separate batch. In principle this should always be the same #my $onHoldkeys = $self->get('onHoldSeqsKeys'); $storageHash->{$lastPkey}->{'seq'} = $sequence; $storageHash->{$lastPkey}->{'outRow'} = $outRow; #$onHoldkeys->{$lastPkey} = ""; #$self->set('onHoldSeqsKeys', $onHoldkeys); $self->set('seqStorageHash', $storageHash); $self->_incrementBatch; } } else { if ($sequence) { $self->_addRow($atable, $outRow, $sequence); } else { $self->_addRow($atable, $outRow, "Sequence unavailable"); } } $outRow = undef; } # End sequence dumping if ($curRow) { # Update the location corresponding to this row my $location = $self->_getLocationFrom($curRow,"pkey","transcript_pkey","chr", "start", "end", "strand", "transcript_count"); $self->_calcSeqOverLocations( $location ); # This IF statement is just there to keep those in sink who have Ensembl 42 or earlier # and this new code. They should not expect transcript_pkey and transcript_count from XML if ($location->{'transcript_pkey'}) { my $storageHash = $self->get('seqStorageHash'); my $prkey = $location->{'pkey'}; my $transcript_key = $location->{'transcript_pkey'}; #print "<BR>TRANS KEY: $transcript_key", " == ", $location->{'transcript_count'}; $storageHash->{$prkey}->{'totalTranscripts'} = $location->{'transcript_count'} || 1; # there are NULL vals in DB if(exists $storageHash->{$prkey}->{'transcriptKeys'}) { if ($storageHash->{$prkey}->{'transcriptKeys'} !~ m/$transcript_key/) { $storageHash->{$prkey}->{'transcriptKeys'} .= $transcript_key.','; $storageHash->{$prkey}->{'transcriptCount'} ++; # print "<BR> I WAS HERE", $storageHash->{$prkey}->{'transcriptCount'}; } } else { $storageHash->{$prkey}->{'transcriptKeys'} = $transcript_key.','; $storageHash->{$prkey}->{'transcriptCount'}++; } #print "<BR>",Dumper ($storageHash); $self->set('seqStorageHash', $storageHash); } } $outRow ||= $self->_initializeReturnRow($curRow); $self->set('lastPkey', $pkey); $self->set('outRow', $outRow); } sub _exonSequences { my ($self, $atable, $curRow) = @_; $curRow || return; # Process row-by-row; can discard last (empty) call return if ($self->_ignoreRow($curRow)); #ignore duplicate exons my $rank = 1; my $locations = {}; $locations->{$rank} = $self->_modFlanks( $self->_getLocationFrom($curRow, "chr", "start", "end", "strand"), 0 ); my $sequence; if ($locations->{1}->{"start"}) { $sequence = $self->_processSequence($locations); $self->_editSequence(\$sequence); } if ($sequence) { $self->_addRow($atable, $self->_initializeReturnRow($curRow), $sequence); } else { $self->_addRow($atable, $self->_initializeReturnRow($curRow), "Sequence unavailable"); } if ($self->get('ignore')) { #will only be true for gene_exon my $ignore = $self->get('ignore'); my $ignore_row = $self->get('ignore_row'); my $ref = $self->_getLocationFrom($curRow, $ignore_row); $ignore->{ $ref->{ $ignore_row } } = 1; #skip duplicate pkeys $self->set('ignore', $ignore); } #else there will be no last entry } sub _rawSequences { my ($self, $atable, $curRow) = @_; 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); } } } sub _utrSequences { my ($self, $atable, $curRow) = @_; my $locations = $self->get('locations'); my $outRow = $self->get('outRow'); if ($curRow) { my $importable_indices = $self->get('importable_indices'); my $pkey = $curRow->[ $importable_indices->{"pkey"} ]; my $lastPkey = $self->get('lastPkey'); if ( $lastPkey && ($pkey ne $lastPkey) ) { my $hasUtr = keys %{$locations}; if ($hasUtr) { my @ranks = sort { $a <=> $b } keys %{$locations}; my $low_rank = shift @ranks; my $hi_rank = ($hasUtr > 1) ? pop @ranks : $low_rank; my $calc_location = $self->get('calc_location'); my $low_loc_strand = $locations->{$low_rank}->{"strand"}; if ($self->get('upstream_flank')) { if ($low_loc_strand < 0) { $calc_location->{"start"} = $calc_location->{"end"} + 1; $calc_location->{"end"} = $calc_location->{"start"} + $self->get('upstream_flank') - 1; } else { $calc_location->{"end"} = $calc_location->{"start"} - 1; $calc_location->{"start"} = $calc_location->{"start"} - $self->get('upstream_flank');# lose a base if include this + 1; $calc_location->{"start"} = 1 if ($calc_location->{"start"} < 1); } #prepend to sequence $locations->{$low_rank - 1} = $calc_location; } elsif ($self->get('downstream_flank')) { if ($low_loc_strand < 0) { $calc_location->{"end"} = $calc_location->{"start"} - 1; $calc_location->{"start"} = $calc_location->{"end"} - $self->get('downstream_flank') + 1; } else { $calc_location->{"start"} = $calc_location->{"end"} + 1; $calc_location->{"end"} = $calc_location->{"start"} + $self->get('downstream_flank') -1; # positive strand & downstream flank, wont lose the base. } #append to sequence $locations->{$hi_rank + 1} = $calc_location; } my $sequence = $self->_processSequence($locations); $self->_editSequence(\$sequence); if ($sequence) { $locations = {}; $hasUtr = 0; $self->_addRow($atable, $outRow, $sequence); $outRow = undef; $self->set('calc_location', undef); } } else { $self->_addRow($atable, $outRow, "Sequence unavailable"); $outRow = undef; } } unless ($outRow) { $outRow = $self->_initializeReturnRow($curRow); } #for this one we build and calc my $rank = $curRow->[ $importable_indices->{"rank"} ]; my $location = $self->_getLocationFrom($curRow, "chr", "start", "end", "strand"); if ($location->{"start"}) { $locations->{$rank} = $location; $self->_calcSeqOverLocations($location); } $self->set('lastPkey', $pkey); } else { # last entry my $hasUtr = keys %{$locations}; if ($hasUtr) { my @ranks = sort { $a <=> $b } keys %{$locations}; my $low_rank = shift @ranks; my $hi_rank = ($hasUtr > 1) ? pop @ranks : $low_rank; my $calc_location = $self->get('calc_location'); my $low_loc_strand = $locations->{$low_rank}->{"strand"}; if ($self->get('upstream_flank')) { if ($low_loc_strand < 0) { $calc_location->{"start"} = $calc_location->{"end"} + 1; $calc_location->{"end"} = $calc_location->{"start"} + $self->get('upstream_flank') - 1; } else { $calc_location->{"end"} = $calc_location->{"start"} - 1; $calc_location->{"start"} = $calc_location->{"start"} - $self->get('upstream_flank') + 1; $calc_location->{"start"} = 1 if ($calc_location->{"start"} < 1); } #prepend to sequence $locations->{$low_rank - 1} = $calc_location; } elsif ($self->get('downstream_flank')) { if ($low_loc_strand < 0) { $calc_location->{"end"} = $calc_location->{"start"} - 1; $calc_location->{"start"} = $calc_location->{"end"} - $self->get('downstream_flank') + 1; } else { $calc_location->{"start"} = $calc_location->{"end"} + 1; $calc_location->{"end"} = $calc_location->{"start"} + $self->get('downstream_flank') - 1; } #append to sequence $locations->{$hi_rank + 1} = $calc_location; } my $sequence = $self->_processSequence($locations); $self->_editSequence(\$sequence); if ($sequence) { $locations = {}; $hasUtr = 0; $self->_addRow($atable, $outRow, $sequence); $outRow = undef; $self->set('calc_location', undef); } } else { $self->_addRow($atable, $outRow, "No UTR is annotated for this transcript"); $outRow = undef; } } $self->set('locations', $locations); $self->set('outRow', $outRow); } sub _snpSequences { my ($self, $atable, $curRow) = @_; my $rank = 1; if ($curRow) { #since locations are just hashes, hack them my $location = $self->_getLocationFrom($curRow, "chr", "pos", "strand", "allele"); if ($location->{"strand"} < 0) { $location->{"start"} = $location->{"pos"} - $self->get('downstream_flank'); $location->{"start"} = 1 if ($location->{"start"} < 1); $location->{"end"} = $location->{"pos"} + $self->get('upstream_flank'); $location->{"off"} = $self->get('upstream_flank'); } else { $location->{"start"} = $location->{"pos"} - $self->get('upstream_flank'); $location->{"end"} = $location->{"pos"} + $self->get('downstream_flank'); $location->{"off"} = $self->get('upstream_flank'); if ($location->{"start"} < 1) { $location->{"off"} = $self->get('upstream_flank') + $location->{"start"} - 1; $location->{"start"} = 1; } } my $locations = {}; $locations->{$rank} = $location; my $sequence = $self->_processSequence($locations); $self->_editSequence(\$sequence); if ($sequence) { substr($sequence, $location->{"off"}, 1) = "%".$location->{"allele"}."%"; $self->_addRow($atable, $self->_initializeReturnRow($curRow), $sequence); } } } =head2 toString Usage : Description: Helper rouinte for Array comparison by converting it to string first. Comparing arrays as "@A" eq "@B" brings flood of warning Returntype : String Exceptions : none Caller : caller =cut sub toString { my ($self, $curRow) = @_; my $string; foreach (@{$curRow}) { $string .= $_ if ($_); } return $string; } 1;