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

use strict;
use warnings;
use Carp;
use Data::Dumper;

use FindBin;
use lib "$FindBin::Bin/../PerlLib";
require "overlapping_nucs.ph";

my $MAX_ALLOWED_PCT_OVERLAP = 50;
my $MIN_PER_ID = 90;

my $DEBUG = 1;

my $usage = "\n\n\tusage: $0 blat.gff3 cdna_targets.fasta\n\n";

my $blat_gff3 = $ARGV[0] or die $usage;
my $cdna_targets_fasta = $ARGV[1] or die $usage;

main: {


    my $headers_file = "$cdna_targets_fasta.headers";
    unless (-s $headers_file) {
        die "Error, cannot locate $headers_file";
    }

    my %trans_id_to_gene_id = &get_trans_to_gene_mapping($headers_file);
     
    my %trans_to_gene_alignments = &get_trans_to_gene_alignments($blat_gff3, \%trans_id_to_gene_id);
    
    &report_fusion_candidates(\%trans_to_gene_alignments);
    



    if (0) {
        
        ## Tier and merge strategy

    
        my %tiered_n_merged_gene_alignments = &tier_and_merge_best_alignments(\%trans_to_gene_alignments);
        
        foreach my $trans_id (keys %tiered_n_merged_gene_alignments) {
            
            my @hits = @{$tiered_n_merged_gene_alignments{$trans_id}};
            @hits = sort {$a->{trans_lend}<=>$b->{trans_lend}} @hits;
            
            if (scalar @hits < 2) {
                # shouldn't happen as filtering is done before this.
                die "Error, have fewer than 2 candidate fusion parts";
            }
            
            ## Report fusion pairs:
            my $geneA = shift @hits;
            while (@hits) {
                my $geneB = shift @hits;
                
                if ($geneA->{gene_orient} eq $geneB->{gene_orient}) {
                    
                    my ($left_gene, $right_gene) = ($geneA->{gene_orient} eq '+') ? ($geneA, $geneB) : ($geneB, $geneA);
                    
                    my $fusion_name = join("--", $left_gene->{gene_name}, $right_gene->{gene_name});
                    
                    print join("\t", $fusion_name, $trans_id,
                               $left_gene->{gene_name}, $left_gene->{trans_lend}, $left_gene->{trans_rend}, $left_gene->{gene_orient}, $left_gene->{per_id} . "\%ID",
                               $right_gene->{gene_name}, $right_gene->{trans_lend}, $right_gene->{trans_rend}, $right_gene->{gene_orient}, $right_gene->{per_id} . "\%ID") . "\n";
                    
                }
                
                $geneA = $geneB;
            }
        }
    }
    
    exit(0);


}


####
sub report_fusion_candidates {
    my ($trans_to_gene_alignments_href) = @_;

    foreach my $trans_id (keys %$trans_to_gene_alignments_href) {
        
        my @aligns = @{$trans_to_gene_alignments_href->{$trans_id}};
        
        @aligns = sort {$a->{trans_lend}<=>$b->{trans_lend}} @aligns;

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

            my $align_i = $aligns[$i];
            
            my $gene_name_i = $align_i->{gene_name};
            my $lend_i = $align_i->{trans_lend};
            my $rend_i = $align_i->{trans_rend};
            my $per_id_i = $align_i->{per_id};

            for (my $j = $i + 1; $j <= $#aligns; $j++) {
                
                my $align_j = $aligns[$j];
                
                my $gene_name_j = $align_j->{gene_name};
                my $lend_j = $align_j->{trans_lend};
                my $rend_j = $align_j->{trans_rend};
                my $per_id_j = $align_j->{per_id};

                my $pct_overlap = &get_percent_overlap([$lend_i, $rend_i], [$lend_j, $rend_j]);

                if ($pct_overlap < $MAX_ALLOWED_PCT_OVERLAP) {

                    my ($left_gene, $right_gene) = ($align_i->{gene_orient} eq '+') ? ($align_i, $align_j) : ($align_j, $align_i);
                    
                    ## report candidate
                    my $fusion_name = join("--", $left_gene->{gene_name}, $right_gene->{gene_name});
                    
                    print join("\t", $fusion_name, $trans_id,
                               $left_gene->{gene_name}, $left_gene->{trans_lend}, $left_gene->{trans_rend}, $left_gene->{gene_orient}, $left_gene->{per_id} . "\%ID",
                               $right_gene->{gene_name}, $right_gene->{trans_lend}, $right_gene->{trans_rend}, $right_gene->{gene_orient}, $right_gene->{per_id} . "\%ID") . "\n";
                    

                }
            }
        }
    }
}


####
sub get_percent_overlap {
    my ($coords_A_aref, $coords_B_aref) = @_;

    if (&coordsets_overlap($coords_A_aref, $coords_B_aref)) {
        
        my $len_A = $coords_A_aref->[1] - $coords_A_aref->[0] + 1;
        my $len_B = $coords_B_aref->[1] - $coords_B_aref->[0] + 1;
        
        my $shorter_len = ($len_A < $len_B) ? $len_A : $len_B;

        my $overlap_len = &nucs_in_common($coords_A_aref->[0], $coords_A_aref->[1], $coords_B_aref->[0], $coords_B_aref->[1]);<commit_msg>reduce redundancy in output, require min anchor 30<commit_after>#!/usr/bin/env perl

use strict;
use warnings;
use Carp;
use Data::Dumper;

use FindBin;
use lib "$FindBin::Bin/../PerlLib";
require "overlapping_nucs.ph";

my $MAX_ALLOWED_PCT_OVERLAP = 10;
my $MIN_PER_ID = 90;
my $MIN_ANCHOR_LEN = 30;

my $DEBUG = 1;

my $usage = "\n\n\tusage: $0 blat.gff3 cdna_targets.fasta\n\n";

my $blat_gff3 = $ARGV[0] or die $usage;
my $cdna_targets_fasta = $ARGV[1] or die $usage;

main: {


    my $headers_file = "$cdna_targets_fasta.headers";
    unless (-s $headers_file) {
        die "Error, cannot locate $headers_file";
    }

    my %trans_id_to_gene_id = &get_trans_to_gene_mapping($headers_file);
     
    my %trans_to_gene_alignments = &get_trans_to_gene_alignments($blat_gff3, \%trans_id_to_gene_id);
    
    &report_fusion_candidates(\%trans_to_gene_alignments);
    



    if (0) {
        
        ## Tier and merge strategy

    
        my %tiered_n_merged_gene_alignments = &tier_and_merge_best_alignments(\%trans_to_gene_alignments);
        
        foreach my $trans_id (keys %tiered_n_merged_gene_alignments) {
            
            my @hits = @{$tiered_n_merged_gene_alignments{$trans_id}};
            @hits = sort {$a->{trans_lend}<=>$b->{trans_lend}} @hits;
            
            if (scalar @hits < 2) {
                # shouldn't happen as filtering is done before this.
                die "Error, have fewer than 2 candidate fusion parts";
            }
            
            ## Report fusion pairs:
            my $geneA = shift @hits;
            while (@hits) {
                my $geneB = shift @hits;
                
                if ($geneA->{gene_orient} eq $geneB->{gene_orient}) {
                    
                    my ($left_gene, $right_gene) = ($geneA->{gene_orient} eq '+') ? ($geneA, $geneB) : ($geneB, $geneA);
                    
                    my $fusion_name = join("--", $left_gene->{gene_name}, $right_gene->{gene_name});
                    
                    print join("\t", $fusion_name, $trans_id,
                               $left_gene->{gene_name}, $left_gene->{trans_lend}, $left_gene->{trans_rend}, $left_gene->{gene_orient}, $left_gene->{per_id} . "\%ID",
                               $right_gene->{gene_name}, $right_gene->{trans_lend}, $right_gene->{trans_rend}, $right_gene->{gene_orient}, $right_gene->{per_id} . "\%ID") . "\n";
                    
                }
                
                $geneA = $geneB;
            }
        }
    }
    
    exit(0);


}


####
sub report_fusion_candidates {
    my ($trans_to_gene_alignments_href) = @_;

    foreach my $trans_id (keys %$trans_to_gene_alignments_href) {
        
        my @aligns = @{$trans_to_gene_alignments_href->{$trans_id}};
        
        @aligns = sort {$a->{trans_lend}<=>$b->{trans_lend}} @aligns;

        my @fusion_preds;

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

            my $align_i = $aligns[$i];
            
            my $gene_name_i = $align_i->{gene_name};
            my $lend_i = $align_i->{trans_lend};
            my $rend_i = $align_i->{trans_rend};
            my $per_id_i = $align_i->{per_id};
            my $seg_len_i = $rend_i - $lend_i + 1;
            
            if ($seg_len_i < $MIN_ANCHOR_LEN) { next; }


            for (my $j = $i + 1; $j <= $#aligns; $j++) {
                
                my $align_j = $aligns[$j];
                
                my $gene_name_j = $align_j->{gene_name};
                my $lend_j = $align_j->{trans_lend};
                my $rend_j = $align_j->{trans_rend};
                my $per_id_j = $align_j->{per_id};
                my $seg_len_j = $rend_j - $lend_j + 1;


                if ($seg_len_j < $MIN_ANCHOR_LEN) { next; }

                my $pct_overlap = &get_percent_overlap([$lend_i, $rend_i], [$lend_j, $rend_j]);

                my $delta = abs($lend_j -1 - $rend_i);
                my ($overlap_len, $gap_len) = ($lend_j >= $rend_i) ? ($delta, 0) : (0, $delta);
                

                if ($pct_overlap < $MAX_ALLOWED_PCT_OVERLAP) {

                    my ($left_gene, $right_gene) = ($align_i->{gene_orient} eq '+') ? ($align_i, $align_j) : ($align_j, $align_i);
                    
                    my $score = ($seg_len_i * $per_id_i) + ($seg_len_j * $per_id_j) 
                        - ($overlap_len * $per_id_i/2) - ($overlap_len * $per_id_j/2)
                        - ($gap_len * 100); 

                    push (@fusion_preds, [$left_gene, $right_gene, $score]);
                    
                }
            }
        }
        
        @fusion_preds = reverse sort {$a->[2]<=>$b->[2]} @fusion_preds;

        my %seen;

        foreach my $fusion_pred (@fusion_preds) {

            my ($left_gene, $right_gene, $score) = @$fusion_pred;
            
            ## report candidate
            my $fusion_name = join("--", $left_gene->{gene_name}, $right_gene->{gene_name});
            
            if ($seen{$fusion_name}) { next; }
            $seen{$fusion_name} = 1;
            
            print join("\t", $fusion_name, $trans_id,
                       $left_gene->{gene_name}, $left_gene->{trans_lend}, $left_gene->{trans_rend}, $left_gene->{gene_orient}, $left_gene->{per_id} . "\%ID",
                       $right_gene->{gene_name}, $right_gene->{trans_lend}, $right_gene->{trans_rend}, $right_gene->{gene_orient}, $right_gene->{per_id} . "\%ID",
                $score) . "\n";
            
        }
    }
}


####
sub get_percent_overlap {
    my ($coords_A_aref, $coords_B_aref) = @_;

    if (&coordsets_overlap($coords_A_aref, $coords_B_aref)) {
        
        my $len_A = $coords_A_aref->[1] - $coords_A_aref->[0] + 1;
        my $len_B = $coords_B_aref->[1] - $coords_B_aref->[0] + 1;
        
        my $shorter_len = ($len_A < $len_B) ? $len_A : $len_B;

        my $overlap_len = &nucs_in_common($coords_A_aref->[0], $coords_A_aref->[1], $coords_B_aref->[0], $coords_B_aref->[1]);
{'fraction_non_alphanumeric': 0.15428361768482848, 'fraction_numerical': 0.006512566501559347, 'mean_word_length': 1.866947146989219, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 109, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3837961', 'n_tokens_mistral': 4248, 'n_tokens_neox': 3875, 'n_words': 664}