view release on metacpan or  search on metacpan

lib/App/Dazz/Command/layout.pm  view on Meta::CPAN

package App::Dazz::Command::layout;
use strict;
use warnings;
use autodie;

use App::Dazz -command;
use App::Dazz::Common;

sub abstract {
    return 'layout anchors within a group';
}

sub opt_spec {
    return (
        [ "outfile|o=s", "output filename", ],
        [ 'border=i', 'length of borders in anchors', { default => 500 }, ],
        [ "max=i",    "max distance",                 { default => 5000 }, ],
        [ 'pa=s',     'prefix of anchors',            { default => "anchor" }, ],
        [ 'oa=s',     'overlaps between anchors', ],
        [ "png",      "write a png file via graphviz", ],
        { show_defaults => 1, },
    );
}

sub usage_desc {
    return "dazz layout [options] <strand.fasta> <.ovlp.tsv> <.relation.tsv>";
}

sub description {
    my $desc;
    $desc .= ucfirst(abstract) . ".\n";
    return $desc;
}

sub validate_args {
    my ( $self, $opt, $args ) = @_;

    if ( @{$args} != 3 ) {
        my $message = "This command need three input files.\n\tIt found";
        $message .= sprintf " [%s]", $_ for @{$args};
        $message .= ".\n";
        $self->usage_error($message);
    }
    for ( @{$args} ) {
        if ( !Path::Tiny::path($_)->is_file ) {
            $self->usage_error("The input file [$_] doesn't exist.");
        }
    }

    if ( $opt->{oa} ) {
        if ( !Path::Tiny::path( $opt->{oa} )->is_file ) {
            $self->usage_error("The overlap file [$opt->{oa}] doesn't exist.\n");
        }
    }

    if ( !exists $opt->{outfile} ) {
        $opt->{outfile} = Path::Tiny::path( $args->[1] )->absolute . ".contig.fasta";
    }
}

sub execute {
    my ( $self, $opt, $args ) = @_;

    #----------------------------#
    # loading sequences
    #----------------------------#
    my $seq_of = App::Fasops::Common::read_fasta( $args->[0] );
    my %is_anchor;
    for my $id ( keys %{$seq_of} ) {
        if ( index( $id, $opt->{pa} . "/" ) == 0 ) {
            $is_anchor{$id} = 1;
        }
    }

    #----------------------------#
    # load overlaps and build graph
    #----------------------------#
    my $graph = Graph->new( directed => 1 );
    my $links_of = {};    # long_id => { anchor_id => overlap_on_long, }
    {

lib/App/Dazz/Command/layout.pm  view on Meta::CPAN

            # we've orient all sequences to the same strand
            next if $info->{g_strand} == 1;

            # skip duplicated overlaps
            my $pair = join( "-", sort ( $info->{f_id}, $info->{g_id} ) );
            next if $seen_pair{$pair};
            $seen_pair{$pair}++;

            if ( $info->{f_B} > 0 ) {
                if ( $info->{f_E} == $info->{f_len} ) {

                    #          f.B        f.E
                    # f ========+---------->
                    # g         -----------+=======>
                    #          g.B        g.E
                    $graph->add_weighted_edge( $info->{f_id}, $info->{g_id},
                        $info->{g_len} - $info->{g_E} );
                }
                else {
                    #          f.B        f.E
                    # f ========+----------+=======>
                    # g         ----------->
                    #          g.B        g.E
                    $graph->add_weighted_edge( $info->{g_id}, $info->{f_id},
                        $info->{f_len} - $info->{f_E} );
                }
            }
            else {
                if ( $info->{g_E} == $info->{g_len} ) {

                    #          f.B        f.E
                    # f         -----------+=======>
                    # g ========+---------->
                    #          g.B        g.E
                    $graph->add_weighted_edge( $info->{g_id}, $info->{f_id},
                        $info->{f_len} - $info->{f_E} );
                }
                else {
                    #          f.B        f.E
                    # f         ----------->
                    # g ========+----------+=======>
                    #          g.B        g.E
                    $graph->add_weighted_edge( $info->{f_id}, $info->{g_id},
                        $info->{g_len} - $info->{g_E} );
                }
            }

            if ( $is_anchor{ $info->{f_id} } and !$is_anchor{ $info->{g_id} } ) {
                my ( $beg, $end ) = App::Dazz::Common::beg_end( $info->{g_B}, $info->{g_E}, );
                $links_of->{ $info->{g_id} }{ $info->{f_id} }
                    = AlignDB::IntSpan->new->add_pair( $beg, $end );
            }
            elsif ( $is_anchor{ $info->{g_id} } and !$is_anchor{ $info->{f_id} } ) {
                my ( $beg, $end ) = App::Dazz::Common::beg_end( $info->{f_B}, $info->{f_E}, );
                $links_of->{ $info->{f_id} }{ $info->{g_id} }
                    = AlignDB::IntSpan->new->add_pair( $beg, $end );
            }
        }
        close $in_fh;
    }
    if ( $opt->{png} ) {
        App::Dazz::Common::g2gv( $graph, $args->[1] . ".linker.png" );
    }

    #----------------------------#
    # Graph of anchors
    #----------------------------#
    my $anchor_graph = Graph->new( directed => 1 );
    {
        my @nodes = $graph->vertices;
        my @linkers = grep { !$is_anchor{$_} } @nodes;

        for my $l (@linkers) {
            my @p = grep { $is_anchor{$_} } $graph->predecessors($l);
            my @s = grep { $is_anchor{$_} } $graph->successors($l);

            for my $p (@p) {
                for my $s (@s) {
                    $anchor_graph->add_edge( $p, $s );
                }
            }

            if ( @p > 1 ) {
                @p = map { $_->[0] }
                    sort { $b->[1] <=> $a->[1] }
                    map { [ $_, $graph->get_edge_weight( $_, $l ) ] } @p;
                for my $i ( 0 .. $#p - 1 ) {
                    $anchor_graph->add_edge( $p[$i], $p[ $i + 1 ] );
                }
            }

            if ( @s > 1 ) {
                @s = map { $_->[0] }
                    sort { $a->[1] <=> $b->[1] }
                    map { [ $_, $graph->get_edge_weight( $l, $_, ) ] } @s;
                for my $i ( 0 .. $#s - 1 ) {
                    $anchor_graph->add_edge( $s[$i], $s[ $i + 1 ] );
                }
            }
        }

        # add anchors missed in $anchor_graph, i.e., not correctly linked to long.fasta
        for my $id ( keys %is_anchor ) {
            if ( !$anchor_graph->has_vertex($id) ) {
                $anchor_graph->add_vertex($id);
            }
        }

        if ( $opt->{png} ) {
            App::Dazz::Common::g2gv( $anchor_graph, $args->[1] . ".png" );
        }
        App::Dazz::Common::transitive_reduction($anchor_graph);
        if ( $opt->{png} ) {
            App::Dazz::Common::g2gv( $anchor_graph, $args->[1] . ".reduced.png" );
        }
    }

    #----------------------------#
    # existing relations
    #----------------------------#
    my $relation_of = {};
    {
        for my $line ( Path::Tiny::path( $args->[2] )->lines( { chomp => 1 } ) ) {
            my @fields = split "\t", $line;

            my $anchor_0   = $fields[0];
            my $anchor_1   = $fields[1];
            my @distances  = split ",", $fields[2];
            my @long_reads = split ",", $fields[3];

            my $pair = join "-", sort ( $anchor_0, $anchor_1, );
            $relation_of->{$pair} = {
                distances  => \@distances,
                long_reads => \@long_reads,
            };
        }
    }

    #----------------------------#
    # link anchors and break branches
    #----------------------------#
    my @paths;
    if ( $anchor_graph->is_dag ) {
        if ( scalar $anchor_graph->exterior_vertices == 2 ) {
            print STDERR "    Linear\n";

            my @ts = $anchor_graph->topological_sort;

            my %idx_of;
            for my $idx ( 0 .. $#ts ) {
                $idx_of{ $ts[$idx] } = $idx;
            }
            for my $wcc ( $anchor_graph->weakly_connected_components ) {
                my @cc = map { $_->[0] }
                    sort { $a->[1] <=> $b->[1] }
                    map { [ $_, $idx_of{$_} ] } @{$wcc};

                push @paths, \@cc;
            }
        }
        else {
            print STDERR "    Branched\n";

            my @ts = $anchor_graph->topological_sort;

            my @branchings;
            my %idx_of;
            for my $idx ( 0 .. $#ts ) {
                $idx_of{ $ts[$idx] } = $idx;

                if (   $anchor_graph->out_degree( $ts[$idx] ) > 1
                    or $anchor_graph->in_degree( $ts[$idx] ) > 1 )
                {
                    push @branchings, $ts[$idx];