PASSAnnotator.pl

#!/usr/bin/perl

use warnings;
use strict;

#use IO::Zlib;
#use Bio::DB::HTS::Tabix;
use FindBin;
use File::Spec;
use Benchmark qw(:all);

use lib "$FindBin::Bin/lib";
use Tests::BedParser;
use Tests::PassParser;


sub parseLineOld();
sub parseLineNew();
sub mysplit($$$);


my $queryBedFile = File::Spec->catfile($FindBin::Bin, 'data', 'refSeq_rn5_Annotated_May2017.bed.gz');
my $queryPassFile = File::Spec->catfile($FindBin::Bin, 'data', 'passData_ClustersRaw_May2017.bed.gz');



#my $obj = BedParser->new($queryBedFile);
#$obj->parse();


my $obj = PassParser->new($queryPassFile, $queryBedFile);
$obj->parse();








=head
MAIN:
{
    my $file_bed = shift;
    my $file_pass = shift;
    
    my $notail = 0;
    my $mito = 0;
    my $intergenic = 0;
    
    my $window_upstream = 20000;
    
    # index bed
    my $obj = Bio::DB::HTS::Tabix->new( filename => $file_bed);
    
    # parse pass
    my $fh = new IO::Zlib;
    $fh->open($file_pass, "rb");
    while(<$fh>)
    {
        # remove new line
        chomp($_);
        
        # current pass record
        my $pass = PassRecord->new($_);
        
        # count no tail
        if ($pass->ispoly() == 0)
        {
            $notail += $pass->{readsSumCoverage};
            next;
        }
        
        # count mito
        if ($pass->ischrom("chrM"))
        {
            $mito += $pass->{readsSumCoverage};
            next;
        }
        
       
        # parse annotations
        my $iter = $obj->query($pass->region($window_upstream));
        my $isintergenic = 1;
        while (my $line = $iter->next)
        {
            # current bed record
            my $bed = BedRecord->new($line);
            
            # skip different strand
            next if($bed->{strand} ne $pass->{strand});
            
            # find closest bed feature
            $bed->closest($pass->{tailsBestBase});
            #print "Tree Result: ",$result,"\n";
            
            
            $isintergenic = 0;
            
        }
        
        # update intergenic
        if ($isintergenic)
        {
            $intergenic += $pass->{readsSumCoverage};
            next;
        }
        
        
        last;
        
        
    }
    $fh->close();
    $obj->close();

    print "NoTail\t",$notail,"\n";
    print "Mito\t",$mito,"\n";
    print "Intergenic\t",$intergenic,"\n";
    
}
=cut