Newer
Older
BIOPZ-Gypas Foivos
committed
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
use Data::Dumper;
my %hash = ();
my %anno = ();
my $readsSum = 0;
open(SAM,$ARGV[0]);
#Annotation file included with beginning and end of every CDS
my $annoFile = (defined($ARGV[1])) ? 1 : 0;
if($annoFile){
open(AN,$ARGV[1]);
while( $line = <AN> ){
chomp $line;
my @F = split /\t/ , $line; #annotation line
my @a = ($F[1],$F[2]);
$anno{$F[0]} = \@a;
}
} else {
print STDERR "Annotation file not specified, using len(5UTR)=200\n";
}
#print Dumper(%anno);
#Parse SAM file and count reads
while( $line = <SAM> ){
chomp $line;
next if $line =~ /^@/; #discard headers/comments
my @F = split /\t/ , $line; #sam line
my $gene_name = $F[2];
if($F[1] == 0) { #only consider hits in plus strand
if($annoFile){
if(($F[3] >= $anno{$gene_name}->[0]-15) & ($F[3] <= $anno{$gene_name}->[1]-15)) { #close to start and stop
$hash{$gene_name}++;
}
} {
if($F[3] >= 186) { #consider that len(5UTR)==200
$hash{$gene_name}++;
}
}
}
}
#print Dumper(\%hash);
print "gene\tcounts\n";
for $key (keys %hash){
print $key,"\t",$hash{$key},"\n";
}