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

use strict;
use warnings;

use Data::Dumper;
use Getopt::Long qw(:config no_ignore_case bundling pass_through);
use List::Util qw (min max);
use Set::IntervalTree;

my $usage = <<__EOUSAGE__;

#####################################################################
#
# Required:
#
#  --gmap_gff3 <string>        gmap gff3 alignment output
#
#  --annot_gtf <string>        transcript structures in gtf file format.
#
######################################################################


__EOUSAGE__

    ;


my $help_flag;
my $gmap_gff3_file;
my $annot_gtf_file;

&GetOptions ( 'h' => \$help_flag,
              'gmap_gff3=s' => \$gmap_gff3_file,
              'annot_gtf=s' => \$annot_gtf_file,
              
              );


if ($help_flag) {
    die $usage;
}

unless ($gmap_gff3_file && $annot_gtf_file) {
    die $usage;
}


my %genes;
my %interval_trees;


main: {

    my %chr_to_gene_coords;
    print STDERR "-parsing $annot_gtf_file\n";
    open (my $fh, $annot_gtf_file) or die "Error, cannot open file $annot_gtf_file";
    while (<$fh>) {
        chomp;
        unless (/\w/) { next; }
        if (/^\#/) { next; }
        s/^>//;
        my @x = split(/\t/);
        
        unless ($x[2] eq "exon") { next; }
        
        my $info = $x[8];
        $info =~ /gene_id \"([^\"]+)/ or die "Error, cannot extract gene_id from $_ [specifically from: $info]";
        my $gene_id = $1 or die "Error, no gene_id from $_";
        
        if ($info =~ /gene_name \"([^\"]+)/) {
            # use gene name instead
            $gene_id = $1;
        }
        
        $info =~ /transcript_id \"([^\"]+)/ or die "Error, cannot extract transcript_id from $_";
        my $transcript_id = $1 or die "Error, no trans id from $_";
        
        my ($lend, $rend) = ($x[3], $x[4]);
        my $chr = $x[0];
        my $orient = $x[6];
        
        
        push (@{$genes{$chr}->{$gene_id}->{$transcript_id}}, { 
            
            gene => $gene_id,
            transcript => $transcript_id,
            chr => $chr,
            lend => $lend,
            rend => $rend,
            orient => $orient,
              }
            );
        
        push (@{$chr_to_gene_coords{$chr}->{$gene_id}}, $lend, $rend);
        
    }
    close $fh;
    
    print STDERR "-building interval tree for fast searching of gene overlaps\n";
    ## Build interval trees
    foreach my $chr (keys %chr_to_gene_coords) {

        my $i_tree = $interval_trees{$chr} = Set::IntervalTree->new;
        
        foreach my $gene_id (keys %{$chr_to_gene_coords{$chr}}) {
            
            my @coords = sort {$a<=>$b} @{$chr_to_gene_coords{$chr}->{$gene_id}};
            my $lend = shift @coords;
            my $rend = pop @coords;

            $i_tree->insert($gene_id, $lend, $rend);
        }
    }
    
    
    my %target_to_aligns;
    print STDERR "-loading alignment data\n";
    open ($fh, $gmap_gff3_file) or die $!;
    while (<$fh>) {
        if (/^\#/) { next; }
        unless (/\w/) { next; }

        chomp;
        my ($chr, $filename, $type, $lend, $rend, $per_id, $orient, $dot, $info) = split(/\t/);
        
        my %info_hash;
        foreach my $keyval (split(/;/, $info)) {
            my ($key, $val) = split(/=/, $keyval);
            $info_hash{$key} = $val;
        }

        my $alignment_ID = $info_hash{ID};
        
        my ($target, $range_lend, $range_rend) = split(/\s+/, $info_hash{Target});

        push (@{$target_to_aligns{$target}->{$alignment_ID}}, { chr => $chr,
                                                                lend => $lend,
                                                                rend => $rend,
                                                                orient => $orient,
                                                                per_id => $per_id,
                                                                
                                                                range_lend => $range_lend,
                                                                range_rend => $range_rend,
                                                            });
        
    }
    close $fh;


    

    print STDERR "-mapping candidate fusion transcripts to gene annotations.\n";
    ## header
    print join("\t", "#transcript", "num_alignments", "align_descr(s)", "[chim_annot_mapping]") . "\n";
    
    foreach my $target (keys %target_to_aligns) {
     

        # examine all the alignments for a given transcript.  A fusion transcript will
        # have at least 2 alignments reported in chimera-alignment mode.
        # For each ordered pair of alignment segments (A,B), map each 
        # to gene annotations and identify candidate fusion gene pairs.
        

        ## create spans
        my @span_ids = keys %{$target_to_aligns{$target}};
        if (scalar @span_ids < 2) { next; } # only want candidate chimeras
        
        my @spans;
        foreach my $span_id (@span_ids) {
            my $exon_hits_aref = $target_to_aligns{$target}->{$span_id};
            my $span_struct = &convert_to_span($exon_hits_aref);
            push (@spans, $span_struct);<commit_msg>restricting to high avg percent identity alignments for fusions (at 98% default)<commit_after>#!/usr/bin/env perl

use strict;
use warnings;

use Data::Dumper;
use Getopt::Long qw(:config no_ignore_case bundling pass_through);
use List::Util qw (min max);
use Set::IntervalTree;

my $MIN_PER_ID = 98;


my $usage = <<__EOUSAGE__;

#####################################################################
#
# Required:
#
#  --gmap_gff3 <string>        gmap gff3 alignment output
#
#  --annot_gtf <string>        transcript structures in gtf file format.
#
#  --min_per_id <float>        minimum percent identity (default: $MIN_PER_ID)
#
######################################################################


__EOUSAGE__

    ;


my $help_flag;
my $gmap_gff3_file;
my $annot_gtf_file;

&GetOptions ( 'h' => \$help_flag,
              'gmap_gff3=s' => \$gmap_gff3_file,
              'annot_gtf=s' => \$annot_gtf_file,
              'min_per_id=f' => \$MIN_PER_ID,
              
              );


if ($help_flag) {
    die $usage;
}

unless ($gmap_gff3_file && $annot_gtf_file) {
    die $usage;
}


my %genes;
my %interval_trees;


main: {

    my %chr_to_gene_coords;
    print STDERR "-parsing $annot_gtf_file\n";
    open (my $fh, $annot_gtf_file) or die "Error, cannot open file $annot_gtf_file";
    while (<$fh>) {
        chomp;
        unless (/\w/) { next; }
        if (/^\#/) { next; }
        s/^>//;
        my @x = split(/\t/);
        
        unless ($x[2] eq "exon") { next; }
        
        my $info = $x[8];
        $info =~ /gene_id \"([^\"]+)/ or die "Error, cannot extract gene_id from $_ [specifically from: $info]";
        my $gene_id = $1 or die "Error, no gene_id from $_";
        
        if ($info =~ /gene_name \"([^\"]+)/) {
            # use gene name instead
            $gene_id = $1;
        }
        
        $info =~ /transcript_id \"([^\"]+)/ or die "Error, cannot extract transcript_id from $_";
        my $transcript_id = $1 or die "Error, no trans id from $_";
        
        my ($lend, $rend) = ($x[3], $x[4]);
        my $chr = $x[0];
        my $orient = $x[6];
        
        
        push (@{$genes{$chr}->{$gene_id}->{$transcript_id}}, { 
            
            gene => $gene_id,
            transcript => $transcript_id,
            chr => $chr,
            lend => $lend,
            rend => $rend,
            orient => $orient,
              }
            );
        
        push (@{$chr_to_gene_coords{$chr}->{$gene_id}}, $lend, $rend);
        
    }
    close $fh;
    
    print STDERR "-building interval tree for fast searching of gene overlaps\n";
    ## Build interval trees
    foreach my $chr (keys %chr_to_gene_coords) {

        my $i_tree = $interval_trees{$chr} = Set::IntervalTree->new;
        
        foreach my $gene_id (keys %{$chr_to_gene_coords{$chr}}) {
            
            my @coords = sort {$a<=>$b} @{$chr_to_gene_coords{$chr}->{$gene_id}};
            my $lend = shift @coords;
            my $rend = pop @coords;

            $i_tree->insert($gene_id, $lend, $rend);
        }
    }
    
    
    my %target_to_aligns;
    print STDERR "-loading alignment data\n";
    open ($fh, $gmap_gff3_file) or die $!;
    while (<$fh>) {
        if (/^\#/) { next; }
        unless (/\w/) { next; }

        chomp;
        my ($chr, $filename, $type, $lend, $rend, $per_id, $orient, $dot, $info) = split(/\t/);
        
        my %info_hash;
        foreach my $keyval (split(/;/, $info)) {
            my ($key, $val) = split(/=/, $keyval);
            $info_hash{$key} = $val;
        }

        my $alignment_ID = $info_hash{ID};
        
        my ($target, $range_lend, $range_rend) = split(/\s+/, $info_hash{Target});

        push (@{$target_to_aligns{$target}->{$alignment_ID}}, { chr => $chr,
                                                                lend => $lend,
                                                                rend => $rend,
                                                                orient => $orient,
                                                                per_id => $per_id,
                                                                
                                                                range_lend => $range_lend,
                                                                range_rend => $range_rend,
                                                            });
        
    }
    close $fh;


    

    print STDERR "-mapping candidate fusion transcripts to gene annotations.\n";
    ## header
    print join("\t", "#transcript", "num_alignments", "align_descr(s)", "[chim_annot_mapping]") . "\n";
    
    foreach my $target (keys %target_to_aligns) {
     

        # examine all the alignments for a given transcript.  A fusion transcript will
        # have at least 2 alignments reported in chimera-alignment mode.
        # For each ordered pair of alignment segments (A,B), map each 
        # to gene annotations and identify candidate fusion gene pairs.
        

        ## create spans
        my @span_ids = keys %{$target_to_aligns{$target}};
        if (scalar @span_ids < 2) { next; } # only want candidate chimeras
        
        my @spans;
        foreach my $span_id (@span_ids) {
            my $exon_hits_aref = $target_to_aligns{$target}->{$span_id};
            my $span_struct = &convert_to_span($exon_hits_aref);
            if ($span_struct->{per_id} >= $MIN_PER_ID) {
                push (@spans, $span_struct);
            }
{'fraction_non_alphanumeric': 0.16272574939489853, 'fraction_numerical': 0.003723701359150996, 'mean_word_length': 1.596181730304495, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 67, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 32, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '660511', 'n_tokens_mistral': 3427, 'n_tokens_neox': 3244, 'n_words': 815}