-
Notifications
You must be signed in to change notification settings - Fork 10
/
dummyai-db-alone.pl
74 lines (66 loc) · 2.01 KB
/
dummyai-db-alone.pl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
=begin
PSI-Sigma: A splicing-detection method for short-read and long-read RNA-seq data
© Kuan-Ting Lin, 2018-2024
PSI-Sigma is free for non-commercial purposes by individuals at an academic or non-profit institution.
For commercial purposes, please contact tech transfer office of CSHL via narayan@cshl.edu
=end
=cut
#!/usr/bin/perl -w
use strict;
use Cwd qw(abs_path);
my ($gtf,$name) = @ARGV;
my $path = abs_path($0);
$path=~s/\/dummyai\-db\-alone\.pl//;
print "Path = $path\n";
my $noveljunctioncriteria = 10;
my %chr;
open(FILE,"$gtf") || die "Aborting.. Can't open $gtf : $!\n";
while(my $line=<FILE>){
chomp $line;
next if($line=~/^\#/);
my @array = split(/\t/,$line);
next if($array[2] ne "transcript");
my ($chr,$cat,$start,$end,$strand,$name) = ($array[0],$array[1],$array[3],$array[4],$array[6],$array[8]);
$chr = "chr" . $chr if($chr!~/chr/);
$chr{$chr}++;
}
close(OUT);
close(FILE);
my $starttime = time;
my $dbname = $name . ".db";
my $bedname = $name . ".bed";
my $chrs;
foreach my $chr(sort keys %chr){
$chrs .= "\t" . $chr;
}
$chrs=~s/\t//;
if(-e $dbname){
if(-z $dbname){
print "Regenerating $dbname...\n";
rundb($noveljunctioncriteria,$gtf,$chrs);
}
}else{
print "Generating $dbname...\n";
rundb($noveljunctioncriteria,$gtf,$chrs);
}
my $stoptime = time;
my $hours = sprintf("%.4f",(($stoptime-$starttime)/3600));
print "===Database spent $hours hours.===\n";
sub rundb{
my $noveljunctioncriteria = shift;
my $gtf = shift;
my $chrs = shift;
my @chromosomes = split(/\t/,$chrs);
foreach my $chromosome(@chromosomes){
next if($chromosome=~/chrGL/);
next if($chromosome=~/chrKI/);
my $commend = "perl " . $path . "/PSIsigma-db-v.1.0.pl $gtf " . $chromosome . " " . $noveljunctioncriteria;
#print "Doing... $commend\n";
print "Doing... $chromosome\n";
system("$commend");
}
system("cat chr*.db > $dbname");
system("cat chr*.bed > $bedname");
system("rm chr*.db");
system("rm chr*.bed");
}