Inheritance diagram for Bio::EnsEMBL::BaseAlignFeature:

Public Member Functions
public Bio::EnsEMBL::BaseAlignFeature	new ()

public String	cigar_string ()

public String	align_type ()

public Int	alignment_length ()

protected Int	_ensembl_cigar_alignment_length ()

public List	ungapped_features ()

public Int	strands_reversed ()

public void	reverse_complement ()

protected void	_ensembl_reverse_complement ()

public Bio::EnsEMBL::BaseAlignFeature	transform ()

protected List	_parse_ensembl_cigar ()

protected	_parse_cigar ()

protected void	_parse_features ()

protected	_parse_ensembl_features ()

protected Int	_hit_unit ()

protected Int	_query_unit ()

protected Int	_mdtag_alignment_length ()

protected Array	_get_mdz_chunks ()

protected Array	_get_mdz_alignment_length ()

protected String	_get_mdz_chunk_type ()

protected Array	_mdz_alignment_string ()

Public Member Functions inherited from Bio::EnsEMBL::FeaturePair
public Bio::EnsEMBL::FeaturePair	new ()

public String	hseqname ()

public Int	hstart ()

public Int	hend ()

public	hstrand ()

public Bio::EnsEMBL::Slice	hslice ()

public String	hseq_region_name ()

public Boolean	hseq_region_strand ()

public Int	hseq_region_start ()

public Int	hseq_region_end ()

public Float	score ()

public Float	percent_id ()

public String	species ()

public String	hspecies ()

public String	coverage ()

public String	hcoverage ()

public String	external_db_id ()

public String	db_name ()

public String	db_display_name ()

public Float	p_value ()

public String	hdescription ()

public String	display_id ()

public Int	identical_matches ()

public Int	positive_matches ()

public Int	group_id ()

public Int	level_id ()

public void	invert ()

public	extra_data ()

public	type ()

Public Member Functions inherited from Bio::EnsEMBL::Feature
public Bio::EnsEMBL::Feature	new ()

public Int	start ()

public Int	end ()

public Int	strand ()

public void	move ()

public Int	length ()

public Bio::EnsEMBL::Analysis	analysis ()

public Bio::EnsEMBL::Slice	slice ()

public Boolean Or Undef	equals ()

public Bio::EnsEMBL::Feature	transform ()

public Bio::EnsEMBL::Feature	transfer ()

public Listref	project_to_slice ()

public Listref	project ()

public String	seqname ()

public String	display_id ()

public String	version ()

public Bio::EnsEMBL::Slice	feature_Slice ()

public String	seq_region_name ()

public Int	seq_region_length ()

public Boolean	seq_region_strand ()

public Int	seq_region_start ()

public Int	seq_region_end ()

public String	coord_system_name ()

public String	seq ()

public Listref	get_all_alt_locations ()

public Boolean	overlaps ()

public Boolean	overlaps_local ()

public List	get_overlapping_Genes ()

public Bio::EnsEMBL::Gene	get_nearest_Gene ()

public String	feature_so_acc ()

public String	feature_so_term ()

public Hashref	summary_as_hash ()

public String	species ()

public	sub_SeqFeature ()

public	add_sub_SeqFeature ()

public	flush_sub_SeqFeature ()

Public Member Functions inherited from Bio::EnsEMBL::Storable
public Bio::EnsEMBL::Storable	new ()

public Instance	new_fast ()

public Int	dbID ()

public Bio::EnsEMBL::DBSQL::BaseAdaptor	adaptor ()

public Boolean	is_stored ()

public	get_all_DAS_Features ()

Detailed Description

Synopsis

my $feat = new Bio::EnsEMBL::DnaPepAlignFeature(
  -slice        => $slice,
  -start        => 100,
  -end          => 120,
  -strand       => 1,
  -hseqname     => 'SP:RF1231',
  -hstart       => 200,
  -hend         => 220,
  -analysis     => $analysis,
  -cigar_string => '10M3D5M2I',
  -align_type   => 'ensembl'
);
 
where $analysis is a Bio::EnsEMBL::Analysis object.
 
Alternatively if you have an array of ungapped features:
 
  my $feat =
    new Bio::EnsEMBL::DnaPepAlignFeature( -features => \@features );
 
where @features is an array of Bio::EnsEMBL::FeaturePair objects.
 
There is a method to (re)create ungapped features from the cigar_string:
 
  my @ungapped_features = $feat->ungapped_features();
 
where @ungapped_features is an array of Bio::EnsEMBL::FeaturePair's.
 
Bio::EnsEMBL::BaseAlignFeature inherits from:
  Bio::EnsEMBL::FeaturePair, which in turn inherits from:
    Bio::EnsEMBL::Feature,
thus methods from both parent classes are available.
 
 
The cigar_string is a condensed representation of the matches and gaps
which make up the gapped alignment (where CIGAR stands for
Concise Idiosyncratic Gapped Alignment Report).
 
CIGAR format is: n <matches> [ x <deletes or inserts> m <matches> ]*
where M = match, D = delete, I = insert; n, m are match lengths;
x is delete or insert length.
 
Spaces are omitted, thus: "23M4I12M2D1M"
as are counts for any lengths of 1, thus 1M becomes M: "23M4I12M2DM"
 
 
To make things clearer this is how a blast HSP would be parsed:
 
>AK014066
       Length = 146
 
  Minus Strand HSPs:
 
  Score = 76 (26.8 bits), Expect = 1.4, P = 0.74
  Identities = 20/71 (28%), Positives = 29/71 (40%), Frame = -1
 
Query:   479 GLQAPPPTPQGCRLIPPPPLGLQAPLPTLRAVGSSHHHP*GRQGSSLSSFRSSLASKASA 300
             G  APPP PQG R   P P G + P   L             + + ++  R  +A   +
Sbjct:     7 GALAPPPAPQG-RWAFPRPTG-KRPATPLHGTARQDRQVRRSEAAKVTGCRGRVAPHVAP 64
 
Query:   299 SSPHNPSPLPS 267
                H P+P P+
Sbjct:    65 PLTHTPTPTPT 75
 
The alignment goes from 479 down to 267 in the query sequence on the reverse
strand, and from 7 to 75 in the subject sequence.
 
The alignment is made up of the following ungapped pieces:
 
query_seq start 447 , sbjct_seq hstart  7 , match length  33 , strand -1
query_seq start 417 , sbjct_seq hstart 18 , match length  27 , strand -1
query_seq start 267 , sbjct_seq hstart 27 , match length 147 , strand -1
 
When assembled into a DnaPepAlignFeature where:
  (seqname, start, end, strand) refer to the query sequence,
  (hseqname, hstart, hend, hstrand) refer to the subject sequence,
these ungapped pieces are represented by the cigar string:
  33M3I27M3I147M
with start 267, end 479, strand -1, and hstart 7, hend 75, hstrand 1.

Definition at line 90 of file BaseAlignFeature.pm.

Member Function Documentation

◆ _ensembl_cigar_alignment_length()

protected Int Bio::EnsEMBL::BaseAlignFeature::_ensembl_cigar_alignment_length ( )

  Arg [1]    : None
  Description: return the alignment length (including indels) based on the cigar_string
  Returntype : int
  Exceptions :
  Caller     :
  Status     : Stable

Code:

click to view

◆ _ensembl_reverse_complement()

protected void Bio::EnsEMBL::BaseAlignFeature::_ensembl_reverse_complement ( )

  Args       : none
  Description: reverse complement the FeaturePair for ensembl cigar string,
               modifing strand, hstrand and cigar_string in consequence
  Returntype : none
  Exceptions : none
  Caller     : general
  Status     : Stable

Code:

click to view

◆ _get_mdz_alignment_length()

protected Array Bio::EnsEMBL::BaseAlignFeature::_get_mdz_alignment_length ( )

  Arg [1]    : array of strings
  Description: calculate the alignment length from the given chunks
  Returntype : array of strings
  Exceptions : none
  Caller     : internal
  Status     : Stable

Code:

click to view

◆ _get_mdz_chunk_type()

protected String Bio::EnsEMBL::BaseAlignFeature::_get_mdz_chunk_type ( )

  Arg [1]    : char
  Description: get the chunk type
  Returntype : string
  Exceptions : none
  Caller     : internal
  Status     : Stable

Code:

click to view

◆ _get_mdz_chunks()

protected Array Bio::EnsEMBL::BaseAlignFeature::_get_mdz_chunks ( )

  Arg [1]    : mdtag string - MD Z String for mismatching positions. Regex : [0-9]+(([A-Z]|\^[A-Z]+)[0-9]+)* (Refer:  SAM/BAM specification)
  Description: parses the mdtag string and group it according the type
               eg: MD:Z:35^VIVALE31^GRPLIQPRRKKAYQLEHTFQGLLGKRSLFTE10 returns ['35', '^', 'VIVALE', '31', '^', 'GRPLIQPRRKKAYQLEHTFQGLLGKRSLFTE', '10']
  Returntype : array of strings
  Exceptions : none
  Caller     : internal
  Status     : Stable

Code:

click to view

◆ _hit_unit()

protected Int Bio::EnsEMBL::BaseAlignFeature::_hit_unit ( )

  Args       : none
  Description: abstract method, overwrite with something that returns
               one or three
  Returntype : int 1,3
  Exceptions : none
  Caller     : internal
  Status     : Stable

Code:

click to view

◆ _mdtag_alignment_length()

protected Int Bio::EnsEMBL::BaseAlignFeature::_mdtag_alignment_length ( )

  Arg [1]    : None
  Description: return the alignment length (including indels) based on the mdtag (mdz) string
  Returntype : int
  Exceptions : none
  Caller     : internal
  Status     : Stable

Code:

click to view

◆ _mdz_alignment_string()

protected Array Bio::EnsEMBL::BaseAlignFeature::_mdz_alignment_string ( )

  Arg [1]    : input sequence
  Arg [2]    : MD Z String for mismatching positions. Regex : [0-9]+(([A-Z]|\^[A-Z]+)[0-9]+)* (Refer:  SAM/BAM specification)
               eg: MD:Z:96^RHKTDSFVGLMGKRALNS0V14
  Example    :

$pf->alignment_strings

  Description: Allows to rebuild the alignment string of both the seq and hseq sequence
  Returntype : array reference containing 2 strings
               the first corresponds to seq
               the second corresponds to hseq
  Exceptions : none
  Caller     : general
  Status     : Stable

Code:

click to view

sub _mdz_alignment_string {
  my ($self, $input_seq, $mdz_string) = @_;
 
  my $chunks = $self->_get_mdz_chunks($mdz_string);
 
  my $target_seq = "";
  my $query_seq = "";
 
   #Handle first chunk, which is always a num
   my $first_chunk = $chunks->[0];
   my $first_chunk_type = $self->_get_mdz_chunk_type($first_chunk);
 
   my $offset = 0;
   if($first_chunk_type eq 'num'){
     if($first_chunk > 0){
       $target_seq = substr($input_seq, $offset, $first_chunk);
 
       #query and target are same at this point
       $query_seq = $target_seq;
       $offset = length($target_seq);
     }
   }else{
     die "First chunk should always be a num...something wrong\n";
   }
 
  for(my $i=1; $i < scalar(@$chunks); $i++){
 
    my $chunk = $chunks->[$i];
    my $chunk_type = $self->_get_mdz_chunk_type($chunk);
 
    if($chunk_type eq 'num'){
      #if 0, what follows next is a INSERT
      if($chunk == 0){
        my $insert_chunk = $chunks->[++$i];
        my $insert_chunk_type = $self->_get_mdz_chunk_type($insert_chunk);
 
        #insert_chunk_type is always alpha
        die if $insert_chunk_type eq "num";
        if($insert_chunk_type eq 'del'){
          $i--;
          next;
        }
        $target_seq = $target_seq . substr($input_seq, $offset, length($insert_chunk));
        $query_seq = $query_seq . $insert_chunk;
 
        $offset = $offset + length($insert_chunk);
 
      }else{
       #if non-0, then it is a match
        my $match_chunk = $chunks->[$i];
        my $match_chunk_type = $self->_get_mdz_chunk_type($match_chunk);
 
         #$match_chunk_type is always num
        die if $match_chunk_type ne "num";
 
        my $match_target_string = substr($input_seq, $offset, $match_chunk);
 
        $target_seq = $target_seq . $match_target_string;
        $query_seq = $query_seq . $match_target_string;
 
        $offset = $offset + $match_chunk;
      }
 
    }elsif($chunk_type eq 'alpha'){
        my $insert_chunk = $chunk;
        my $insert_chunk_type = $self->_get_mdz_chunk_type($insert_chunk);
        die if $insert_chunk_type ne "alpha";
 
        $target_seq = $target_seq . substr($input_seq, $offset, length($insert_chunk));
        $query_seq = $query_seq . $insert_chunk;
 
        $offset = $offset + length($insert_chunk);
 
    }elsif($chunk_type eq 'del'){
      #get next chunk which contains the deleted sequence
      my $del_chunk = $chunks->[++$i];
      my $del_chunk_type = $self->_get_mdz_chunk_type($del_chunk);
 
      #del_chunk_type is always alpha
      die if $del_chunk_type ne 'alpha';
 
      $target_seq = $target_seq . "-" x length($del_chunk);
      $query_seq = $query_seq . $del_chunk;
 
      #no change in offset
      $offset = $offset;
 
 
    }
 
  }
 
  return [$target_seq, $query_seq];
 
}

◆ _parse_cigar()

protected Bio::EnsEMBL::BaseAlignFeature::_parse_cigar ( )

Undocumented method

Code:

click to view

◆ _parse_ensembl_cigar()

protected List Bio::EnsEMBL::BaseAlignFeature::_parse_ensembl_cigar ( )

  Args       : none
  Description: PRIVATE (internal) method - creates ungapped features from 
               internally stored cigar line in ensembl format
  Returntype : list of Bio::EnsEMBL::FeaturePair
  Exceptions : none
  Caller     : ungapped_features
  Status     : Stable

Code:

click to view

sub _parse_ensembl_cigar {
  my ( $self ) = @_;
 
  my $query_unit = $self->_query_unit();
  my $hit_unit = $self->_hit_unit();
 
  my $string = $self->{'cigar_string'};
 
  throw("No cigar string defined in object") if(!defined($string));
 
  my @pieces = ( $string =~ /(\d*[MDI])/g );
 
  my @features;
  my $strand1 = $self->{'strand'} || 1;
  my $strand2 = $self->{'hstrand'}|| 1;
 
  my ( $start1, $start2 );
 
  if( $strand1 == 1 ) {
    $start1 = $self->{'start'};
  } else {
    $start1 = $self->{'end'};
  }
 
  if( $strand2 == 1 ) {
    $start2 = $self->{'hstart'};
  } else {
    $start2 = $self->{'hend'};
  }
 
  #
  # Construct ungapped blocks as FeaturePairs objects for each MATCH
  #
  foreach my $piece (@pieces) {
 
    my ($length) = ( $piece =~ /^(\d*)/ );
    if( $length eq "" ) { $length = 1 }
    my $mapped_length;
 
    # explicit if statements to avoid rounding problems
    # and make sure we have sane coordinate systems
    if( $query_unit == 1 && $hit_unit == 3 ) {
      $mapped_length = $length*3;
    } elsif( $query_unit == 3 && $hit_unit == 1 ) {
      $mapped_length = $length / 3;
    } elsif ( $query_unit == 1 && $hit_unit == 1 ) {
      $mapped_length = $length;
    } else {
      throw("Internal error $query_unit $hit_unit, currently only " .
            "allowing 1 or 3 ");
    }
 
    if( int($mapped_length) != $mapped_length and
        ($piece =~ /M$/ or $piece =~ /D$/)) {
      throw("Internal error with mismapped length of hit, query " .
            "$query_unit, hit $hit_unit, length $length");
    }
 
    if( $piece =~ /M$/ ) {
      #
      # MATCH
      #
      my ( $qstart, $qend);
      if( $strand1 == 1 ) {
        $qstart = $start1;
        $qend = $start1 + $length - 1;
        $start1 = $qend + 1;
      } else {
        $qend = $start1;
        $qstart = $start1 - $length + 1;
        $start1 = $qstart - 1;
      }
 
      my ($hstart, $hend);
      if( $strand2 == 1 ) {
        $hstart = $start2;
        $hend = $start2 + $mapped_length - 1;
        $start2 = $hend + 1;
      } else {
        $hend = $start2;
        $hstart = $start2 - $mapped_length + 1;
        $start2 = $hstart - 1;
      }
 
 
      push @features, Bio::EnsEMBL::FeaturePair->new
        (-SLICE      => $self->{'slice'},
         -SEQNAME   => $self->{'seqname'},
         -START      => $qstart,
         -END        => $qend,
         -STRAND     => $strand1,
         -HSLICE      => $self->{'hslice'},
         -HSEQNAME   => $self->{'hseqname'},
         -HSTART     => $hstart,
         -HEND       => $hend,
         -HSTRAND    => $strand2,
         -SCORE      => $self->{'score'},
         -PERCENT_ID => $self->{'percent_id'},
         -ANALYSIS   => $self->{'analysis'},
         -P_VALUE    => $self->{'p_value'},
         -EXTERNAL_DB_ID => $self->{'external_db_id'},
         -HCOVERAGE   => $self->{'hcoverage'},
         -GROUP_ID    => $self->{'group_id'},
         -LEVEL_ID    => $self->{'level_id'});
      
 
      # end M cigar bits 
    } elsif( $piece =~ /I$/ ) {
      #
      # INSERT
      #
      if( $strand1 == 1 ) {
        $start1 += $length;
      } else {
        $start1 -= $length;
      }
    } elsif( $piece =~ /D$/ ) {
      #
      # DELETION
      #
      if( $strand2 == 1 ) {
        $start2 += $mapped_length;
      } else {
        $start2 -= $mapped_length;
      }
    } else {
      throw( "Illegal cigar line $string!" );
    }
  }
 
  return \@features;
}

◆ _parse_ensembl_features()

protected Bio::EnsEMBL::BaseAlignFeature::_parse_ensembl_features ( )

Undocumented method

Code:

click to view

sub _parse_ensembl_features {
  my ($self,$features ) = @_;
 
  if (ref($features) ne "ARRAY") {
    throw("features must be an array reference not a [".ref($features)."]");
  } elsif (scalar(@$features) == 0) {
    throw("features array must not be empty");
  }
 
  my $query_unit = $self->_query_unit();
  my $hit_unit = $self->_hit_unit();
 
  my $strand  = $features->[0]->strand;
 
  throw ('FeaturePair needs to have strand == 1 or strand == -1') if(!$strand);
 
  my @f;
 
  #
  # Sort the features on their start position
  # Ascending order on positive strand, descending on negative strand
  #
  if( $strand == 1 ) {
    @f = sort {$a->start <=> $b->start} @$features;
  } else {
    @f = sort { $b->start <=> $a->start} @$features;
  }
 
  my $hstrand     = $f[0]->hstrand;
  my $slice       = $f[0]->slice();
  my $hslice       = $f[0]->hslice();
  my $name        = $slice ? $slice->name() : undef;
  my $hname       = $f[0]->hseqname;
  my $score       = $f[0]->score;
  my $percent     = $f[0]->percent_id;
  my $analysis    = $f[0]->analysis;
  my $pvalue      = $f[0]->p_value();
  my $external_db_id = $f[0]->external_db_id;
  my $hcoverage   = $f[0]->hcoverage;
  my $group_id    = $f[0]->group_id;
  my $level_id    = $f[0]->level_id;
 
  my $seqname = $f[0]->seqname;
  # implicit strand 1 for peptide sequences
  $strand  ||= 1;
  $hstrand ||= 1;
  my $ori = $strand * $hstrand;
 
  throw("No features in the array to parse") if(scalar(@f) == 0);
 
  my $prev1; # where last feature q part ended
  my $prev2; # where last feature s part ended
 
  my $string;
 
  # Use strandedness info of query and hit to make sure both sets of 
  # start and end  coordinates are oriented the right way around.
  my $f1start;
  my $f1end;
  my $f2end;
  my $f2start;
 
  if ($strand == 1) {
    $f1start = $f[0]->start;
    $f1end = $f[-1]->end;
  } else {
    $f1start = $f[-1]->start;
    $f1end = $f[0]->end;
  }
 
  if ($hstrand == 1) {
    $f2start = $f[0]->hstart;
    $f2end = $f[-1]->hend;
  } else {
    $f2start = $f[-1]->hstart;
    $f2end = $f[0]->hend;
  }
 
  #
  # Loop through each portion of alignment and construct cigar string
  #
 
  foreach my $f (@f) {
    #
    # Sanity checks
    #
 
    if (!$f->isa("Bio::EnsEMBL::FeaturePair")) {
      throw("Array element [$f] is not a Bio::EnsEMBL::FeaturePair");
    }
    if( defined($f->hstrand()) && $f->hstrand() != $hstrand ) {
      throw("Inconsistent hstrands in feature array");
    }
    if( defined($f->strand()) && ($f->strand != $strand)) {
      throw("Inconsistent strands in feature array");
    }
    if ( defined($name) && $name ne $f->slice->name()) {
      throw("Inconsistent names in feature array [$name - ".
            $f->slice->name()."]");
    }
    if ( defined($hname) && $hname ne $f->hseqname) {
      throw("Inconsistent hit names in feature array [$hname - ".
            $f->hseqname . "]");
    }
    if ( defined($score) && $score ne $f->score) {
      throw("Inconsisent scores in feature array [$score - " .
            $f->score . "]");
    }
    if (defined($f->percent_id) && $percent ne $f->percent_id) {
      throw("Inconsistent pids in feature array [$percent - " .
            $f->percent_id . "]");
    }
    if(defined($pvalue) && $pvalue != $f->p_value()) {
      throw("Inconsistant p_values in feature arraw [$pvalue " .
            $f->p_value() . "]");
    }
    if($seqname && $seqname ne $f->seqname){
      throw("Inconsistent seqname in feature array [$seqname - ".
            $f->seqname . "]");
    }
    my $start1 = $f->start;      #source sequence alignment start
    my $start2 = $f->hstart();   #hit sequence alignment start
 
    #
    # More sanity checking
    #
    if (defined($prev1)) {
      if ( $strand == 1 ) {
        if ($f->start < $prev1) {
          throw("Inconsistent coords in feature array (forward strand).\n" .
               "Start [".$f->start()."] in current feature should be greater " .
           "than previous feature end [$prev1].");
        }
      } else {
        if ($f->end > $prev1) {
          throw("Inconsistent coords in feature array (reverse strand).\n" .
                "End [".$f->end() ."] should be less than previous feature " .
                "start [$prev1].");
        }
      }
    }
 
    my $length = ($f->end - $f->start + 1);    #length of source seq alignment
    my $hlength = ($f->hend - $f->hstart + 1); #length of hit seq alignment
 
    # using multiplication to avoid rounding errors, hence the
    # switch from query to hit for the ratios
 
    #
    # Yet more sanity checking
    #
    if($query_unit > $hit_unit){
      # I am going to make the assumption here that this situation will 
      # only occur with DnaPepAlignFeatures, this may not be true
      my $query_p_length = sprintf "%.0f", ($length/$query_unit);
      my $hit_p_length = sprintf "%.0f", ($hlength * $hit_unit);
      if( $query_p_length != $hit_p_length) {
        throw( "Feature lengths not comparable Lengths:" .$length .
               " " . $hlength . " Ratios:" . $query_unit . " " .
               $hit_unit );
      }
    } else{
      my $query_d_length = sprintf "%.0f", ($length*$hit_unit);
      my $hit_d_length = sprintf "%.0f", ($hlength * $query_unit);
      if( $length * $hit_unit != $hlength * $query_unit ) {
        throw( "Feature lengths not comparable Lengths:" . $length .
               " " . $hlength . " Ratios:" . $query_unit . " " .
               $hit_unit );
      }
    }
 
    my $hlengthfactor = ($query_unit/$hit_unit);
 
    #
    # Check to see if there is an I type (insertion) gap:
    #   If there is a space between the end of the last source sequence 
    #   alignment and the start of this one, then this is an insertion
    #
 
    my $insertion_flag = 0;
    if( $strand == 1 ) {
      if( ( defined $prev1 ) && ( $f->start > $prev1 + 1  )) {
 
        #there is an insertion
        $insertion_flag = 1;
        my $gap = $f->start - $prev1 - 1;
        if( $gap == 1 ) {
          $gap = ""; # no need for a number if gap length is 1
        }
        $string .= "$gap"."I";
 
      }
 
      #shift our position in the source seq alignment
      $prev1 = $f->end();
    } else {
 
      if(( defined $prev1 ) && ($f->end + 1 < $prev1 )) {
 
        #there is an insertion
        $insertion_flag = 1;
        my $gap = $prev1 - $f->end() - 1;
        if( $gap == 1 ) {
          $gap = ""; # no need for a number if gap length is 1
        }
        $string .= "$gap"."I";
      }
 
      #shift our position in the source seq alignment
      $prev1 = $f->start();
    }
 
    #
    # Check to see if there is a D type (deletion) gap
    #   There is a deletion gap if there is a space between the end of the
    #   last portion of the hit sequence alignment and this one
    #
    if( $hstrand == 1 ) {
      if((  defined $prev2 ) && ( $f->hstart() > $prev2 + 1 )) {
 
        #there is a deletion
        my $gap = $f->hstart - $prev2 - 1;
        my $gap2 = int( $gap * $hlengthfactor + 0.5 );
    
        if( $gap2 == 1 ) {
          $gap2 = "";  # no need for a number if gap length is 1
        }
        $string .= "$gap2"."D";
 
        #sanity check,  Should not be an insertion and deletion
        if($insertion_flag) {
          if ($message_only_once) {
            warning("Should not be an deletion and insertion on the " .
                    "same alignment region. cigar_line=$string\n");
            $message_only_once = 0;
          }
        }
      }
      #shift our position in the hit seq alignment
      $prev2 = $f->hend();
 
     } else {
      if( ( defined $prev2 ) && ( $f->hend() + 1 < $prev2 )) {
 
        #there is a deletion
        my $gap = $prev2 - $f->hend - 1;
        my $gap2 = int( $gap * $hlengthfactor + 0.5 );
    
        if( $gap2 == 1 ) {
          $gap2 = "";  # no need for a number if gap length is 1
        }
        $string .= "$gap2"."D";
 
        #sanity check,  Should not be an insertion and deletion
        if($insertion_flag) {
          if ($message_only_once) {
            warning("Should not be an deletion and insertion on the " .
                    "same alignment region. prev2 = $prev2; f->hend() = " .
                    $f->hend() . "; cigar_line = $string;\n");
            $message_only_once = 0;
          }
        }
      }
      #shift our position in the hit seq alignment
      $prev2 = $f->hstart();
    }
 
    my $matchlength = $f->end() - $f->start() + 1;
    if( $matchlength == 1 ) {
      $matchlength = "";
    }
    $string .= $matchlength."M";
  }
 
  $self->{'start'}      = $f1start;
  $self->{'end'}        = $f1end;
  $self->{'seqname'}    = $seqname;
  $self->{'strand'}     = $strand;
  $self->{'score'}      = $score;
  $self->{'percent_id'} = $percent;
  $self->{'analysis'}   = $analysis;
  $self->{'slice'}      = $slice;
  $self->{'hslice'}     = $hslice;
  $self->{'hstart'}     = $f2start;
  $self->{'hend'}       = $f2end;
  $self->{'hstrand'}    = $hstrand;
  $self->{'hseqname'}   = $hname;
  $self->{'cigar_string'} = $string;
  $self->{'p_value'}      = $pvalue;
  $self->{'external_db_id'} = $external_db_id; 
  $self->{'hcoverage'}    = $hcoverage; 
  $self->{'group_id'}     = $group_id;
  $self->{'level_id'}     = $level_id;
}

◆ _parse_features()

protected void Bio::EnsEMBL::BaseAlignFeature::_parse_features ( )

  Arg  [1]   : listref Bio::EnsEMBL::FeaturePair $ungapped_features
  Description: creates internal cigar_string and start,end hstart,hend
               entries.
  Returntype : none, fills in values of self
  Exceptions : argument list undergoes many sanity checks - throws under many
               invalid conditions
  Caller     : new
  Status     : Stable

Code:

click to view

◆ _query_unit()

protected Int Bio::EnsEMBL::BaseAlignFeature::_query_unit ( )

  Args       : none
  Description: abstract method, overwrite with something that returns
               one or three
  Returntype : int 1,3
  Exceptions : none
  Caller     : internal
  Status     : Stable

Code:

click to view

◆ align_type()

public String Bio::EnsEMBL::BaseAlignFeature::align_type ( )

  Arg [1]    : type $align_type
  Example    :

$feature->align_type( "ensembl" );

  Description: get/set for attribute align_type.
               align_type specifies which cigar string 
               is used to describe the alignment:
               The default is 'ensembl'
  Returntype : string
  Exceptions : none
  Caller     : general
  Status     : Stable

Code:

click to view

◆ alignment_length()

public Int Bio::EnsEMBL::BaseAlignFeature::alignment_length ( )

  Arg [1]    : None
  Description: return the alignment length (including indels) based on the alignment_type ('ensembl', 'mdtag')
  Returntype : int
  Exceptions : 
  Caller     : 
  Status     : Stable

Code:

click to view

◆ cigar_string()

public String Bio::EnsEMBL::BaseAlignFeature::cigar_string ( )

  Arg [1]    : string $cigar_string
  Example    :

$feature->cigar_string( "12MI3M" );

  Description: get/set for attribute cigar_string.
               cigar_string describes the alignment:
                 "xM" stands for x matches (or mismatches),
                 "xI" for x inserts into the query sequence,
                 "xD" for x deletions from the query sequence
                 where the query sequence is specified by (seqname, start, ...)
                 and the subject sequence by (hseqname, hstart, ...).
               An "x" that is 1 can be omitted.
               See the SYNOPSIS for an example.
  Returntype : string
  Exceptions : none
  Caller     : general
  Status     : Stable

Code:

click to view

◆ new()

public Bio::EnsEMBL::BaseAlignFeature Bio::EnsEMBL::BaseAlignFeature::new ( )

  Arg [..]   : List of named arguments. (-cigar_string , -features, -align_type) defined
               in this constructor, others defined in FeaturePair and 
               SeqFeature superclasses.  Either cigar_string or a list
               of ungapped features should be provided - not both.
  Example    :

$baf = new BaseAlignFeatureSubclass(-cigar_string => '3M3I12M', -align_type => 'ensembl');

  Description: Creates a new BaseAlignFeature using either a cigar string or
               a list of ungapped features.  BaseAlignFeature is an abstract
               baseclass and should not actually be instantiated - rather its
               subclasses should be.
  Returntype : Bio::EnsEMBL::BaseAlignFeature
  Exceptions : thrown if both feature and cigar string args are provided
               thrown if neither feature nor cigar string args are provided
               warn if cigar string is provided without cigar type
  Caller     : general
  Status     : Stable

Code:

click to view

◆ reverse_complement()

public void Bio::EnsEMBL::BaseAlignFeature::reverse_complement ( )

  Args       : none
  Description: reverse complement the FeaturePair based on the cigar type
               modifing strand, hstrand and cigar_string in consequence
  Returntype : none
  Exceptions : none
  Caller     : general
  Status     : Stable

Code:

click to view

◆ strands_reversed()

public Int Bio::EnsEMBL::BaseAlignFeature::strands_reversed ( )

 
  Arg [1]    : int $strands_reversed
  Description: get/set for attribute strands_reversed
               0 means that strand and hstrand are the original strands obtained
                 from the alignment program used
               1 means that strand and hstrand have been flipped as compared to
                 the original result provided by the alignment program used.
                 You may want to use the reverse_complement method to restore the
                 original strandness.
  Returntype : int
  Exceptions : none
  Caller     : general
  Status     : Stable

Code:

click to view

◆ transform()

public Bio::EnsEMBL::BaseAlignFeature Bio::EnsEMBL::BaseAlignFeature::transform ( )

  Arg  1     : String $coordinate_system_name
  Arg [2]    : String $coordinate_system_version
  Example    :

$feature = $feature->transform('contig');

$feature = $feature->transform('chromosome', 'NCBI33');

  Description: Moves this AlignFeature to the given coordinate system.
               If the feature cannot be transformed to the destination 
               coordinate system undef is returned instead.
  Returntype : Bio::EnsEMBL::BaseAlignFeature;
  Exceptions : wrong parameters
  Caller     : general
  Status     : Medium Risk

Code:

click to view

◆ ungapped_features()

public List Bio::EnsEMBL::BaseAlignFeature::ungapped_features ( )

  Args       : none
  Example    :

@ungapped_features = $align_feature->get_feature

  Description: converts the internal cigar_string into an array of
               ungapped feature pairs
  Returntype : list of Bio::EnsEMBL::FeaturePair
  Exceptions : cigar_string not set internally
  Caller     : general
  Status     : Stable

Code:

click to view

The documentation for this class was generated from the following file:

ensembl/modules/Bio/EnsEMBL/BaseAlignFeature.pm

Public Member Functions

Detailed Description

Synopsis

Member Function Documentation

◆ _ensembl_cigar_alignment_length()

◆ _ensembl_reverse_complement()

◆ _get_mdz_alignment_length()

◆ _get_mdz_chunk_type()

◆ _get_mdz_chunks()

◆ _hit_unit()

◆ _mdtag_alignment_length()

◆ _mdz_alignment_string()

◆ _parse_cigar()

◆ _parse_ensembl_cigar()

◆ _parse_ensembl_features()

◆ _parse_features()

◆ _query_unit()

◆ align_type()

◆ alignment_length()

◆ cigar_string()

◆ new()

◆ reverse_complement()

◆ strands_reversed()

◆ transform()

◆ ungapped_features()