PassToMap.pl

#!/usr/bin/perl


use warnings;
use strict;
use Bio::DB::HTS::Tabix;

sub parseChromSizes($);

MAIN:
{
    my $file_annotation = shift;
    my $file_clusters = shift;
    my $file_chromSizes = shift;
    
    my $window_fp = 5000;
    my $window_tp = 20000;
    
    # parse chrom sizes
    my $chromSizes = parseChromSizes($file_chromSizes);
    
    # index annotation
    my $obj = Bio::DB::HTS::Tabix->new( filename => $file_annotation);
    
    # parse clusters
    my $mito = 0;
    my $notail = 0;
    my $intergenic = 0;
    open (my $fh, "gunzip -c $file_clusters|") or die $!;
    while(<$fh>)
    {
        # remove new line
        chomp($_);
        
        # split pass line
        my @line_pass = split("\t", $_, 15);
        
        # state no polyA tail
        if ($line_pass[11] == 0)
        {
            $notail += $line_pass[8];
            next;
        }
        
        # state mitochondrial chromosome
        if ($line_pass[0] eq "chrM")
        {
            $mito += $line_pass[8];
            next;
        }
        
        
        # state check closest
        my $regionLeft = ($line_pass[5] eq "+") ? ($line_pass[14] - $window_fp) : ($line_pass[14] - $window_tp);
        $regionLeft = 0 if($regionLeft < 0);
        
        my $regionRight = ($line_pass[5] eq "+") ? ($line_pass[14] + $window_tp) : ($line_pass[14] + $window_fp);
        $regionLeft = $chromSizes->{$line_pass[0]} if($regionLeft > $chromSizes->{$line_pass[0]});
        
        my $region = $line_pass[0] . ":" . $regionLeft . "-" . $regionRight;
        
        
        my $iter = $obj->query($region);
        my $flag = 0;
        while (my $bed = $iter->next)
        {
            # split reference bed line
            my @line_bed = split("\t", $bed, 12);
            
            # skip strand
            next if($line_bed[5] ne $line_pass[5]);
            
            # check 5'UTR_extended
            print "5pUTR_extended\n" if($line_pass[14] < $line_bed[1]);
            print "3pUTR_extended\n" if($line_pass[14] > $line_bed[2]);
            print "5pUTR\n" if(($line_pass[14] >= $line_bed[1]) && ($line_pass[14] < $line_bed[6]));
            print "3pUTR\n" if(($line_pass[14] <= $line_bed[2]) && ($line_pass[14] > $line_bed[7]));
            print "CDS/intron\n" if(($line_bed[6] <= $line_pass[14]) && ($line_pass[14] <= $line_bed[7]));
            
            
            $flag = 1;
            #print $bed,"\n";
            
        }
    
        if($flag == 1)
        {
            print $line_pass[3],"\t",$region,";",$line_pass[5],"\t",$line_pass[14],"\n";
            
            $intergenic += $line_pass[8];
            
            last;
            
        }
        
    }
    close($fh);
    $obj->close;
    
    print "NoTail: $notail\n";
    print "Mito: $mito\n";
    print "Intergenic: $intergenic\n";
    
}


sub parseChromSizes($)
{
    my $file_chromSizes = $_[0];
    my %hash_chromSizes = ();
    
    open(my $fh, "<", $file_chromSizes) or die $!;
    while(<$fh>)
    {
        chomp($_);
        my ($chrom, $span) = split("\t", $_, 2);
        $hash_chromSizes{$chrom} = $span;
    }
    close($fh);
    return \%hash_chromSizes;
}