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
|
#!/usr/bin/env perl
use strict;
use warnings;
use lib ($ENV{EUK_MODULES});
use WigParser;
use Gene_obj;
my $usage = "usage: $0 clip.wig\n\n";
my $clip_wig = $ARGV[0] or die $usage;
main: {
my %mol_to_clips = &get_clip_pts($clip_wig);
foreach my $mol (keys %mol_to_clips) {
my @clip_pos = @{$mol_to_clips{$mol}};
foreach my $pos (@clip_pos) {
my ($end5, $end3) = ($pos-1, $pos+1);
my $gene_obj = new Gene_obj();
$gene_obj->populate_gene_object( {$end5=>$end3}, {$end5=>$end3});
$gene_obj->{asmbl_id} = $mol;
$gene_obj->{com_name} = "$mol-C:$pos";
print $gene_obj->to_BED_format();
}
}
exit(0);
}
####
sub get_clip_pts {
my ($jaccard_clips) = @_;
my %trans_to_clips;
my $trans_acc = "";
open (my $fh, $jaccard_clips) or die "Error, cannot open file $jaccard_clips";
while (<$fh>) {
chomp;
if (/^variableStep chrom=(\S+)/) {
$trans_acc = $1;
}
elsif (/^(\d+)/) {
my ($coord, $val) = split(/\t/);
if ($val) {
push (@{$trans_to_clips{$trans_acc}}, $coord);
}
}
}
close $fh;
return(%trans_to_clips);
}
|
