comparison GALAXY_FILES/tools/EMBER/Expression_Pattern_Distance.pl @ 1:e62b2ba92070 default tip

Uploaded

1:e62b2ba92070

#!/usr/bin/perl
# computes difference distances between motif models

if( $#ARGV != 1 ){
	print "Usage: ./Expression_Pattern_Distance.pl model1 model2\n";
	exit;
}

@model1 = &read_model( $ARGV[0] );
@model2 = &read_model( $ARGV[1] );
if( $#model1 != $#model2 ){
	print "Error: model dimensions do not match\n";
	exit;
}


#
# second distance: sum of absolute difference, divided by absolute sum
# bounded between 0 and 1; equals 1 if all signs are different
#
$dist = 0;
$denom = 0;
for($i=0; $i< $classes; $i++){
	for($j=0; $j< $dims; $j++){
		if( $model1[$i][$j] ne "nan" && $model2[$i][$j] ne "nan" ){
			$dist+= abs($model1[$i][$j] - $model2[$i][$j]);
			$denom+= abs($model1[$i][$j]) + abs($model2[$i][$j]);
		}
	}
}
$dist/=$denom;

printf("%0.3f\n", $dist);


exit;




################################

# function to read a motif model from file $_[0]
sub read_model{
	$dims = 0;
	$classes = 0;
	open(IN,"$_[0]") || die "Error: can't open file $_[0]\n";
	@modpars = ();
	# burn the first two lines
	$line = <IN>;
	$line = <IN>;
	while( $line = <IN> ){
		chomp($line);
		@parts = split(' ',$line);
		if($parts[0] !~ /#/ ){
			@parts2 = ();
			for($i=1; $i<= $#parts; $i++){
				if( $parts[$i] ne "NA" ){
					push(@parts2, $parts[$i]);
				}
				else{
					push(@parts2, 0);
				}
			}
			push(@modpars, [@parts2] );
			$classes++;
			$dims = $#parts;
		}
	}
	close(IN);
	return @modpars;
}