Mercurial > repos > portiahollyoak > temp

comparison scripts/make.bp.bed.pl @ 0:28d1a6f8143f draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 9672fe07a232
comparison
equal deleted inserted replaced
-1:000000000000 0:28d1a6f8143f
1 #! /usr/bin/perl
2
3 use strict;
4
5 my @sample=();
6 open (in, "<$ARGV[0]") or die "Can't open $ARGV[0] since $!\n";
7 my $line=<in>;
8 close in;
9
10 my %chrs=();
11 my @a=split(/\t/, $line);
12 for my $i (0..$#a) {
13 if ($a[$i] =~ /_class$/) {
14 my $name=$a[$i];
15 $name =~ s/_class//;
16 my $j=$i+1;
17 my $k=$i+2;
18 my $l=$i+3;
19 system("cut -f7,4,6,$j,$k,$l $ARGV[0] > temp");
20 open (input, "<temp") or die "Can't open temp since $!\n";
21 open (output, ">>$name.insertion.bp.bed") or die "Can't open $name.insertion.bp.bed since $!\n";
22 my $header=<input>;
23 while (my $line=<input>) {
24 chomp($line);
25 my @b=split(/\t/, $line);
26 if (($b[4] ne "0")||($b[5] ne "0")) {
27 my @c=split(/\:/, $b[2]);
28 my @d=split(/\./, $c[1]);
29 if ($d[0] > $d[1]) {
30 my $temp=$d[0];
31 $d[0]=$d[1];
32 $d[1]=$temp;
33 }
34 my $lower=$d[0];
35 my $upper=$d[1];
36 if (($lower >= 0) && ($upper >= 0)) {
37 print output "$c[0]\t$lower\t$upper\t$b[0]\t$b[1]\t$b[3]\t$b[4]\t$b[5]\n";
38 }
39 $chrs{$c[0]}=1;
40 }
41 }
42 close input;
43 close output;
44 system("rm temp");
45
46 if ($ARGV[1] ne "") {
47 open (input, "<$name.insertion.bp.bed") or die "Can't open $name.insertion.bp.bed since $!\n";
48 open (output, ">tmp") or die "Can't tmp since $!\n";
49 while (my $line=<input>) {
50 chomp($line);
51 my @a=split(/\t/, $line);
52 if (($a[0] =~ /^\d{1,2}$/) || ($a[0] eq "X") || ($a[0] eq "Y")) {$a[0]="chr$a[0]";}
53 my $strand="+";
54 if ($a[4] eq "antisense") {$strand="-";}
55 print output "$a[0]\t$a[1]\t$a[2]\t$a[3]\t\.\t$strand\t$a[5]\t$a[6]\t$a[7]\n";
56 }
57 close input;
58 close output;
59
60 system("bedtools intersect -a tmp -b $ARGV[1] -f 0.1 -wo -s > tmp1");
61 if ($ARGV[2] eq "") {
62 system("awk -F \"\\t\" '{OFS=\"\\t\"; if ((\$4==\$13)&&(\$6==\$15)) print \$1,\$2,\$3,\$4,\$5,\$6}' tmp1 > tmp2");
63 }
64 else {
65 my %family=();
66 open (input, "<$ARGV[2]") or die "Can't open $ARGV[2] since $!\n";
67 while (my $line=<input>) {
68 chomp($line);
69 my @a=split(/\t/, $line);
70 $family{$a[0]}=$a[1];
71 }
72 close input;
73
74 open (input, "<tmp1") or die "Can't open tmp1 since $!\n";
75 open (output, ">>tmp2") or die "Can't open tmp2 since $!\n";
76 while (my $line=<input>) {
77 chomp($line);
78 my @a=split(/\t/, $line);
79 if (($family{$a[3]} eq $family{$a[12]}) && ($a[5] eq $a[14])) {
80 print output "$a[0]\t$a[1]\t$a[2]\t$a[3]\t$a[4]\t$a[5]\n";
81 }
82 }
83 close input;
84 close output;
85 }
86
87 if (-s "tmp2") {
88 system("bedtools subtract -a tmp -b tmp2 -f 1.0 > tmp3");
89 open (input, "<tmp3") or die "Can't open tmp3 since $!\n";
90 open (output, ">$name.insertion.bp.bed") or die "Can't open $name.insertion.bp.bed since $!\n";
91 while (my $line=<input>) {
92 chomp($line);
93 my @a=split(/\t/, $line);
94 my $direction="sense";
95 if ($a[5] eq "-") {$direction="antisense";}
96 my $chr_num=$a[0];
97 $chr_num =~ s/chr//;
98 if (($chrs{$a[0]} == 1) && (! defined $chrs{$chr_num})) {$chr_num=$a[0];}
99 print output "$chr_num\t$a[1]\t$a[2]\t$a[3]\t$direction\t$a[6]\t$a[7]\t$a[8]\n";
100 }
101 close input;
102 close output;
103 }
104 }
105
106 system("rm tmp*");
107
108 }
109 }
110