Mercurial > repos > portiahollyoak > temp
comparison scripts/refine_breakpoint.in.pl @ 0:28d1a6f8143f draft
planemo upload for repository https://github.com/portiahollyoak/Tools commit 132bb96bba8e7aed66a102ed93b7744f36d10d37-dirty
author | portiahollyoak |
---|---|
date | Mon, 25 Apr 2016 13:08:56 -0400 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:28d1a6f8143f |
---|---|
1 #! /usr/bin/perl | |
2 | |
3 use strict; | |
4 | |
5 #my @files=<*.bp.bed.sfcp>; | |
6 my @files=<*.clipped.reads.aln>; | |
7 for my $file (@files) { | |
8 | |
9 my $title=$file; | |
10 $title =~ s/clipped.reads.aln/insertion.refined.bp/; | |
11 my $title2=$file; | |
12 $title2 =~ s/clipped.reads.aln/insertion.bp.bed/; | |
13 | |
14 open (input, "<$file") or die "Can't open $file since $!\n"; | |
15 open (input2, "<$title2") or die "Can't open $title2 since $!\n"; | |
16 open (output, ">>$title") or die "Can't open $title since $!\n"; | |
17 while (my $line=<input>) { | |
18 chomp($line); | |
19 my @a=split(/\t/, $line); | |
20 my @b=split(/\;/, $a[4]); | |
21 my $plusmax=""; | |
22 my $minusmax=""; | |
23 my $plus=0; | |
24 my $minus=0; | |
25 my $bp=""; | |
26 | |
27 my $line2=<input2>; | |
28 my @z=split(/\t/, $line2); | |
29 my $psup=abs($z[6]); | |
30 my $msup=abs($z[7]); | |
31 my $strand=$z[4]; | |
32 my $class=$z[5]; | |
33 | |
34 for my $element (@b) { | |
35 my @c=split(/\:/, $element); | |
36 chop($c[0]); | |
37 my @d=split(/\(/, $c[0]); | |
38 if (($d[1] eq "+") && ($c[1] > $plus)) { | |
39 $plusmax=$d[0]; | |
40 $plus=$c[1]; | |
41 } | |
42 elsif (($d[1] eq "-") && ($c[1] > $minus)) { | |
43 $minusmax=$d[0]; | |
44 $minus=$c[1]; | |
45 } | |
46 } | |
47 | |
48 if ($a[1] > 0) { | |
49 $a[1] += 15; | |
50 $a[2] -= 15; | |
51 print output "$a[0]\t$a[1]\t$a[2]\t$a[3]\t$strand\t$class\t"; | |
52 if (($minus >= 1)&&($plus >= 1)&&(abs($plusmax-$minusmax) <= 25)) { | |
53 print output "$plusmax\(\+\)\t$minusmax\(\-\)\t$plus\t$minus\t"; | |
54 } | |
55 elsif (($plus >= $minus)&&($plus >= 2)&&($plusmax >= $a[1])&&($plusmax <= $a[2])) { | |
56 print output "$plusmax\(\+\)\t$plusmax\(\-\)\t$plus\t0\t"; | |
57 } | |
58 elsif (($minus >= 2)&&($minusmax >= $a[1])&&($minusmax <= $a[2])) { | |
59 print output "$minusmax\(\+\)\t$minusmax\(\-\)\t0\t$minus\t"; | |
60 } | |
61 else { | |
62 my $mid=int(($a[1] + $a[2])/2); | |
63 print output "$mid\(\+\)\t$mid\(\-\)\t0\t0\t"; | |
64 } | |
65 print output "$psup\t$msup\n"; | |
66 } | |
67 } | |
68 close input; | |
69 close input2; | |
70 close output; | |
71 system("uniq $title > temp"); | |
72 system("mv temp $title"); | |
73 } |