Public Member Functions
public A	new ()

public String	create_tempdir ()

public void	remove_tempdir ()

public void	write_sequence ()

public void	run_blastz ()

public void	lav_to_axt ()

public void	find_best_alignment ()

public void	parse_blastz_output ()

public void	cleanup_tmpfiles ()

public void	found_match ()

public void	adjust_coords ()

public	filter_overlaps ()

public void	write_assembly ()

public Int	stats_incr ()

public void	init_stats ()

public Int	get_stats ()

public void	log_block_stats ()

public void	log_overall_stats ()

protected Values	_by_chr_num ()

public String	AUTOLOAD ()

Detailed Description

Synopsis

my $support = new Bio::EnsEMBL::Utils::ConversionSupport($SERVERROOT);
my $aligner = BlastzAligner->new(-SUPPORT => $support);
 
# create a tempdir for storing input and output
$aligner->create_tempdir;
 
# write sequences to fasta and nib files
$aligner->write_sequence(
    $alt_slice,
    $support->param('altassembly'), 
    "alt_sequence.1"
);
$aligner->write_sequence(
    $ref_slice,
    $support->param('assembly'), 
    "ref_sequence.1"
);
 
# run blastz
$aligner->run_blastz("alt_sequence.1", "ref_sequence.1");

Description

This modules contains helper functions to generate a whole genome alignment
between two closely related assemblies using blastz from scratch. Alignments
are then stored in an Ensembl assembly table.

The module is part of a series of scripts to create a mapping between two
assemblies. See "Related scripts" below for an overview of the whole process.

Definition at line 43 of file BlastzAligner.pm.

Member Function Documentation

◆ _by_chr_num()

protected Values AssemblyMapper::BlastzAligner::_by_chr_num ( )

  Example     :

my @sorted = sort _by_chr_num qw(X, 6-COX, 14, 7);

  Description : Subroutine to use in sort for sorting chromosomes. Sorts
                numerically, then alphabetically
  Return type : values to be used by sort
  Exceptions  : none
  Caller      : internal

Code:

click to view

◆ adjust_coords()

public void AssemblyMapper::BlastzAligner::adjust_coords ( )

  Arg[1]      : Int $A_start - start of alternatvie block in chromosomal coords
  Arg[2]      : Int $A_end - end of alternative block in chromosomal coords
  Arg[3]      : Arrayref $R_coords - list of start/end pairs of reference
                blocks in chromosomal coords
  Example     :

my $R_coords = [ [1, 1000], [3000, 5000] ];

$aligner->adjust_coords(1, 30000000, $R_coords);

  Description : Adjusts coordinates of blastz alignments to chromosomal coords.
  Return type : none
  Exceptions  : none
  Caller      : general

Code:

click to view

sub adjust_coords {
    my ($self, $A_start, $A_end, $R_coords) = @_;
    my $R_chr = $self->seq_region_name;
    my $id = $self->id;
 
    for (my $align = 0; $align < scalar(@{ $self->{'_match'}->{$R_chr}->{$id} || [] }); $align++) {
        for (my $c = 0; $c < scalar(@{ $self->{'_match'}->{$R_chr}->{$id}->[$align] || []}); $c++) {
            $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[0] += $A_start - 1;
            $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[1] += $A_start - 1;
 
            # forward strand match
            if ($self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[4] == 1) {
                $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[2] += $R_coords->{$self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[5]}->[0] - 1;
                $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[3] += $R_coords->{$self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[5]}->[0] - 1;
            
            # reverse strand match
            } else {
                my $tmp_start = $R_coords->{$self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[5]}->[1] - $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[3] + 1;
 
                $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[3] = $R_coords->{$self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[5]}->[1] - $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[2] + 1;
 
                $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[2] = $tmp_start;
            }
 
            # sanity check: aligned region pairs must have same length
            my $A_len = $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[1] - $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[0];
            my $R_len = $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[3] - $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[2];
            $self->support->log_warning("Length mismatch: $A_len <> $R_len in block $id, alignment $align, stretch $c\n", 2) unless ($A_len == $R_len);
        }
    }
}

◆ AUTOLOAD()

public String AssemblyMapper::BlastzAligner::AUTOLOAD ( )

  Arg[1]      : (optional) String/Object - attribute to set
  Example     :

 # setting a attribute
                $self->attr($val);
@section autotoc_md0 getting the attribute
                $self->attr;
@section autotoc_md1 undefining an attribute
                $self->attr(undef);

  Description : lazy function generator for getters/setters
  Return type : String/Object
  Exceptions  : none
  Caller      : general

Code:

click to view

◆ cleanup_tmpfiles()

public void AssemblyMapper::BlastzAligner::cleanup_tmpfiles ( )

  Arg[1-N]    : (optional) list @files - additional tmp files to delete
  Example     :

$self->cleanup_tmpfiles('e_seq.fa', 'v_seq.fa');

  Description : deletes temporary files
  Return type : none
  Exceptions  : Warning if file cannot be deleted
  Caller      : general
  Status      : stable

Code:

click to view

◆ create_tempdir()

public String AssemblyMapper::BlastzAligner::create_tempdir ( )

  Arg[1]      : String $tmpdir - temporary directory name
  Example     :

$aligner->create_tempdir;

  Description : Creates a temporary directory in /tmp with a semi-random name
                (username.timestamp.randomnumber). Alternatively, you can pass
                it the name of a directory to use.
  Return type : String - name of the tempdir created
  Exceptions  : Thrown if tempdir can't be created
  Caller      : general

Code:

click to view

◆ filter_overlaps()

public AssemblyMapper::BlastzAligner::filter_overlaps ( )

  Description : Deprecated:
. This filtering algorithm didn't work well. Please
                run the separate script fix_overlaps.pl.

Code:

click to view

sub filter_overlaps {
    my $self = shift;
    my $id = $self->id;
    
    foreach my $R_chr (sort keys %{ $self->{'_match'} }) {
        my $coord_check = [];
        # rearrange the datastructure so that we can find overlaps
        foreach my $id (keys %{ $self->{'_match'}->{$R_chr} }) {
            for (my $align = 0; $align < scalar(@{ $self->{'_match'}->{$R_chr}->{$id} }); $align++) {
                for (my $c = 0; $c < scalar(@{ $self->{'_match'}->{$R_chr}->{$id}->[$align] || []}); $c++) {
                    push @{ $coord_check }, [
                        $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[0],
                        $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[1],
                        $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[2],
                        $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[3],
                        $id,
                        $align,
                        $c,
                    ];
                }
            }
        }
 
        my @A_sort = sort { $a->[0] <=> $b->[0] } @{ $coord_check };
        my @R_sort = sort { $a->[2] <=> $b->[2] } @{ $coord_check };
 
        # sanity check: alternative alignments must not overlap (axtBest should
        # guarantee that)
        my $last;
        foreach my $c (@A_sort) {
            $self->support->log_warning("Overlapping alternative alignment at ".join(':', $R_chr, $c->[0], $c->[1])." (last_end ".$last->[1].")\n", 1) if ($last and $c->[0] <= $last->[1]);
            $last = $c;
        }
 
        # now filter reference overlaps
        my @seen;
        $last = undef;
        foreach my $c (@R_sort) {
            if ($last and $c->[2] <= $last->[3]) {
                $self->support->log_verbose("Overlapping reference alignment at ".join(':', $R_chr, $c->[2], $c->[3])." (last_end ".$last->[3].")\n", 1);
 
                # if last alignment was longer, delete this one
                if ($last->[3]-$last->[2] > $c->[3]-$c->[2]) {
                    undef $self->{'_match'}->{$R_chr}->{$c->[4]}->[$c->[5]]->[$c->[6]];
 
                # if last alignment was shorter, trace back and delete all
                # overlapping shorter alignments
                } else {
                    foreach my $s (@seen) {
                        # earlier alignment still overlapping
                        if ($c->[2] <= $s->[3]) {
                            # earlier alignment shorter -> delete it
                            if ($s->[3]-$s->[2] < $c->[3]-$c->[2]) {
                                undef $self->{'_match'}->{$R_chr}->{$s->[4]}->[$s->[5]]->[$s->[6]];
 
                            # this alignment shorter -> delete it
                            } else {
                                undef $self->{'_match'}->{$R_chr}->{$c->[4]}->[$c->[5]]->[$c->[6]];
                                $last = $s;
                                last;
                            }
                        } else {
                            $last = $s;
                            last;
                        }
                    }                
                    
                    $last = $c;
                }
            }
            unshift @seen, $c;
            $last = $c unless ($last);
        }
    }
}

◆ find_best_alignment()

public void AssemblyMapper::BlastzAligner::find_best_alignment ( )

  Example     :

$aligner->find_best_alignment;

  Description : Finds the best set of non-overlapping alignments by running
                axtBest.
  Return type : none
  Exceptions  : thrown if axtBest fails
  Caller      : general

Code:

click to view

◆ found_match()

public void AssemblyMapper::BlastzAligner::found_match ( )

  Arg[1]      : Boolean $match_flag - flag indicating if last bp was a match
  Arg[2]      : Int $j - current bp position in the alignment
  Arg[3]      : Hashref $coords - alignment coordinates and strand from blastz
                output
  Description : Populates a datastructure describing blocks of alignment.
  Return type : none
  Exceptions  : none
  Caller      : internal

Code:

click to view

◆ get_stats()

public Int AssemblyMapper::BlastzAligner::get_stats ( )

  Arg[1]      : String $code - stats code
  Example     :

my $num_mismatches = $aligner->get_stats('mismatch');

  Description : Stats getter.
  Return type : Int - stats value
  Exceptions  : none
  Caller      : general

Code:

click to view

◆ init_stats()

public void AssemblyMapper::BlastzAligner::init_stats ( )

  Arg[1]      : Array @codes - stats codes to initialise
  Example     :

$aligner->init_stats('match', 'mismatch', 'alignments');

  Description : Initialises stats (i.e. sets to 0).
  Return type : none
  Exceptions  : none
  Caller      : general

Code:

click to view

◆ lav_to_axt()

public void AssemblyMapper::BlastzAligner::lav_to_axt ( )

  Example     :

$aligner->lav_to_axt;

  Description : Converts blastz output from lav to axt format. Target and query
                sequences must be present in nib format in the temporary
                directory for this to work.
  Return type : none
  Exceptions  : thrown if lavToAxt fails
  Caller      : general

Code:

click to view

◆ log_block_stats()

public void AssemblyMapper::BlastzAligner::log_block_stats ( )

  Arg[1]      : Int $indent - indentation level
  Example     :

$aligner->log_block_stats(3);

  Description : Logs stats for an alignment block.
  Return type : none
  Exceptions  : none
  Caller      : general

Code:

click to view

◆ log_overall_stats()

public void AssemblyMapper::BlastzAligner::log_overall_stats ( )

  Example     :

$aligner->log_overall_stats;

  Description : Logs overall alignment stats.
  Return type : none
  Exceptions  : none
  Caller      : general

Code:

click to view

sub log_overall_stats {
    my $self = shift;
    
    # blastz
    $self->support->log("\nOverall blastz alignment stats:\n");
 
    unless ($self->get_stats('alignments')) {
      $self->support->log("No alignments found.\n", 1);
      return;
    }
    
    $self->support->log(sprintf(FMT1,
        "Matches:",
        $self->get_stats('match_total'),
        $self->get_stats('match_total')/$self->get_stats('bp_total')*100),
    1);
    $self->support->log(sprintf(FMT1,
        "Mismatches:",
        $self->get_stats('mismatch_total'),
        $self->get_stats('mismatch_total')/$self->get_stats('bp_total')*100),
    1);
    $self->support->log(sprintf(FMT1,
        "Gaps:",
        $self->get_stats('gap_total'),
        $self->get_stats('gap_total')/$self->get_stats('bp_total')*100),
    1);
 
    # alignments to be written to assembly table
    $self->support->log_verbose("\nAlignments that will be written to assembly table:\n");
    $self->support->log_verbose(sprintf(FMT3,
        qw(CHR BLOCK ALIGNMENT ALT_START ALT_END REF_START REF_END)),
    1);
    $self->support->log_verbose(('-'x63)."\n", 1);
    foreach my $R_chr (sort _by_chr_num keys %{ $self->{'_match'} }) {
        foreach my $id (sort { $a <=> $b } keys %{ $self->{'_match'}->{$R_chr} }) {
            for (my $align = 0; $align < scalar(@{ $self->{'_match'}->{$R_chr}->{$id} }); $align++) {
                for (my $c = 0; $c < scalar(@{ $self->{'_match'}->{$R_chr}->{$id}->[$align] }); $c++) {
                    if ($self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]) {
                        $self->support->log_verbose(sprintf(FMT4,
                            $R_chr,
                            $id,
                            $align+1,
                            @{ $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c] }),
                        1);
                        my $l = $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[1] - $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[0];
                        $self->stats_incr('alignments_total', 1);
                        $self->stats_incr('short1_10_total', 1) if ($l < 11);
                        $self->stats_incr('short11_100_total', 1) if ($l > 10 and $l < 101);
                    }
                }
            }
        }
    }
    $self->support->log("\nAssembly entry stats:\n");
    $self->support->log(sprintf(FMT2,
        "Total alignment blocks:",
        $self->get_stats('alignments_total')),
    1);
    $self->support->log(sprintf(FMT2,
        "Alignments 1-10 bp:",
        $self->get_stats('short1_10_total')),
    1);
    $self->support->log(sprintf(FMT2,
        "Alignments 11-100 bp:",
        $self->get_stats('short11_100_total')),
    1);
}

◆ new()

public A AssemblyMapper::BlastzAligner::new ( )

  Arg[-SUPPORT] : a Bio::EnsEMBL::Utils::ConversionSupport object
  Example       :

my $support = new Bio::EnsEMBL::Utils::ConversionSupport(
                   $SERVERROOT);
                 my $aligner = BlastzAligner->new(-SUPPORT => $support);

  Description   : object constructor method
  Return type   : a BlastzAligner object
  Exceptions    : none
  Caller        : general

Code:

click to view

◆ parse_blastz_output()

public void AssemblyMapper::BlastzAligner::parse_blastz_output ( )

  Example     :

$aligner->parse_blastz_output;

  Description : Reads a blastz alignment result from an axt file and creates
                a datastructure containing ungapped alignments from it. Note
                that mismatches are allowed, but separate stats will be
                collected for them.
  Return type : none
  Exceptions  : none
  Caller      : general

Code:

click to view

sub parse_blastz_output {
    my $self = shift;
 
    # read file
    my $tmpdir = $self->tempdir;
    my $id = $self->id;
    my $fh = $self->support->filehandle('<', "$tmpdir/blastz.$id.best.axt");
 
    # initialize stats
    $self->init_stats(qw(match mismatch gap alignments bp));
    
    my $i = 1;
    my ($header, $A_seq, $R_seq);
 
    while (my $line = <$fh>) {
        # there are blocks of 4 lines, where line 1 is the header, line 2 is
        # A_seq, line3 is R_seq
        $header = $line unless (($i-1) % 4);
        $A_seq = $line unless (($i-2) % 4);
        chomp $A_seq;
        my @A_arr = split(//, $A_seq);
        $R_seq = $line unless (($i-3) % 4);
        chomp $R_seq;
        my @R_arr = split(//, $R_seq);
 
        # compare sequences letter by letter
        if ($i % 4 == 0) {
            my $match_flag = 0;
            $self->init_stats(qw(A_gap R_gap));
            my %coords;
            @coords{'R_id', 'A_start', 'A_end', 'R_start', 'R_end', 'strand'} =
                (split(/ /, $header))[4, 2, 3, 5, 6, 7];
            $coords{'R_id'} =~ s/ref_seq\.(.*)/$1/;
            ($coords{'strand'} eq '+') ? ($coords{'strand'} = 1) :
                                         ($coords{'strand'} = -1);
            for (my $j = 0; $j < scalar(@A_arr); $j++) {
                # gap
                if ($A_arr[$j] eq '-' or $R_arr[$j] eq '-') {
                    $self->stats_incr('gap', 1);
                    $self->stats_incr('A_gap', 1) if ($A_arr[$j] eq '-');
                    $self->stats_incr('R_gap', 1) if ($R_arr[$j] eq '-');
                    $match_flag = 0;
                } else {
                    $self->found_match($match_flag, $j, \%coords);
                    $match_flag = 1;
 
                    # match
                    if ($A_arr[$j] eq $R_arr[$j]) {
                        $self->stats_incr('match', 1);
                    # mismatch
                    } else {
                        $self->stats_incr('mismatch', 1);
                    }
                }
            }
            $self->stats_incr('bp', scalar(@A_arr));
            $self->stats_incr('alignments', 1);
        }
        
        $i++;
    }
}

◆ remove_tempdir()

public void AssemblyMapper::BlastzAligner::remove_tempdir ( )

  Example     :

$aligner->remove_tempdir;

  Description : Removes a temporary directory created by $self->create_tempdir.
  Return type : none
  Exceptions  : none
  Caller      : general

Code:

click to view

◆ run_blastz()

public void AssemblyMapper::BlastzAligner::run_blastz ( )

  Arg[1]      : String $A_basename - basename of alternative fasta file
  Arg[2]      : String $R_basename - basename of reference fasta file
  Example     :

$aligner->run_blastz('alt_seq.1', 'ref_seq.1');

  Description : Runs blastz between an alternative and multiple reference
                sequences.
  Return type : none
  Exceptions  : thrown if blastz fails
  Caller      : general

Code:

click to view

◆ stats_incr()

public Int AssemblyMapper::BlastzAligner::stats_incr ( )

  Arg[1]      : String $code - stats code
  Arg[2]      : Int $incr - number by which stats should be incremented
  Example     :

$aligner->stats_incr('total_alignments', 1);

  Description : Increments stats.
  Return type : Int - stat value after increment
  Exceptions  : none
  Caller      : general

Code:

click to view

◆ write_assembly()

public void AssemblyMapper::BlastzAligner::write_assembly ( )

  Arg[1]      : Bio::*::DBSQL::DBAdaptor $R_dba - reference database adaptor
  Example     :

my $R_dba = $support->get_database('ensembl');

$aligner->write_assembly($R_dba);

  Description : Writes assembly entries for blastz alignments to the database.
  Return type : none
  Exceptions  : none
  Caller      : general

Code:

click to view

sub write_assembly {
    my ($self, $R_dba) = @_;
 
    my $R_dbh = $R_dba->dbc->db_handle;
    my $R_sa = $R_dba->get_SliceAdaptor;
 
    my $sth = $R_dbh->prepare(qq(
        INSERT IGNORE INTO assembly (asm_seq_region_id, cmp_seq_region_id,
            asm_start,asm_end, cmp_start, cmp_end, ori)
        VALUES (?, ?, ?, ?, ?, ?, ?)
    ));
 
    $self->support->log("Adding assembly entries for alignments...\n");
 
    my $i;
    foreach my $R_chr (sort _by_chr_num keys %{ $self->{'_match'} }) {
        # get seq_region_id for alternative and reference chromosome
        my $A_chr = $R_chr;
        my $R_slice = $R_sa->fetch_by_region('toplevel', $R_chr, undef, undef, undef, $self->support->param('assembly'));
        my $R_sid = $R_sa->get_seq_region_id($R_slice);
        # we need to fetch the alternative slice from the reference db
        # explicitely by coord_system, since toplevel attribute is not set
        # there
        my $cs_name = $R_slice->coord_system_name;
        my $A_sid = $R_sa->get_seq_region_id($R_sa->fetch_by_region($cs_name, $A_chr, undef, undef, undef, $self->support->param('altassembly')));
 
        foreach my $id (sort { $a <=> $b } keys %{ $self->{'_match'}->{$R_chr} }) {
            for (my $align = 0; $align < scalar(@{ $self->{'_match'}->{$R_chr}->{$id} }); $align++) {
                for (my $c = 0; $c < scalar(@{ $self->{'_match'}->{$R_chr}->{$id}->[$align] }); $c++) {
                    if ($self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]) {
                        $sth->execute(
                            $R_sid,
                            $A_sid,
                            $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[2],
                            $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[3],
                            $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[0],
                            $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[1],
                            $self->{'_match'}->{$R_chr}->{$id}->[$align]->[$c]->[4],
                        );
                        $i++;
                    }
                }
            }
        }
    }
    $self->support->log("Done inserting $i entries.\n");
}

◆ write_sequence()

public void AssemblyMapper::BlastzAligner::write_sequence ( )

  Arg[1]      : Bio::EnsEMBL::Slice $slice - slice for which to write sequence
  Arg[2]      : String $assembly - assembly name
  Arg[3]      : String $basename1 - basename of single sequence file
  Arg[4]      : (optional) String $basename2 - basename of multiple sequence
                file
  Example     :

$aligner->write_sequence($slice, 'NCBI35', 'e_seq.1');

  Description : Writes a slice's sequence to a fasta file and converts it to nib
                format. Optionally appends the sequence to another
                multi-sequence fasta file.
  Return type : none
  Exceptions  : thrown if faToNib or file appending fails
  Caller      : general

Code:

click to view

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

ensembl/misc-scripts/assembly/AssemblyMapper/BlastzAligner.pm

Public Member Functions

Detailed Description

Synopsis

Description

Member Function Documentation

◆ _by_chr_num()

◆ adjust_coords()

◆ AUTOLOAD()

◆ cleanup_tmpfiles()

◆ create_tempdir()

◆ filter_overlaps()

◆ find_best_alignment()

◆ found_match()

◆ get_stats()

◆ init_stats()

◆ lav_to_axt()

◆ log_block_stats()

◆ log_overall_stats()

◆ new()

◆ parse_blastz_output()

◆ remove_tempdir()

◆ run_blastz()

◆ stats_incr()

◆ write_assembly()

◆ write_sequence()