|
0
|
1 #!/usr/bin/perl
|
|
|
2 # This script takes a matrix as input and generates an R script to look for
|
|
|
3 # simularity between rows.
|
|
|
4 #
|
|
|
5 # 'tab2list' is a complement to list2tab. It reads a 2-dimensional table and
|
|
|
6 # produces a list in the form: row label - column label - value.
|
|
|
7 #
|
|
|
8 use Getopt::Long;
|
|
|
9
|
|
|
10 $keycols = 1;
|
|
|
11 $usage = 0;
|
|
|
12 @cis;
|
|
|
13 %hash;
|
|
|
14 open OUT, "> FISH/Fisher_pre.txt";
|
|
|
15 open OUT2, "> FISH/Fisher_pre_named.txt";
|
|
|
16 GetOptions("keycols=i" => \$keycols, "usage" => \$usage, "help"=>\$usage);
|
|
|
17
|
|
|
18 if ($usage) {
|
|
|
19 print STDERR "usage: tab2list < infile > outfile
|
|
|
20 Options:
|
|
|
21 --keycols=X ... How many cols from left are to be kept (default: 1)\n";
|
|
|
22 exit 0;
|
|
|
23 }
|
|
|
24
|
|
|
25 $keycols--;
|
|
|
26
|
|
|
27 $_ = <>;
|
|
|
28 chomp;
|
|
|
29 @headline = split /\t/, $_;
|
|
|
30 @anim = @headline;
|
|
|
31 shift(@anim);
|
|
|
32 while (<>) {
|
|
|
33 chomp;
|
|
|
34 @line = split /\t/, $_;
|
|
|
35 $head = "";
|
|
|
36 foreach my $i (0 .. $keycols) {
|
|
|
37 $head .= "$line[$i]\t";
|
|
|
38 }
|
|
|
39 chop $head;
|
|
|
40 push (@cis,$head);
|
|
|
41
|
|
|
42 foreach my $i ($keycols+1 .. $#line) {
|
|
|
43 # print "$head\t$headline[$i]\t$line[$i]\n";
|
|
|
44 $value=$head."|".$headline[$i];
|
|
|
45 $hash{$value}=$line[$i];
|
|
|
46
|
|
|
47 }
|
|
|
48 }
|
|
|
49 foreach (@cis){
|
|
|
50 #print "$_";
|
|
|
51 }
|
|
|
52
|
|
|
53 print "\n";
|
|
|
54 foreach (@anim){
|
|
|
55 #print "$_";
|
|
|
56 }
|
|
|
57
|
|
|
58 print "\n";
|
|
|
59 foreach (keys %hash) {
|
|
|
60 # print "The key $_ $hash{$_}\n";
|
|
|
61 }
|
|
|
62
|
|
|
63 %seen;
|
|
|
64 $first=0;
|
|
|
65 $second=0;
|
|
|
66 $third=0;
|
|
|
67 $fourth=0;
|
|
|
68
|
|
|
69 foreach $x(@cis){
|
|
|
70 foreach $y(@cis){
|
|
|
71 $first=0;
|
|
|
72 $second=0;
|
|
|
73 $third=0;
|
|
|
74 $fourth=0;
|
|
|
75 if ($x eq $y){
|
|
|
76 next;}
|
|
|
77 $z = $x."|".$y;
|
|
|
78 if (exists $seen{$z}){
|
|
|
79 next;}
|
|
|
80 $seen{$z}= '1';
|
|
|
81 $z2 = $y."|".$x;
|
|
|
82 $seen{$z2}= '1';
|
|
|
83 foreach $a(@anim) {
|
|
|
84 $b =$x."|".$a;
|
|
|
85 $c =$y."|".$a;
|
|
|
86 if (($hash{$b}==0) and ($hash{$c}==0)){
|
|
|
87 $fourth = $fourth + 1;
|
|
|
88 }
|
|
|
89 elsif (($hash{$b}==0) and ($hash{$c}==1)){
|
|
|
90 $third = $third + 1;
|
|
|
91 }
|
|
|
92 elsif (($hash{$b}==1) and ($hash{$c}==0)){
|
|
|
93 $second = $second + 1;
|
|
|
94 }
|
|
|
95 elsif (($hash{$b}==1) and ($hash{$c}==1)){
|
|
|
96 $first = $first + 1;
|
|
|
97 }
|
|
|
98 #print "$a|$x|$y\n";
|
|
|
99 #print "$hash{$b}|$hash{$c}\n";
|
|
|
100 }
|
|
|
101 print OUT "object__object\t$first\t$second\t$third\t$fourth\n";
|
|
|
102 print OUT2 "$x"."__"."$y\t$first\t$second\t$third\t$fourth\n";
|
|
|
103 }
|
|
|
104 }
|
|
|
105 close OUT;
|
|
|
106 close OUT2;
|
|
|
107 open SOURCE, "< FISH/Fisher_pre.txt";
|
|
|
108 open OUT, "> FISH/Fisher.R";
|
|
|
109 while (defined($line = <SOURCE>)) {
|
|
|
110 chomp $line;
|
|
|
111 @field= split(/\t/, $line);
|
|
|
112 print OUT "$field[0] <- matrix(c($field[1],$field[2],$field[3],$field[4]), nr=2, dimnames=list(c(\"cis\", \"nocis\"), c(\"pheno\", \"notpheno\")))\n$field[0]\nfisher.test($field[0]".",alternative = \"greater\")\n";
|
|
|
113 }
|
|
|
114 #system "module load R";
|
|
|
115 #system "R < FISH/Fisher.R --vanilla > FISH/fisher_dump.txt";
|
|
|
116 open (my $fisher_dump, ">", "FISH/fisher_dump.txt") || die "Unable to open FISH/fisher_dump.txt. $!\n";
|
|
|
117 open (my $out, "Rscript --vanilla FISH/Fisher.R | ");
|
|
|
118 while (<$out>) {
|
|
|
119 print $fisher_dump $_;
|
|
|
120 }
|
|
|
121 close($out);
|
|
|
122 close($fisher_dump);
|
|
|
123 close OUT;
|
|
|
124
|
|
|
125 open SOURCE, "< FISH/fisher_dump.txt";
|
|
|
126 open OUT, "> FISH/R_result.txt";
|
|
|
127 while (defined($line = <SOURCE>)) {
|
|
|
128 chomp ($line);
|
|
|
129 @field= split(/\t/, $line);
|
|
|
130 if ($field[0] =~ m/data/){
|
|
|
131 print OUT "$field[0]\t";
|
|
|
132 }
|
|
|
133 if ($field[0] =~ m/p-value/){
|
|
|
134 print OUT "$field[0]\n";
|
|
|
135 }
|
|
|
136 }
|
|
|
137 close OUT;
|
|
|
138 print "cooperation analyses complete";
|
|
|
139
|
