Raw content of Bio::EnsEMBL::Analysis::RunnableDB::UTR_Builder # # Ensembl module for Bio::EnsEMBL::Analysis::RunnableDB::UTR_Builder # # Cared for by EnsEMBL <ensembl-dev@ebi.ac.uk> # # Copyright GRL and EBI # # You may distribute this module under the same terms as perl itself # POD documentation - main docs before the code =head1 NAME Bio::EnsEMBL::Analysis::RunnableDB::UTR_Builder =head1 SYNOPSIS my $utrbuilder_runnable = new Bio::EnsEMBL::Analysis::RunnableDB::UTR_Builder( -db => $db, -input_id => $input_id ); $utrbuilder_runnable->fetch_input(); $utrbuilder_runnable->run(); $utrbuilder_runnable->output(); $utrbuilder_runnable->write_output(); #writes to DB =head1 DESCRIPTION This is the new version of the UTR-addition procedure. It combines predictions made from proteins with cDNA alignments to add UTR regions to the gene models. It can also inlcude ESTs and ditags. It uses code from Coalescer/Consensus to produce score for the alternative models and chose the best option. It also includes ("look-for-both") code to correct the phases of the transcripts unless they are "blessed" and inculdes the option to check for predefined protein/cDNA pairing as a first step, looking for NM/NPentries in a GeneBank file. Config files to set-up are Bio::EnsEMBL::Analysis::Config::GeneBuild::UTR_Builder Bio::EnsEMBL::Analysis::Config::Databases Bio::EnsEMBL::Analysis::Config::GeneBuild::TranscriptConsensus (just a copy of the example file) Bio::EnsEMBL::Analysis::Config::GeneBuild::KillListFilter =head1 CONTACT Ensembl - ensembl-dev@ebi.ac.uk =head1 APPENDIX The rest of the documentation details each of the object methods. Internal methods are usually preceded with a '_' =cut package Bio::EnsEMBL::Analysis::RunnableDB::UTR_Builder; use vars qw(@ISA); use strict; use Bio::SeqIO; use Bio::EnsEMBL::Utils::Exception qw( throw warning verbose); use Bio::EnsEMBL::Analysis::RunnableDB; use Bio::EnsEMBL::Gene; use Bio::EnsEMBL::Transcript; use Bio::EnsEMBL::Translation; use Bio::EnsEMBL::Exon; use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::ExonUtils qw( clone_Exon transfer_supporting_evidence validate_Exon_coords ); use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::TranscriptUtils qw( is_Transcript_sane all_exons_are_valid intron_lengths_all_less_than_maximum set_start_codon set_stop_codon clone_Transcript has_no_unwanted_evidence ); use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::TranslationUtils qw( validate_Translation_coords compute_translation contains_internal_stops print_Translation ); use Bio::EnsEMBL::Analysis::Tools::Algorithms::ClusterUtils; use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::TranscriptExtended; use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::ExonExtended; use Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::CollapsedCluster; use Bio::EnsEMBL::Analysis::Tools::Utilities; use Bio::EnsEMBL::Analysis::Runnable::TranscriptConsensus; use Bio::EnsEMBL::Analysis::RunnableDB::BaseGeneBuild; use Bio::EnsEMBL::KillList::KillList; use Bio::EnsEMBL::KillList::DBSQL::DBAdaptor; #CONFIG FILES use Bio::EnsEMBL::Analysis::Config::GeneBuild::UTR_Builder; use Bio::EnsEMBL::Analysis::Config::GeneBuild::TranscriptConsensus; use Bio::EnsEMBL::Analysis::Config::GeneBuild::KillListFilter; @ISA = qw ( Bio::EnsEMBL::Analysis::RunnableDB::BaseGeneBuild ); =head2 new Arg [1] : Bio::EnsEMBL::Analysis::RunnableDB::UTR_Builder Function : instatiates object & check config file Returntype: Bio::EnsEMBL::Analysis::RunnableDB::UTR_Builder =cut sub new { my ($class, @args) = @_; my $self = $class->SUPER::new(@args); $self->read_and_check_config($UTR_BUILDER_CONFIG_BY_LOGIC); return $self; } my $totalgenes = 0; =head2 fetch_input Arg [1] : none Description: Get all raw data needed from the databases Returntype : none =cut sub fetch_input{ my ($self) = @_; my $slice = $self->fetch_sequence(undef, $self->INPUT_DB); if(!$slice){ throw "can't fetch input slice!\n" } $self->query($slice); # fetch all general input genes foreach my $input_genetype (@{$self->INPUT_GENETYPES}) { my $input_genes = $self->query->get_all_Genes_by_type($input_genetype); print STDERR "got " . scalar(@{$input_genes}) . " $input_genetype genes [ ". $self->INPUT_DB->dbname() . "@" . $self->INPUT_DB->host ." ]\n"; $self->gw_genes( $input_genes ); } # get blessed genes my $blessed_slice; my @blessed_genes; my $blessed_type; if($self->BLESSED_DB and scalar(@{$self->BLESSED_GENETYPES})){ if ($self->BLESSED_DB){ #fetch blessed genes from blessed db $blessed_slice = $self->BLESSED_DB->get_SliceAdaptor->fetch_by_name($self->input_id); } else{ #fetch blessed genes from gw db $blessed_slice = $self->query; } foreach my $bgt( @{$self->BLESSED_GENETYPES} ){ my $blessed_genes = $blessed_slice->get_all_Genes_by_type($bgt); print STDERR "got " . scalar(@{$blessed_genes}) . " $bgt genes [ ". $self->BLESSED_DB->dbname() . "@" . $self->BLESSED_DB->host ." ]\n"; $self->blessed_genes( $blessed_slice, $blessed_genes ); $blessed_type .= $bgt.""; } # store all blessed type names for VIP treatment $self->{'blessed_type'} = ($blessed_type.$self->BLESSED_UTR_GENETYPE); } # are there any genes here at all? if(!scalar @{$self->gw_genes} and !scalar @{$self->blessed_genes}){ print STDERR "\nNo genes found here.\n\n"; return 0; } # get cdnas if(!(defined($self->cDNA_GENETYPE) and $self->CDNA_DB)){ warn("NOT using any cDNAs!?\n"); } else{ my $cdna_vc = $self->CDNA_DB->get_SliceAdaptor->fetch_by_name($self->input_id); $self->_cdna_slice($cdna_vc); foreach my $cdna_type (@{$self->cDNA_GENETYPE}){ #my @cdna_genes = @{$cdna_vc->get_all_Genes_by_type($cdna_type)}; my @cdna_genes = @{$self->_cdna_slice->get_all_Genes_by_type($cdna_type)}; print STDERR "got ".scalar(@cdna_genes)." ".$cdna_type." cDNAs.\n" if $self->VERBOSE; $self->CDNA_DB->dbc->disconnect_when_inactive(1); # filter cdnas my $filtered_cdna = $self->_filter_cdnas(\@cdna_genes, 0); $self->cdna_genes($filtered_cdna); print STDERR "got " . scalar(@{$filtered_cdna}) . " cdnas after filtering.\n" if $self->VERBOSE; } } # get ESTs if(defined($self->EST_GENETYPE) && $self->EST_DB){ my $est_vc = $self->EST_DB->get_SliceAdaptor->fetch_by_name($self->input_id); my @est_genes = @{$est_vc->get_all_Genes_by_type($self->EST_GENETYPE)}; print STDERR "got " . scalar(@est_genes) . " ".$self->EST_GENETYPE." ESTs.\n"; if($self->FILTER_ESTS){ my $filtered_ests = $self->_filter_cdnas(\@est_genes, 1); $self->ests($filtered_ests); } else{ $self->ests(\@est_genes); } print STDERR "got " . scalar(@{$self->ests}) . " ESTs after filtering.\n" if $self->VERBOSE; } # get ditags my ($dfa, $ditag_slice); my @ditags; if((scalar @{$self->DITAG_TYPE_NAMES}) && $self->DITAG_DB){ $dfa = $self->DITAG_DB->get_DitagFeatureAdaptor; $ditag_slice = $self->DITAG_DB->get_SliceAdaptor->fetch_by_name($self->input_id); foreach my $ditag_type (@{$self->DITAG_TYPE_NAMES}) { my @type_ditags = @{$dfa->fetch_pairs_by_Slice($ditag_slice, $ditag_type)}; print STDERR "got " . scalar(@type_ditags) . " ".$ditag_type." ditags.\n" if $self->VERBOSE; push(@ditags, @type_ditags); } } if(scalar @ditags){ @ditags = sort {($a->{'start'} <=> $b->{'start'}) or ($a->{'end'} <=> $b->{'end'})} @ditags; $self->ditags(\@ditags); print STDERR "got " . scalar(@ditags) ." ditags.\n"; } else{ print STDERR "not using Ditags.\n"; } # db disconnections $self->INPUT_DB->dbc->disconnect_when_inactive(1); $self->CDNA_DB->dbc->disconnect_when_inactive(1) if($self->CDNA_DB); $self->EST_DB->dbc->disconnect_when_inactive(1) if($self->EST_DB); $self->BLESSED_DB->dbc->disconnect_when_inactive(1) if($self->BLESSED_DB); $self->DITAG_DB->dbc->disconnect_when_inactive(1) if($self->DITAG_DB); # set evidence sets for Coalescer code my (@est_biotypes, @cdna_logicnames); push(@est_biotypes, $self->EST_GENETYPE); my @simgw_biotypes = @{ $self->INPUT_GENETYPES }; push(@cdna_logicnames, @{$self->cDNA_GENETYPE}); $self->{evidence_sets} = { 'est' => \@est_biotypes, 'simgw' => \@simgw_biotypes, 'cdna' => \@cdna_logicnames, }; # prune genes during filtering? # (usually not used, as too many good things are thrown out.) $self->prune($self->PRUNE_GENES); # get rid of identical models? # (usually not used, as we want the same number of output as input genes) $self->{'remove_redundant'} = 0; #get killlist entries for cDNAs to ignore my $kill_list_type = "cDNA"; my $kill_list = $self->populate_kill_list($kill_list_type); $self->kill_list($kill_list); if($self->LOOK_FOR_KNOWN){ #read KnowUTR pairing into hash $self->create_predefined_pairing($self->KNOWNUTR_FILE); } } =head2 run Description: general run method Returntype : none =cut sub run { my ($self) = @_; # filter the genewise predictions to remove fragments. # Make sure we keep the fragments to let the # genebuilder do the final sorting out, but don't add UTRs to them. print "TOTAL GENES IN: ".(scalar @{$self->gw_genes} + scalar @{$self->blessed_genes})."\n"; my $ininumber = scalar @{$self->gw_genes}; my $filtered_genes = $self->filter_genes($self->gw_genes); print STDERR "filtered from $ininumber to ".(scalar @$filtered_genes). " genes.\n" if $self->VERBOSE && $filtered_genes; #for blessed genes, only check region and existing UTR $ininumber = scalar @{$self->blessed_genes}; my $filtered_blessed_genes = $self->filter_genes($self->blessed_genes, 1); print STDERR "filtered blessed from $ininumber to ".(scalar @$filtered_blessed_genes). " genes.\n" if $self->VERBOSE && $filtered_blessed_genes; # match gene-model to cDNAs/ESTs for blessed & other genes $self->run_matching($filtered_blessed_genes, $self->BLESSED_UTR_GENETYPE, 1); $self->run_matching($filtered_genes, $self->UTR_GENETYPE, 0); print STDERR "Have ".(scalar @{$self->combined_genes})." combined_genes & ". (scalar @{$self->unmatched_genes})." unmatched_genes.\n" if $self->VERBOSE; # check translation, start and stop my @remapped; push( @remapped, @{$self->remap_genes( $self->combined_genes, undef )} ); push( @remapped, @{$self->remap_genes( $self->unmatched_genes, $self->EXTEND_BIOTYPE_OF_UNCHANGED_GENES )} ); #results are stored in "output" $self->output(@remapped); print "TOTAL GENES OUT: ".(scalar @remapped)."\n"; if($self->VERBOSE){ foreach ( @remapped ) { print "final biotype : " . $_->biotype . "\n"; } } } =head2 filter_genes Arg [1] : ref to array of Bio::EnsEMBL::Gene Description: We do not want to add UTRs to fragments of genes otherwise they get boosted in importance during the final gene build and produce rubbish alternative transcripts and in the worst cases can misjoin clusters especially if there is some misassembled/duplicated sequence. One solution is to cluster the genewises at this stage in much the same way as we do during the GeneBuilder run, and attempt to add UTRs ony to the best ones eg fullest genomic extent. Also checks for existing UTR, these genes are just passed on Returntype : ref to array of Bio::EnsEMBL::Gene =cut sub filter_genes { my ($self, $genesref, $blessed) = @_; my @tested_genes; my $pruned_CDS; #run tests: get rid of long-intron genes, etc. GENES: foreach my $test_gene (@$genesref){ if( $test_gene ){ my $test_transcript = $test_gene->get_all_Transcripts->[0]; #would like to avoid slow sorting call - but the "get_all_Transcripts" does not do this job properly... $test_transcript->sort; if($test_transcript->three_prime_utr or $test_transcript->five_prime_utr){ print STDERR "Gene has UTR already. Skipping.\n" if $self->VERBOSE; $self->unmatched_genes($test_gene); next GENES; } if($blessed){ #for blessed genes only check the region if($test_transcript->start < 1 && $test_transcript->end > 1){ print STDERR "ignoring blessed gene: falls off the slice by its lower end\n" if $self->VERBOSE; } else{ push(@tested_genes, $test_gene); } } else{ if( is_Transcript_sane($test_transcript) && all_exons_are_valid($test_transcript, $self->MAX_EXON_LENGTH) && intron_lengths_all_less_than_maximum($test_transcript, $self->MAX_INTRON_LENGTH) && ($test_transcript->start < $self->query->length) && ($test_transcript->end > 1) && has_no_unwanted_evidence($test_transcript) ) { push(@tested_genes, $test_gene); } else{ #test if it falls off the slice on the left, these will be picked up by another job if($test_transcript->start < 1 && $test_transcript->end > 1){ print STDERR "ignoring gene: falls off the slice by its lower end\n" if $self->VERBOSE; } else{ #keep them in the set $self->unmatched_genes($test_gene); } } } } else{ throw "\nERROR: Cant check genes!\n\n"; } } $genesref = \@tested_genes; #cluster genes (stores internally) if(!$blessed){ # all the genes are single transcript at this stage, cluster them my $clustered_genes = $self->cluster_CDS($genesref); #pruning if desired if($self->prune){ $genesref = $self->prune_CDS($clustered_genes); } } return $genesref; } =head2 cluster_CDS Arg [1] : ref to array of Bio::EnsEMBL::Gene Description: Separates CDSs according to strand and then clusters each set by calling Bio::EnsEMBL::Analysis::Tools::Algorithms::ClusterUtils:cluster_Genes Returntype : ref to array of Bio::EnsEMBL::Analysis::Tools::Algorithms::GeneCluster =cut sub cluster_CDS { my ($self, $CDSref) = @_; my @forward_CDS; my @reverse_CDS; foreach my $gene (@$CDSref){ my @exons = @{ $gene->get_all_Exons }; if ( $exons[0]->strand == 1 ){ push( @forward_CDS, $gene ); } else{ push( @reverse_CDS, $gene ); } } my (@clusters, @forward_clusters, @reverse_clusters); my ($forward_clusters, $forward_clusters_rest); my ($reverse_clusters, $reverse_clusters_rest); if ( @forward_CDS ){ ($forward_clusters, $forward_clusters_rest) = cluster_Genes( \@forward_CDS, $self->{evidence_sets} ); push(@forward_clusters, @$forward_clusters); push(@forward_clusters, @$forward_clusters_rest); } if ( @reverse_CDS ){ ($reverse_clusters, $reverse_clusters_rest) = cluster_Genes( \@reverse_CDS, $self->{evidence_sets} ); push(@reverse_clusters, @$reverse_clusters); push(@reverse_clusters, @$reverse_clusters_rest); } #store the two cluster sets and return combined set if ( scalar @forward_clusters ){ #store for later $self->forward_genewise_clusters(\@forward_clusters); push( @clusters, @forward_clusters); } if ( scalar @reverse_clusters ){ #store for later $self->reverse_genewise_clusters(\@reverse_clusters); push( @clusters, @reverse_clusters); } return \@clusters; } =head2 prune_CDS Arg [1] : ref to array of GeneClusters Description: rejects duplicate CDS, makes sure they are kept on unmatched_genes or they (and their supporting evidence) will be lost for the rest of the build Note: pruning will not be used unless specified otherwise Returntype : ref to array of Bio::EnsEMBL::Gene =cut sub prune_CDS { my ($self, $gene_cluster_ref) = @_; my $cluster_count = 0; my @pruned_transcripts; my @pruned_genes; CLUSTER: foreach my $gene_cluster ( @$gene_cluster_ref ){ $cluster_count++; my @unpruned_genes; #check for multi-transcript genes foreach my $unpruned_gene ($gene_cluster->get_Genes){ if(scalar(@{$unpruned_gene->get_all_Transcripts}) > 1){ throw($unpruned_gene->dbID . " has >1 transcript - can't handle it yet \n"); } push(@unpruned_genes, $unpruned_gene); } # sort the unpruned_genes by total exon length of their single transcripts - there are no UTRs to worry about yet @unpruned_genes = sort{$b->get_all_Transcripts->[0]->length <=> $a->get_all_Transcripts->[0]->length} @unpruned_genes; # do we really just want to take the first transcript only? # If there's a very long single exon gene we will lose any underlying multi-exon transcripts # this may increase problems with the loss of valid single exon genes as mentioned below. # it's a balance between keeping multi exon transcripts and losing single exon ones my $maxexon_number = 0; foreach my $gene (@unpruned_genes){ my $exon_number = scalar(@{$gene->get_all_Transcripts->[0]->get_all_Exons}); if ( $exon_number > $maxexon_number ){ $maxexon_number = $exon_number; } } if ($maxexon_number == 1){ # ie the longest transcript is a single exon one # take the longest: push (@pruned_genes, $unpruned_genes[0]); print STDERR "VAC: found single_exon_transcript: " .$unpruned_genes[0]->dbID. "\n" if $self->VERBOSE; shift @unpruned_genes; $self->unmatched_genes(@unpruned_genes); next CLUSTER; } # otherwise we need to deal with multi exon transcripts and reject duplicates. # links each exon in the transcripts of this cluster with a hash of other exons it is paired with my %pairhash; # allows retrieval of exon objects by exon->id - convenience my %exonhash; # prune redundant transcripts GENE: foreach my $gene (@unpruned_genes) { my @exons = @{$gene->get_all_Transcripts->[0]->get_all_Exons}; my $i = 0; my $found = 1; # 10.1.2002 VAC we know there's a potential problem here - single exon transcripts which are in a # cluster where the longest transcript has > 1 exon are not going to be considered in # this loop, so they'll always be marked "transcript already seen" # How to sort them out? If the single exon overlaps an exon in a multi exon transcript then # by our rules it probably ought to be rejected the same way transcripts with shared exon-pairs are. # Tough one. EXONS: for ($i = 0; $i < $#exons; $i++) { my $foundpair = 0; my $exon1 = $exons[$i]; my $exon2 = $exons[$i+1]; # go through the exon pairs already stored in %pairhash. # If there is a pair whose exon1 overlaps this exon1, and # whose exon2 overlaps this exon2, then these two transcripts are paired EXONHASH: foreach my $first_exon_id (keys %pairhash) { my $first_exon = $exonhash{$first_exon_id}; foreach my $second_exon_id (keys %{$pairhash{$first_exon}}) { my $second_exon = $exonhash{$second_exon_id}; if ( $exon1->overlaps($first_exon) && $exon2->overlaps($second_exon) ) { $foundpair = 1; last EXONHASH; # this method allows a transcript to be covered by exon pairs # from different transcripts, rejecting possible # splicing variants. Needs rethinking } } } if ($foundpair == 0) { # ie this exon pair does not overlap with a pair yet found in another transcript $found = 0; # ie currently this transcript is not paired with another # store the exons so they can be retrieved by id $exonhash{$exon1} = $exon1; $exonhash{$exon2} = $exon2; # store the pairing between these 2 exons $pairhash{$exon1}{$exon2} = 1; } } # end of EXONS # decide whether this is a new transcript or whether it has already been seen if ($found == 0) { push(@pruned_genes, $gene); } elsif ($found == 1 && $#exons == 0){ $self->unmatched_genes($gene); } else { # dont know what this is about... $self->unmatched_genes($gene); if ( $gene == $unpruned_genes[0] ){ } } } # end of this gene } #end CLUSTER return \@pruned_genes; } =head2 run_matching Arg [1] : ref to array of Bio::EnsEMBL::Gene Arg [2] : string representing genetype to be associated with UTR-modified genes Arg [3] : bool indicating that we are NOT trying to map to predefined cDNA (blessed genes) Description: main function which matches CDS (genewise) predictions to cDNA alignments considering only terminal CDS exons, looking for "knowns" first, optionally using ESTs and ditags. Returntype : none =cut sub run_matching{ my ($self, $genesref, $combined_genetype, $blessed) = @_; # merge exons with frameshifts into a big exon my @merged_genes = (); #maybe we should not do merging with blessed genes at all?! @merged_genes = $self->_merge_genes($genesref); print STDERR "\n --- \ngot " . scalar(@merged_genes) . " merged " . $combined_genetype . " genes\n" if $self->VERBOSE; # sort genewises by exonic length and genomic length @merged_genes = sort { my $result = ( ($b->get_all_Transcripts->[0]->end - $b->get_all_Transcripts->[0]->start + 1) <=> ($a->get_all_Transcripts->[0]->end - $b->get_all_Transcripts->[0]->start + 1) ); unless ($result){ return ( ($b->get_all_Transcripts->[0]->length) <=> ($a->get_all_Transcripts->[0]->length) ); } return $result; } @merged_genes; # find one-2-one matching between proteins and cdnas CDS: foreach my $cds (@merged_genes){ print STDERR "--next cds: ".$cds->seq_region_start."-".$cds->seq_region_end." --\n" if $self->VERBOSE; # should be only 1 transcript my @cds_tran = @{$cds->get_all_Transcripts}; my @cds_exons = @{$cds_tran[0]->get_all_Exons}; # ordered array of exons my $strand = $cds_exons[0]->strand; my $cdna_match; my $usingKnown = 0; if($cds_exons[$#cds_exons]->strand != $strand){ warn("first and last cds exons have different strands - can't make a sensible combined gene\n"); # get and store unmerged version of cds my $unmerged_cds = $self->retrieve_unmerged_gene($cds); $self->unmatched_genes($unmerged_cds); next CDS; } #look for a pre-defined pairing between a protein and a cDNA the gene was build from my $predef_match = undef; my $combined_transcript = undef; if($self->LOOK_FOR_KNOWN){ $predef_match = $self->check_for_predefined_pairing($cds, $blessed); if(defined $predef_match){ $combined_transcript = $self->combine_genes($cds, $predef_match); $combined_transcript->sort; # just check combined transcript works before throwing away the original transcript if ( $combined_transcript && is_Transcript_sane($combined_transcript) && all_exons_are_valid($combined_transcript, $self->MAX_EXON_LENGTH, 1) && intron_lengths_all_less_than_maximum($combined_transcript, $self->MAX_INTRON_LENGTH) && validate_Translation_coords($combined_transcript, 1) && !contains_internal_stops($combined_transcript) && $combined_transcript->translate && has_no_unwanted_evidence($combined_transcript) ){ # make sure combined transcript doesn't misjoin any genewise clusters if($self->find_cluster_joiners($combined_transcript)){ print STDERR "Found a cluster_joiner!\n" if $self->VERBOSE; print STDERR $combined_transcript->seq_region_start."/".$combined_transcript->seq_region_end."\n" if $self->VERBOSE; #dont use this one $combined_transcript = undef; } } else{ print SDERR "Did not pass tests.\n" if $self->VERBOSE; $combined_transcript = undef; } } } if($predef_match && $combined_transcript){ #use this cDNA print STDERR "Using predefined cDNA!\n" if $self->VERBOSE; $cdna_match = $predef_match; $usingKnown = 1; } else{ # find matching cdnas using scoring of all evidence available my ($matching_cdnas, $utr_length_hash, $UTR_side_indicator_hash) = $self->match_protein_to_cdna($cds, 0); #we probably dont need any of this UTR stuff anymore!? my %utr_length_hash = %$utr_length_hash; if(!(scalar @$matching_cdnas)){ warn("Could not identify any matching cDNA!\n"); #try to use ESTs instead #we could use coalescer code to produce joined ESTs (find longest path, etc.)? ($matching_cdnas, $utr_length_hash, $UTR_side_indicator_hash) = $self->match_protein_to_cdna($cds, 1); %utr_length_hash = %$utr_length_hash; if(!(scalar @$matching_cdnas)){ my $unmerged_cds = $self->retrieve_unmerged_gene($cds); warn("Could not identify any matching ESTs either!\n"); $self->unmatched_genes($unmerged_cds); next CDS; } } ## scoring code... #convert genes to extended objects my $matching_extended_cdnas = $self->convert_to_extended_genes($matching_cdnas); #cluster matching cDNA or EST genes #call cluster method from ClusterUtils # print "EXTENDED : ",join(' ',@{$matching_extended_cdnas}),"\n"; #foreach my $keys (keys %{$self->{evidence_sets}}){ # print "Evidence: ",@{$self->{evidence_sets}{$keys}},"\n"; #} my ($clusters, $non_clusters) = cluster_Genes( $matching_extended_cdnas, $self->{evidence_sets} ) ; #store genes seperated by strand #USAGE?? my $genes_by_strand; my @possible_transcripts; foreach my $gene (@$matching_extended_cdnas){ push @{$genes_by_strand->{$gene->strand}}, $gene; } #apply collapsing method from TranscriptConsensus (creates the scores) my @cluster_to_use = (); if(scalar @$clusters){ push(@cluster_to_use, @$clusters); } if(scalar @$non_clusters){ push(@cluster_to_use, @$non_clusters); } foreach my $cluster (@cluster_to_use){ print STDERR "CLUSTER: ".$cluster->start." ".$cluster->end." ".$cluster->strand."\n" if $self->VERBOSE; my $collapsed_cluster = Bio::EnsEMBL::Analysis::Runnable::TranscriptConsensus->collapse_cluster($cluster, $genes_by_strand); my @potential_genes = $cluster->get_Genes; foreach my $potential_gene (@potential_genes){ foreach my $potential_trans (@{$potential_gene->get_all_Transcripts}) { #create an extended transcript (with scores) my @exon_array = ( @{$potential_trans->get_all_Exons} ); my $new_transcript = Bio::EnsEMBL::Analysis::Runnable::TranscriptConsensus->transcript_from_exons( \@exon_array, undef, $collapsed_cluster); #add supporting features $new_transcript->add_supporting_features(@{$potential_trans->get_all_supporting_features}); #add a ditag score print STDERR "new_transcript score :".$new_transcript->score()."\n" if $self->VERBOSE; my ($ditag_score, $index) = $self->score_ditags($self->ditags, $new_transcript, 0); $new_transcript->score($new_transcript->score() + $ditag_score); #add a UTR length score my $UTR_score = $self->calculate_UTR_score($cds_tran[0], $new_transcript); $new_transcript->score($new_transcript->score() + $UTR_score); push(@possible_transcripts, $new_transcript); } } } #clusters #get the highest scoring transcript, that survives tests @possible_transcripts = sort { $a->score <=> $b->score } @possible_transcripts if @possible_transcripts; my ($cdnaname, $genename); my $round = 0; POS: while(my $chosen_transcript = pop @possible_transcripts){ my $chosen_feats = $chosen_transcript->get_all_supporting_features; #make a gene from this $cdna_match = Bio::EnsEMBL::Gene->new; $cdna_match->slice($chosen_transcript->slice); $cdna_match->add_Transcript($chosen_transcript); $cdna_match->analysis($chosen_transcript->analysis); $cdna_match->get_all_Transcripts->[0]->add_supporting_features(@$chosen_feats); #combine it with the cds gene $combined_transcript = $self->combine_genes($cds, $cdna_match); # just check combined transcript works before throwing away the original transcript if ( defined($combined_transcript)){ $combined_transcript->sort; if(is_Transcript_sane($combined_transcript) && all_exons_are_valid($combined_transcript, $self->MAX_EXON_LENGTH, 1) && intron_lengths_all_less_than_maximum($combined_transcript, $self->MAX_INTRON_LENGTH) && has_no_unwanted_evidence($combined_transcript) ){ # make sure combined transcript doesn't misjoin any genewise clusters if($self->find_cluster_joiners($combined_transcript)){ print STDERR "Found a cluster_joiner!\n" if $self->VERBOSE; print STDERR $combined_transcript->seq_region_start."/".$combined_transcript->seq_region_end."\n" if $self->VERBOSE; #dont use this one $combined_transcript = undef; } } else{ print STDERR "Didn't pass checks.\n" if $self->VERBOSE; $combined_transcript = undef; } } else{ $combined_transcript = undef; } if(defined $combined_transcript){ #leave the loop last POS; } } } #find match if (defined $combined_transcript){ #transfer evidence $combined_transcript = $self->_transfer_evidence($combined_transcript, $cdna_match); # #make the new gene with UTR # my $genetype; # if($usingKnown){ # $genetype = $self->KNOWN_UTR_GENETYPE; # } else{ # if ( $self->EXTEND_ORIGINAL_BIOTYPE && (length($self->EXTEND_ORIGINAL_BIOTYPE)>0)) { # my $sep = ""; # $sep = "_" unless ( $self->UTR_GENETYPE =~ m/^_/ ); # $genetype = $combined_transcript->biotype. $sep. $self->UTR_GENETYPE ; # } else { # $genetype = $combined_genetype; # } # } #set biotype my $genetype; if ( $self->EXTEND_ORIGINAL_BIOTYPE && (length($self->EXTEND_ORIGINAL_BIOTYPE)>0)) { my $sep = ""; if($usingKnown){ $sep = "_" unless ( $self->KNOWN_UTR_GENETYPE =~ m/^_/ ); $genetype = $combined_transcript->biotype. $sep. $self->KNOWN_UTR_GENETYPE; } else{ $sep = "_" unless ( $self->UTR_GENETYPE =~ m/^_/ ); $genetype = $combined_transcript->biotype. $sep. $self->UTR_GENETYPE; } } else { if($usingKnown){ $genetype = $self->KNOWN_UTR_GENETYPE; } else{ $genetype = $combined_genetype; } } #print STDERR "MAKING_GENE FROM ".." AND ".$cdna_match->hit_name."\n"; my $combined_genes = $self->make_gene($genetype, $combined_transcript); #check phases, etc.- if it is not a blessed gene my $combined_gene; if(!$blessed){ $combined_gene = $self->look_for_both($combined_genes->[0]); my $combined_transcript = $combined_gene->get_all_Transcripts->[0]; $combined_transcript->sort; if(! (is_Transcript_sane($combined_transcript) && all_exons_are_valid($combined_transcript, $self->MAX_EXON_LENGTH, 1) && intron_lengths_all_less_than_maximum($combined_transcript, $self->MAX_INTRON_LENGTH) && validate_Translation_coords($combined_transcript, 1) && !contains_internal_stops($combined_transcript) && $combined_transcript->translate && has_no_unwanted_evidence($combined_transcript) ) ){ #revert to previous version if tests failed $combined_gene = $combined_genes->[0]; } } else{ #use blessed genes without modifications $combined_gene = $combined_genes->[0]; } #store as combined $self->combined_genes($combined_gene); print STDERR "RESULT-CombinedGene: ".$combined_gene->seq_region_name." ".$combined_gene->seq_region_start. "-".$combined_gene->seq_region_end."\n"; } else{ #retrieving unmerged if no UTR could be defined my $unmerged_cds = $self->retrieve_unmerged_gene($cds); $self->unmatched_genes($unmerged_cds); print STDERR "RESULT-UnmergedGene: ".$unmerged_cds->seq_region_name." ".$unmerged_cds->seq_region_start. "-".$unmerged_cds->seq_region_end."\n"; next CDS; } } print STDERR "At the end of the matching I have ".(scalar @{$self->combined_genes})." combined_genes genes". " and ".(scalar @{$self->unmatched_genes})." unmatched_genes\n" if $self->VERBOSE; if(!$blessed && ($self->{'remove_redundant'})){ #remove redundatant models from the umatched group my $unique_unmatched_genes = $self->remove_redundant_models($self->combined_genes, $self->unmatched_genes); #replace all unmatched genes $self->{'_unmatched_genes'} = []; $self->unmatched_genes(@{$unique_unmatched_genes}); print STDERR "without redundant models, we have ".(scalar @{$self->combined_genes})." combined_genes genes". " and ".(scalar @{$self->unmatched_genes})." unmatched_genes\n" if $self->VERBOSE; } } =head2 calculate_UTR_score Arg [1] : genewise transcript Arg [2] : transcript adding UTR region Description: calculate a score favouring transcripts that add UTR on both sides Could also favour longer UTRs if desired Returntype : int score =cut sub calculate_UTR_score { my ($self, $gw_gene, $matching_gene) = @_; my $BONUS_FOR_BOTH_SIDES = 1; my $UTR_score = 0; if(($gw_gene->start - $matching_gene->start) && ($matching_gene->end - $gw_gene->end)){ $UTR_score += $BONUS_FOR_BOTH_SIDES; } #...to be extended if needed return $UTR_score; } =head2 score_ditags Arg [1] : array ref with mapped ditags to analyse Arg [2] : transcript to analyse Arg [3] : (optional) int index, where to start looking in the ditag-array Description: score ditags in the region of a exon, transcripts or est that might support it Give a higher score to those features that match well to ditags and have a high tag_count. Returntype : int score: high value indicates many/well matching ditags int index: updated array index Exceptions : none =cut sub score_ditags{ my ($self, $ditags, $feature, $index) = @_; my $ditag_score = 0; if(!$index){ $index = 0 } #check start & end seperately, use best matching for(my $i = $index; $i < (scalar @$ditags); $i++){ my $ditag = $ditags->[$i]; if(($ditag->{'start'} < $feature->end) && ($ditag->{'end'} > $feature->start)){ my $start_distance = abs($ditag->{'start'} - $feature->start); my $end_distance = abs($ditag->{'end'} - $feature->end); if(($start_distance < $self->DITAG_WINDOW) && ($end_distance < $self->DITAG_WINDOW)){ #matching ditag; produce a score favoring those transcripts, #that have a high ditag count &/| perfectly positioned ditags. #magic equation to generate a score my $score = 0; $score += ($ditag->{'tag_count'} / (($start_distance +1) * 0.9)); $score += ($ditag->{'tag_count'} / (($end_distance +1) * 0.9)); if($score > 0){ $ditag_score += $score; } } } if($ditag->{'end'} < $feature->start){ $index++; } elsif($ditag->{'start'} > $feature->start){ last; } } print STDERR " returning ditagscore $ditag_score.\n" if $self->VERBOSE; return($ditag_score, $index); } =head2 check_for_predefined_pairing Arg [1] : Bio::EnsEMBL::Gene CDS gene Arg [2] : Boolean to indicate blessed gene Description: check whether there is a specific cDNA assigned to the given gene model Returntype : Bio::EnsEMBL::Gene cDNA gene =cut sub check_for_predefined_pairing { my ($self, $gene, $blessed) = @_; my $cdna = undef; my $protein_id; my $cdna_id; # get hash with cDNAs: $cdna_evidence{$evidence->hseqname()} = $cdna-gene; my $cdna_evidence = $self->_cdna_evidence(); # get the protein id, the gene was built from; # do this for blessed genes also, as this should be most reliable EXON: foreach my $exon(@{$gene->get_all_Exons}){ my @feat = @{$exon->get_all_supporting_features}; foreach my $feat(@feat){ $protein_id = $feat->hseqname; last EXON if (defined $protein_id); } } if ((!defined $protein_id || $protein_id eq '') && $blessed){ #try to use RefSeq xrefs for the blessed genes #The CCDS genes are stored with NM-entries in the xref table, #so these can be used directly as cDNA-ids. #Might have to be adjusted if this changes. $protein_id = 'blessed'; my @xrefs = @{$gene->get_all_DBLinks()}; if(scalar @xrefs){ foreach my $xref (@xrefs){ if($xref->display_id =~ "^NM_"){ $cdna_id = $xref->display_id; $cdna_id =~ s/\.\S+$//; print STDERR "have xref $cdna_id\n" if $self->VERBOSE; last; } } } else{ print STDERR "no suitable xref\n" if $self->VERBOSE; } } else{ if (!defined $protein_id || $protein_id eq ''){ print STDERR "no protein_id for gene.\n" if $self->VERBOSE; return undef; } # Using the protein id of the targetted gene, determine the # corresponding cDNA id from the pre-loaded hash. #remove version info $protein_id =~ s/\.\S+$//; $cdna_id = $self->get_cdna_id_from_protein_id($protein_id); } if (!defined $cdna_id || $cdna_id eq ''){ print STDERR "no predefined cDNA found.\n" if $self->VERBOSE; return undef; } else{ print STDERR "found predefined cDNA $cdna_id for $protein_id\n" if $self->VERBOSE; } #make sure it's not on the kill list if(defined ($self->kill_list()->{$cdna_id})){ print STDERR "skipping " . $cdna_id . " present in kill list\n"; return undef; } #get the gene for this cDNA my $cdna_gene = $cdna_evidence->{$cdna_id}; #check if they really are overlapping to avoid disappointment when combining if($cdna_gene){ if($self->VERBOSE){ print STDERR "found cdna gene $cdna_gene / $cdna_id.\n" } if( !(($gene->seq_region_name eq $cdna_gene->seq_region_name) && ($gene->strand == $cdna_gene->strand) && ($gene->seq_region_start < $cdna_gene->seq_region_end) && ($cdna_gene->seq_region_start < $gene->seq_region_end) && ($gene->get_all_Transcripts->[0]->get_all_Exons->[0]->strand == $cdna_gene->get_all_Transcripts->[0]->get_all_Exons->[0]->strand)) ){ print STDERR "not overlapping properly, not using." if $self->VERBOSE; $cdna_gene = undef; } } else { print STDERR "Unable to fetch cDNA gene [$cdna_id].\n"; return undef; } return $cdna_gene; } =head2 get_cdna_id_from_protein_id Arg [1] : String protein-id Description: Find the cDNA-id a given protein-id is linked to as a "known" evidence Returntype : String cDNA-id =cut sub get_cdna_id_from_protein_id { my ($self, $protein_id) = @_; my $cdna_id = undef; if(defined($self->_known_pairs()->{$protein_id})){ $cdna_id = $self->_known_pairs()->{$protein_id}; } return $cdna_id; } =head2 create_predefined_pairing Arg [1] : none Description: Read GenBank sequence file linking NP proteins to a specific NM cDNA to use as "known" UTR evidence Returns : nothing =cut sub create_predefined_pairing { my ($self, $known_utr_file) = @_; print STDERR "Parsing GenBank file $known_utr_file for KnowUTR pairing.\n" if $self->VERBOSE; open(REFSEQ, "<$known_utr_file") or die "Can't open ".$known_utr_file.": $@\n"; my $cdna_id; my $protein_id; my %predefined_pairing = (); while(<REFSEQ>){ #use only these 2 fields of the file: next unless /^VERSION|DBSOURCE/; if(/VERSION/){ next unless /(NP\S+)/; if(defined $protein_id){ die("previous protein_id [$protein_id] has not been cleared out\n"); } if(defined $cdna_id){ die("previous cdna_id [$cdna_id] has not been cleared out ...\n"); } $protein_id = $1; } if(/DBSOURCE/){ # don't want NCs or NGs, etc. if(/(NC\_\S+)|(NG\_\S+)|(XM\S+)|(AC\_\S+)/){ $cdna_id = undef; $protein_id = undef; next; } if(!defined $protein_id){ die("something very wrong - no protein_id for $_\n"); } if (defined $cdna_id){ die("previous cdna_id [$cdna_id] has not been cleared out ...\n"); } next unless /(NM\S+)/; $cdna_id = $1; #remove version info $cdna_id =~ s/\.\S+$//; $protein_id =~ s/\.\S+$//; if(exists($predefined_pairing{$protein_id})){ print STDERR $protein_id." allready exists!\n"; } $predefined_pairing{$protein_id} = $cdna_id; $cdna_id = undef; $protein_id = undef; } } close REFSEQ; $self->_known_pairs(\%predefined_pairing); } =head2 find_cluster_joiners Arg [1] : ref to array of GeneClusters Description: screens UTR modified transcripts for those that potentially misjoin genewise clusters, and remove them Returntype : ref to array of Bio::EnsEMBL::Gene =cut sub find_cluster_joiners{ my ($self, $transcript) = @_; my @clusters; my $transcript_start; my $transcript_end; my $overlaps_previous_cluster = 0; my $matching_clusters = 0; if ($transcript->get_all_Exons->[0]->strand == 1){ @clusters = @{$self->forward_genewise_clusters}; $transcript_start = $transcript->start_Exon->start; $transcript_end = $transcript->end_Exon->end; } else{ @clusters = @{$self->reverse_genewise_clusters}; $transcript_start = $transcript->end_Exon->start; $transcript_end = $transcript->start_Exon->end; } print STDERR "transcript spans $transcript_start - $transcript_end ". $transcript->seq_region_start." - ".$transcript->seq_region_end."\n". "clusters span: \n" if $self->VERBOSE; if(!scalar(@clusters)){ print STDERR "Odd, no clusters\n" if $self->VERBOSE; return 0; } #NEW CODE: CLUSTERING WITH EXON COORDINATES foreach my $cluster(@clusters){ if( _overlapping_genes($transcript, $cluster) ){ $matching_clusters++; } if($matching_clusters>1){ print STDERRR "transcript joins clusters - discard it\n" if $self->VERBOSE; return 1; } } return 0; } =head2 _overlapping_genes Arg [1] : Bio::EnsEMBL::Transcript Arg [2] : Bio::EnsEMBL::Analysis::Tools::Algorithms::GeneCluster Description: Check whether any of the exons of the transcript overlap any of the exons of the cluster. Returntype : boolean, true if overlap exists =cut sub _overlapping_genes { my ($gene1, $cluster) = @_; # quit if genes do not have genomic overlap if ($gene1->end < $cluster->start || $gene1->start > $cluster->end) { return 0; } # overlap check based on all (noncoding + coding) Exons foreach my $exon1 (@{$gene1->get_all_Exons}){ foreach my $exon2 (@{$cluster->get_all_Exons}){ if ( ($exon1->overlaps($exon2)) && ($exon1->strand == $exon2->strand) ){ return 1; } } } return 0; } =head2 _filter_cdnas Arg [1] : ref to array of Bio::EnsEMBL::Gene Arg [2] : flag whether these genes are EST genes rather than cDNA genes Description: This method checks the cDNA and EST genes translation is not checked as these genescome without translation Returntype : ref to array of filtered Bio::EnsEMBL::Gene =cut sub _filter_cdnas{ my ($self, $cdna_arrayref, $ests) = @_; my @newcdna; my %cdna_evidence; cDNA_GENE: foreach my $cdna (@{$cdna_arrayref}) { cDNA_TRANSCRIPT: foreach my $tran (@{$cdna->get_all_Transcripts}) { $tran->sort; if(!$ests){ #store cDNA gene evidence for later foreach my $evidence (@{ $tran->get_all_supporting_features }){ #remove version number?! my $evidence_name = $evidence->hseqname(); $evidence_name =~ s/\.\S+$//; #print STDERR "evidence: ".$evidence->hseqname()." => ".$evidence_name."\n"; $cdna_evidence{$evidence_name} = $cdna; } } # rejecting on basis of intron length may not be valid here # - it may not be that simple in the same way as it isn;t that simple in Targetted & Similarity builds next cDNA_TRANSCRIPT unless ( is_Transcript_sane($tran) && all_exons_are_valid($tran, $self->MAX_EXON_LENGTH, 1) && intron_lengths_all_less_than_maximum($tran, $self->MAX_INTRON_LENGTH) && has_no_unwanted_evidence($tran) ); push(@newcdna, $cdna); } } if(!$ests){ #store cDNAs for later use $self->_cdna_evidence(\%cdna_evidence); } return \@newcdna; } =head2 write_output Arg [1] : none Description: Do some cleaning up uf the result and store the genes to the database Returntype : none =cut sub write_output { my($self) = @_; # write genes in the output database my $gene_adaptor = $self->OUTPUT_DB->get_GeneAdaptor; print STDERR "Have ".scalar (@{$self->output})."(".$totalgenes.") genes to write\n"; GENE: foreach my $gene (@{$self->output}) { if(!$gene->analysis || $gene->analysis->logic_name ne $self->analysis->logic_name){ $gene->analysis($self->analysis); } # double check gene coordinates $gene->recalculate_coordinates; #As all genes/exons are stored as new in the target db #it's save to remove the old adaptor & old dbID here, #to avoid warnings from the store function. foreach my $exon (@{$gene->get_all_Exons}) { $exon->adaptor(undef); $exon->dbID(undef); } $gene->dbID(undef) ; $gene->adaptor(undef); for ( @{$gene->get_all_Transcripts} ) { for ( @{$_->get_all_supporting_features}){ $_->dbID(undef); $_->adaptor(undef); } $_->dbID(undef); $_->adaptor(undef); } eval { $gene_adaptor->store($gene); print STDERR "wrote gene dbID " . $gene->dbID . "\t".$gene->biotype . "\n" ; }; if( $@ ) { print STDERR "UNABLE TO WRITE GENE " . $gene->dbID. "of type " . $gene->type . "\n\n$@\n\nSkipping this gene\n"; } } } =head2 convert_to_extended_genes Arg [1] : ref to array of Bio::EnsEMBL::Gene Description: converts transcripts to Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::TranscriptExtended Returntype : ref to array of TranscriptExtended Exceptions : throws if gene has multiple transcripts =cut sub convert_to_extended_genes { my ($self, $genes) =@_; my @new_genes; for my $g (@$genes){ #my $st = _get_evidence_set( $g->biotype ); #my $set = $g->biotype; my ( $set ) = $self->_get_evidence_set ( $g->biotype ) ; for my $t (@{$g->get_all_Transcripts}){ throw ("gene has more than one trancript - only processing 1-gene-1-transcript-genes") if (@{$g->get_all_Transcripts}>1); # conversion my $gene_from_pt = Bio::EnsEMBL::Gene->new( -start => $t->start , -end => $t->end , -strand => $t->strand , -slice => $t->slice , -biotype => $g->biotype, -analysis => $t->analysis, ); my $new_tr = Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::TranscriptExtended->new( -BIOTYPE => $g->biotype , -ANALYSIS => $t->analysis , ); $new_tr->ev_set($set); $new_tr->add_supporting_features(@{$t->get_all_supporting_features}); my @pt_exons = @{$t->get_all_Exons} ; for (my $i=0 ; $i<scalar(@pt_exons) ; $i++) { # converting Bio::EnsEMBL::PredictionExon into ExonExtened (ISA Bio::EnsEMBL::Exon) my $pte = $pt_exons[$i]; bless $pte,"Bio::EnsEMBL::Analysis::Tools::GeneBuildUtils::ExonExtended"; $pte->biotype($g->biotype); $pte->ev_set($set); $pte->end_phase(0); $pte->phase(0); $pte->next_exon($pt_exons[$i+1]) ; if ($i == 0 ) { $pte->prev_exon(0); } else { $pte->prev_exon($pt_exons[$i-1]); } $pte->transcript($new_tr) ; $pte->analysis($t->analysis) ; }; # Extending the Bio::EnsEMBL::Transcript object by ev_set methods for (@pt_exons) { $new_tr->add_Exon($_); } $new_tr->sort; $gene_from_pt->add_Transcript($new_tr); push @new_genes , $gene_from_pt; } } return \@new_genes; } =head2 make_gene Arg [1] : string representing genetyoe to be associated with genes Arg [2] : array of Bio::EnsEMBL::Transcript Description: Constructs Bio::EnsEMBL::Gene objects from UTR-modified transcripts Returntype : none; new genes are stored in self->combined_genes Exceptions : throws when missing genetype or analysis =cut sub make_gene{ my ($self, $genetype, @transcripts) = @_; unless ( $genetype ){ throw("You must define UTR_GENETYPE in Bio::EnsEMBL::Analysis::Conf::GeneBuild::UTR_Builder"); } # an analysis should be passed in via the RunnableDB.m parent class: my $analysis = $self->analysis; unless ($analysis){ throw("You have to pass an analysis to this RunnableDB through new()"); } my @genes; my $count=0; foreach my $trans(@transcripts){ $trans->sort; unless ( is_Transcript_sane($trans) && intron_lengths_all_less_than_maximum($trans, $self->MAX_INTRON_LENGTH) && all_exons_are_valid($trans, $self->MAX_EXON_LENGTH, 1) && has_no_unwanted_evidence($trans) ){ print STDERR "\nrejecting transcript\n"; return; } my $gene = new Bio::EnsEMBL::Gene; $gene->biotype($genetype); $trans->biotype($genetype); $trans->analysis($analysis); $gene->add_Transcript($trans); $gene->analysis($analysis); # do not modify the analysis of the supporting features # they should be the original ones: cdna, targetted_genewise or similarity_genewise if($self->validate_gene($gene)){ push (@genes,$gene); $count++; } else{ print STDERR "\nCould not validate gene!\n"; } } return(\@genes); } =head2 match_protein_to_cdna Arg [1] : Bio::EnsEMBL::Gene Arg [2] : int is_est Description: this method tries to find the cdnas that can be merged with the genewise genes. Basically it checks for exact internal splice matching in the 5' and 3' exons of the genewise gene. In order to match a starting genewise exon with an cDNA exon, we need to have a. exactly coinciding exon boundaries b. either cdna exon has start <= genewise exon start, OR cdna exon start lies within $exon_slop bp of genewise exon start AND the cdna transcript will add exra UTR exons. Substitute "end" for "start" for 3prime ends of transcripts BUT do not allow through any cDNA that will result just in a shortened peptide without additional UTR exons. Returntype : ref to array of Bio::EnsEMBL::Gene, ref to hash relating Bio::EnsEMBL::Gene to UTR length Exceptions : =cut sub match_protein_to_cdna{ my ($self, $gw, $is_est) = @_; print STDERR "Looking at ESTs.\n" if $is_est; my (%UTR_hash, %UTR_side_indicator_hash); my @other_genes; my @matching_e2g; my @gw_tran = @{$gw->get_all_Transcripts}; my @gw_exons = @{$gw_tran[0]->get_all_Exons}; my $strand = $gw_exons[0]->strand; if($gw_exons[$#gw_exons]->strand != $strand){ warn("first and last gw exons have different strands ". "- can't make a sensible combined gene\n with ".$gw_tran[0]->dbId ); return undef; } if ( @gw_exons ){ if ($strand == 1 ){ @gw_exons = sort { $a->start <=> $b->start } @gw_exons; } else{ @gw_exons = sort { $b->start <=> $a->start } @gw_exons; } } else{ warn("gw gene without exons: ".$gw->dbID.", skipping it"); return undef; } my $exon_slop = 20; my $cds_length = length($gw_tran[0]->translateable_seq()); my @genes; if($is_est){ @genes = @{$self->ests}; } else{ @genes = @{$self->cdna_genes}; } cDNA: foreach my $e2g(@genes){ my @egtran = @{$e2g->get_all_Transcripts}; my @eg_exons = @{$egtran[0]->get_all_Exons}; $strand = $eg_exons[0]->strand; if($eg_exons[$#eg_exons]->strand != $strand){ warn("first and last e2g exons have different strands - skipping transcript ".$egtran[0]->dbID); next cDNA; } if ($strand == 1 ){ @eg_exons = sort { $a->start <=> $b->start } @eg_exons; } else{ @eg_exons = sort { $b->start <=> $a->start } @eg_exons; } my $fiveprime_match = 0; my $threeprime_match = 0; my $left_exon; my $right_exon; my $left_diff = 0; my $right_diff = 0; # Lets deal with single exon genes first if ($#gw_exons == 0) { foreach my $current_exon (@eg_exons) { if($current_exon->strand != $gw_exons[0]->strand){ next cDNA; } # don't yet deal with genewise leakage for single exon genes if ($gw_exons[0]->end <= $current_exon->end && $gw_exons[0]->start >= $current_exon->start){ $fiveprime_match = 1; $threeprime_match = 1; $left_exon = $current_exon; $right_exon = $current_exon; $left_diff = $gw_exons[0]->start - $current_exon->start; $right_diff = $current_exon->end - $gw_exons[0]->end; } } } else { # Now the multi exon genewises cDNA_EXONS: foreach my $current_exon (@eg_exons) { if($current_exon->strand != $gw_exons[0]->strand){ next cDNA; } if($gw_exons[0]->strand == 1){ #FORWARD: # 5prime if ($gw_exons[0]->end == $current_exon->end && # either e2g exon starts before genewise exon ($current_exon->start <= $gw_exons[0]->start || # or e2g exon is a bit shorter but there are spliced UTR exons as well (abs($current_exon->start - $gw_exons[0]->start) <= $exon_slop && $current_exon != $eg_exons[0]))){ $fiveprime_match = 1; $left_exon = $current_exon; $left_diff = $gw_exons[0]->start - $current_exon->start; } # 3prime elsif($gw_exons[$#gw_exons]->start == $current_exon->start && # either e2g exon ends after genewise exon ($current_exon->end >= $gw_exons[$#gw_exons]->end || # or there are UTR exons to be added (abs($current_exon->end - $gw_exons[$#gw_exons]->end) <= $exon_slop && $current_exon != $eg_exons[$#eg_exons]))){ $threeprime_match = 1; $right_exon = $current_exon; $right_diff = $current_exon->end - $gw_exons[0]->end; } } elsif($gw_exons[0]->strand == -1){ #REVERSE: # 5prime if ($gw_exons[0]->start == $current_exon->start && # either e2g exon ends after gw exon ($current_exon->end >= $gw_exons[0]->end || # or there are UTR exons to be added (abs($current_exon->end - $gw_exons[0]->end) <= $exon_slop && $current_exon != $eg_exons[0]))){ #print STDERR "fiveprime reverse match\n"; $fiveprime_match = 1; $right_exon = $current_exon; $right_diff = $current_exon->end - $gw_exons[0]->end; } #3prime elsif ($gw_exons[$#gw_exons]->end == $current_exon->end && # either e2g exon starts before gw exon ($current_exon->start <= $gw_exons[$#gw_exons]->start || # or there are UTR exons to be added (abs($current_exon->start - $gw_exons[$#gw_exons]->start) <= $exon_slop && $current_exon != $eg_exons[$#eg_exons]))){ #print STDERR "threeprime reverse match\n"; $threeprime_match = 1; $left_exon = $current_exon; $left_diff = $gw_exons[0]->start - $current_exon->start; } } } } # can match either end, or both if($fiveprime_match || $threeprime_match){ my ($UTR_length, $left_UTR_length, $right_UTR_length) = $self->_compute_UTRlength($egtran[0], $left_exon, $left_diff, $right_exon, $right_diff); my $UTR_diff = $egtran[0]->length; #make sure CDS is not much smaller than UTR if(($cds_length * 10) > $UTR_diff){ $UTR_hash{$e2g} = $UTR_length; $UTR_side_indicator_hash{$e2g} = 1; print STDERR "considering cDNA ".$e2g->seq_region_start."-".$e2g->seq_region_end."[".($cds_length * 10)." > ".$UTR_diff."]\n"; push(@matching_e2g, $e2g); } else{ print STDERR "didnt pass UTR length check [".($cds_length * 10)." > ".$UTR_diff."]: ".$e2g->seq_region_start."-".$e2g->seq_region_end."\n" if $self->VERBOSE; } } } print STDERR "returning ".(scalar @matching_e2g)." candidates.\n"; return (\@matching_e2g,\%UTR_hash, \%UTR_side_indicator_hash); } =head2 _compute_UTRlength Arg [1] : eg-transcript Arg [2] : first matching exon Arg [3] : (int) UTR bases of first matching exon Arg [4] : last matching exon Arg [5] : (int) UTR bases of first matching exon Description: add up genomic extend of UTR regions Returntype : int (basepairs) total UTR-length, left UTR-length, right UTR-lenght Exceptions : none =cut sub _compute_UTRlength{ my ($self, $transcript, $left_exon, $left_diff, $right_exon, $right_diff) = @_; my $strand = $transcript->start_Exon->strand; my @exons = sort { $a->start <=> $b->start } @{ $transcript->get_all_Exons }; my $UTRlength = 0; my $in_UTR = 1; my $start_flag = 0; my $left_UTR = $left_diff; my $right_UTR = $right_diff; foreach my $exon ( @exons ){ if ( defined $left_exon && $exon == $left_exon ){ $UTRlength += $left_diff; $in_UTR = 0; } elsif( defined $right_exon && $exon == $right_exon ){ $UTRlength += $right_diff; $in_UTR = 1; } elsif( $in_UTR == 1 ){ $UTRlength += $exon->length; if(!$start_flag){ $left_UTR += $UTRlength; $start_flag = 1; } else{ $right_UTR += $UTRlength; } } } return($UTRlength, $left_UTR, $right_UTR); } =head2 _merge_genes Arg [1] : ref to array of Bio::EnsEMBL::Gene Description: merges adjacent exons if they are frameshifted; stores component exons as subSeqFeatures of the merged exon Returntype : ref to array of Bio::EnsEMBL::Gene =cut sub _merge_genes { my ($self, $genesref) = @_; my @merged; my $count = 1; UNMERGED_GENE: foreach my $unmerged (@{$genesref}){ my $gene = new Bio::EnsEMBL::Gene; $gene->dbID($unmerged->dbID); my @pred_exons; my $ecount = 0; # order is crucial my @trans = @{$unmerged->get_all_Transcripts}; if(scalar(@trans) != 1) { throw("Gene with dbID " . $unmerged->dbID . " has NO or more than one related transcript, where a 1-gene-to-1-transcript-relation is assumed.". " Check preceding analysis \n"); } $trans[0]->sort; # check the sanity of the transcript if( !is_Transcript_sane($trans[0]) || !all_exons_are_valid($trans[0], $self->MAX_EXON_LENGTH, 1) || !has_no_unwanted_evidence($trans[0]) || !intron_lengths_all_less_than_maximum($trans[0], $self->MAX_INTRON_LENGTH) ){ print STDERR "transcript ".$unmerged->dbID."did NOT pass sanity check!\n"; print STDERR "found at ".$unmerged->seq_region_name.", ".$unmerged->seq_region_start.", ".$unmerged->seq_region_end."\n"; next UNMERGED_GENE; } ### we follow here 5' -> 3' orientation ### # $trans[0]->sort; my $cloned_translation = new Bio::EnsEMBL::Translation; my @unmerged_exons = @{$trans[0]->get_all_Exons}; my $strand = $unmerged_exons[0]->strand; my $previous_exon; EXON: foreach my $exon(@unmerged_exons){ ## genewise frameshift? we merge here two exons separated by max 10 bases into a single exon #if ($ecount && $pred_exons[$ecount-1]){ # $previous_exon = $pred_exons[$ecount-1]; #} $ecount++; my $separation = 0; my $merge_it = 0; ## we put every exon following a frameshift into the first exon before the frameshift ## following the ordering 5' -> 3' if (defined($previous_exon)){ #print STDERR "previous exon: ".$previous_exon->start."-".$previous_exon->end."\n"; #print STDERR "current exon : ".$exon->start."-".$exon->end."\n"; if ($strand == 1){ $separation = $exon->start - $previous_exon->end - 1; } elsif( $strand == -1 ){ $separation = $previous_exon->end - $exon->start - 1; } if ($separation <=10){ $merge_it = 1; } } if ( defined($previous_exon) && $merge_it == 1){ # combine the two # the first exon (5'->3' orientation always) is the containing exon, # which gets expanded and the other exons are added into it #print STDERR "merging $exon into $previous_exon\n"; #print STDERR $exon->start."-".$exon->end." into ".$previous_exon->start."-".$previous_exon->end."\n"; if ($strand == 1) { $previous_exon->end($exon->end); } else { $previous_exon->start($exon->start); } $previous_exon->add_sub_SeqFeature($exon,''); # if this is end of translation, keep that inf: if ( $exon == $trans[0]->translation->end_Exon ){ $cloned_translation->end_Exon( $previous_exon ); $cloned_translation->end($trans[0]->translation->end); } my %evidence_hash; foreach my $sf( @{$exon->get_all_supporting_features}){ if ( $evidence_hash{$sf->hseqname}{$sf->hstart}{$sf->hend}{$sf->start}{$sf->end} ){ next; } #print STDERR $sf->start."-".$sf->end." ".$sf->hstart."-".$sf->hend." ".$sf->hseqname."\n"; $evidence_hash{$sf->hseqname}{$sf->hstart}{$sf->hend}{$sf->start}{$sf->end} = 1; $previous_exon->add_supporting_features($sf); } next EXON; } else{ # make a new Exon - clone using Bio::EnsEMBL::Analysis::Tools::GeneB uildUtils::ExonUtils my $cloned_exon = clone_Exon($exon); $cloned_exon->add_sub_SeqFeature($exon,''); # if this is start/end of translation, keep that info: if ( $exon == $trans[0]->translation->start_Exon ){ $cloned_translation->start_Exon( $cloned_exon ); $cloned_translation->start($trans[0]->translation->start); } if ( $exon == $trans[0]->translation->end_Exon ){ $cloned_translation->end_Exon( $cloned_exon ); $cloned_translation->end($trans[0]->translation->end); } $previous_exon = $cloned_exon; push(@pred_exons, $cloned_exon); } } # transcript my $merged_transcript = new Bio::EnsEMBL::Transcript; $merged_transcript->dbID($trans[0]->dbID); foreach my $pe(@pred_exons){ $merged_transcript->add_Exon($pe); } $merged_transcript->sort; $merged_transcript->translation($cloned_translation); my @seqeds = @{$trans[0]->translation->get_all_SeqEdits}; if (scalar(@seqeds)) { print "Copying sequence edits\n" if $self->VERBOSE; foreach my $se (@seqeds) { my $new_se = Bio::EnsEMBL::SeqEdit->new( -CODE => $se->code, -NAME => $se->name, -DESC => $se->description, -START => $se->start, -END => $se->end, -ALT_SEQ => $se->alt_seq ); my $attribute = $new_se->get_Attribute(); $cloned_translation->add_Attributes($attribute); } } my @support = @{$trans[0]->get_all_supporting_features}; if (scalar(@support)) { $merged_transcript->add_supporting_features(@support); } $merged_transcript->add_Attributes(@{$trans[0]->get_all_Attributes}); # and gene $gene->add_Transcript($merged_transcript); $gene->biotype($unmerged->biotype); push(@merged, $gene); $count++; # store match between merged and original gene so we can easily retrieve the latter if we need to $self->merged_unmerged_pairs($gene,$unmerged); } # end UNMERGED_GENE return @merged; } =head2 combine_genes Arg [1] : genewise gene Arg [2] : cDNA gene Description: combine gene with matching cDNA Returntype : Bio::EnsEMBL::Transcript =cut sub combine_genes{ my ($self, $gw, $e2g) = @_; my $modified_peptide = 0; my @combined_transcripts = (); # should be only 1 transcript my @gw_tran = @{$gw->get_all_Transcripts}; $gw_tran[0]->sort; my @gw_exons = @{$gw_tran[0]->get_all_Exons}; # ordered array of exons my @egtran = @{$e2g->get_all_Transcripts}; $egtran[0]->sort; my @e2g_exons = @{$egtran[0]->get_all_Exons}; # ordered array of exons # clone transcript my $newtranscript = new Bio::EnsEMBL::Transcript; if ( $self->EXTEND_ORIGINAL_BIOTYPE && (length($self->EXTEND_ORIGINAL_BIOTYPE)>0) ) { # original biotype of gene / transcript will only be extended, not overwritten $newtranscript->biotype($gw->biotype) ; } foreach my $exon(@gw_exons){ my $new_exon = clone_Exon($exon); $newtranscript->add_Exon($new_exon); } my @support = @{$gw_tran[0]->get_all_supporting_features}; if (scalar(@support)) { $newtranscript->add_supporting_features(@support); } $newtranscript->add_Attributes(@{$gw_tran[0]->get_all_Attributes}); my $translation = new Bio::EnsEMBL::Translation; $translation->start($gw_tran[0]->translation->start); $translation->end($gw_tran[0]->translation->end); $translation->start_Exon($gw_tran[0]->translation->start_Exon); $translation->end_Exon($gw_tran[0]->translation->end_Exon); my @seqeds = @{$gw_tran[0]->translation->get_all_SeqEdits}; if (scalar(@seqeds)) { print "Copying sequence edits\n" if $self->VERBOSE; foreach my $se (@seqeds) { my $new_se = Bio::EnsEMBL::SeqEdit->new( -CODE => $se->code, -NAME => $se->name, -DESC => $se->description, -START => $se->start, -END => $se->end, -ALT_SEQ => $se->alt_seq ); my $attribute = $new_se->get_Attribute(); $translation->add_Attributes($attribute); } } $newtranscript->translation($translation); $newtranscript->translation->start_Exon($newtranscript->start_Exon); $newtranscript->translation->end_Exon($newtranscript->end_Exon); my $eecount = 0; my $modified_peptide_flag = 0; EACH_E2G_EXON: foreach my $ee (@e2g_exons){ # check strands are consistent if ($ee->strand != $gw_exons[0]->strand){ warn("gw and e2g exons have different strands - can't combine genes\n") ; return undef; } # single exon genewise prediction? if(scalar(@gw_exons) == 1) { ($newtranscript, $modified_peptide_flag) = $self->transcript_from_single_exon_genewise( $ee, $gw_exons[0], $newtranscript, $translation, $eecount, @e2g_exons); } else { ($newtranscript, $modified_peptide_flag) = $self->transcript_from_multi_exon_genewise($ee, $newtranscript, $translation, $eecount, $gw, $e2g) } if ( $modified_peptide_flag ){ $modified_peptide = 1; } # increment the exon $eecount++; } # end of EACH_E2G_EXON #don't modify translation of blessed genes my $biotype = $gw->biotype; if(defined($newtranscript) && $modified_peptide && (($self->{'blessed_type'}) =~ m/$biotype/)){ print STDERR "translation of blessed gene would need to be modified - not using combined gene.\n"; return undef; } ############################## # expand merged exons ############################## # the new transcript is made from a merged genewise gene # check the transcript and expand frameshifts in all but original 3' gw_exon # (the sub_SeqFeatures have been flushed for this exon) if (defined($newtranscript)){ $newtranscript->sort; foreach my $ex (@{$newtranscript->get_all_Exons}){ if($ex->sub_SeqFeature && scalar($ex->sub_SeqFeature) > 1 ){ my @sf = sort {$a->start <=> $b->start} $ex->sub_SeqFeature; my $first = shift(@sf); $ex->end($first->end); # add back the remaining component exons foreach my $s(@sf){ $newtranscript->add_Exon($s); $newtranscript->sort; } # flush the sub_SeqFeatures $ex->flush_sub_SeqFeature; } } # check that the resulting transcript unless( is_Transcript_sane($newtranscript) && all_exons_are_valid($newtranscript, $self->MAX_EXON_LENGTH, 1) && intron_lengths_all_less_than_maximum($newtranscript, $self->MAX_INTRON_LENGTH) && has_no_unwanted_evidence($newtranscript) ){ print STDERR "problems with this combined transcript, return undef\n"; return undef; } #check the translation unless( validate_Translation_coords($newtranscript, 1) && validate_Translation_coords($newtranscript, 1) && !contains_internal_stops($newtranscript) && $newtranscript->translate){ print STDERR "problems with this combined translation, return undef\n"; return undef; } # check translation is the same as for the genewise gene we built from my $foundtrans = 0; # the genewise translation can be modified due to a disagreement in a # splice site with cdnas. This can happen as neither blast nor genewise can # always find very tiny exons. # we then recalculate the translation using genomewise: my $newtrans; if ( $modified_peptide ){ my $strand = $newtranscript->start_Exon->strand; $newtrans = $self->_recalculate_translation($newtranscript, $strand); unless( validate_Translation_coords($newtrans, 1) && is_Transcript_sane($newtrans) && all_exons_are_valid($newtrans, $self->MAX_EXON_LENGTH, 1) && intron_lengths_all_less_than_maximum($newtrans, $self->MAX_INTRON_LENGTH) ){ print STDERR "problems with this genomewise alternative model, returning original transript.\n"; $newtrans = $newtranscript } } else{ $newtrans = $newtranscript; } return $newtrans; } else{ warn("No combination could be built\n"); return undef; } } =head2 transcript_from_single_exon_genewise Arg [1] : Description: This method will actually do the combination of both cdna and genewise gene. Note that if there is a match on one end but not on the other, the code will extend one end, but will leave the other as it is in the genewise genes. This will explit cdna matches that look fine on one end and we disregard the mismatching part. Returntype : =cut sub transcript_from_single_exon_genewise { my ($self, $eg_exon, $gw_exon, $transcript, $translation, $exoncount, @e2g_exons) = @_; # save out current translation end - we will need this if we have to unmerge frameshifted exons later my $orig_tend = $translation->end; # stay with being strict about gw vs e2g coords - may change this later ... # the overlapping e2g exon must at least cover the entire gw_exon if ($gw_exon->start >= $eg_exon->start && $gw_exon->end <= $eg_exon->end){ my $egstart = $eg_exon->start; my $egend = $eg_exon->end; my $gwstart = $gw_exon->start; my $gwend = $gw_exon->end; # modify the coordinates of the first exon in $newtranscript my $ex = $transcript->start_Exon; $ex->start($eg_exon->start); $ex->end($eg_exon->end); # need to explicitly set the translation start & end exons here. $translation->start_Exon($ex); # end_exon may be adjusted by 3' coding exon frameshift expansion. Ouch. $translation->end_Exon($ex); # need to deal with translation start and end this time - varies depending on strand #FORWARD: if($gw_exon->strand == 1){ my $diff = $gwstart - $egstart; my $tstart = $translation->start; my $tend = $translation->end; #print STDERR "diff: ".$diff." translation start : ".$tstart." end: ".$tend."\n"; #print STDERR "setting new translation to start: ".($tstart+$diff)." end: ".($tend+$diff)."\n"; $translation->start($tstart + $diff); $translation->end($tend + $diff); if($translation->start < 0){ warn("Forward strand: setting very dodgy translation start: " . $translation->start. "\n"); return(undef, 1); } if($translation->end > $translation->end_Exon->length){ warn("Forward strand: setting dodgy translation end: " . $translation->end . " exon_length: " . $translation->end_Exon->length . "\n"); return(undef, 1); } } #REVERSE: elsif($gw_exon->strand == -1){ my $diff = $egend - $gwend; my $tstart = $translation->start; my $tend = $translation->end; $translation->start($tstart+$diff); $translation->end($tend + $diff); if($translation->start < 0){ warn("Forward strand: setting very dodgy translation start: " . $translation->start. "\n"); return(undef, 1); } if($translation->end > $translation->end_Exon->length){ warn("Forward strand: setting dodgy translation end: " . $translation->end . " exon_length: " . $translation->end_Exon->length . "\n"); return(undef, 1); } } # expand frameshifted single exon genewises back from one exon to multiple exons if(defined($ex->sub_SeqFeature) && (scalar($ex->sub_SeqFeature) > 1)){ #print STDERR "frameshift in a single exon genewise\n"; my @sf = $ex->sub_SeqFeature; # save current start and end of modified exon my $cstart = $ex->start; my $cend = $ex->end; my $exlength = $ex->length; # get first exon - this has same id as $ex my $first = shift(@sf); $ex->end($first->end); # NB end has changed!!! # don't touch start. # no need to modify translation start # get last exon my $last = pop(@sf); $last->end($cend); $transcript->add_Exon($last); # and adjust translation end - the end is still relative to the merged gw exon $translation->end_Exon($last); # put back the original end translation $translation->end($orig_tend); # get any remaining exons foreach my $s(@sf){ $transcript->add_Exon($s); } $transcript->sort; # flush the sub_SeqFeatures $ex->flush_sub_SeqFeature; } # need to add back exons, both 5' and 3' $self->add_5prime_exons($transcript, $exoncount, @e2g_exons); $self->add_3prime_exons($transcript, $exoncount, @e2g_exons); } return ($transcript,0); } =head2 transcript_from_multi_exon_genewise Arg [1] : Description: This method will actually do the combination of both cdna and genewise gene. Note that if there is a match on one end but not on the other, the code will extend one end, but will leave the other as it is in the genewise genes. This will explit cdna matches that look fine on one end and we disregard the mismatching part. Returntype : =cut sub transcript_from_multi_exon_genewise { my ($self, $current_exon, $transcript, $translation, $exoncount, $gw_gene, $eg_gene) = @_; # $current_exon is the exon one the e2g_transcript we are in at the moment # $exoncount is the position of the e2g exon in the array # save out current translation->end - we'll need it if we have to expand 3prime exon later # my $orig_tend = $translation->end; my @gwtran = @{$gw_gene->get_all_Transcripts}; $gwtran[0]->sort; my @gwexons = @{$gwtran[0]->get_all_Exons}; my @egtran = @{$eg_gene->get_all_Transcripts}; $egtran[0]->sort; my @egexons = @{$egtran[0]->get_all_Exons}; # in order to match a starting genewise exon with an e2g exon, we need to have # a. exactly coinciding exon ends # b. exon starts lying within $exon_slop bp of each other. # previously we had required e2g start to be strictly <= gw start, but this will lose us some valid UTRs # substitute "end" for "start" for 3' ends of transcripts # compare to the first genewise exon if($gwexons[0]->strand == 1){ return $self->transcript_from_multi_exon_genewise_forward($current_exon, $transcript, $translation, $exoncount, $gw_gene, $eg_gene); } elsif( $gwexons[0]->strand == -1 ){ return $self->transcript_from_multi_exon_genewise_reverse($current_exon, $transcript, $translation, $exoncount, $gw_gene, $eg_gene); } } =head2 transcript_from_multi_exon_genewise_forward Arg [1] : Description: Returntype : =cut sub transcript_from_multi_exon_genewise_forward { my ($self, $current_exon, $transcript, $translation, $exoncount, $gw_gene, $eg_gene) = @_; my $modified_peptide = 0; my @gwtran = @{$gw_gene->get_all_Transcripts}; $gwtran[0]->sort; my @gwexons = @{$gwtran[0]->get_all_Exons}; my @egtran = @{$eg_gene->get_all_Transcripts}; $egtran[0]->sort; my @egexons = @{$egtran[0]->get_all_Exons}; # save out current translation->end - we'll need it if we have to expand 3prime exon later my $orig_tend = $translation->end; my $exon_slop = 20; #5_PRIME: if (#they have a coincident end $gwexons[0]->end == $current_exon->end && # either e2g exon starts before genewise exon ($current_exon->start <= $gwexons[0]->start || # or e2g exon is a bit shorter but there are spliced UTR exons as well (abs($current_exon->start - $gwexons[0]->start) <= $exon_slop && $current_exon != $egexons[0]))){ my $current_start = $current_exon->start; my $gwstart = $gwexons[0]->start; # this exon will be the start of translation, convention: phase = -1 my $ex = $transcript->start_Exon; $ex->phase(-1); # modify the coordinates of the first exon in $newtranscript if # e2g is larger on this end than gw. if ( $current_exon->start < $gwexons[0]->start ){ $ex->start($current_exon->start); } elsif( $current_exon->start == $gwexons[0]->start ){ $ex->start($gwstart); $ex->phase($gwexons[0]->phase); } # if the e2g exon starts after the gw exon, # modify the start only if this e2g exon is not the first of the transcript elsif( $current_start > $gwstart && $exoncount != 0 ) { $ex->start($current_exon->start); } # add all the exons from the est2genome transcript, previous to this one transfer_supporting_evidence($current_exon, $ex); $self->add_5prime_exons($transcript, $exoncount, @egexons); # fix translation start if($gwstart >= $current_start){ # take what it was for the gw gene, and add on the extra my $tstart = $translation->start; print STDERR "Forward 5': original translation start: $tstart "; ## $tstart += ($gwstart - $current_start); $translation->start($tstart); print STDERR "re-setting translation start to: $tstart\n"; ## } # only trust a smaller cdna exon if it is not the first of the transcript # (it could be a truncated cdna) elsif($gwstart < $current_start && $exoncount != 0){ $modified_peptide = 1; #print STDERR "SHORTENING GENEWISE TRANSLATION\n"; # genewise has leaked over the start. Tougher call - we need to take into account the # frame here as well #print STDERR "gw exon starts: $gwstart < new start: $current_start\n"; #print STDERR "modifying exon, as cdna exon is not the first of transcript-> exoncount = $exoncount\n"; # $diff is the number of bases we chop from the genewise exon my $diff = $current_start - $gwstart; my $tstart = $translation->start; warn("this is a case where gw translation starts at $tstart > 1") if ($tstart>1); print STDERR "gw translation start: ".$tstart."\n"; #print STDERR "start_exon: ".$translation->start_Exon->start. #"-".$translation->start_Exon->end. #" length: ".($translation->start_Exon->end - $translation->start_Exon->start + 1). #" phase: ".$translation->start_Exon->phase. #" end_phase: ".$translation->start_Exon->end_phase."\n"; if($diff % 3 == 0) { # we chop exactily N codons from the beginning of translation $translation->start(1); } elsif ($diff % 3 == 1) { # we chop N codons plus one base $translation->start(3); } elsif ($diff % 3 == 2) { # we chop N codons plus 2 bases $translation->start(2); } else { $translation->start(1); warn("very odd - $diff mod 3 = " . $diff % 3 . "\n"); } } else{ #print STDERR "gw exon starts: $gwstart > new start: $current_start"; #print STDERR "but cdna exon is the first of transcript-> exoncount = $exoncount, so we don't modify it\n"; } throw("setting very dodgy translation start: " . $translation->start. "\n") unless $translation->start > 0; } # end 5' exon # 3_PRIME: elsif (# they have coincident start $gwexons[$#gwexons]->start == $current_exon->start && # either e2g exon ends after genewise exon ($current_exon->end >= $gwexons[$#gwexons]->end || # or we allow to end before if there are UTR exons to be added (abs($current_exon->end - $gwexons[$#gwexons]->end) <= $exon_slop && $current_exon != $egexons[$#egexons]))){ my $end_translation_shift = 0; # modify the coordinates of the last exon in $newtranscript # e2g is larger on this end than gw. my $ex = $transcript->end_Exon; # this exon is the end of translation, convention: end_phase = -1 $ex->end_phase(-1); if ( $current_exon->end > $gwexons[$#gwexons]->end ){ $ex->end($current_exon->end); } elsif( $current_exon->end == $gwexons[$#gwexons]->end ){ $ex->end($gwexons[$#gwexons]->end); $ex->end_phase($gwexons[$#gwexons]->end_phase); } # if the e2g exon ends before the gw exon, # modify the end only if this e2g exon is not the last of the transcript elsif ( $current_exon->end < $gwexons[$#gwexons]->end && $exoncount != $#egexons ){ $modified_peptide = 1; #print STDERR "SHORTENING GENEWISE TRANSLATION\n"; ## fix translation end iff genewise has leaked over - will need truncating my $diff = $gwexons[$#gwexons]->end - $current_exon->end; #print STDERR "diff: $diff\n"; my $tend = $translation->end; my $gw_exon_length = $gwexons[$#gwexons]->end - $gwexons[$#gwexons]->start + 1; my $cdna_exon_length = $current_exon->end - $current_exon->start + 1; #print STDERR "gw exon length : $gw_exon_length\n"; #print STDERR "cdna exon length: $cdna_exon_length\n"; my $length_diff = $gw_exon_length - $cdna_exon_length; #print STDERR "length diff: ".$length_diff."\n"; # should be == diff $ex->end($current_exon->end); if($diff % 3 == 0) { # we chop exactily N codons from the end of the translation # so it can end where the cdna exon ends $translation->end($cdna_exon_length); $end_translation_shift = $length_diff; } elsif ($diff % 3 == 1) { # we chop N codons plus one base # it should end on a full codon, so we need to end translation 2 bases earlier: $translation->end($cdna_exon_length - 2); $end_translation_shift = $length_diff + 2; } elsif ($diff % 3 == 2) { # we chop N codons plus 2 bases # it should end on a full codon, so we need to end translation 1 bases earlier: $translation->end($cdna_exon_length - 1); $end_translation_shift = $length_diff + 1; } else { # absolute genebuild paranoia 8-) $translation->end($cdna_exon_length); warn("very odd - $diff mod 3 = " . $diff % 3 . "\n"); } #print STDERR "Forward: translation end set to : ".$translation->end."\n"; } # need to explicitly set the translation end exon for translation to work out my $end_ex = $transcript->end_Exon; $translation->end_Exon($end_ex); # strand = 1 my $expanded = $self->expand_3prime_exon($ex, $transcript, 1); if($expanded){ # set translation end to what it originally was in the unmerged genewise gene # taking into account the diff #print STDERR "Forward: expanded 3' exon, re-setting end of translation from ".$translation->end." to orig_end ($orig_tend)- ( length_diff + shift_due_to_phases ) ($end_translation_shift)".($orig_tend - $end_translation_shift)."\n"; $translation->end($orig_tend - $end_translation_shift); } # finally add any 3 prime e2g exons transfer_supporting_evidence($current_exon, $ex); $self->add_3prime_exons($transcript, $exoncount, @egexons); } # end 3' exon return ($transcript,$modified_peptide); } =head2 transcript_from_multi_exon_genewise_reverse Arg [1] : Description: Returntype : Exceptions : Example : =cut sub transcript_from_multi_exon_genewise_reverse { my ($self, $current_exon, $transcript, $translation, $exoncount, $gw_gene, $eg_gene) = @_; my $modified_peptide = 0; my @gwtran = @{$gw_gene->get_all_Transcripts}; $gwtran[0]->sort; my @gwexons = @{$gwtran[0]->get_all_Exons}; my @egtran = @{$eg_gene->get_all_Transcripts}; $egtran[0]->sort; my @egexons = @{$egtran[0]->get_all_Exons}; # save out current translation->end - we'll need it if we have to expand 3prime exon later my $orig_tend = $translation->end; my $exon_slop = 20; # 5_PRIME: if ($gwexons[0]->start == $current_exon->start && # either e2g exon ends after gw exon ($current_exon->end >= $gwexons[0]->end || # or there are UTR exons to be added (abs($current_exon->end - $gwexons[0]->end) <= $exon_slop && $current_exon != $egexons[0]))){ # sort out translation start my $tstart = $translation->start; if($current_exon->end >= $gwexons[0]->end){ # take what it was for the gw gene, and add on the extra $tstart += $current_exon->end - $gwexons[0]->end; $translation->start($tstart); } elsif( $current_exon->end < $gwexons[0]->end && $current_exon != $egexons[0] ){ # genewise has leaked over the start. Tougher call - we need to take into account the # frame here as well $modified_peptide = 1; #print STDERR "SHORTENING GENEWISE TRANSLATION\n"; #print STDERR "In Reverse strand. gw exon ends: ".$gwexons[0]->end." > cdna exon end: ".$current_exon->end."\n"; #print STDERR "modifying exon, as cdna exon is not the first of transcript-> exoncount = $exoncount\n"; my $diff = $gwexons[0]->end - $current_exon->end; my $gwstart = $gwexons[0]->end; my $current_start = $current_exon->end; my $tstart = $translation->start; #print STDERR "start_exon: ".$translation->start_Exon->start. # "-".$translation->start_Exon->end. # " length: ".($translation->start_Exon->end - $translation->start_Exon->start + 1). # " phase: ".$translation->start_Exon->phase. # " end_phase: ".$translation->start_Exon->end_phase."\n"; if ($diff % 3 == 0) { $translation->start(1); } elsif ($diff % 3 == 1) { $translation->start(3); } elsif ($diff % 3 == 2) { $translation->start(2); } else { $translation->start(1); warn("very odd - $diff mod 3 = " . $diff % 3 . "\n");} } throw("setting very dodgy translation start: " . $translation->start. "\n") unless $translation->start > 0; # this exon is the start of translation, convention: phase = -1 my $ex = $transcript->start_Exon; $ex->phase(-1); # modify the coordinates of the first exon in $newtranscript if ( $current_exon->end > $gwexons[0]->end){ ## HERE WAS THE PROBLEM WHICH CHANGED THE SOURCE COORDINATES! ## $ex->end($current_exon->end); $ex->phase(-1); } elsif ( $current_exon->end == $gwexons[0]->end){ $ex->end($gwexons[0]->end); $ex->phase($gwexons[0]->phase); } elsif ( ($current_exon->end < $gwexons[0]->end) && ($current_exon != $egexons[0]) ){ $ex->end($current_exon->end); } # need to explicitly set the translation start exon for translation to work out $translation->start_Exon($ex); transfer_supporting_evidence($current_exon, $ex); $self->add_5prime_exons($transcript, $exoncount, @egexons); } # end 5' exon # 3_PRIME: elsif ($gwexons[$#gwexons]->end == $current_exon->end && # either e2g exon starts before gw exon ($current_exon->start <= $gwexons[$#gwexons]->start || # or there are UTR exons to be added (abs($current_exon->start - $gwexons[$#gwexons]->start) <= $exon_slop && $current_exon != $egexons[$#egexons]))){ my $end_translation_shift = 0; # this exon is the end of translation, convention: end_phase = -1 my $ex = $transcript->end_Exon; $ex->end_phase(-1); # modify the coordinates of the last exon in $newtranscript if ( $current_exon->start < $gwexons[$#gwexons]->start ){ # no need to modify translation->end as the 'end' of this exon has not changed $ex->start($current_exon->start); $ex->end_phase(-1); } elsif( $current_exon->start == $gwexons[$#gwexons]->start){ $ex->start($gwexons[$#gwexons]->start); $ex->end_phase($gwexons[$#gwexons]->end_phase); } # if the e2g exon starts after the gw exon, # modify the end only if this e2g exon is not the last of the transcript elsif ( $current_exon->start > $gwexons[$#gwexons]->start && $exoncount != $#egexons ){ $modified_peptide = 1; #print STDERR "SHORTENING GENEWISE TRANSLATION\n"; #print STDERR "In Reverse strand: gw exon start: ".$gwexons[$#gwexons]->start." < cdna exon start: ".$current_exon->start."\n"; #print STDERR "modifying exon, as cdna exon is not the last of transcript-> exoncount = $exoncount, and #egexons = $#egexons\n"; ## adjust translation my $diff = $current_exon->start - $gwexons[$#gwexons]->start; #print STDERR "diff: $diff\n"; my $tend = $translation->end; my $gw_exon_length = $gwexons[$#gwexons]->end - $gwexons[$#gwexons]->start + 1; my $cdna_exon_length = $current_exon->end - $current_exon->start + 1; #print STDERR "gw exon length : $gw_exon_length\n"; #print STDERR "cdna exon length: $cdna_exon_length\n"; my $length_diff = $gw_exon_length - $cdna_exon_length; # modify the combined exon coordinate to be that of the cdna $ex->start($current_exon->start); if($diff % 3 == 0) { # we chop exactily N codons from the end of the translation # so it can end where the cdna exon ends $translation->end($cdna_exon_length); $end_translation_shift = $length_diff; } elsif ($diff % 3 == 1) { # we chop N codons plus one base # it should end on a full codon, so we need to end translation 2 bases earlier: $translation->end($cdna_exon_length - 2); $end_translation_shift = $length_diff + 2; } elsif ($diff % 3 == 2) { # we chop N codons plus 2 bases # it should end on a full codon, so we need to end translation 1 bases earlier: $translation->end($cdna_exon_length - 1); $end_translation_shift = $length_diff + 1; } else { # absolute genebuild paranoia 8-) $translation->end($cdna_exon_length); warn("very odd - $diff mod 3 = " . $diff % 3 . "\n"); } } # strand = -1 my $expanded = $self->expand_3prime_exon($ex, $transcript, -1); # need to explicitly set the translation end exon for translation to work out my $end_ex = $transcript->end_Exon; $translation->end_Exon($end_ex); if($expanded){ # set translation end to what it originally was in the unmerged genewise gene #print STDERR "Reverse: expanded 3' exon, re-setting translation exon ".$translation->end. # " to original end( $orig_tend ) - shifts_due_to_phases_etc ( $end_translation_shift ) :". # ($orig_tend - $end_translation_shift)."\n"; $translation->end($orig_tend - $end_translation_shift); } transfer_supporting_evidence($current_exon, $ex); $self->add_3prime_exons($transcript, $exoncount, @egexons); } # end 3' exon return ($transcript, $modified_peptide); } =head2 add_5prime_exons Arg [1] : Description: Returntype : =cut sub add_5prime_exons { my ($self, $transcript, $exoncount, @e2g_exons) = @_; # add all the exons from the est2genome transcript, previous to this one # db handle will be screwed up, need to mak new exons from these my $c = 0; my $modified = 0; while($c < $exoncount){ my $newexon = new Bio::EnsEMBL::Exon; my $oldexon = $e2g_exons[$c]; $newexon->start($oldexon->start); $newexon->end($oldexon->end); $newexon->strand($oldexon->strand); # these are all 5prime UTR exons $newexon->phase(-1); $newexon->end_phase(-1); $newexon->slice($oldexon->slice); my %evidence_hash; #print STDERR "adding evidence at 5':\n"; foreach my $sf( @{$oldexon->get_all_supporting_features} ){ if ( $evidence_hash{$sf->hseqname}{$sf->hstart}{$sf->hend}{$sf->start}{$sf->end} ){ next; } $evidence_hash{$sf->hseqname}{$sf->hstart}{$sf->hend}{$sf->start}{$sf->end} = 1; #print STDERR $sf->start."-".$sf->end." ".$sf->hstart."-".$sf->hend." ".$sf->hseqname."\n"; $newexon->add_supporting_features($sf); } # print STDERR "Adding 5prime exon " . $newexon->start . " " . $newexon->end . "\n"; $transcript->add_Exon($newexon); $modified = 1; $c++; } if ($modified == 1){ $transcript->translation->start_Exon->phase(-1); } } # $exon is the terminal exon in the genewise transcript, $transcript. We need # to expand any frameshifts we merged in the terminal genewise exon. # The expansion is made by putting $exon to be the last (3' end) component, so we modify its # start but not its end. The rest of the components are added. The translation end will have to be modified, # this happens in the method _transcript_from_multi_exon.... =head2 expand_3prime_exon Arg [1] : Description: $exon is the terminal exon in the genewise transcript, $transcript. We need to expand any frameshifts we merged in the terminal genewise exon. The expansion is made by putting $exon to be the last (3 end) component, so we modify its start but not its end. The rest of the components are added. The translation end will have to be modified, this happens in the method _transcript_from_multi_exon.... Returntype : =cut sub expand_3prime_exon{ my ($self, $exon, $transcript, $strand) = @_; if(defined($exon->sub_SeqFeature) && (scalar($exon->sub_SeqFeature) > 1)){ #print STDERR "expanding 3'prime frameshifted exon $exon in strand $strand: ". #$exon->start."-".$exon->end." phase: ".$exon->phase." end_phase: ".$exon->end_phase."\n"; my @sf = $exon->sub_SeqFeature; my $last = pop(@sf); #print STDERR "last component: ".$last->start."-".$last->end." phase ".$last->phase." end_phase ".$last->end_phase."\n"; #print STDERR "setting exon $exon start: ".$last->start." phase: ".$last->phase."\n"; $exon->start($last->start); # but don't you dare touch the end! $exon->dbID($last->dbID); $exon->phase($last->phase); # add back the remaining component exons foreach my $s(@sf){ #print STDERR "adding exon: ".$s->start."-".$s->end."\n"; $transcript->add_Exon($s); $transcript->sort; } # flush the sub_SeqFeatures so we don't try to re-expand later $exon->flush_sub_SeqFeature; return 1; } # else, no expansion return 0; } =head2 add_3prime_exons Arg [1] : Description: $exoncount tells us which position in the array of e2g exons corresponds to the end of the genewise transcript so we add back exons 3 to that position. $exon and $transcript are references to Exon and Transcript objects. Returntype : =cut sub add_3prime_exons { my ($self, $transcript, $exoncount, @e2g_exons) = @_; # need to deal with frameshifts - 3' exon is a special case as its end might have changed # add all the exons from the est2genome transcript, subsequent to this one my $c = $#e2g_exons; my $modified = 0; while($c > $exoncount){ my $newexon = new Bio::EnsEMBL::Exon; my $oldexon = $e2g_exons[$c]; $newexon->start($oldexon->start); $newexon->end($oldexon->end); $newexon->strand($oldexon->strand); # these are all exons with UTR: $newexon->phase(-1); $newexon->end_phase(-1); $newexon->slice($oldexon->slice); #print STDERR "adding evidence in 3':\n"; my %evidence_hash; foreach my $sf( @{$oldexon->get_all_supporting_features }){ if ( $evidence_hash{$sf->hseqname}{$sf->hstart}{$sf->hend}{$sf->start}{$sf->end} ){ next; } $evidence_hash{$sf->hseqname}{$sf->hstart}{$sf->hend}{$sf->start}{$sf->end} = 1; #print STDERR $sf->start."-".$sf->end." ".$sf->hstart."-".$sf->hend." ".$sf->hseqname."\n"; $newexon->add_supporting_features($sf); } # print STDERR "Adding 3prime exon " . $newexon->start . " " . $newexon->end . "\n"; $transcript->add_Exon($newexon); $modified = 1; $c--; } if ($modified == 1){ $transcript->translation->end_Exon->end_phase(-1); } } =head2 remap_genes Description: strictly speaking this is no longer remapping anything... checks translation, set start/stop, can set biotype Returntype : ref to array of Bio::EnsEMBL::Gene =cut sub remap_genes { my ($self, $genes, $biotype_suffix) = @_; my @remapped_genes = (); my $blessed_type = $self->{'blessed_type'}; print STDERR "remapping ".scalar @$genes." genes.\n" if $self->VERBOSE; GENE: foreach my $gene (@$genes) { #leave the blessed genes alone my $biotype = $gene->biotype; #force a centain biotype? if($biotype_suffix && (length($biotype_suffix) > 0)){ $gene->biotype($gene->biotype().$biotype_suffix); } if(defined $blessed_type && $blessed_type =~ m/$biotype/){ print STDERR "not remapping ".$biotype."\n" if $self->VERBOSE; push(@remapped_genes, $gene); next GENE; } print STDERR "remapping ".$gene->biotype."\n" if $self->VERBOSE; my @t = @{$gene->get_all_Transcripts}; my $tran = $t[0]; # check that it translates unless(validate_Translation_coords($tran, 1) && !contains_internal_stops($tran) && $tran->translate){ print STDERR "\nERROR: Gene at ".$gene->seq_region_name." ".$gene->seq_region_start."-".$gene->seq_region_end." ". $gene->seq_region_strand." doesn't translate!\n\n"; push(@remapped_genes, $gene); ##added next GENE; } foreach my $transcript ( @{$gene->get_all_Transcripts} ){ if($biotype){ $transcript->biotype($gene->biotype); $transcript->analysis($gene->analysis); } # set start and stop codons set_start_codon($transcript); set_stop_codon($transcript); } push(@remapped_genes, $gene); } return \@remapped_genes; } =head2 validate_gene Arg [1] : Bio::EnsEMBL::Gene Description: checks start and end coordinates of each exon of each transcript are sane Returntype : 1 if gene is valid, otherwise zero =cut sub validate_gene { my ($self, $gene) = @_; # should be only a single transcript my @transcripts = @{$gene->get_all_Transcripts}; if(scalar(@transcripts) != 1) { print STDERR "Rejecting gene - should have one transcript, not " . scalar(@transcripts) . "\n"; return 0; } foreach my $transcript(@transcripts){ foreach my $exon(@{$transcript->get_all_Exons}){ unless( validate_Exon_coords($exon, 1) ){ print STDERR "Rejecting gene because of invalid exon\n"; return 0; } } } return 1; } =head2 _recalculate_translation Arg[1] : Bio::EnsEMBL::Transcript Arg[2] : strand of the transcript Description : a transcript is used as evidence for genomewise to recalculate the ORF. The idea is to use this when the peptide has been shortened, due to a genewise model being incompatible with the cdna splicing. This can happen when genewise cannot find very short exons and attaches them to one of the flanking exons. We tell genomewise to keep the splice boundaries pretty much static, so that we preserve the original splicing structure. Returntype : Bio::EnsEMBL::Transcript =cut sub _recalculate_translation { my ($self, $mytranscript, $strand) = @_; my $this_is_my_transcript = clone_Transcript($mytranscript); compute_translation($mytranscript); # check that everything is sane: unless(validate_Translation_coords($mytranscript, 1) && contains_internal_stops($mytranscript) && $mytranscript->translate){ print STDERR "problem with the translation. Returning the original transcript\n"; return $this_is_my_transcript; } return $mytranscript; } =head2 _transfer_evidence Arg [1] : reference to Bio::EnsEMBL::Tanscript $combined_transcript Arg [2] : reference to Bio::EnsEMBL::Transcript $cdna_transcript Description: transfers cdna evidence to combined transcript Returntype: Bio::EnsEMBL::Transcript =cut sub _transfer_evidence { my ($self, $combined_transcript, $cdna_transcript) = @_; my $first_support_id; foreach my $combined_exon(@{$combined_transcript->get_all_Exons}){ foreach my $cdna_exon(@{$cdna_transcript->get_all_Exons}){ # exact match or overlap? # exact match # if($combined_exon->start == $cdna_exon->start && # $combined_exon->end == $cdna_exon->end && # $combined_exon->strand == $cdna_exon->strand){ # print STDERR "exact match " . $combined_exon->dbID . " with " . $cdna_exon->dbID . "; transferring evidence\n"; # transfer_supporting_evidence($cdna_exon, $combined_exon); # } # overlap - feature boundaries may well be wonky if($combined_exon->overlaps($cdna_exon)){ if($combined_exon->strand != $cdna_exon->strand){ print STDERR "OVERLAPPING-BUT-DIFFERENT_STRANDS!\n"; } else{ transfer_supporting_evidence($cdna_exon, $combined_exon); } } } } my $cdna_trans = $cdna_transcript->get_all_Transcripts()->[0]; foreach my $tsf (@{$cdna_trans->get_all_supporting_features}) { print STDERR "adding supporting feature: ".$tsf->hseqname."\n" if $self->VERBOSE; $combined_transcript->add_supporting_features($tsf); } return $combined_transcript; } =head2 look_for_both Arg [1] : Bio:EnsEMBL:Gene Description: a copy of Steve's look_for_both-script, checks phases, etc. Returntype : Bio:EnsEMBL:Gene =cut sub look_for_both { my ($self, $gene) = @_; my $time = time; my $nupdated_start = 0; my $nupdated_end = 0; my $metcnt = 1; my $maxterdist = 150; foreach my $trans (@{$gene->get_all_Transcripts}) { if ($trans->translation) { my $tln = $trans->translation; my $coding_start = $trans->cdna_coding_start; my $orig_coding_start = $coding_start; $trans->sort; my $cdna_seq = uc($trans->spliced_seq); my @pepgencoords = $trans->pep2genomic(1,1); if(scalar(@pepgencoords) > 2) { print STDERR "pep start does not map cleanly\n"; goto TLNEND; # I swore I'd never use this - this code desperately needs a rewrite } my $pepgenstart = $pepgencoords[0]->start; my $pepgenend = $pepgencoords[$#pepgencoords]->end; unless($pepgencoords[0]->isa('Bio::EnsEMBL::Mapper::Coordinate')) { print STDERR "pep start maps to gap\n"; goto TLNEND; # I swore I'd never use this - this code desperately needs a rewrite } unless($pepgencoords[$#pepgencoords]->isa('Bio::EnsEMBL::Mapper::Coordinate')) { print STDERR "pep start (end of) maps to gap\n"; goto TLNEND; # I swore I'd never use this - this code desperately needs a rewrite } print STDERR "Pep genomic location = " . $pepgenstart . " " . $pepgenend . "\n" if $self->VERBOSE; my $startseq= substr($cdna_seq,$coding_start-1,3); print "cdna seq for pep start = " . $startseq . "\n"; if ($startseq ne "ATG") { if ($coding_start > 3) { my $had_stop = 0; while ($coding_start > 3 && !$had_stop) { my $testseq = substr($cdna_seq,$coding_start-4,3); if ($testseq eq "ATG") { print_Translation($trans) if $self->VERBOSE; my @coords = $trans->cdna2genomic($coding_start-3,$coding_start-1,$gene->strand); my $new_start; my $new_end; if(scalar(@coords) > 2) { print STDERR "Shouldn't happen - new start does not map cleanly - I'm out of here\n"; next; } elsif (scalar(@coords) == 2) { print STDERR "WOW ISN'T NATURE HORRIBLE: new start crosses intron\n"; print STDERR "coord[0] = " . $coords[0]->start . " " . $coords[0]->end ."\n"; print STDERR "coord[1] = " . $coords[1]->start . " " . $coords[1]->end ."\n"; if ($gene->strand == 1) { $new_start = $coords[0]->start; $new_end = $coords[$#coords]->end; } else { $new_start = $coords[0]->end; $new_end = $coords[$#coords]->start; } } else { $new_start = $coords[0]->start; $new_end = $coords[0]->end; } unless($coords[0]->isa('Bio::EnsEMBL::Mapper::Coordinate')) { print STDERR "Shouldn't happen - new start maps to gap - I'm out of here\n"; next; } unless($coords[$#coords]->isa('Bio::EnsEMBL::Mapper::Coordinate')) { print STDERR "Shouldn't happen - new start (end of) maps to gap - I'm out of here\n"; next; } print STDERR "genomic pos for new start = " . $new_start . " " . $new_end . "\n" if $self->VERBOSE; if ($new_end - $new_start == 2) { $nupdated_start++; my $newstartexon; foreach my $exon (@{$trans->get_all_Exons}) { if ($exon->end >= $new_start && $exon->start <= $new_start) { $newstartexon = $exon; last; } } if ($newstartexon == $tln->start_Exon) { if ($tln->start_Exon->strand == 1) { $tln->start($new_start - $tln->start_Exon->start + 1); } else { $tln->start($tln->start_Exon->end - $new_end + 1); } # NAUGHTY, but hey I should have to do this - I've changed the translation after all $trans->{'transcript_mapper'} = undef; $trans->{'coding_region_start'} = undef; $trans->{'coding_region_end'} = undef; $trans->{'cdna_coding_start'} = undef; $trans->{'cdna_coding_end'} = undef; } else { # find exon if (!defined($newstartexon)) { print STDERR "Failed finding new start exon - how can this be?\n"; next; } # create a copy of if and of current start exon (because of phase change) my $copyexon = new Bio::EnsEMBL::Exon( -start => $tln->start_Exon->start, -end => $tln->start_Exon->end, -strand => $gene->strand, ); my $copynewstartexon = new Bio::EnsEMBL::Exon( -start => $newstartexon->start, -end => $newstartexon->end, -strand => $gene->strand, ); # $copyexon->phase(0); $copyexon->end_phase($tln->start_Exon->end_phase); $copyexon->contig($tln->start_Exon->contig); if ($tln->start_Exon->stable_id) { $copyexon->stable_id($tln->start_Exon->stable_id . "MET" . $metcnt++); $copyexon->created($time); $copyexon->modified($time); $copyexon->version(1); } $copynewstartexon->phase($newstartexon->phase); # $copynewstartexon->end_phase(0); $copynewstartexon->contig($newstartexon->contig); if ($newstartexon->stable_id) { $copynewstartexon->stable_id($newstartexon->stable_id . "MET" . $metcnt++); $copynewstartexon->created($time); $copynewstartexon->modified($time); $copynewstartexon->version(1); } # TODO evidence if ($copynewstartexon->strand == 1) { $tln->start($new_start - $copynewstartexon->start + 1); } else { $tln->start($copynewstartexon->end - $new_end + 1); } # Replace exons in transcript, and fix phases my @newexons; my $inrange = 0; foreach my $exon (@{$trans->get_all_Exons}) { if ($inrange) { $exon->phase( $newexons[$#newexons]->end_phase ); $exon->end_phase(($exon->length + $exon->phase) % 3); } if ($exon == $tln->start_Exon) { $copyexon->phase( $newexons[$#newexons]->end_phase ); push @newexons,$copyexon; $inrange = 0; } elsif ($exon == $newstartexon) { push @newexons,$copynewstartexon; $copynewstartexon->end_phase(($exon->length - $tln->start + 1)%3); print STDERR "Setting end_phase on new start exon to " . $copynewstartexon->end_phase . " l = " . $exon->length . " ts = " . $tln->start . "\n" if $self->VERBOSE; $inrange = 1; } else { push @newexons,$exon; } } $trans->flush_Exons; foreach my $exon (@newexons) { $trans->add_Exon($exon); } # Reset translation start exon if ($tln->end_Exon == $tln->start_Exon) { $tln->end_Exon($copyexon); } $tln->start_Exon($copynewstartexon); } print_Translation($trans); } else { print STDERR "Across exons - not handling this\n"; } last; } else { if ($testseq =~ /TAA/ or $testseq =~ /TGA/ or $testseq =~ /TAG/) { $had_stop = 1; } else { $coding_start -= 3; } } } } else { print STDERR "Not enough bases upstream - NOT looking into genomic\n"if $self->VERBOSE; } } TLNEND: { my $coding_end = $trans->cdna_coding_end; my $orig_coding_end = $coding_end; $trans->sort; my $peplen = $trans->translate->length; my @pepgencoords = $trans->pep2genomic($peplen,$peplen); if(scalar(@pepgencoords) > 2) { print STDERR "pep end does not map cleanly\n"; next; } my $pepgenstart = $pepgencoords[0]->start; my $pepgenend = $pepgencoords[$#pepgencoords]->end; unless($pepgencoords[0]->isa('Bio::EnsEMBL::Mapper::Coordinate')) { print STDERR "pep end maps to gap\n"; next; } unless($pepgencoords[$#pepgencoords]->isa('Bio::EnsEMBL::Mapper::Coordinate')) { print STDERR "pep start (end of) maps to gap\n"; next; } #print "End Pep genomic location = " . $pepgenstart . " " . $pepgenend . "\n"; my $endseq= substr($cdna_seq,$coding_end-3,3); my $cdnalen = length($cdna_seq); #print "cdna seq for pep end = " . $endseq . "\n"; my $longendseq= substr($cdna_seq,$coding_end-6,12); #print "long end seq (-3 to len 12) = $longendseq\n"; # if (!($endseq ne "TGA" and $endseq ne "TAA" and $endseq ne "TAG")) { # print "Has end " . $trans->translateable_seq . "\n"; # } if ($endseq ne "TGA" and $endseq ne "TAA" and $endseq ne "TAG") { if (($cdnalen-$coding_end) > 3) { while (($cdnalen-$coding_end) > 3 && ($coding_end-$orig_coding_end) <= $maxterdist) { my $testseq = substr($cdna_seq,$coding_end,3); #print "Test seq = $testseq\n" if $self->VERBOSE ; if ($testseq eq "TGA" or $testseq eq "TAA" or $testseq eq "TAG") { my @coords = $trans->cdna2genomic($coding_end+1,$coding_end+3,$gene->strand); my $new_start; my $new_end; if(scalar(@coords) > 2) { throw("new end does not map cleanly\n"); next; } elsif (scalar(@coords) == 2) { print STDERR "WOW ISN'T NATURE HORRIBLE: new end crosses intron\n"; print STDERR "coord[0] = " . $coords[0]->start . " " . $coords[0]->end ."\n"; print STDERR "coord[1] = " . $coords[1]->start . " " . $coords[1]->end ."\n"; if ($gene->strand == 1) { $new_start = $coords[0]->start; $new_end = $coords[$#coords]->end; } else { $new_start = $coords[0]->end; $new_end = $coords[$#coords]->start; } } else { $new_start = $coords[0]->start; $new_end = $coords[0]->end; } unless($coords[0]->isa('Bio::EnsEMBL::Mapper::Coordinate')) { print STDERR "new start maps to gap\n"; next; } unless($coords[$#coords]->isa('Bio::EnsEMBL::Mapper::Coordinate')) { print STDERR "new start (end of) maps to gap\n"; next; } if ($new_end - $new_start == 2) { #print "Sequence of genomic pos of new end = " . $slice->subseq($new_start,$new_end,$gene->strand) . "\n"; $nupdated_end++; my $newendexon; foreach my $exon (@{$trans->get_all_Exons}) { if ($exon->end >= $new_start && $exon->start <= $new_start) { $newendexon = $exon; last; } } if ($newendexon == $tln->end_Exon) { if ($tln->end_Exon->strand == 1) { $tln->end($new_end - $tln->end_Exon->start + 1); } else { $tln->end($tln->end_Exon->end - $new_start + 1); } # NAUGHTY, but hey I should have to do this - I've changed the translation after all $trans->{'transcript_mapper'} = undef; $trans->{'coding_region_start'} = undef; $trans->{'coding_region_end'} = undef; $trans->{'cdna_coding_start'} = undef; $trans->{'cdna_coding_end'} = undef; } else { # find exon if (!defined($newendexon)) { print STDERR "Failed finding new end exon - how can this be?\n"; next; } # create a copy of if and of current end exon (because of phase change) my $copyexon = new Bio::EnsEMBL::Exon( -start => $tln->end_Exon->start, -end => $tln->end_Exon->end, -strand => $gene->strand, ); my $copynewendexon = new Bio::EnsEMBL::Exon( -start => $newendexon->start, -end => $newendexon->end, -strand => $gene->strand, ); $copyexon->phase($tln->end_Exon->phase); $copyexon->end_phase($tln->end_Exon->end_phase); $copyexon->contig($tln->end_Exon->contig); if ($tln->end_Exon->stable_id) { $copyexon->stable_id($tln->end_Exon->stable_id . "TER" . $metcnt++); $copyexon->created($time); $copyexon->modified($time); $copyexon->version(1); } $copynewendexon->phase($newendexon->phase); # $copynewendexon->end_phase(0); $copynewendexon->contig($newendexon->contig); if ($newendexon->stable_id) { $copynewendexon->stable_id($newendexon->stable_id . "TER" . $metcnt++); $copynewendexon->created($time); $copynewendexon->modified($time); $copynewendexon->version(1); } # TODO evidence if ($copynewendexon->strand == 1) { $tln->end($new_end - $copynewendexon->start + 1); } else { $tln->end($copynewendexon->end - $new_start + 1 ); my $tercodon = $copynewendexon->seq->subseq($copynewendexon->end - $new_start-1, $copynewendexon->end - $new_start +1); #reverse($tercodon); #$tercodon =~ tr /ACGT/TGCA/; } # Replace exons in transcript, and fix phases my @newexons; my $inrange = 0; foreach my $exon (@{$trans->get_all_Exons}) { if ($inrange) { print STDERR "in range exon before phase = " . $exon->phase . " endphase " . $exon->end_phase . "\n" if $self->VERBOSE; $exon->phase( $newexons[$#newexons]->end_phase ); $exon->end_phase(($exon->length + $exon->phase) % 3); print STDERR "in range exon after phase = " . $exon->phase . " endphase " . $exon->end_phase . "\n" if $self->VERBOSE; } if ($exon == $tln->end_Exon) { my $phase = $exon->phase; if ($phase == -1) { $phase = 0; } if ($exon == $tln->start_Exon) { $copyexon->end_phase(($exon->length - $tln->start + 1)%3); } else { $copyexon->end_phase(($exon->length + $exon->phase)%3); } print STDERR "Setting end_phase on old end exon to " . $copyexon->end_phase . " l = " . $exon->length . "\n" if $self->VERBOSE; push @newexons,$copyexon; $inrange = 1; } elsif ($exon == $newendexon) { $copynewendexon->phase( $newexons[$#newexons]->end_phase ); $copynewendexon->end_phase( -1); push @newexons,$copynewendexon; $inrange = 0; } else { push @newexons,$exon; } } $trans->flush_Exons; foreach my $exon (@newexons) { $trans->add_Exon($exon); } # Reset translation start exon if ($tln->end_Exon == $tln->start_Exon) { $tln->start_Exon($copyexon); } $tln->end_Exon($copynewendexon); } print_Translation($trans); # print "translateable seq = \n"; # print $trans->translateable_seq . "\n"; } else { print STDERR "Across exons - not handling this\n" if $self->VERBOSE; } last; } $coding_end += 3; } } else { print STDERR "Not enough bases downstream - NOT looking into genomic\n" if $self->VERBOSE; } } } } # These lines force loads from the database to stop attempted lazy # loading during the write (which fail because they are to the wrong # db) $trans->get_all_supporting_features(); my @exons= @{$trans->get_all_Exons}; my $get = $trans->translation; $trans->_translation_id(undef); foreach my $exon (@exons) { $exon->stable_id; $exon->contig($gene->slice); $exon->get_all_supporting_features; } } return($gene); } =head2 forward_genewise_clusters Arg [1] : Description: get/set for genewise clusters Returntype : =cut sub forward_genewise_clusters{ my ($self, $cluster_ref) = @_; if (!defined($self->{'_forward_genewise_clusters'})) { $self->{'_forward_genewise_clusters'} = []; } if (defined $cluster_ref && scalar(@{$cluster_ref})) { push(@{$self->{'_forward_genewise_clusters'}},@{$cluster_ref}); # store them sorted @{$self->{'_forward_genewise_clusters'}} = sort { $a->start <=> $b->start } @{$self->{'_forward_genewise_clusters'}}; } return $self->{'_forward_genewise_clusters'}; } =head2 reverse_genewise_clusters Arg [1] : Description: get/set for genewise clusters Returntype : Exceptions : Example : =cut sub reverse_genewise_clusters{ my ($self, $cluster_ref) = @_; if (!defined($self->{'_reverse_genewise_clusters'})) { $self->{'_reverse_genewise_clusters'} = []; } if (defined $cluster_ref && scalar(@{$cluster_ref})) { push(@{$self->{'_reverse_genewise_clusters'}},@{$cluster_ref}); # store them sorted @{$self->{'_reverse_genewise_clusters'}} = sort { $a->start <=> $b->start } @{$self->{'_reverse_genewise_clusters'}}; } return $self->{'_reverse_genewise_clusters'}; } =head2 cdna_genes Arg [1] : Description: get/set for e2g gene array Returntype : =cut sub cdna_genes { my ($self, $genes) = @_; if (!defined($self->{'_cdna_genes'})) { $self->{'_cdna_genes'} = []; } if (defined $genes && scalar(@{$genes})) { push(@{$self->{'_cdna_genes'}},@{$genes}); } return ($self->{'_cdna_genes'}); } =head2 prune Arg [1] : (optional) bool Description: get/set for option to prune genes =cut sub prune { my ($self, $bool) = @_; if (!defined($self->{'_prune_genes'})) { $self->{'_prune_genes'} = undef; } if (defined $bool) { $self->{'_prune_genes'} = $bool; } return ($self->{'_prune_genes'}); } =head2 ests Arg [1] : (optional) ref to array with ests Description: get/set for ests Returntype : array ref with ST objects =cut sub ests { my ($self, $ests) = @_; if (!defined($self->{'_ests'})) { $self->{'_ests'} = []; } if (defined $ests && scalar(@{$ests})) { push(@{$self->{'_ests'}},@{$ests}); } return($self->{'_ests'}); } =head2 ditags Arg [1] : (optional) ref to array with ditags Description: get/set ditags Returntype : array ref with ditag objects Exceptions : none =cut sub ditags { my ($self, $ditags) = @_; if (!defined($self->{'_ditags'})) { $self->{'_ditags'} = []; } if (defined $ditags && scalar(@{$ditags})) { push(@{$self->{'_ditags'}}, @{$ditags}); } return($self->{'_ditags'}); } =head2 gw_genes Arg [1] : Description: get/set for genewise gene array Returntype : Exceptions : Example : =cut sub gw_genes { my ($self, $genes) = @_; if (!defined($self->{'_gw_genes'})) { $self->{'_gw_genes'} = []; } if (defined $genes && scalar(@{$genes})) { push(@{$self->{'_gw_genes'}},@{$genes}); } return $self->{'_gw_genes'}; } =head2 blessed_genes Arg [1] : Description: get/set for blessed gene array clones the genes Returntype : =cut sub blessed_genes { my ($self, $slice, $genes) = @_; if (!defined($self->{'_blessed_genes'})) { $self->{'_blessed_genes'} = []; } if (defined $slice && defined $genes && scalar(@{$genes})) { # split input genes into one transcript per gene; keep type the same OLDGENE: foreach my $gene(@{$genes}){ foreach my $transcript(@{$gene->get_all_Transcripts}){ # make a new gene my $newgene = new Bio::EnsEMBL::Gene; $newgene->biotype($gene->biotype); #preserve the xref of the blessed genes! foreach my $dblink (@{$gene->get_all_DBLinks()}){ #print STDERR "Adding xref ".$dblink->display_id."\n"; $newgene->add_DBEntry($dblink); } # clone transcript my $newtranscript = new Bio::EnsEMBL::Transcript; $newtranscript->slice($slice); # clone translation my $newtranslation = new Bio::EnsEMBL::Translation; my @seqeds = @{$transcript->translation->get_all_SeqEdits}; if (scalar(@seqeds)) { print "Copying sequence edits\n" if $self->VERBOSE; foreach my $se (@seqeds) { my $new_se = Bio::EnsEMBL::SeqEdit->new( -CODE => $se->code, -NAME => $se->name, -DESC => $se->description, -START => $se->start, -END => $se->end, -ALT_SEQ => $se->alt_seq ); my $attribute = $new_se->get_Attribute(); $newtranslation->add_Attributes($attribute); } } my @support = @{$transcript->get_all_supporting_features}; if (scalar(@support)) { $newtranscript->add_supporting_features(@support); } $newtranscript->add_Attributes(@{$transcript->get_all_Attributes}); $newtranscript->translation($newtranslation); foreach my $exon(@{$transcript->get_all_Exons}){ # clone the exon my $newexon = clone_Exon($exon); # get rid of stable ids $newexon->stable_id(''); # if this is start/end of translation, keep that info: if ( $exon == $transcript->translation->start_Exon ){ $newtranslation->start_Exon( $newexon ); $newtranslation->start($transcript->translation->start); } if ( $exon == $transcript->translation->end_Exon ){ $newtranslation->end_Exon( $newexon ); $newtranslation->end($transcript->translation->end); } $newtranscript->add_Exon($newexon); #add sf $newexon->add_supporting_features(@{$exon->get_all_supporting_features}); } $newgene->add_Transcript($newtranscript); push(@{$self->{'_blessed_genes'}}, $newgene); } } } return $self->{'_blessed_genes'}; } =head2 combined_genes Arg [1] : ref to Bio::EnsEMBL::Gene Description: get/set for combined gene array Returntype : =cut sub combined_genes { my ($self, $genesref) = @_; if (!defined($self->{'_combined_genes'})) { $self->{'_combined_genes'} = []; } if (defined $genesref) { push(@{$self->{'_combined_genes'}}, $genesref); } # if (defined $genesref && scalar(@{$genesref})) { # push(@{$self->{'_combined_genes'}},@{$genesref}); # } return $self->{'_combined_genes'}; } =head2 unmatched_genes Arg [1] : array of genes Description: get/set for unmatched gene array Returntype : none =cut sub unmatched_genes { my ($self, @genes) = @_; if (!defined($self->{'_unmatched_genes'})) { $self->{'_unmatched_genes'} = []; } if(@genes){ push(@{$self->{'_unmatched_genes'}},@genes); } return $self->{'_unmatched_genes'}; } =head2 populate_kill_list Arg [1] : String mol-type Description: read ids to ignore from KillList db Returntype : ref to hash =cut sub populate_kill_list { my ($self, $type) = @_; my $kill_list_object = Bio::EnsEMBL::KillList::KillList->new(-TYPE => $type); my %kill_list = %{$kill_list_object->get_kill_list()}; return \%kill_list; } =head2 kill_list Arg [1] : optional hash ref with kill list Description: get/set kill list object Returntype : ref to hash =cut sub kill_list { my ($self, $kill_list) = @_; if($kill_list){ $self->{'_kill_list'} = $kill_list; } return $self->{'_kill_list'}; } =head2 retrieve_unmerged_gene Arg [1] : Bio::EnsEMBL::Gene Description: Returns unmeregd, frameshifted version of input gene Returntype : Bio::EnsEMBL::Gene =cut sub retrieve_unmerged_gene{ my ($self, $merged_gene) = @_; my %pairs = %{$self->merged_unmerged_pairs()}; if(!exists $pairs{$merged_gene}){ print STDERR "Can't retrieve unmerged \n"; print STDERR "for gene ".$merged_gene->dbID." \n"; } return $pairs{$merged_gene}; } =head2 merged_unmerged_pairs Arg [1] : Description: get/set for pairs of frameshift merged and unmerged genes. Key is merged gene, value is unmerged Returntype : =cut sub merged_unmerged_pairs { my ($self, $merged_gene, $unmerged_gene) = @_; if (!defined($self->{'_merged_unmerged_pairs'})) { $self->{'_merged_unmerged_pairs'} = {}; } if ($unmerged_gene && $merged_gene) { $self->{'_merged_unmerged_pairs'}{$merged_gene}= $unmerged_gene; } # hash ref return $self->{'_merged_unmerged_pairs'}; } =head2 modified_unmodified_pairs Arg [1] : Description: get/set for pairs of UTR modified and unmodified genes. Key is modified gene, value is unmodified Returntype : =cut sub modified_unmodified_pairs { my ($self, $modified_gene, $unmodified_gene) = @_; if (!defined($self->{'_modified_unmodified_pairs'})) { $self->{'_modified_unmodified_pairs'} = {}; } if ($unmodified_gene && $modified_gene) { $self->{'_modified_unmodified_pairs'}{$modified_gene}= $unmodified_gene; } # hash ref return $self->{'_modified_unmodified_pairs'}; } =head2 output Arg [1] : optional RunnableDB output Description: get/set for result of Analysis Returntype : array ref with genes =cut sub output{ my ($self,@genes) = @_; if (!defined($self->{'_output'})) { $self->{'_output'} = []; } if(@genes){ push(@{$self->{'_output'}},@genes); } return $self->{'_output'}; } ############################################### =head2 _get_evidence_set ($logic_name_or_biotype) Name : get_evidence_set( $logic_name_or_biotype ) Arg : String Func : returns the name of the evidence_set of a genee / PredictionTranscript Returnval: String describing evidence_set_name =cut sub _get_evidence_set { my ($self, $logic_name_or_biotype) = @_ ; my %ev_sets = %{ $self->{evidence_sets} } ; my $result_set_name; for my $set_name (keys %ev_sets){ my @logic_names = @{$ev_sets{$set_name}} ; for my $ln (@logic_names ) { if ($ln eq $logic_name_or_biotype){ $result_set_name = $set_name ; } } } return $result_set_name; } #=head2 make_seqfetcher # Title : make_seqfetcher # Usage : # Description: for the KnownUTR module # Example : # Returns : Bio::DB::RandomAccessI # Args : #=cut #sub make_seqfetcher{ # my ( $self, $index, $seqfetcher_class ) = @_; # my $seqfetcher; # (my $class = $seqfetcher_class) =~ s/::/\//g; # eval{ # require "$class.pm"; # }; # if($@){ warn("Could not require seqfetcher class.\n") } # print "using index $index \n" if $self->VERBOSE; # if(defined $index && $index ne ''){ # my @db = ( $index ); # # make sure that your class is compatible with the index type # $seqfetcher = "$seqfetcher_class"->new('-db' => \@db, ); # } # else{ # throw("can't make seqfetcher\n"); # } # return $seqfetcher; #} #=head2 _cdna_seqfetcher # Arg [1] : # Description: # Returntype : #=cut #sub _cdna_seqfetcher{ # my ($self, $seqfetcher) = @_; # if($seqfetcher){ # $self->{'_cdna_seqfetcher'} = $seqfetcher; # } # return $self->{'_cdna_seqfetcher'}; #} =head2 _cdna_slice Arg [1] : Description: Returntype : =cut sub _cdna_slice { my ($self, $slice) = @_; if($slice){ $self->{'_cdna_slice'} = $slice; } return $self->{'_cdna_slice'}; } =head2 _cdna_evidence Arg [1] : Description: Returntype : =cut sub _cdna_evidence { my ($self, $cdna_evidence) = @_; if($cdna_evidence){ $self->{'_cdna_evidence'} = $cdna_evidence; } return $self->{'_cdna_evidence'}; } =head2 _known_pairs Arg [1] : optional hash ref Description: store the links between NM and NP entries Returntype : has ref =cut sub _known_pairs { my ($self, $hashref) = @_; if (defined $hashref) { $self->{_known_pairs} = $hashref; } return $self->{_known_pairs}; } ######## database access functions ########### =head2 CDNA_DB Arg [1] : optional Bio::EnsEMBL::DBSQL::DBAdaptor Description: get/set for db storing exonerate alignments of cDNAs Returntype : Bio::EnsEMBL::DBSQL::DBAdaptor Exceptions : Throws if input isn't DBAdaptor or db parameters haven't been defined =cut sub CDNA_DB { my( $self, $cdna_db ) = @_; if ($cdna_db){ $self->{_cdna_db} = get_db_adaptor_by_string($cdna_db,1); } return $self->{_cdna_db}; } =head2 EST_DB Arg [1] : optional Bio::EnsEMBL::DBSQL::DBAdaptor Description: get/set for db storing exonerate alignments of ESTs Returntype : Bio::EnsEMBL::DBSQL::DBAdaptor Exceptions : Throws if input isn't DBAdaptor or db parameters haven't been defined =cut sub EST_DB { my( $self, $est_db ) = @_; if ($est_db){ $self->{_est_db} = get_db_adaptor_by_string($est_db,1); } return $self->{_est_db}; } =head2 DITAG_DB Arg [1] : optional Bio::EnsEMBL::DBSQL::DBAdaptor Description: get/set for db storing exonerate alignments of ditags Returntype : Bio::EnsEMBL::DBSQL::DBAdaptor Exceptions : Throws if input isn't DBAdaptor or db parameters haven't been defined =cut sub DITAG_DB { my( $self, $ditag_db ) = @_; if ($ditag_db){ $self->{_ditag_db} = get_db_adaptor_by_string($ditag_db,1); } if((defined $self->{_ditag_db}) && (!$self->{_ditag_db}) && (scalar @{$self->DITAG_TYPE_NAMES})){ throw("You have defined DITAG_TYPE_NAMES, but no DITAG_DB parameters.\n"); } return $self->{_ditag_db}; } =head2 INPUT_DB Arg [1] : optional Bio::EnsEMBL::DBSQL::DBAdaptor Description: get/set for db to read input genes from Returntype : Bio::EnsEMBL::DBSQL::DBAdaptor Exeptions : Throws if input isn't DBAdaptor or db parameters haven't been defined =cut sub INPUT_DB { my( $self, $input_db ) = @_; if ($input_db){ $self->{_input_db} = get_db_adaptor_by_string($input_db,1); } if(!$self->{_input_db}){ throw("Please define database parameters for input db.\n"); } return $self->{_input_db}; } =head2 blessed_db Arg [1] : optional Bio::EnsEMBL::DBSQL::DBAdaptor Description: get/set for db storing blessed gene structures Returntype : Bio::EnsEMBL::DBSQL::DBAdaptor Exceptions : Throws if input isn't DBAdaptor or db parameters haven't been defined =cut sub BLESSED_DB { my( $self, $blessed_db ) = @_; if ($blessed_db){ $self->{_blessed_db} = get_db_adaptor_by_string($blessed_db,1); } return $self->{_blessed_db}; } =head2 OUTPUT_DB Arg [1] : optional Bio::EnsEMBL::DBSQL::DBAdaptor Description: get/set for db in which UTR modified genes should be stored Returntype : Bio::EnsEMBL::DBSQL::DBAdaptor Exceptions : Throws if input isn't DBAdaptor or db parameters haven't been defined =cut sub OUTPUT_DB { my( $self, $output_db ) = @_; if ($output_db){ $self->{_output_db} = get_db_adaptor_by_string($output_db,1); } if(!$self->{_output_db}){ throw("Please define database parameters for output db.\n"); } return $self->{_output_db}; } #### other config variable get/set methods ###### =head2 various_config_variable_methods Arg [1] : optional parameter Description: setter / getter for config vars Returntype : config value =cut sub INPUT_GENETYPES { my ($self, $input_genetypes) = @_; if($input_genetypes){ $self->{'_input_genetypes'} = $input_genetypes; } return $self->{'_input_genetypes'}; } sub BLESSED_GENETYPES { my ($self,$value) = @_; if (defined $value) { $self->{_blessed_genetypes} = $value; } return $self->{_blessed_genetypes}; } sub VERBOSE { my ($self,$value) = @_; if (defined $value) { $self->{_verbose} = $value; } return $self->{_verbose}; } sub MAX_EXON_LENGTH { my ($self,$value) = @_; if (defined $value) { $self->{_max_exon_length} = $value; } return $self->{_max_exon_length}; } sub DITAG_TYPE_NAMES { my ($self,$value) = @_; if (defined $value) { $self->{_ditag_type_names} = $value; } return $self->{_ditag_type_names}; } sub cDNA_GENETYPE { my ($self,$value) = @_; if (defined $value) { $self->{_cdna_genetype} = $value; } return $self->{_cdna_genetype}; } sub EXTEND_ORIGINAL_BIOTYPE { my ($self,$value) = @_; if (defined $value) { $self->{_extend_original_biotype} = $value; } return $self->{_extend_original_biotype}; } sub EXTEND_BIOTYPE_OF_UNCHANGED_GENES { my ($self,$value) = @_; if (defined $value) { $self->{_extend_biotype_of_unchanged_genes} = $value; } return $self->{_extend_biotype_of_unchanged_genes}; } sub KNOWN_UTR_GENETYPE { my ($self,$value) = @_; if (defined $value) { $self->{_known_utr_genetype} = $value; } return $self->{_known_utr_genetype}; } sub UTR_GENETYPE { my ($self,$value) = @_; if (defined $value) { $self->{_utr_genetype} = $value; } return $self->{_utr_genetype}; } sub PRUNE_GENES { my ($self,$value) = @_; if (defined $value) { $self->{_prune_genes} = $value; } return $self->{_prune_genes}; } sub FILTER_ESTS { my ($self,$value) = @_; if (defined $value) { $self->{_filter_ests} = $value; } return $self->{_filter_ests}; } sub MAX_INTRON_LENGTH { my ($self,$value) = @_; if (defined $value) { $self->{_max_intron} = $value; } return $self->{_max_intron}; } sub EST_GENETYPE { my ($self,$value) = @_; if (defined $value) { $self->{_est_genetype} = $value; } return $self->{_est_genetype}; } sub BLESSED_UTR_GENETYPE { my ($self,$value) = @_; if (defined $value) { $self->{_blessed_utr_genetype} = $value; } return $self->{_blessed_utr_genetype}; } sub LOOK_FOR_KNOWN { my ($self,$value) = @_; if (defined $value) { $self->{_look_for_known} = $value; } return $self->{_look_for_known}; } sub DITAG_WINDOW { my ($self,$value) = @_; if (defined $value) { $self->{_ditag_window} = $value; } return $self->{_ditag_window}; } sub KNOWNUTR_FILE { my ($self,$value) = @_; if (defined $value) { $self->{_knownutr_file} = $value; } return $self->{_knownutr_file}; } 1;