-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathbedCluster.pl
executable file
·51 lines (51 loc) · 2.18 KB
/
bedCluster.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
#!/bin/env perl
use strict;
my %junction;
use 5.010;
while(<>){
chomp;
my @fields=split "\t";
my ($chr,$start,$end,$strand)=@fields[0,1,2,5];
my @blockSizes=split ",",$fields[10];
my $junctionStart=$start+$blockSizes[0];
my $junctionEnd=$end-$blockSizes[1];
if(exists $junction{$chr}{$strand}{$junctionStart}{$junctionEnd}){
$junction{$chr}{$strand}{$junctionStart}{$junctionEnd}{NO}++;
if ($start<$junction{$chr}{$strand}{$junctionStart}{$junctionEnd}{leftBoundary}){
$junction{$chr}{$strand}{$junctionStart}{$junctionEnd}{leftBoundary}=$start;
}
if ($start<$junction{$chr}{$strand}{$junctionStart}{$junctionEnd}{rightBoundary}){
$junction{$chr}{$strand}{$junctionStart}{$junctionEnd}{rightBoundary}=$end;
}
}else{
$junction{$chr}{$strand}{$junctionStart}{$junctionEnd}{NO}=1;
$junction{$chr}{$strand}{$junctionStart}{$junctionEnd}{leftBoundary}=$start;
$junction{$chr}{$strand}{$junctionStart}{$junctionEnd}{rightBoundary}=$end;
}
}
my $counter=1;
while( my ($chr,$iStrand)=each %junction){
while( my ($strand,$iJunctionStart)=each %$iStrand){
while( my ($junctionStart,$iJunctionEnd)=each %$iJunctionStart){
while( my ($junctionEnd,$junction)=each %{$iJunctionEnd}){
my $start=$junction->{leftBoundary};
my $end=$junction->{rightBoundary};
my $blockSizes=($junctionStart-$start).','.($end-$junctionEnd);
my $blockStarts='0,'.($junctionEnd-$start);
say join "\t",($chr,
$start,
$end,
"Junction".$counter++,#name
$junction->{NO},#read Cover No.
$strand,
$start,#thickStart
$end,#thickEnd
'255,0,0',#itemRgb
2,#blockCount
$blockSizes,
$blockStarts
);
}
}
}
}