view release on metacpan or  search on metacpan

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

                    $full_set_of->{$g_id} = AlignDB::IntSpan->new->add_pair( 1, $g_len );
                }
            }
            elsif ( $anchor_range->contains($g_id) and !$anchor_range->contains($f_id) ) {
                my ( $beg, $end ) = App::Anchr::Common::beg_end( $f_B, $f_E, );
                $links_of->{$f_id}{$g_id} = AlignDB::IntSpan->new->add_pair( $beg, $end );

                if ( $contained eq "contained" ) {
                    $contained_long->add($f_id);
                }

                if ( !exists $full_set_of->{$f_id} ) {
                    $full_set_of->{$f_id} = AlignDB::IntSpan->new->add_pair( 1, $f_len );
                }
            }
        }
        close $in_fh;
    }

    for my $long_id ( $multi_matched->as_array ) {
        if ( exists $links_of->{$long_id} ) {
            delete $links_of->{$long_id};
        }
    }

    for my $long_id ( $contained_long->as_array ) {
        if ( exists $links_of->{$long_id} ) {
            delete $links_of->{$long_id};
        }
    }

    my $region_of = {};    # long_id => IntSpan
    for my $long_id ( keys %{$links_of} ) {

        # covered parts of long reads
        my $covered = AlignDB::IntSpan->new;
        for my $anchor_id ( keys %{ $links_of->{$long_id} } ) {

            #@type AlignDB::IntSpan
            my $set = $links_of->{$long_id}{$anchor_id};
            $set = $set->trim( $opt->{border} );
            $covered->add($set);    # avoid overlapped anchors
        }

        my $multiplied = int( $opt->{len} * $opt->{power} );
        my $rest_set
            = $full_set_of->{$long_id}->diff($covered)->pad($multiplied);
        $region_of->{$long_id} = $full_set_of->{$long_id}->intersect($rest_set);
    }

    {
        $tempdir->child("break.fasta")->remove;
        for my $serial ( sort { $a <=> $b } keys %{$region_of} ) {

            #@type AlignDB::IntSpan
            my $region = $region_of->{$serial};

            for my $set ( $region->sets ) {
                my $cmd;
                $cmd .= "DBshow -U $fn_dazz $serial";
                $cmd .= " | faops frag -l 0 stdin @{[$set->min]} @{[$set->max]} stdout";
                $cmd .= " >> break.fasta";
                App::Anchr::Common::exec_cmd( $cmd, { verbose => $opt->{verbose}, } );
            }
        }

        if ( !$tempdir->child("break.fasta")->is_file ) {
            Carp::croak "Failed: create break.fasta\n";
        }

        YAML::Syck::DumpFile(
            "break.yml",
            {   "Contained"           => $contained_long->runlist,
                "Contained count"     => $contained_long->size,
                "Multi-matched"       => $multi_matched->runlist,
                "Multi-matched count" => $multi_matched->size,
                "region_of" => { map { $_ => $region_of->{$_}->runlist } keys %{$region_of} },
            }
        );
    }

    {
        # Outputs. stdout is handeld by faops
        my $cmd;
        $cmd .= "faops filter -l 0 break.fasta";
        $cmd .= " $opt->{outfile}";
        App::Anchr::Common::exec_cmd( $cmd, { verbose => $opt->{verbose}, } );

        $tempdir->child("break.yml")->copy("$opt->{outfile}.break.yml");
    }

    chdir $cwd;
}

1;