Mercurial > repos > mmaiensc > ember

comparison GALAXY_FILES/tools/EMBER/Expression_Pattern_Distance.pl @ 3:037c3edda16e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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