<commit_before>#!/usr/bin/env perl

use strict;
use warnings;

my $usage = "\n\n\tusage: $0 LR_blastn.outfmt6 iworm.fasta [glue_setting=2]\n\n";

my $blast_outfile = $ARGV[0] or die $usage;
my $iworm_fasta = $ARGV[1];
my $glue = $ARGV[2] || 2;

main: {

    my %LR_to_iworm_hits;
    my %iworm_to_LR_hits;

    {
        open (my $fh, $blast_outfile) or die "Error, cannot open file $blast_outfile";
        while (<$fh>) {
            chomp;
            my @x = split(/\t/);
            my ($iworm_acc, $LR_acc, $bitscore) = ($x[0], $x[1], $x[11]);
            
            my $score = $iworm_to_LR_hits{$iworm_acc}->{$LR_acc};
            if ( (! defined $score) || $score < $bitscore) {
                $LR_to_iworm_hits{$LR_acc}->{$iworm_acc} = $bitscore;
                $iworm_to_LR_hits{$iworm_acc}->{$LR_acc} = $bitscore;
            }
        }
        close $fh;
    }

    my %LR_to_ordered_iworms;
    {
        foreach my $LR (keys %LR_to_iworm_hits) {

            my $LR_hits_href = $LR_to_iworm_hits{$LR};
            
            my @iworm_accs = reverse sort {$LR_hits_href->{$a} <=> $LR_hits_href->{$b}} keys %$LR_hits_href;

            $LR_to_ordered_iworms{$LR} = \@iworm_accs;
        }
    }

    
    ## get fasta index values for the iworm contigs:
    my %iworm_to_fasta_index;
    {
        my $index = 0;
        open (my $fh, $iworm_fasta) or die "Error, cannot open file $iworm_fasta";
        while (<$fh>) {
            chomp;
            if (/^>(\S+)/) {
                my $acc = $1;
                $iworm_to_fasta_index{$acc} = $index;
                $index++;
            }
        }
        close $fh;
    }

    ## Assign each iworm contig to the best-matching iworm contig linked to its best matching LR sequence.

    foreach my $iworm_acc (keys %iworm_to_LR_hits) {
        my $LR_hits_href = $iworm_to_LR_hits{$iworm_acc};
        
        my @LR_accs = reverse sort {$LR_hits_href->{$a} <=> $LR_hits_href->{$b} } keys %$LR_hits_href;
        
        my $top_LR = shift @LR_accs;
        
        my $ordered_iworms_aref = $LR_to_ordered_iworms{$top_LR};
        
        my $top_iworm_acc = $ordered_iworms_aref->[0];
        if ($iworm_acc eq $top_iworm_acc) {
            if ($#$ordered_iworms_aref > 0) {
                $top_iworm_acc = $ordered_iworms_aref->[1];
            }
            else {
                next; # no other hits
            }
        }
        
        # output record
        my $iworm_acc_A = $iworm_acc;
        my $iworm_index_A = $iworm_to_fasta_index{$iworm_acc_A};
        
        my $iworm_acc_B = $top_iworm_acc;
        my $iworm_index_B = $iworm_to_fasta_index{$iworm_acc_B};

        print join("\t", 
                   $iworm_acc_A, $iworm_index_A, 
                   $iworm_acc_B, $iworm_index_B,
                   $glue) . "\n";
    }
 

    exit(0);
}
<commit_msg>update to LR based scaffolding<commit_after>#!/usr/bin/env perl

use strict;
use warnings;
use List::Util qw(min);

my $usage = "\n\n\tusage: $0 LR_blastn.outfmt6 iworm.fasta min_containment=98 KMER_LENGTH=25\n\n";

my $blast_outfile = $ARGV[0] or die $usage;
my $iworm_fasta = $ARGV[1];
my $MIN_CONTAINMENT = $ARGV[2] || 98;
my $KMER_LENGTH = $ARGV[3] || 25;

main: {

    my %LR_to_iworm_hits;

    {
        open (my $fh, $blast_outfile) or die "Error, cannot open file $blast_outfile";
        while (<$fh>) {
            chomp;
            my @x = split(/\t/);
            my ($iworm_acc, $LR_acc, 
                $iworm_end5, $iworm_end3,
                $LR_end5, $LR_end3,
                $bitscore,
                $max_either_containment,
                ) = ($x[0], $x[1], 
                     $x[6], $x[7], 
                     $x[8], $x[9],
                     $x[11],
                     $x[16],
                    );
            
            unless ($max_either_containment >= $MIN_CONTAINMENT) { next; }
            

            my ($LR_lend, $LR_rend) = sort {$a<=>$b} ($LR_end5, $LR_end3);

            push (@{$LR_to_iworm_hits{$LR_acc}}, { iworm_acc => $iworm_acc,
                                                   lend => $LR_lend,
                                                   rend => $LR_rend,
                  });
            
            
            
        }
        close $fh;
    }

    
    my %iworm_pairs;
    
    foreach my $LR (keys %LR_to_iworm_hits) {

        my @iworm_hits = sort {$a->{lend}<=>$b->{lend}
                               ||
                                   $a->{rend} <=> $b->{rend}} @{$LR_to_iworm_hits{$LR}};

        
        for (my $i = 0; $i < $#iworm_hits; $i++) {

            my $iworm_i = $iworm_hits[$i];

            for (my $j = $i + 1; $j <= $#iworm_hits; $j++) {

                my $iworm_j = $iworm_hits[$j];

                if ($iworm_j->{lend} > $iworm_i->{rend}) {
                    last; 
                }

                if (&overlap_by_Kminus1($iworm_i, $iworm_j, $KMER_LENGTH)) {
                    
                    my $pair_token = join("$;", sort ($iworm_i->{iworm_acc}, $iworm_j->{iworm_acc}));
                    $iworm_pairs{$pair_token}++;
                    
                }
            }
        }

    }
    

    ## get fasta index values for the iworm contigs:
    my %iworm_to_fasta_index;
    {
        my $index = 0;
        open (my $fh, $iworm_fasta) or die "Error, cannot open file $iworm_fasta";
        while (<$fh>) {
            chomp;
            if (/^>(\S+)/) {
                my $acc = $1;
                $iworm_to_fasta_index{$acc} = $index;
                $index++;
            }
        }
        close $fh;
    }

    
    foreach my $pair (keys %iworm_pairs) {

        my ($iworm_acc_A, $iworm_acc_B) = split(/$;/, $pair);
        
        # a179572;19
        
        # assign glue as the lower of the coverage values:

        my ($coreA, $covA) = split(/;/, $iworm_acc_A);
        my ($coreB, $covB) = split(/;/, $iworm_acc_B);
        
        my $glue = min($covA, $covB);
        
        # output record
        my $iworm_index_A = $iworm_to_fasta_index{$iworm_acc_A};
        
        my $iworm_index_B = $iworm_to_fasta_index{$iworm_acc_B};

                
        print join("\t", 
                   $iworm_acc_A, $iworm_index_A, 
                   $iworm_acc_B, $iworm_index_B,
                   $glue) . "\n";
    }



    exit(0);
    
}


####
sub overlap_by_Kminus1 {
    my ($iworm_i, $iworm_j, $K) = @_;

    ## require:
    ##
    ##   -----------------   iworm_i
    ##             -------------------   iworm_j
    ##
    ## ----------------------------------  LR
    ##
    ##  overlap by at least K-1
    ##  staggered overlaps along LR
    ##  iworm_i and iworm_j alignment lengths >= 2 * (K-1)

    
    my $iworm_i_align_len = $iworm_i->{rend} - $iworm_i->{lend} + 1;
    my $iworm_j_align_len = $iworm_j->{rend} - $iworm_j->{lend} + 1;

    my $min_align_len = 2 * ($K-1);
    
    if ($iworm_i_align_len < $min_align_len 
        ||
        $iworm_j_align_len < $min_align_len) { 
       
        return(0);
    }

    if ($iworm_i->{rend} > $iworm_j->{rend}) {
        return(0);
    }
    
    
    my $overlap_len = $iworm_i->{rend} - $iworm_j->{lend} + 1;
    
    if ($overlap_len >= $K - 1) {
        
        return(1);
    }
    else {
        return(0);
    }
    
}

{'fraction_non_alphanumeric': 0.15292025714676516, 'fraction_numerical': 0.010258514567090685, 'mean_word_length': 1.4495812395309882, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6337986', 'n_tokens_mistral': 2900, 'n_tokens_neox': 2676, 'n_words': 505}