view release on metacpan or  search on metacpan

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

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

use App::Anchr -command;
use App::Anchr::Common;

use constant abstract => "layout anthor group";

sub opt_spec {
    return (
        [ "outfile|o=s", "output filename", ],
        [ 'border=i', 'length of borders in anchors', { default => 100 }, ],
        [ "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 "anchr layout [options] <.ovlp.tsv> <.relation.tsv> <strand.fasta>";
}

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 one or more 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->[0] )->absolute . ".contig.fasta";
    }
}

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

    #----------------------------#
    # load overlaps and build graph
    #----------------------------#
    my $graph = Graph->new( directed => 1 );
    my %is_anchor;
    my $links_of = {};    # long_id => { anchor_id => overlap_on_long, }
    {
        open my $in_fh, "<", $args->[0];

        my %seen_pair;
        while ( my $line = <$in_fh> ) {
            my $info = App::Anchr::Common::parse_ovlp_line($line);

            # ignore self overlapping
            next if $info->{f_id} eq $info->{g_id};

            # we've orient all sequences to the same strand

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

                    # 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::Anchr::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::Anchr::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;
    }

    #----------------------------#
    # 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 ] );
                }
            }
        }

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

    #----------------------------#
    # existing relations
    #----------------------------#
    my $relation_of = {};
    {
        for my $line ( Path::Tiny::path( $args->[1] )->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,
            };
        }
    }

    #----------------------------#
    # loading sequences
    #----------------------------#
    my $seq_of = App::Fasops::Common::read_fasta( $args->[2] );

    #----------------------------#
    # 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;
            push @paths, \@ts;
        }
        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];
                }
            }

            $anchor_graph->delete_vertex($_) for @branchings;

            for my $wcc ( $anchor_graph->weakly_connected_components ) {