Inheritance diagram for Bio::EnsEMBL::Compara::AlignedMemberSet:

Class Summary

Description

AlignedMemberSet is a specialized version of MemberSet for members
that have a sequence and that are aligned (i.e. AlignedMember).

It is the actual base class for Family, Homology and GeneTree.
It only adds methods to deal with the alignment.

Since SeqMember can hold several sequences per member (seq_type),
AlignedMemberSet can also do so. An AlignedMemberSet defaults on a
given seq_type, but can also be mapped to another one. For instance,
the alignments of protein sequences can be transformed into alignments
of the initial CDS sequences.

Synopsis

Global properties of the alignment:
 - seq_type()
 - aln_length()
 - aln_method()
 - consensus_cigar_line()

I/O:
 - load_cigars_from_file()
 - print_alignment_to_file()
 - get_SimpleAlign()
 - get_4D_SimpleAlign()

Definition at line 45 of file AlignedMemberSet.pm.

Available Methods
protected Array	_attr_to_copy_list ()
protected Array	_get_Member ()
protected String	_toString ()
public Bio::EnsEMBL::DBSQL::BaseAdaptor	adaptor ()
public void	add_Member ()
public Int	aln_length ()
public String	aln_method ()
public	clear ()
public String	consensus_cigar_line ()
public Int	dbID ()
public Bio::EnsEMBL::Compara::MemberSet	deep_copy ()
public String	description ()
public Array	gene_list ()
public Bio::SimpleAlign	get_4D_SimpleAlign ()
public	get_all_DAS_Features ()
public Array	get_all_GeneMembers ()
public Array	get_all_GenomeDBs_by_member_source_name ()
public Array	get_all_Members ()
public Array	get_all_taxa_by_member_source_name ()
public Array	get_Member_by_GenomeDB ()
public Array	get_Member_by_source ()
public Array	get_Member_by_source_GenomeDB ()
public Array	get_Member_by_source_taxon ()
public Bio::SimpleAlign	get_SimpleAlign ()
public Boolean	has_species_by_name ()
public Boolean	is_stored ()
public void	load_cigars_from_file ()
public Int	Member_count_by_GenomeDB ()
public Int	Member_count_by_source ()
public Int	Member_count_by_source_GenomeDB ()
public Int	Member_count_by_source_taxon ()
public Bio::EnsEMBL::Compara::MethodLinkSpeciesSet	method_link_species_set ()
public Int	method_link_species_set_id ()
public Bio::EnsEMBL::Compara::AlignedMemberSet	new ()
public Instance	new_fast ()
public Bio::SimpleAlign	print_alignment_to_file ()
public Number	print_sequences_to_file ()
public void	remove_Member ()
public String	seq_type ()
public String	stable_id ()
public String	toString ()
public String	update_alignment_stats ()
public Int	version ()
Public Attributes
public	String: member_class( )

Method Documentation

protected Array Bio::EnsEMBL::Compara::AlignedMemberSet::_attr_to_copy_list ( )

  Description: Returns the list of all the attributes to be copied by deep_copy()
  Returntype : Array of String
  Caller     : General

Code:

click to view

Reimplemented from Bio::EnsEMBL::Compara::MemberSet.

Reimplemented in Bio::EnsEMBL::Compara::Family, and Bio::EnsEMBL::Compara::GeneTree.

protected Array Bio::EnsEMBL::Compara::MemberSet::_get_Member ( ) [inherited]

  Arg [1]    : string $scope
  Arg [2]    : string $key
  Example    :

 $family->_get_Member('_members_by_source', 'ENSEMBLPEP')

  Description: Used as the generic reference point for all 
               get_Memeber_by* methods. The method searches the given
               scope & if the values cannot be found will initialize that value
               to an empty array reference.
  Returntype : array reference of Bio::EnsEMBL::Compara::Member
  Exceptions : None.
  Caller     : internal

Code:

click to view

protected String Bio::EnsEMBL::Compara::AlignedMemberSet::_toString ( )

  Description : Helper method for MemberSet::toString that provides class-specific information
  Returntype  : String
  Exceptions  : none
  Caller      : internal

Code:

click to view

Reimplemented from Bio::EnsEMBL::Compara::MemberSet.

Reimplemented in Bio::EnsEMBL::Compara::GeneTree.

public void Bio::EnsEMBL::Compara::MemberSet::add_Member ( ) [inherited]

  Arg [1]    : Bio::EnsEMBL::Compara::Member $member
  Example    :

 $gene_tree->add_Member($member);

  Description: Add a new Member to this set
  Returntype : none
  Exceptions : Throws if input objects don't check
  Caller     : general
  Status     : Stable

Code:

click to view

Reimplemented in Bio::EnsEMBL::Compara::GeneTree.

public Int Bio::EnsEMBL::Compara::AlignedMemberSet::aln_length ( )

  Description : Getter/Setter for the aln_length field. This field contains
                the length of the alignment
  Returntype  : Integer
  Example     :

 my $len = $tree->aln_length();

  Caller      : General

Code:

click to view

public String Bio::EnsEMBL::Compara::AlignedMemberSet::aln_method ( )

  Description : Getter/Setter for the aln_method field. This field should
                represent the method used for the alignment
  Returntype  : String
  Example     :

 my $method = $tree->aln_method();

  Caller      : General

Code:

click to view

public Bio::EnsEMBL::Compara::MemberSet::clear ( ) [inherited]

  Arg [1]    : None
  Description: Clears the list of members

Code:

click to view

public String Bio::EnsEMBL::Compara::AlignedMemberSet::consensus_cigar_line ( )

  Example    :

 my $consensus_cigar = $gene_tree->consensus_cigar_line();

  Description: Creates a consensus cigar line for all the members of the
               set. See Bio::EnsEMBL::Compara::AlignedMemberSet
  Returntype : string
  Caller     : general

Code:

click to view

public Bio::EnsEMBL::Compara::MemberSet Bio::EnsEMBL::Compara::MemberSet::deep_copy ( ) [inherited]

  Description: Returns a copy of $self. All the members are themselves copied
  Returntype : Bio::EnsEMBL::Compara::MemberSet
  Caller     : general
  Status     : Stable

Code:

click to view

public String Bio::EnsEMBL::Compara::MemberSet::description ( ) [inherited]

  Arg [1]    : string $description (optional)
  Example    :

  Description: Getter/setter for the description corresponding to this relation
  Returntype : string
  Exceptions : none
  Caller     : general
  Status     : Stable

Code:

click to view

public Array Bio::EnsEMBL::Compara::MemberSet::gene_list ( ) [inherited]

  Example    :

 my $pair = $homology->gene_list

  Description: return the pair of members for the homology
  Returntype : array ref of (2) Bio::EnsEMBL::Compara::Member objects
  Caller     : general

Code:

click to view

public Bio::SimpleAlign Bio::EnsEMBL::Compara::AlignedMemberSet::get_4D_SimpleAlign ( )

  Example    :

 $4d_align = $homology->get_4D_SimpleAlign();

  Description: get 4 times degenerate positions pairwise simple alignment
  Returntype : Bio::SimpleAlign
  Caller     : general

Code:

click to view

sub get_4D_SimpleAlign {
    my $self = shift;
    my $keep_gaps = shift;
    my $aa_must_be_identical = shift;

    Bio::EnsEMBL::Compara::Utils::Preloader::load_all_sequences($self->adaptor->db->get_SequenceAdaptor, 'cds', $self) if $self->adaptor;

    my $sa = Bio::SimpleAlign->new();

    my %member_seqstr;
    foreach my $member (@{$self->get_all_Members}) {
        # Only peptides can have a CDS sequence
        next unless $member->source_name =~ /PEP$/;
        my $seqstr = $member->alignment_string('cds');
        next if(!$seqstr);
        #print STDERR $seqstr,"\n";
        $seqstr =~ s/(.{3})/$1 /g;
        my @tmp_tab =  split /\s+/, $seqstr;
        #print STDERR "tnp_tab 0: ", $tmp_tab[0],"\n";
        $member_seqstr{$member->stable_id} = \@tmp_tab;
    }

    my $seqstr_length;
    foreach my $seqid (keys %member_seqstr) {
        unless (defined $seqstr_length) {
            #print STDERR $member_seqstr{$seqid}->[0],"\n";
            $seqstr_length = scalar @{$member_seqstr{$seqid}};
            next;
        }
        unless ($seqstr_length == scalar @{$member_seqstr{$seqid}}) {
            die "Length of dna alignment are not the same, $seqstr_length and " . scalar @{$member_seqstr{$seqid}} ." respectively for homology_id " . $self->dbID . "\n";
        }
    }

    my %FourD_member_seqstr;
    for (my $i=0; $i < $seqstr_length; $i++) {
        my $FourD_codon = 1;
        my $FourD_aminoacid;
        foreach my $seqid (keys %member_seqstr) {
            if (defined $FOURD_CODONS{$member_seqstr{$seqid}->[$i]} or ($keep_gaps and $member_seqstr{$seqid}->[$i] eq '---')) {
                #print STDERR "YES ",$FOURD_CODONS{$member_seqstr{$seqid}->[$i]}," ",$member_seqstr{$seqid}->[$i],"\n";
                next unless $aa_must_be_identical;

                if (defined $FourD_aminoacid && $FourD_aminoacid eq $FOURD_CODONS{$member_seqstr{$seqid}->[$i]}) {
                    #print STDERR "YES ",$FOURD_CODONS{$member_seqstr{$seqid}->[$i]}," ",$member_seqstr{$seqid}->[$i],"\n";
                    next;
                } elsif (defined $FourD_aminoacid) {
                    #print STDERR "NO ",$FOURD_CODONS{$member_seqstr{$seqid}->[$i]}," ",$member_seqstr{$seqid}->[$i],"\n";
                    $FourD_codon = 0;
                    last;
                } else {
                    $FourD_aminoacid = $FOURD_CODONS{$member_seqstr{$seqid}->[$i]};
                    #print STDERR $FOURD_CODONS{$member_seqstr{$seqid}->[$i]}," ",$member_seqstr{$seqid}->[$i]," ";
                }
                next;
            } else {
                #print STDERR "NO ",$CODONS{$member_seqstr{$seqid}->[$i]}," ",$member_seqstr{$seqid}->[$i],"\n";
                $FourD_codon = 0;
                last;
            }
        }
        next unless ($FourD_codon);
        foreach my $seqid (keys %member_seqstr) {
            $FourD_member_seqstr{$seqid} .= substr($member_seqstr{$seqid}->[$i],2,1);
        }
    }

    foreach my $seqid (keys %FourD_member_seqstr) {

        my $seq = Bio::LocatableSeq->new(
                -SEQ    => $FourD_member_seqstr{$seqid},
                -START  => 1,
                -END    => length($FourD_member_seqstr{$seqid}),
                -ID     => $seqid,
                -STRAND => 0
        );

        $sa->add_seq($seq);
    }

    $sa = $sa->remove_gaps(undef, 1);

    return $sa;
}

public Array Bio::EnsEMBL::Compara::MemberSet::get_all_GeneMembers ( ) [inherited]

  Arg [1]    : (optional) genome db id 
  Example    :

 my $gene_members = $ortholog->get_all_GeneMembers($genome_db->dbID')

  Example    :

 my $gene_members = $ortholog->get_all_GeneMembers($genome_db->dbID)

  Description: 
  Returntype : array reference of Bio::EnsEMBL::Compara::Member
  Exceptions : 
  Caller     : public

Code:

click to view

public Array Bio::EnsEMBL::Compara::MemberSet::get_all_GenomeDBs_by_member_source_name ( ) [inherited]

  Arg [1]    : string $source_name
               e.g. "ENSEMBLPEP"
  Example    :

  Description: Returns the distinct GenomeDBs found in this family. Please note
               that if you specify a source other than an EnsEMBL based one
               the chances of getting back GenomeDBs are very low.
  Returntype : array reference of distinct Bio::EnsEMBL::Compara::GenomeDB 
               objects found in this family
  Exceptions : 
  Caller     : public

Code:

click to view

public Array Bio::EnsEMBL::Compara::MemberSet::get_all_Members ( ) [inherited]

  Example    :

  Description: Returns all the members in this set
  Returntype : array reference of Bio::EnsEMBL::Compara::Member
  Exceptions : 
  Caller     :

Code:

click to view

Reimplemented in Bio::EnsEMBL::Compara::GeneTree.

public Array Bio::EnsEMBL::Compara::MemberSet::get_all_taxa_by_member_source_name ( ) [inherited]

  Arg [1]    : string $source_name
               e.g. "ENSEMBLPEP"
  Example    :

  Description: Returns the distinct taxons found in this family across
               the specified source. If you do not specify a source then
               the code will return all taxons in this family.
  Returntype : array reference of distinct Bio::EnsEMBL::Compara::NCBITaxon 
               objects found in this family
  Exceptions : 
  Caller     : public

Code:

click to view

public Array Bio::EnsEMBL::Compara::MemberSet::get_Member_by_GenomeDB ( ) [inherited]

  Arg [1]    : Bio::EnsEMBL::Compara::GenomeDB $genome_db
  Example    :

 $family->get_Member_by_GenomeDB($genome_db)

  Description: Returns all [Member] entries linked to this GenomeDB. 
               This will only return EnsEMBL based entries since UniProtKB 
               entries are not linked to a GenomeDB.
  Returntype : array reference of Bio::EnsEMBL::Compara::Member
  Exceptions : If input is undefined & genome db is not of expected type
  Caller     : public

Code:

click to view

public Array Bio::EnsEMBL::Compara::MemberSet::get_Member_by_source ( ) [inherited]

  Arg [1]    : string $source_name
               e.g. "ENSEMBLPEP"
  Example    :

  Description: 
  Returntype : array reference of Bio::EnsEMBL::Compara::Member
  Exceptions : 
  Caller     : public

Code:

click to view

public Array Bio::EnsEMBL::Compara::MemberSet::get_Member_by_source_GenomeDB ( ) [inherited]

  Arg [1]    : string $source_name
  Arg [2]    : Bio::EnsEMBL::Compara::GenomeDB $genome_db
  Example    :

 $family->get_Member_by_source_taxon('ENSEMBLPEP', $genome_db)

  Description: Returns all [Member] entries linked to this GenomeDB
               and the given source_name. This will only return EnsEMBL based 
               entries since UniProtKB entries are not linked to a GenomeDB.
  Returntype : array reference of Bio::EnsEMBL::Compara::Member
  Exceptions : If input is undefined & genome db is not of expected type
  Caller     : public

Code:

click to view

public Array Bio::EnsEMBL::Compara::MemberSet::get_Member_by_source_taxon ( ) [inherited]

  Arg [1]    : string $source_name
  Arg [2]    : int $taxon_id
  Example    :

 $family->get_Member_by_source_taxon('ENSEMBLPEP',9606)

  Description: 
  Returntype : array reference of Bio::EnsEMBL::Compara::Member
  Exceptions : 
  Caller     : public

Code:

click to view

public Bio::SimpleAlign Bio::EnsEMBL::Compara::AlignedMemberSet::get_SimpleAlign ( )

    Arg [-UNIQ_SEQ] : (opt) boolean (default: false)
        : whether redundant sequences should be discarded
    Arg [-ID_TYPE] (opt) string (one of 'STABLE'*, 'SEQ', 'MEMBER', 'DISPLAY')
        : which identifier should be used as sequence names: the stable_id, the sequence_id, the seq_member_id, or the display_label
    Arg [-STOP2X] (opt) boolean (default: false)
        : whether the stop codons (character '*') should be replaced with gaps (character 'X')
          other unusual aminoacids (U and O) are also replaced by their closest match (C and K)
    Arg [-APPEND_TAXON_ID] (opt) boolean (default: false)
        : whether the taxon_ids should be added to the sequence names
    Arg [-APPEND_SP_SHORT_NAME] (opt) boolean (default: false)
        : whether the species (in short name format) should be added to the sequence names
    Arg [-APPEND_GENOMEDB_ID] (opt) boolean (default: false)
        : whether the genome_db_id should be added to the sequence names
    Arg [-APPEND_SPECIES_TREE_NODE_ID] (opt) hashref (default:  undef)
        : a hashref {genome_db_id => SpeciesTreeNode} that is used to append species_tree_node_ids
          to the sequence names
    Arg [-REMOVE_GAPS] (opt) boolean (default: false)
        : whether columns that only contain gaps should be removed from the alignment
    Arg [-SEQ_TYPE] (opt) string
        : which sequence should be used instead of the default one.
        : Can be 'exon_cased' for proteins and ncRNAs, and 'cds' for proteins only
    Arg [-REMOVED_COLUMNS] (opt) arrayref of arrayrefs of 2 integers
        : intervals of columns that have to be filtered out from the alignment
        : Used to remove low-information sites

  Example    :

 $tree->get_SimpleAlign(-SEQ_TYPE => 'cds');

  Description: Returns the alignment as a BioPerl object
  Returntype : Bio::SimpleAlign
  Exceptions : none
  Caller     : general

Code:

click to view

sub get_SimpleAlign {

    my ($self, @args) = @_;

    my $id_type = 0;
    my $unique_seqs = 0;
    my $stop2x = 0;
    my $append_taxon_id = 0;
    my $append_sp_short_name = 0;
    my $append_genomedb_id = 0;
    my $append_stn_id = undef;
    my $remove_gaps = 0;
    my $seq_type = undef;
    my $removed_columns = undef;
    my $removed_members = undef;
    my $map_long_seq_names = undef;
    my $annotations = undef;
    my $member_sources = undef;
    if (scalar @args) {
        ($unique_seqs,  $id_type, $stop2x, $append_taxon_id, $append_sp_short_name, $append_genomedb_id, $append_stn_id, $remove_gaps, $seq_type, $removed_columns, $removed_members, $map_long_seq_names, $annotations, $member_sources) =
            rearrange([qw(UNIQ_SEQ ID_TYPE STOP2X APPEND_TAXON_ID APPEND_SP_SHORT_NAME APPEND_GENOMEDB_ID APPEND_SPECIES_TREE_NODE_ID REMOVE_GAPS SEQ_TYPE REMOVED_COLUMNS REMOVED_MEMBERS MAP_LONG_SEQ_NAMES ANNOTATIONS MEMBER_SOURCES)], @args);
    }

    die "-SEQ_TYPE and \$self->seq_type cannot have different (defined) values" if $seq_type and $self->seq_type and ($seq_type ne $self->seq_type);
    $seq_type = $self->seq_type unless $seq_type;

    my $sa = Bio::SimpleAlign->new();
    if ($self->stable_id) {
        $sa->id($self->stable_id);
    } elsif ($self->dbID) {
        my $type = ref($self);
        $type =~ s/^.*://;
        $sa->id($type."/".($self->dbID));
    } else {
        $sa->id("$self");
    }

    Bio::EnsEMBL::Compara::Utils::Preloader::load_all_sequences($self->adaptor->db->get_SequenceAdaptor, $seq_type, $self) if $self->adaptor;

    my $seq_hash = {};

    #Counter for the unique temporary sequence name.
    my $countSeq = 0;

    #This is to assure that we have always the same order when dumping and mapping an alignment. 
    my @all_members; 
    if ($map_long_seq_names){
        @all_members = sort {$b->dbID <=> $a->dbID} @{$self->get_all_Members};
    }else{
        @all_members = @{$self->get_all_Members};
    }

    # Filter by source_name
    if ($member_sources) {
        my %sources_ok = map {$_ => 1} @$member_sources;
        @all_members = grep {$sources_ok{$_->source_name}} @all_members;
    }

    foreach my $member (@all_members) {

        next if $member->source_name =~ m/^Uniprot/i and $seq_type;

        next if $removed_members->{$member->dbID};

        # When seq_types are consistent (there is a die above that ensures
        # that both can't be set to different values) prevent the x3 expansion
        if ($seq_type and $self->seq_type) {
            $member->{_expansion_factor} = 1;
        }

        my $seqstr = $member->alignment_string($seq_type);
        next unless $seqstr;

        # Print unique sequences only ?
        next if $unique_seqs and $seq_hash->{$seqstr};
        $seq_hash->{$seqstr} = 1;

        my $alphabet = $member->source_name =~ /TRANS$/ ? 'dna' : 'protein';
        $alphabet = 'dna' if $seq_type and ($seq_type eq 'cds');

        # Sequence name
        my $prefix;
        $prefix = $member->stable_id;
        $prefix = $member->sequence_id if $id_type and $id_type =~ m/^SEQ/i;
        $prefix = $member->seq_member_id if $id_type and $id_type =~ m/^MEM/i;
        $prefix = $member->display_label if $id_type and $id_type =~ m/^DIS/i;

        my $suffix = '';
        $suffix = $member->taxon_id if($append_taxon_id);
        $suffix = $member->genome_db_id if ($append_genomedb_id);
        $suffix = $append_stn_id->{$member->genome_db_id}->node_id if ($append_stn_id);

        my $seqID;

        # Many phylogenetic methods take as input the sequence format phylip, which has a limit of 10 characters for the sequence names.
        # For those cases we need to create a mapping of the original names to be re-mapped whenever needed.
        # If map_long_seq_names is defined we need to store the mapping 
        if ($map_long_seq_names){
            #Instead, the sequence counter needs to be incremented for each sequence.
            $countSeq++;

            #Define the new sequence naming
            #Flanking with chars to assure that the \b regex will work just fine.
            $seqID = "x".$countSeq."x";

            #Store the mapping
            #Permanent sequence name
            $map_long_seq_names->{$seqID}->{'seq'} = $prefix;

            #Permanent suffix that will be used in the DB, it needs to be here to be replaced by the derived classes
            $map_long_seq_names->{$seqID}->{'suf'} = $suffix;

        }else{
            $seqID = $prefix;
            $seqID .= '_' . $suffix if $suffix;
        }

        ## Append $seqID with species short name, if required
        if ($append_sp_short_name and $member->genome_db_id) {
            my $species = $member->genome_db->get_short_name;
            $species =~ s/\s/_/g;
            $seqID .= "_" . $species;
        }

        if ($stop2x) {
            $seqstr =~ s/\*/X/g;
            if ($alphabet eq 'protein') {
                $seqstr =~ s/U/C/g;
                $seqstr =~ s/O/K/g;
            }
        }

        my $n_gaps = ($seqstr =~ tr/-/-/);

        my $seq = Bio::LocatableSeq->new(
                -SEQ        => $seqstr,
                -ALPHABET   => $alphabet,
                -START      => 1,
                -END        => length($seqstr)-$n_gaps,
                -ID         => $seqID,
                -STRAND     => 0
        );

        $sa->add_seq($seq);
    }
    $sa = $sa->remove_gaps(undef, 1) if $remove_gaps;
    $sa = $sa->remove_columns(@$removed_columns) if $removed_columns and scalar(@$removed_columns);

    #Add cutoffs annotations to the HMM stockholm file, which will be caught by hmmbuild
    if ($annotations) {
        # AnnotationCollection from the SimpleAlign object
        my $factory = Bio::Annotation::AnnotationFactory->new( -type => 'Bio::Annotation::SimpleValue' );

        foreach my $annot_tag (keys %{$annotations}) {
            my $annot_obj = $factory->create_object( -value => $annotations->{$annot_tag}, -tagname => $annot_tag );
            $sa->annotation->add_Annotation( 'custom', $annot_obj );
        }
    }

    return $sa;
}

public Boolean Bio::EnsEMBL::Compara::MemberSet::has_species_by_name ( ) [inherited]

  Arg [1]    : string $species_name
  Example    :

 my $ret = $homology->has_species_by_name("Homo sapiens");

  Description: return TRUE or FALSE whether one of the members in the homology is from the given species
  Returntype : 1 or 0
  Exceptions :
  Caller     :

Code:

click to view

public void Bio::EnsEMBL::Compara::AlignedMemberSet::load_cigars_from_file ( )

  Arg [1]    : string $file 
               The name of the file containing the multiple alignment
  Arg [-IMPORT_SEQ] : (opt) boolean (default: false)
               Whether the sequences of the members should be reassigned
               using the alignment file
  Arg [-CHECK_SEQ] : (opt) boolean (default: false)
               Checks that all the sequences recovered from the alignment file
               match the data in the Members
  Arg [-FORMAT]     : (opt) string (default: undef)
               The format of the alignment. By default, BioPerl will try to
               guess it from the file extension.  Refer to
               https://bioperl.org/howtos/AlignIO_and_SimpleAlign_HOWTO.html
               for a list of the supported formats.
  Arg [-ID_TYPE] (opt) string (one of 'STABLE', 'SEQ', 'MEMBER'*)
                : which identifier should be used as sequence names:
                  the stable_id, the sequence_id, or the seq_member_id (default)
  Example    :

 $family->load_cigars_from_file('/tmp/clustalw.aln');

  Description: Parses the multiple alignment fileand sets the cigar lines
               of each of the members of this AlignedMemberSet
  Returntype : none
  Exceptions : thrown if file cannot be parsed
               dies if a sequence identifier cannot be found in the set
               dies if there is a sequence mismatch with the set

Code:

click to view

sub load_cigars_from_file {
    my ($self, $file, @args) = @_;

    my $format;
    my $import_seq;
    my $id_type;
    my $check_seq;
    if (scalar @args) {
        ($import_seq, $format, $id_type, $check_seq) =
            rearrange([qw(IMPORT_SEQ FORMAT ID_TYPE CHECK_SEQ)], @args);
    }

    ## First read the alignment file and put the data in a hash
    my $alignio = Bio::AlignIO->new(-file => $file, -format => $format);
    my $aln = $alignio->next_aln or die "Bio::AlignIO could not get next_aln() from file '$file'";

    my $aln_length;
    my %aln_string_hash;
    foreach my $seq ($aln->each_seq) {

        my $id = $seq->display_id;
        $id =~ s/_.*$//;

        my $seqseq = $seq->seq();
        throw("No sequence for alignment portion: [$id]") unless $seqseq;

        $aln_length = length($seqseq) unless $aln_length;
        throw("While parsing the alignment, [$id] did not return the expected alignment length $aln_length\n") if $aln_length != length($seqseq);

        my $cigar_line = Bio::EnsEMBL::Compara::Utils::Cigars::cigar_from_alignment_string($seqseq);
        throw("empty cigar_line for [$id]\n") unless length($cigar_line);

        $seqseq =~ s/-//g;

        # The following code is disabled because it slows down the method by 10%
        # We should rather have unit-testing on Bio::EnsEMBL::Compara::Utils::Cigars

        ## Check that the cigar length (Ms) matches the sequence length
        ## Take the M lengths into an array
        #my @cigar_match_lengths = map { $_ eq '' ? 1 : $_ } map { $_ =~ /^(\d*)/ } ( $cigar_line =~ /(\d*[M])/g );
        ## Sum up the M lengths
        #my $seq_cigar_length; map { $seq_cigar_length += $_ } @cigar_match_lengths;
        #if ($seq_cigar_length != length($seqseq)) {
            #warn "$id:$seq_cigar_length:".length($seqseq)."\n$seqseq\n$cigar_line\n".$seq->seq."\n";
            #$self->throw("While storing the cigar line, the returned cigar length did not match the sequence length\n");
        #}

        $aln_string_hash{$id} = [$seqseq, $cigar_line];
    }
    $self->aln_length($aln_length);

    my %unseen = (map {$_ => 1} (keys %aln_string_hash));

    ## Then we associate each member with its alignment, based on the required $id_type
    if ($self->isa('Bio::EnsEMBL::Compara::GeneTree') && $self->root->is_supertree()) {
        $self->expand_subtrees; # TODO: remove this additional call once ENSCOMPARASW-4276 is resolved
    }
    foreach my $member (@{$self->get_all_Members}) {

        my $seqID = $member->seq_member_id;
        $seqID = $member->sequence_id if $id_type and $id_type =~ m/^SEQ/i;
        $seqID = $member->stable_id if $id_type and $id_type =~ m/^STA/i;

        throw($member->stable_id." ($seqID) cannot be found in the alignment\n") unless exists $aln_string_hash{$seqID};

        my $seqseq = $aln_string_hash{$seqID}->[0];
        my $cigar_line = $aln_string_hash{$seqID}->[1];

        $member->cigar_line($cigar_line);

        if ($import_seq) {
            $member->sequence($seqseq);
        } elsif ($check_seq) {
            my $member_sequence = $member->sequence;
            # Stop-codons are often misrepresented so we just ignore them
            # (this is in line with SeqMember::bioseq)
            $seqseq =~ tr/*/X/;
            $member_sequence =~ tr/*/X/;
            if ($member_sequence ne $seqseq) {
                throw($member->stable_id." ($seqID) has a different sequence in the alignment file '$file'");
            }
        }

        delete $unseen{$seqID} if exists $unseen{$seqID};
    }

    if (scalar(keys %unseen)) {
        throw("No member for alignment ids: ".join(", ", keys %unseen));
    }
}

public Int Bio::EnsEMBL::Compara::MemberSet::Member_count_by_GenomeDB ( ) [inherited]

  Arg [1]    : Bio::EnsEMBL::Compara::GenomeDB $genome_db
  Example    :

 Member_count_by_GenomeDB($genome_db);

  Description: Convenience wrapper for member counts by a GenomeDB
  Returntype : int
  Exceptions : Thrown by subroutines this call. See get_Member
               equivalent
  Caller     : public

Code:

click to view

public Int Bio::EnsEMBL::Compara::MemberSet::Member_count_by_source ( ) [inherited]

  Arg [1]    : string $source_name
               e.g. "ENSEMBLPEP"
  Example    :

 $family->Member_count_by_source('ENSEMBLPEP');

  Description: 
  Returntype : int
  Exceptions : 
  Caller     : public

Code:

click to view

public Int Bio::EnsEMBL::Compara::MemberSet::Member_count_by_source_GenomeDB ( ) [inherited]

  Arg [1]    : string $source_name
  Arg [2]    : Bio::EnsEMBL::Compara::GenomeDB $genome_db
  Example    :

 Member_count_by_source_GenomeDB('ENSEMBLPEP', $genome_db);

  Description: Convenience wrapper for member counts by a GenomeDB
  Returntype : int
  Exceptions : Thrown by subroutines this call. See get_Member
               equivalent
  Caller     : public

Code:

click to view

public Int Bio::EnsEMBL::Compara::MemberSet::Member_count_by_source_taxon ( ) [inherited]

  Arg [1]    : string $source_name
  Arg [2]    : int $taxon_id
  Example    :

 Member_count_by_source_taxon('ENSEMBLPEP',9606);

  Description: 
  Returntype : int
  Exceptions : 
  Caller     : public

Code:

click to view

public Bio::EnsEMBL::Compara::MethodLinkSpeciesSet Bio::EnsEMBL::Compara::MemberSet::method_link_species_set ( ) [inherited]

  Arg [1]    : MethodLinkSpeciesSet object (optional)
  Example    :

  Description: getter/setter method for the MethodLinkSpeciesSet for this relation.
               Can lazy-load the method_link_species_set from the method_link_species_set_id
               if that one is set and the adaptor is set.
  Returntype : Bio::EnsEMBL::Compara::MethodLinkSpeciesSet
  Exceptions : throws if setting to an unsuitable object 
  Caller     : general
  Status     : Stable

Code:

click to view

public Int Bio::EnsEMBL::Compara::MemberSet::method_link_species_set_id ( ) [inherited]

  Arg [1]    : integer (optional)
  Example    :

  Description: getter/setter method for the internal ID of the MethodLinkSpeciesSet
               for this relation.
  Returntype : integer
  Exceptions : none
  Caller     : general
  Status     : Stable

Code:

click to view

public Bio::EnsEMBL::Compara::AlignedMemberSet Bio::EnsEMBL::Compara::AlignedMemberSet::new ( )

  Example    :

  Description:
  Returntype : Bio::EnsEMBL::Compara::AlignedMemberSet
  Exceptions :
  Caller     :

Code:

click to view

Reimplemented from Bio::EnsEMBL::Compara::MemberSet.

Reimplemented in Bio::EnsEMBL::Compara::Family, and Bio::EnsEMBL::Compara::GeneTree.

public Bio::SimpleAlign Bio::EnsEMBL::Compara::AlignedMemberSet::print_alignment_to_file ( )

  Arg [1]     : scalar (string or file handle) - output file
  Arg [-FORMAT] : string - format of the output (cf BioPerl capabilities) (example: 'fasta')
  Arg [-FLAT_DISPLAYNAME] : boolean - whether the sequence names in the alignment should be displayed "flat".
                            "flat" is a BioPerl option that means: without alignment start/end (e.g "/1-2542").
  Arg [...]   : all the other arguments of get_SimpleAlign
  Example     :

 $family->print_alignment_to_file(-APPEND_TAXON_ID => 1);

  Description : Wrapper around get_SimpleAlign to print the alignment to a file (stdout by default)
  Returntype  : Bio::SimpleAlign
  Exceptions  : none
  Caller      : general
  Status      : Stable

Code:

click to view

public Number Bio::EnsEMBL::Compara::MemberSet::print_sequences_to_file ( ) [inherited]

  Arg [1]     : scalar (string or file handle) - output file
  Arg [-FORMAT]   : string - format of the output (cf BioPerl capabilities) (example: 'fasta')
  Arg [-UNIQ_SEQ] : boolean - whether only 1 copy of each sequence should be printed
                    (when multiple proteins share the same sequence)
  Arg [...]  : all the other arguments of SeqMember::bioseq()
  Example    :

 $family->print_sequences_to_file(-file => 'output.fa', -format => 'fasta', -id_type => 'MEMBER');

  Description: Prints the sequences of the members into a file
  Returntype : number of unique sequences in the set
  Exceptions : none
  Caller     : general
  Status     : Stable

Code:

click to view

public void Bio::EnsEMBL::Compara::MemberSet::remove_Member ( ) [inherited]

  Arg [1]    : Bio::EnsEMBL::Compara::Member $member
  Example    :

 $gene_tree->remove_Member($member);

  Description: Remove a new Member from this set
  Returntype : none
  Exceptions : Throws if input objects don't check
  Caller     : general
  Status     : Stable

Code:

click to view

public String Bio::EnsEMBL::Compara::AlignedMemberSet::seq_type ( )

  Description : Getter/Setter for the seq_type field. This field contains
                the type of sequence used for the members. If undefined,
                the usual sequence is used. Otherwise, there must be a
                matching sequence in the other_member_sequence table.
  Returntype  : String
  Example     :

 my $type = $tree->seq_type();

  Caller      : General

Code:

click to view

public String Bio::EnsEMBL::Compara::MemberSet::stable_id ( ) [inherited]

  Arg [1]    : string $stable_id (optional)
  Example    :

  Description: Getter/setter for the stable ID of this relation
  Returntype : string
  Exceptions : none
  Caller     : general
  Status     : Stable

Code:

click to view

public String Bio::EnsEMBL::Compara::MemberSet::toString ( ) [inherited]

  Example     :

 $member_set->toString();

  Description : Returns a description of this object as a string
  Returntype  : String
  Exceptions  : none
  Caller      : general

Code:

click to view

Reimplemented in Bio::EnsEMBL::Compara::Homology.

public String Bio::EnsEMBL::Compara::AlignedMemberSet::update_alignment_stats ( )

  Arg [1]    : (optional) $seq_type
  Example    :

 my $consensus_cigar = $gene_tree->consensus_cigar_line();

  Description: Update the AlignedMemberSet properties (% identity, % coverage, etc)
               for the given seq_type (the default sequence is used otherwise)
  Returntype : string
  Caller     : internal

Code:

click to view

sub update_alignment_stats {
    my $self = shift;
    my $seq_type = shift;

    Bio::EnsEMBL::Compara::Utils::Preloader::load_all_sequences($self->adaptor->db->get_SequenceAdaptor, $seq_type, $self) if $self->adaptor;

    my $genes = $self->get_all_Members;
    my $ngenes = scalar(@$genes);

    # New cigars
    my @new_cigars = Bio::EnsEMBL::Compara::Utils::Cigars::minimize_cigars(map {$_->cigar_line} @$genes);
    for (my $j=0; $j<$ngenes; $j++) {
        $genes->[$j]->cigar_line($new_cigars[$j]);
    }

    my @aln = map {$_->alignment_string($seq_type)} @$genes;
    my $aln_length = length($aln[0]);

    my @seq_length   = (0) x $ngenes;
    my @nmatch_id    = (0) x $ngenes;
    my @nmatch_pos   = (0) x $ngenes;
    my @nmismatch    = (0) x $ngenes;

    if ($ngenes == 2) {
        # This code is >4 times faster with pairs of genes

        my $gene1 = $genes->[0];
        my $gene2 = $genes->[1];

        my @aln1 = split(//, $aln[0]);
        my @aln2 = split(//, $aln[1]);

        for (my $i=0; $i <= $#aln1; $i++) {
            my $gap1 = ($aln1[$i] eq '-');
            my $gap2 = ($aln2[$i] eq '-');
            $seq_length[0]++ unless $gap1;
            $seq_length[1]++ unless $gap2;
            if (not $gap1 and not $gap2) {
                if ($aln1[$i] eq $aln2[$i]) {
                    $nmatch_id[0]++;
                    $nmatch_pos[0]++;
                } elsif ($matrix_hash{uc $aln1[$i]}{uc $aln2[$i]} > 0) {
                    $nmatch_pos[0]++;
                    $nmismatch[0]++;
                } else {
                    $nmismatch[0]++;
                }
            }
        }

        $nmatch_id[1] = $nmatch_id[0];
        $nmatch_pos[1] = $nmatch_pos[0];
        $nmismatch[1] = $nmismatch[0];

    } else {

        my $min_seq = shift;
        $min_seq = int($min_seq * $ngenes) if $min_seq <= 1;

        for (my $i=0; $i < $aln_length; $i++) {

            my @char_i =  map {substr($_, $i, 1)} @aln;
            #print "pos $i: ", join('/', @char_i), "\n";

            my %seen;
            map {$seen{$_}++} @char_i;
            next if $seen{'-'} == $ngenes;
            my $is_cov_match = ($seen{'-'} <= $ngenes-2);
            delete $seen{'-'};

            my %pos_chars = ();
            my @chars = keys %seen;
            while (my $c1 = shift @chars) {
                foreach my $c2 (@chars) {
                    if (($matrix_hash{uc $c1}{uc $c2} > 0) and ($seen{$c1}+$seen{$c2} >= $min_seq)) {
                        $pos_chars{$c1} = 1;
                        $pos_chars{$c2} = 1;
                    }
                }
            }

            for (my $j=0; $j<$ngenes; $j++) {
                if ($char_i[$j] ne '-') {
                    $seq_length[$j]++;
                    if ($seen{$char_i[$j]} >= $min_seq) {
                        $nmatch_id[$j]++;
                        $nmatch_pos[$j]++;
                    } elsif (exists $pos_chars{$char_i[$j]}) {
                        $nmatch_pos[$j]++;
                        $nmismatch[$j]++;
                    } else {
                        $nmismatch[$j]++;
                    }
                }
            }
        }
    }

    for (my $j=0; $j<$ngenes; $j++) {
        if ($seq_length[$j]) {
            $genes->[$j]->num_matches( $nmatch_id[$j] );
            $genes->[$j]->num_pos_matches( $nmatch_pos[$j] );
            $genes->[$j]->num_mismatches( $nmismatch[$j] );
            $genes->[$j]->perc_id( 100.0 * $nmatch_id[$j] / $seq_length[$j] );
            $genes->[$j]->perc_pos( 100.0 * $nmatch_pos[$j] / $seq_length[$j] );
            $genes->[$j]->perc_cov( 100.0 * ($nmatch_id[$j]+$nmismatch[$j]) / $seq_length[$j] );
        }
    }
}

public Int Bio::EnsEMBL::Compara::MemberSet::version ( ) [inherited]

  Arg [1]    : string $version (optional)
  Example    :

  Description: Getter/setter for the version number of the stable ID
  Returntype : int
  Exceptions : none
  Caller     : general
  Status     : Stable

Code:

click to view

Member Data Documentation

public Bio::EnsEMBL::Compara::AlignedMemberSet::String

Reimplemented from Bio::EnsEMBL::Compara::MemberSet.

Reimplemented in Bio::EnsEMBL::Compara::GeneTree.

Definition at line 191 of file AlignedMemberSet.pm.

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

Bio/EnsEMBL/Compara/AlignedMemberSet.pm

Upcoming Ensembl Platform Transition

Class Summary

Description

Synopsis

Available Methods

Public Attributes

Method Documentation

Member Data Documentation

About Us

Get help

Our sister sites

Follow us