Archive Ensembl HomeArchive Ensembl Home
DnaAlignFeatureAdaptor.pm
Go to the documentation of this file.
00001 # Copyright (c) 1999-2012
00002 #
00003 # Ensembl module for Bio::EnsEMBL::Compara::DBSQL::DnaAlignFeatureAdaptor
00004 #
00005 # You may distribute this module under the same terms as perl itself
00006 
00007 =head1 NAME
00008 
00009 Bio::EnsEMBL::Compara::DBSQL::DnaAlignFeatureAdaptor
00010 
00011 =head1 SYNOPSIS
00012 
00013 $dafa = $compara_dbadaptor->get_DnaAlignFeatureAdaptor;
00014 @align_features = @{$dafa->fetch_by_Slice_species($slice, $qy_species)};
00015 
00016 =head1 DESCRIPTION
00017 
00018 Retrieves alignments from a compara database in the form of DnaDnaAlignFeatures
00019 
00020 =head1 CONTACT
00021 
00022 Post questions to the EnsEMBL developer list: <dev@ensembl.org>
00023 
00024 =cut
00025 
00026 
00027 package Bio::EnsEMBL::Compara::DBSQL::DnaAlignFeatureAdaptor;
00028 use strict;
00029 use vars qw(@ISA);
00030 
00031 use Bio::EnsEMBL::DBSQL::BaseAdaptor;
00032 use Bio::EnsEMBL::Utils::Cache; #CPAN LRU cache
00033 use Bio::EnsEMBL::DnaDnaAlignFeature;
00034 
00035 use Bio::EnsEMBL::Utils::Exception;
00036 
00037 
00038 @ISA = qw(Bio::EnsEMBL::DBSQL::BaseAdaptor);
00039 
00040 my $CACHE_SIZE = 4;
00041 
00042 =head2 new
00043 
00044   Arg [1]    : list of args to super class constructor
00045   Example    : $dafa = new Bio::EnsEMBL::Compara::Genomi
00046   Description: Creates a new DnaAlignFeatureAdaptor.  The superclass
00047                constructor is extended to initialise an internal cache.  This
00048                class should be instantiated through the get method on the
00049                DBAdaptor rather than calling this method directory.
00050   Returntype : none
00051   Exceptions : none
00052   Caller     : Bio::EnsEMBL::DBSQL::DBConnection
00053 
00054 =cut
00055 
00056 sub new {
00057   my ($class, @args) = @_;
00058 
00059   my $self = $class->SUPER::new(@args);
00060 
00061   #initialize internal LRU cache
00062   tie(%{$self->{'_cache'}}, 'Bio::EnsEMBL::Utils::Cache', $CACHE_SIZE);
00063 
00064   return $self;
00065 }
00066 
00067 
00068 
00069 =head2 fetch_all_by_species_region
00070 
00071  Arg [1]    : string $cs_species
00072               e.g. "Homo sapiens"
00073  Arg [2]    : string $cs_assembly (can be undef)
00074               e.g. "NCBI_31" if undef assembly_default will be taken
00075  Arg [3]    : string $qy_species
00076               e.g. "Mus musculus"
00077  Arg [4]    : string $qy_assembly (can be undef)
00078               e.g. "MGSC_3", if undef assembly_default will be taken
00079  Arg [5]    : string $chromosome_name
00080               the name of the chromosome to retrieve alignments from (e.g. 'X')
00081  Arg [6]    : int start
00082  Arg [7]    : int end
00083  Arg [8]    : string $alignment_type
00084               The type of alignments to be retrieved
00085               e.g. WGA or WGA_HCR
00086  Example    : $gaa->fetch_all_by_species_region("Homo sapiens", "NCBI_31",
00087                         "Mus musculus", "MGSC_3",
00088                                                 "X", 250_000, 750_000,"WGA");
00089  Description: Retrieves alignments between the consensus and query species
00090               from a specified region of the consensus genome.
00091  Returntype : an array reference of Bio::EnsEMBL::DnaDnaAlignFeature objects
00092  Exceptions : returns a ref to an empty list if requested DnaFrag or MethodLinkSpeciesSet
00093               are not in the compara DB.
00094  Caller     : general
00095 
00096 =cut
00097 
00098 sub fetch_all_by_species_region {
00099   my ($self, $consensus_species, $consensus_assembly,
00100       $query_species, $query_assembly,
00101       $chromosome_name, $start, $end, $alignment_type, $limit,$dnafrag_type) = @_;
00102 
00103   $limit = 0 unless (defined $limit);
00104 
00105   #get the genome database for each species
00106   my $genome_db_adaptor = $self->db->get_GenomeDBAdaptor;
00107   my $consensus_genome_db = $genome_db_adaptor->fetch_by_name_assembly($consensus_species, $consensus_assembly);
00108   my $query_genome_db = $genome_db_adaptor->fetch_by_name_assembly($query_species, $query_assembly);
00109 
00110   #retrieve dna fragments from the subjects species region of interest
00111   my $dna_frag_adaptor = $self->db->get_DnaFragAdaptor;
00112   my $this_dnafrag = $dna_frag_adaptor->fetch_by_GenomeDB_and_name(
00113           $consensus_genome_db,
00114           $chromosome_name,
00115       );
00116   return [] if (!$this_dnafrag);
00117 
00118   # Get the Bio::EnsEMBL::Compara::MethodLinkSpeciesSet object corresponding to the alignment_type and
00119   # the couple of genomes
00120   my $method_link_species_set_adaptor = $self->db->get_MethodLinkSpeciesSetAdaptor;
00121   my $method_link_species_set;
00122   if ($consensus_genome_db->dbID == $query_genome_db->dbID) {
00123     # Allow to fetch the right method_link_species_set for self-comparisons!
00124     $method_link_species_set = $method_link_species_set_adaptor->fetch_by_method_link_type_GenomeDBs(
00125         $alignment_type, [$consensus_genome_db]);
00126   } else {
00127     # Normal, pairwise comparisons...
00128     $method_link_species_set = $method_link_species_set_adaptor->fetch_by_method_link_type_GenomeDBs(
00129           $alignment_type,
00130           [$consensus_genome_db, $query_genome_db]
00131       );
00132   }
00133   return [] if (!$method_link_species_set);
00134 
00135 #   my $gaa = $self->db->get_GenomicAlignAdaptor;
00136   my $genomic_align_block_adaptor = $self->db->get_GenomicAlignBlockAdaptor;
00137 
00138   my @out = ();
00139 
00140   my $consensus_slice_adaptor = $consensus_genome_db->db_adaptor->get_SliceAdaptor;
00141   my $query_slice_adaptor;
00142   eval {
00143     $query_slice_adaptor = $query_genome_db->db_adaptor->get_SliceAdaptor;
00144   };
00145   #caclulate coords relative to start of dnafrag
00146   ## Bio::EnsEMBL::Compara::Dnafrag::start is always 1
00147   #     my $this_dnafrag_start = $start - $this_dnafrag->start + 1;
00148   #     my $this_dnafrag_end   = $end   - $this_dnafrag->start + 1;
00149   my $this_dnafrag_start = $start;
00150   my $this_dnafrag_end   = $end;
00151 
00152   #constrain coordinates so they are completely within the dna frag
00153   my $this_dnafrag_length = $this_dnafrag->length;
00154   $this_dnafrag_start = 1 unless (defined $this_dnafrag_start);
00155   $this_dnafrag_start = ($this_dnafrag_start < 1)  ? 1 : $this_dnafrag_start;
00156 
00157   $this_dnafrag_end = $this_dnafrag_length unless (defined $this_dnafrag_end);
00158   $this_dnafrag_end   = ($this_dnafrag_end > $this_dnafrag_length) ? $this_dnafrag_length : $this_dnafrag_end;
00159 
00160   #fetch all alignments in the region we are interested in
00161   my $genomic_align_blocks = $genomic_align_block_adaptor->fetch_all_by_MethodLinkSpeciesSet_DnaFrag(
00162           $method_link_species_set,
00163           $this_dnafrag,
00164           $this_dnafrag_start,
00165           $this_dnafrag_end,
00166           $limit
00167       );
00168 
00169   #convert genomic align blocks to dna align features
00170   my $dafs = _convert_GenomicAlignBlocks_into_DnaDnaAlignFeatures($genomic_align_blocks);
00171 
00172   # We need to attach slices of the entire seq region to the features.
00173   # The features come without any slices at all, but their coords are
00174   # relative to the beginning of the seq region.
00175   my $top_slice = $consensus_slice_adaptor->fetch_by_region($dnafrag_type, $chromosome_name);
00176 
00177   map {$_->slice($top_slice)} @$dafs;
00178 
00179   return $dafs;
00180 }
00181 
00182 
00183 
00184 
00185 =head2 fetch_all_by_Slice
00186 
00187  Arg [1]    : Bio::EnsEMBL::Slice
00188  Arg [2]    : string $qy_species
00189               The query species to retrieve alignments against
00190  Arg [3]    : string $qy_assembly
00191  Arg [4]    : string $$alignment_type
00192               The type of alignments to be retrieved
00193               e.g. WGA or WGA_HCR
00194  Example    : $gaa->fetch_all_by_Slice($slice, "Mus musculus","WGA");
00195  Description: find matches of query_species in the region of a slice of a
00196               subject species
00197  Returntype : an array reference of Bio::EnsEMBL::DnaDnaAlignFeature objects
00198  Exceptions : none
00199  Caller     : general
00200 
00201 =cut
00202 
00203 sub fetch_all_by_Slice {
00204   my ($self, $orig_slice, $qy_species, $qy_assembly, $alignment_type,
00205       $limit) = @_;
00206 
00207   unless($orig_slice && ref $orig_slice &&
00208          $orig_slice->isa('Bio::EnsEMBL::Slice')) {
00209     throw("Invalid slice argument [$orig_slice]\n");
00210   }
00211 
00212   unless($qy_species) {
00213     throw("Query species argument is required");
00214   }
00215 
00216   $limit = 0 unless (defined $limit);
00217 
00218   my $genome_db_adaptor = $self->db->get_GenomeDBAdaptor();
00219   my $cs_genome_db = $genome_db_adaptor->fetch_by_Slice($orig_slice);
00220   my $qy_genome_db = $genome_db_adaptor->fetch_by_name_assembly($qy_species, $qy_assembly);
00221   if(! $qy_genome_db && ! $qy_assembly) {
00222     #Only use the registry if we cannot find one by name & assembly was not defined.
00223     $qy_genome_db = $genome_db_adaptor->fetch_by_registry_name($qy_species);
00224   }
00225   return [] if (!$cs_genome_db or !$qy_genome_db);
00226 
00227   my $method_link_species_set_adaptor = $self->db->get_MethodLinkSpeciesSetAdaptor();
00228   my $method_link_species_set;
00229   if ($cs_genome_db->dbID == $qy_genome_db->dbID) {
00230     $method_link_species_set = $method_link_species_set_adaptor->
00231         fetch_by_method_link_type_GenomeDBs($alignment_type, [$cs_genome_db]);
00232   } else {
00233     $method_link_species_set = $method_link_species_set_adaptor->
00234         fetch_by_method_link_type_GenomeDBs($alignment_type,
00235             [$cs_genome_db, $qy_genome_db]);
00236   }
00237 
00238   my $key = uc(join(':', $orig_slice->name,
00239                     $cs_genome_db->name, $qy_genome_db->name, $qy_genome_db->assembly, $alignment_type));
00240 
00241   if(exists $self->{'_cache'}->{$key}) {
00242     return $self->{'_cache'}->{$key};
00243   }
00244 
00245   my $genomic_align_block_adaptor = $self->db->get_GenomicAlignBlockAdaptor();
00246   my $genomic_align_blocks = $genomic_align_block_adaptor->fetch_all_by_MethodLinkSpeciesSet_Slice(
00247       $method_link_species_set, $orig_slice, $limit);
00248 
00249   my $dafs = _convert_GenomicAlignBlocks_into_DnaDnaAlignFeatures($genomic_align_blocks);
00250 
00251   #update the cache
00252   $self->{'_cache'}->{$key} = $dafs;
00253   return $dafs;
00254 }
00255 
00256 
00257 =head2 interpolate_best_location
00258 
00259   Arg [1]    : Bio::EnsEMBL::Slice $slice
00260   Arg [2]    : string $species
00261                e.g. "Mus musculus"
00262   Arg [3]    : string $alignment_type
00263                e.g. "BLASTZ_NET"
00264   Arg [4]    : string $seq_region_name
00265                e.g. "6-COX"
00266   Example    :
00267   Description:
00268   Returntype : array with 3 elements
00269   Exceptions :
00270   Caller     :
00271 
00272 =cut
00273 
00274 sub interpolate_best_location {
00275   my ($self,$slice,$species,$alignment_type,$seq_region_name) = @_;
00276 
00277 #warn $slice->name,"\t$species\t$alignment_type\t$seq_region_name";
00278 
00279   $| =1 ;
00280   my $max_distance_for_clustering = 10000;
00281   my $dafs = $self->fetch_all_by_Slice($slice, $species, undef, $alignment_type);
00282   my %name_strand_clusters;
00283   my $based_on_group_id = 1;
00284   foreach my $daf (@{$dafs}) {
00285     next if ($seq_region_name && $daf->hseqname ne $seq_region_name);
00286     if (defined $daf->group_id && $daf->group_id > 0 && $alignment_type ne "TRANSLATED_BLAT") {
00287       push @{$name_strand_clusters{$daf->group_id}}, $daf;
00288     } else {
00289       $based_on_group_id = 0 if ($based_on_group_id);
00290       push @{$name_strand_clusters{$daf->hseqname. "_" .$daf->hstrand}}, $daf;
00291     }
00292   }
00293 
00294   if ($based_on_group_id) {
00295     my @ordered_name_strands = sort {scalar @{$name_strand_clusters{$b}} <=> scalar @{$name_strand_clusters{$a}}} keys %name_strand_clusters;
00296 
00297     my @best_blocks = sort {$a->hstart <=> $b->hend} @{$name_strand_clusters{$ordered_name_strands[0]}||[]};
00298 
00299     return undef if( !@best_blocks );
00300     return ($best_blocks[0]->hseqname,
00301             $best_blocks[0]->hstart
00302             + int(($best_blocks[-1]->hend - $best_blocks[0]->hstart)/2),
00303             $best_blocks[0]->hstrand * $slice->strand,
00304             $best_blocks[0]->hstart,
00305             $best_blocks[-1]->hend);
00306 
00307   } else {
00308 
00309     my @refined_clusters;
00310     foreach my $name_strand (keys %name_strand_clusters) {
00311       # an array of arrayrefs
00312       # name, strand, start, end, nb of blocks
00313       my @sub_clusters;
00314       foreach my $block (sort {$a->hstart <=> $b->hstart} @{$name_strand_clusters{$name_strand}||[]}) {
00315         unless (scalar @sub_clusters) {
00316           push @sub_clusters, [$block->hseqname,$block->hstrand, $block->hstart, $block->hend, 1];
00317           next;
00318         }
00319         my $block_clustered = 0;
00320         foreach my $arrayref (@sub_clusters) {
00321           my ($n,$st,$s,$e,$c) = @{$arrayref};
00322           if ($block->hstart<=$e &&
00323               $block->hend>=$s) {
00324             # then overlaps.
00325             $arrayref->[2] = $block->hstart if ($block->hstart < $s);
00326             $arrayref->[3] = $block->hend if ($block->hend > $e);
00327             $arrayref->[4]++;
00328             $block_clustered = 1;
00329           } elsif ($block->hstart <= $e + $max_distance_for_clustering &&
00330                    $block->hstart > $e) {
00331             # then is downstream
00332             $arrayref->[3] = $block->hend;
00333             $arrayref->[4]++;
00334             $block_clustered = 1;
00335           } elsif ($block->hend >= $s - $max_distance_for_clustering &&
00336                    $block->hend < $s) {
00337             # then is upstream
00338             $arrayref->[2] = $block->hstart;
00339             $arrayref->[4]++;
00340             $block_clustered = 1;
00341           }
00342         }
00343         unless ($block_clustered) {
00344           # do not overlap anything already seen, so adding as new seeding cluster
00345           push @sub_clusters, [$block->hseqname,$block->hstrand, $block->hstart, $block->hend, 1];
00346         }
00347       }
00348       push @refined_clusters, @sub_clusters;
00349     }
00350 
00351     # sort by the max number of blocks desc
00352     @refined_clusters = sort {$b->[-1] <=> $a->[-1]} @refined_clusters;
00353 
00354     return undef if(!@refined_clusters);
00355     return ($refined_clusters[0]->[0], #hseqname,
00356             $refined_clusters[0]->[2]
00357             + int(($refined_clusters[0]->[3] - $refined_clusters[0]->[2])/2),
00358             $refined_clusters[0]->[1] * $slice->strand,
00359             $refined_clusters[0]->[2],
00360             $refined_clusters[0]->[3]);
00361 
00362   }
00363 }
00364 
00365 
00366 =head2 deleteObj
00367 
00368   Arg [1]    : none
00369   Example    : none
00370   Description: Called automatically by DBConnection during object destruction
00371                phase. Clears the cache to avoid memory leaks.
00372   Returntype : none
00373   Exceptions : none
00374   Caller     : none
00375 
00376 =cut
00377 
00378 sub deleteObj {
00379   my $self = shift;
00380 
00381   $self->SUPER::deleteObj;
00382 
00383   #clear the cache, removing references
00384   %{$self->{'_cache'}} = ();
00385 }
00386 
00387 
00388 =head2 deleteObj
00389 
00390   Arg [1]    : none
00391   Example    : none
00392   Description: Called automatically by DBConnection during object destruction
00393                phase. Clears the cache to avoid memory leaks.
00394   Returntype : none
00395   Exceptions : none
00396   Caller     : none
00397 
00398 =cut
00399 
00400 sub _convert_GenomicAlignBlocks_into_DnaDnaAlignFeatures {
00401   my ($genomic_align_blocks) = @_;
00402   my $dna_dna_align_features = [];
00403 
00404   my $query_slice_adaptor;
00405   foreach my $this_genomic_align_block (@$genomic_align_blocks) {
00406 
00407     ## KNOWN BUG: This will ignore third and following parts of a multiple alignment...
00408     ## This adaptor cannot deal with multiple alignments. Use the new
00409     ## Bio::EnsEMBL::Compara::DBSQL::GenomicAlignBlockAdaptor instead.
00410     my $consensus_genomic_align = $this_genomic_align_block->reference_genomic_align;
00411     my $query_genomic_align = $this_genomic_align_block->get_all_non_reference_genomic_aligns->[0];
00412 
00413     my $top_slice;
00414     if (!defined($query_slice_adaptor)) {
00415       $query_slice_adaptor = $query_genomic_align->get_Slice->adaptor;
00416     }
00417     if ($query_slice_adaptor) {
00418       $top_slice = $query_slice_adaptor->fetch_by_region(
00419               $query_genomic_align->dnafrag->coord_system_name,
00420               $query_genomic_align->dnafrag->name
00421           );
00422     } else {
00423       $top_slice = undef;
00424     }
00425 
00426     ## The code for transforming GenomicAlignBlocks into DnaDnaAlignFeatures assumes that
00427     ## reference_genomic_align is on the forward strand!
00428     if ($consensus_genomic_align->dnafrag_strand == -1) {
00429       $this_genomic_align_block->reverse_complement;
00430     }
00431     my $cstart = $consensus_genomic_align->dnafrag_start;
00432     my $cend   = $consensus_genomic_align->dnafrag_end;
00433 
00434     #skip features which do not overlap the requested region
00435     #next if ($cstart > $end || $cend < $start);
00436 
00437     my $ga_cigar_line;
00438     do {
00439       my @consensus_cigar_pieces = split(/(\d*[DIMG])/, $consensus_genomic_align->cigar_line);
00440       my @query_cigar_pieces = split(/(\d*[DIMG])/, $query_genomic_align->cigar_line);
00441 
00442       my @consensus_gapped_pieces;
00443       foreach my $piece (@consensus_cigar_pieces) {
00444         next if ($piece !~ /^(\d*)([MDIG])$/);
00445         my $num = $1;
00446         my $type = $2;
00447         $num = 1 if ($num !~ /^\d+$/);
00448         if( $type eq "M" ) {
00449           for (my $i=0; $i<$num; $i++) {push(@consensus_gapped_pieces, "N")}
00450         } else {
00451           for (my $i=0; $i<$num; $i++) {push(@consensus_gapped_pieces, '-')}
00452         }
00453       }
00454       my @query_gapped_pieces;
00455       foreach my $piece (@query_cigar_pieces) {
00456         next if ($piece !~ /^(\d*)([MDIG])$/);
00457         my $num = $1;
00458         my $type = $2;
00459         $num = 1 if ($num !~ /^\d+$/);
00460         if( $type eq "M" ) {
00461           for (my $i=0; $i<$num; $i++) {push(@query_gapped_pieces, "N")}
00462         } else {
00463           for (my $i=0; $i<$num; $i++) {push(@query_gapped_pieces, '-')}
00464         }
00465       }
00466       throw if (scalar(@consensus_gapped_pieces) != scalar(@query_gapped_pieces));
00467       my $type = "";
00468       my $num = 0;
00469       for (my $i=0; $i<@consensus_gapped_pieces; $i++) {
00470         if ($consensus_gapped_pieces[$i] eq "N" and $query_gapped_pieces[$i] eq "N") {
00471           if ($type ne "M") {
00472             $ga_cigar_line .= (($num==1)?"":$num).$type if ($num);
00473             $num = 0;
00474             $type = "M";
00475           }
00476         } elsif ($consensus_gapped_pieces[$i] eq "N" and $query_gapped_pieces[$i] eq "-") {
00477           if ($type ne "I") {
00478             $ga_cigar_line .= (($num==1)?"":$num).$type if ($num);
00479             $num = 0;
00480             $type = "I";
00481           }
00482         } elsif ($consensus_gapped_pieces[$i] eq "-" and $query_gapped_pieces[$i] eq "N") {
00483           if ($type ne "D") {
00484             $ga_cigar_line .= (($num==1)?"":$num).$type if ($num);
00485             $num = 0;
00486             $type = "D";
00487           }
00488         } else {
00489           throw "no double gaps can occur in a pairwise alignment!";
00490         }
00491         $num++;
00492       }
00493       $ga_cigar_line .= (($num==1)?"":$num).$type;
00494     };
00495     my $df_name = $consensus_genomic_align->dnafrag->name;
00496     my $score = $this_genomic_align_block->score;
00497     my $perc_id = $this_genomic_align_block->perc_id;
00498     my $qdf_start = 1;
00499     my $ga_query_start = $query_genomic_align->dnafrag_start;
00500     my $ga_query_end = $query_genomic_align->dnafrag_end;
00501     my $ga_query_strand = $query_genomic_align->dnafrag_strand;
00502     my $qdf_name = $query_genomic_align->dnafrag->name;
00503     my $ga_strands_reversed = 0;
00504     if ($consensus_genomic_align->dnafrag_strand == -1) {
00505       $ga_strands_reversed = 1;
00506       $ga_query_strand = -$ga_query_strand;
00507     }
00508 
00509     #Stored on the genomic_align_block
00510     my $ga_group_id = $this_genomic_align_block->group_id();
00511     my $ga_level_id = $this_genomic_align_block->level_id();
00512 
00513     my ($slice, $seq_name, $start, $end, $strand);
00514     if ($this_genomic_align_block->reference_slice) {
00515       $slice = $this_genomic_align_block->reference_slice;
00516       $seq_name = $this_genomic_align_block->reference_slice->name;
00517       $start = $this_genomic_align_block->reference_slice_start;
00518       $end = $this_genomic_align_block->reference_slice_end;
00519       $strand = $this_genomic_align_block->reference_slice_strand;
00520     } else {
00521       $slice = $consensus_genomic_align->get_Slice();
00522       $seq_name = $consensus_genomic_align->dnafrag->name;
00523       $start = $consensus_genomic_align->dnafrag_start;
00524       $end = $consensus_genomic_align->dnafrag_end;
00525       $strand = $consensus_genomic_align->dnafrag_strand;
00526     }
00527 
00528     my $f = Bio::EnsEMBL::DnaDnaAlignFeature->new_fast
00529       ({'cigar_string' => $ga_cigar_line,
00530         'slice'        => $slice,
00531         'seqname'      => $seq_name,
00532         'start'        => $start,
00533         'end'          => $end,
00534         'strand'       => $strand,
00535         'species'      => $consensus_genomic_align->genome_db->name,
00536         'score'        => $score,
00537         'percent_id'   => $perc_id,
00538         'hstart'       => $qdf_start + $ga_query_start - 1,
00539         'hend'         => $qdf_start + $ga_query_end -1,
00540         'hstrand'      => $ga_query_strand,
00541         'hseqname'     => $qdf_name,
00542         'hspecies'     => $query_genomic_align->genome_db->name,
00543         'hslice'       => $top_slice,
00544         'alignment_type' => $this_genomic_align_block->method_link_species_set->method_link_type,
00545         'group_id'     => $ga_group_id,
00546         'level_id'     => $ga_level_id,
00547         'strands_reversed' => $ga_strands_reversed});
00548 
00549     push @$dna_dna_align_features, $f;
00550   }
00551 
00552   return $dna_dna_align_features;
00553 }
00554 
00555 
00556 1;
00557 
00558