view release on metacpan or  search on metacpan

lib/App/Rangeops/Command/clean.pm  view on Meta::CPAN

                if ( !$graph->has_vertex($line_j) ) {
                    $graph->add_vertex($line_j);
                }

                my $intspan0_i = $info_of->{$range0_i}{intspan};
                my $intspan1_i = $info_of->{$range1_i}{intspan};

                my $intspan0_j = $info_of->{$range0_j}{intspan};
                my $intspan1_j = $info_of->{$range1_j}{intspan};

                if (    $intspan0_i->overlap($intspan0_j) >= $opt->{bundle}
                    and $intspan1_i->overlap($intspan1_j) >= $opt->{bundle} )
                {
                    $graph->add_edge( $line_i, $line_j );
                }
            }
        }

        # bundle connected lines
        my @cc = grep { scalar @{$_} > 1 } $graph->connected_components;
        for my $c (@cc) {
            printf STDERR "\n" . " " x 4 . "Merge %s lines\n", scalar @{$c}
                if $opt->{verbose};

            my @merged_range;

            for my $i ( 0, 1 ) {
                my ( $chr, $strand );
                my $merge_intspan = AlignDB::IntSpan->new;

                for my $line ( @{$c} ) {
                    @lines = grep { $_ ne $line } @lines;
                    my $range = ( split /\t/, $line )[$i];
                    $chr    = $info_of->{$range}{chr};
                    $strand = $info_of->{$range}{strand};
                    $merge_intspan->merge( $info_of->{$range}{intspan} );
                }

                $merged_range[$i] = App::RL::Common::encode_header(
                    {   chr    => $chr,
                        strand => $strand,
                        start  => $merge_intspan->min,
                        end    => $merge_intspan->max,
                    }
                );
            }

            my $line = join "\t", @merged_range;
            push @lines, $line;
            print STDERR " " x 8 . "$line\n" if $opt->{verbose};

            # new range introduced, update $info_of
            $info_of = App::Rangeops::Common::build_info_intspan( [$line],
                $info_of );
        }

        @lines = @{ App::Rangeops::Common::sort_links( \@lines ) };
    }

    #----------------------------#
    # Links of nearly identical ranges escaped from merging
    #----------------------------#
    if ( $opt->{replace} ) {
        print STDERR "==> Remove self links\n" if $opt->{verbose};

        my @same_pair_lines = map {
            my ( $r0, $r1 ) = split /\t/;
            $info_of->{$r0}{chr} eq $info_of->{$r1}{chr} ? ($_) : ()
        } @lines;

        for my $line (@same_pair_lines) {
            my ( $range0, $range1 ) = split /\t/, $line;

            my $intspan0 = $info_of->{$range0}{intspan};
            my $intspan1 = $info_of->{$range1}{intspan};

            my $intspan_i = $intspan0->intersect($intspan1);
            if ( $intspan_i->is_not_empty ) {
                if (    $intspan_i->size / $intspan0->size > 0.5
                    and $intspan_i->size / $intspan1->size > 0.5 )
                {
                    @lines = grep { $_ ne $line } @lines;
                }
            }
        }
    }

    #----------------------------#
    # Output
    #----------------------------#
    my $out_fh;
    if ( lc( $opt->{outfile} ) eq "stdout" ) {
        $out_fh = \*STDOUT;
    }
    else {
        open $out_fh, ">", $opt->{outfile};
    }

    print {$out_fh} "$_\n" for @lines;

    close $out_fh;
}

1;