File Coverage

blib/lib/App/Fasops/Command/links.pm
Criterion Covered Total %
statement 98 101 97.0
branch 18 24 75.0
condition 6 9 66.6
subroutine 13 13 100.0
pod 6 6 100.0
total 141 153 92.1


line stmt bran cond sub pod time code
1             package App::Fasops::Command::links;
2 21     21   14767 use strict;
  21         46  
  21         643  
3 21     21   151 use warnings;
  21         43  
  21         529  
4 21     21   107 use autodie;
  21         41  
  21         113  
5              
6 21     21   106618 use App::Fasops -command;
  21         57  
  21         228  
7 21     21   7454 use App::RL::Common;
  21         50  
  21         574  
8 21     21   122 use App::Fasops::Common;
  21         64  
  21         26946  
9              
10             sub abstract {
11 2     2 1 53 return 'scan blocked fasta files and output bi/multi-lateral range links';
12             }
13              
14             sub opt_spec {
15             return (
16 6     6 1 57 [ "outfile|o=s", "Output filename. [stdout] for screen." ],
17             [ "pair|p", "pairwise links" ],
18             [ "best|b", "best-to-best pairwise links" ],
19             { show_defaults => 1, }
20             );
21             }
22              
23             sub usage_desc {
24 6     6 1 59172 return "fasops links [options] [more infiles]";
25             }
26              
27             sub description {
28 1     1 1 1001 my $desc;
29 1         5 $desc .= ucfirst(abstract) . ".\n";
30 1         3 $desc .= <<'MARKDOWN';
31              
32             * are paths to axt files, .axt.gz is supported
33             * infile == stdin means reading from STDIN
34              
35             MARKDOWN
36              
37 1         4 return $desc;
38             }
39              
40             sub validate_args {
41 5     5 1 5388 my ( $self, $opt, $args ) = @_;
42              
43 5 100       12 if ( !@{$args} ) {
  5         25  
44 1         3 my $message = "This command need one or more input files.\n\tIt found";
45 1         2 $message .= sprintf " [%s]", $_ for @{$args};
  1         3  
46 1         3 $message .= ".\n";
47 1         10 $self->usage_error($message);
48             }
49 4         9 for ( @{$args} ) {
  4         12  
50 4 50       15 next if lc $_ eq "stdin";
51 4 100       30 if ( !Path::Tiny::path($_)->is_file ) {
52 1         134 $self->usage_error("The input file [$_] doesn't exist.");
53             }
54             }
55              
56 3 50       301 if ( !exists $opt->{outfile} ) {
57 0         0 $opt->{outfile} = Path::Tiny::path( $args->[0] )->absolute . ".tsv";
58             }
59             }
60              
61             sub execute {
62 3     3 1 23 my ( $self, $opt, $args ) = @_;
63              
64 3         7 my @links;
65 3         7 for my $infile ( @{$args} ) {
  3         10  
66 3         6 my $in_fh;
67 3 50       13 if ( lc $infile eq "stdin" ) {
68 0         0 $in_fh = *STDIN{IO};
69             }
70             else {
71 3         25 $in_fh = IO::Zlib->new( $infile, "rb" );
72             }
73              
74 3         5631 my $content = ''; # content of one block
75 3         7 while (1) {
76 84 100 66     344 last if $in_fh->eof and $content eq '';
77 81         3527 my $line = '';
78 81 50       251 if ( !$in_fh->eof ) {
79 81         3180 $line = $in_fh->getline;
80             }
81 81 50       9682 next if substr( $line, 0, 1 ) eq "#";
82              
83 81 100 66     480 if ( ( $line eq '' or $line =~ /^\s+$/ ) and $content ne '' ) {
      66        
84 9         47 my $info_of = App::Fasops::Common::parse_block( $content, 1 );
85 9         16 $content = '';
86              
87 9         16 my @headers = keys %{$info_of};
  9         48  
88              
89 9 100       266 if ( $opt->{best} ) {
    100          
90 3         9 my @matrix = map { [ (undef) x ( scalar @headers ) ] } 0 .. $#headers;
  12         30  
91              
92             # distance is 0 for same sequence
93 3         12 for my $i ( 0 .. $#headers ) {
94 12         20 $matrix[$i][$i] = 0;
95             }
96              
97             # compute a triangle, fill full matrix
98 3         12 for ( my $i = 0; $i <= $#headers; $i++ ) {
99 12         39 for ( my $j = $i + 1; $j <= $#headers; $j++ ) {
100             my $D = App::Fasops::Common::pair_D(
101             [ $info_of->{ $headers[$i] }{seq},
102             $info_of->{ $headers[$j] }{seq},
103 18         82 ]
104             );
105 18         72 $matrix[$i][$j] = $D;
106 18         58 $matrix[$j][$i] = $D;
107             }
108             }
109              
110             # print YAML::Syck::Dump \@matrix;
111              
112             # best_pairwise
113 3         7 my @pair_ary;
114 3         7 for my $i ( 0 .. $#headers ) {
115 12         21 my @row = @{ $matrix[$i] };
  12         28  
116 12         16 $row[$i] = 999; # remove the score (zero) of this item
117 12         31 my $min = List::Util::min(@row);
118 12     19   64 my $min_idx = App::Fasops::Common::firstidx { $_ == $min } @row;
  19         68  
119              
120             # to remove duplications of a:b and b:a
121 12         53 push @pair_ary, join ":", sort { $a <=> $b } ( $i, $min_idx );
  12         50  
122             }
123 3         13 @pair_ary = App::Fasops::Common::uniq(@pair_ary);
124              
125 3         7 for (@pair_ary) {
126 9         25 my ( $i, $j ) = split ":";
127 9         66 push @links, [ $headers[$i], $headers[$j] ];
128             }
129             }
130             elsif ( $opt->{pair} ) {
131 3         22 for ( my $i = 0; $i <= $#headers; $i++ ) {
132 12         66 for ( my $j = $i + 1; $j <= $#headers; $j++ ) {
133 18         56 push @links, [ $headers[$i], $headers[$j] ];
134             }
135             }
136             }
137             else {
138 3         39 push @links, \@headers;
139             }
140             }
141             else {
142 72         149 $content .= $line;
143             }
144             }
145              
146 3         526 $in_fh->close;
147             }
148              
149 3         605 my $out_fh;
150 3 50       19 if ( lc( $opt->{outfile} ) eq "stdout" ) {
151 3         11 $out_fh = *STDOUT{IO};
152             }
153             else {
154 0         0 open $out_fh, ">", $opt->{outfile};
155             }
156              
157 3         10 for my $link (@links) {
158 30         429 printf {$out_fh} join( "\t", @{$link} );
  30         69  
  30         105  
159 30         549 printf {$out_fh} "\n";
  30         82  
160             }
161 3         62 close $out_fh;
162             }
163              
164             1;