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 }