-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathgffcmp_stats.pl
executable file
·91 lines (82 loc) · 2.3 KB
/
gffcmp_stats.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
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
#!/usr/bin/perl
use strict;
use Getopt::Std;
use FindBin;use lib $FindBin::Bin;
my $usage = q/Usage:
gffcmp_stats.pl [-H] [-s <set>] <gffcmp.stats>
Collect the relevant stats from the gffcompare *.stats output
in a tab-delimited table:
qfname rtnum rlocnum qtnum qmet qlocnum matchICnum matchTnum iSn iPr icSn icPr tSn tPr
Use -H option to print the header like the above
/;
umask 0002;
getopts('Hs:o:') || die($usage."\n");
my $outfile=$Getopt::Std::opt_o;
my $set=$Getopt::Std::opt_s;
if ($outfile) {
open(OUTF, '>'.$outfile) || die("Error creating output file $outfile\n");
select(OUTF);
}
# --
if ($Getopt::Std::opt_H) {
print "set\t" if $set;
print join("\t",
split(/\s+/, q/query rtnum rlocnum qtnum qmet qlocnum matchIC matchT iSn iPr icSn icPr tSn tPr/))."\n";
}
my ($qf, $rtnum, $rlocnum, $qtnum, $qmet, $qlocnum, $icMatch, $tMatch, $iSn, $iPr, $icSn, $icPr,
$tSn, $tPr);
while (<>) {
if (m/^#= Summary for dataset:\s+(\S+)/) {
my $f=$1;
if ($qf) {
print "$set\t" if $set;
print join("\t", ($qf, $rtnum, $rlocnum, $qtnum, $qmet, $qlocnum, $icMatch, $tMatch, $iSn, $iPr,
$icSn, $icPr, $tSn, $tPr))."\n";
}
$tMatch='';
#$qf=$f;
($qf)=($f=~m{([^/]+)$});
$qf=~s/\.\w+$//;
$qf=~s/\.\w+$//;
next;
}
if (m/^#\s+Query mRNAs\s*:\s*(\d+)\s+in\s+(\d+)\s+loci\s+\((\d+)\s+multi\-exon/) {
($qtnum, $qlocnum, $qmet)=($1, $2, $3);
next;
}
if (m/^#\s+Reference mRNAs\s*:\s*(\d+)\s+in\s+(\d+)\s+loci/) {
($rtnum, $rlocnum)=($1, $2);
next;
}
if (m/^\s+Matching intron chains\s*:\s*(\d+)/) {
$icMatch=$1;
next;
}
if (m/^\s+Matching transcripts\s*:\s*(\d+)/) {
$tMatch=$1;
next;
}
if (m/^\s*Intron level\s*:\s*([\d\.]+)[\s\|]+([\d\.]+)/) {
($iSn, $iPr)=($1, $2);
next;
}
if (m/^\s*Intron chain level\s*:\s*([\d\.]+)[\s\|]+([\d\.]+)/) {
($icSn, $icPr)=($1, $2);
next;
}
if (m/^\s*Transcript level\s*:\s*([\d\.]+)[\s\|]+([\d\.]+)/) {
($tSn, $tPr)=($1, $2);
next;
}
} #while <>
if ($qf && $tMatch) {
print "$set\t" if $set;
print join("\t", ($qf, $rtnum, $rlocnum, $qtnum, $qmet, $qlocnum, $icMatch, $tMatch,
$iSn, $iPr, $icSn, $icPr, $tSn, $tPr))."\n";
}
# --
if ($outfile) {
select(STDOUT);
close(OUTF);
}
#************ Subroutines **************