view release on metacpan or  search on metacpan

lib/App/Anchr/Common.pm  view on Meta::CPAN

    my $cmd = sprintf "DBshow -n %s %s ", $dazz_db, join( " ", @{$serials} );
    my @lines = map { $_ =~ s/^>//; $_; } grep {defined} split /\n/, `$cmd`;

    my %name_of;
    for my $i ( 0 .. $#lines ) {
        $name_of{ $serials->[$i] } = $lines[$i];
    }

    return \%name_of;
}

sub judge_distance {
    my $d_ref = shift;
    my $max_dis = shift || 5000;

    return 0 unless defined $d_ref;

    my $sum = 0;
    my $min = $d_ref->[0];
    my $max = $min;
    for my $d ( @{$d_ref} ) {
        $sum += $d;
        if ( $d < $min ) { $min = $d; }
        if ( $d > $max ) { $max = $d; }
    }
    my $avg = $sum / scalar( @{$d_ref} );
    return 0 if $avg > $max_dis;

    # max k-mer is 127.
    # For k-unitigs, overlaps are less than k-mer
    my $v = $max - $min;
    if ( $v < 20 or abs( $v / $avg ) < 0.2 ) {
        return 1;
    }
    else {
        return 0;
    }
}

sub g2gv {

    #@type Graph
    my $g  = shift;
    my $fn = shift;

    my $gv = GraphViz->new( directed => 1 );

    for my $v ( $g->vertices ) {
        $gv->add_node($v);
    }

    for my $e ( $g->edges ) {
        if ( $g->has_edge_weight( @{$e} ) ) {
            $gv->add_edge( @{$e}, label => $g->get_edge_weight( @{$e} ) );
        }
        else {
            $gv->add_edge( @{$e} );
        }
    }

    Path::Tiny::path($fn)->spew_raw( $gv->as_png );
}

sub g2gv0 {

    #@type Graph
    my $g  = shift;
    my $fn = shift;

    my $gv = GraphViz->new( directed => 0 );

    for my $v ( $g->vertices ) {
        $gv->add_node($v);
    }

    for my $e ( $g->edges ) {
        $gv->add_edge( @{$e} );
    }

    Path::Tiny::path($fn)->spew_raw( $gv->as_png );
}

sub transitive_reduction {

    #@type Graph
    my $g = shift;

    my $count = 0;
    my $prev_count;
    while (1) {
        last if defined $prev_count and $prev_count == $count;
        $prev_count = $count;

        for my $v ( $g->vertices ) {
            next if $g->out_degree($v) < 2;

            my @s = sort { $a cmp $b } $g->successors($v);
            for my $i ( 0 .. $#s ) {
                for my $j ( 0 .. $#s ) {
                    next if $i == $j;
                    if ( $g->is_reachable( $s[$i], $s[$j] ) ) {
                        $g->delete_edge( $v, $s[$j] );

                        $count++;
                    }
                }
            }
        }
    }

    return $count;
}

sub poa_consensus {
    my $seq_refs = shift;

    my $aln_prog = "poa";

    # temp in and out
    my $temp_in  = Path::Tiny->tempfile("seq_in_XXXXXXXX");
    my $temp_out = Path::Tiny->tempfile("seq_out_XXXXXXXX");

    # msa may change the order of sequences
    my @indexes = 0 .. scalar( @{$seq_refs} - 1 );
    {
        my $fh = $temp_in->openw();
        for my $i (@indexes) {
            printf {$fh} ">seq_%d\n", $i;
            printf {$fh} "%s\n",      $seq_refs->[$i];
        }
        close $fh;
    }

    my @args;
    {
        push @args, "-hb";
        push @args, "-read_fasta " . $temp_in->absolute->stringify;
        push @args, "-clustal " . $temp_out->absolute->stringify;
        push @args, File::ShareDir::dist_file( 'App-Anchr', 'poa-blosum80.mat' );
    }