Mercurial > repos > mmaiensc > ember
comparison GALAXY_FILES/tools/EMBER/Expression_Pattern_Distance.pl @ 3:037c3edda16e
Uploaded
author | mmaiensc |
---|---|
date | Thu, 22 Mar 2012 13:49:52 -0400 |
parents | 003f802d4c7d |
children |
comparison
equal
deleted
inserted
replaced
2:1a84b8178b45 | 3:037c3edda16e |
---|---|
1 #!/usr/bin/perl | |
2 # computes difference distances between motif models | |
3 | |
4 if( $#ARGV != 1 ){ | |
5 print "Usage: ./Expression_Pattern_Distance.pl model1 model2\n"; | |
6 exit; | |
7 } | |
8 | |
9 @model1 = &read_model( $ARGV[0] ); | |
10 @model2 = &read_model( $ARGV[1] ); | |
11 if( $#model1 != $#model2 ){ | |
12 print "Error: model dimensions do not match\n"; | |
13 exit; | |
14 } | |
15 | |
16 | |
17 # | |
18 # second distance: sum of absolute difference, divided by absolute sum | |
19 # bounded between 0 and 1; equals 1 if all signs are different | |
20 # | |
21 $dist = 0; | |
22 $denom = 0; | |
23 for($i=0; $i< $classes; $i++){ | |
24 for($j=0; $j< $dims; $j++){ | |
25 if( $model1[$i][$j] ne "nan" && $model2[$i][$j] ne "nan" ){ | |
26 $dist+= abs($model1[$i][$j] - $model2[$i][$j]); | |
27 $denom+= abs($model1[$i][$j]) + abs($model2[$i][$j]); | |
28 } | |
29 } | |
30 } | |
31 $dist/=$denom; | |
32 | |
33 printf("%0.3f\n", $dist); | |
34 | |
35 | |
36 exit; | |
37 | |
38 | |
39 | |
40 | |
41 ################################ | |
42 | |
43 # function to read a motif model from file $_[0] | |
44 sub read_model{ | |
45 $dims = 0; | |
46 $classes = 0; | |
47 open(IN,"$_[0]") || die "Error: can't open file $_[0]\n"; | |
48 @modpars = (); | |
49 # burn the first two lines | |
50 $line = <IN>; | |
51 $line = <IN>; | |
52 while( $line = <IN> ){ | |
53 chomp($line); | |
54 @parts = split(' ',$line); | |
55 if($parts[0] !~ /#/ ){ | |
56 @parts2 = (); | |
57 for($i=1; $i<= $#parts; $i++){ | |
58 if( $parts[$i] ne "NA" ){ | |
59 push(@parts2, $parts[$i]); | |
60 } | |
61 else{ | |
62 push(@parts2, 0); | |
63 } | |
64 } | |
65 push(@modpars, [@parts2] ); | |
66 $classes++; | |
67 $dims = $#parts; | |
68 } | |
69 } | |
70 close(IN); | |
71 return @modpars; | |
72 } | |
73 |