スピアマンの順位相関係数を perl で計算

September 30, 2006

大枚はたいて導入している matlab には、 スピアマンの順位相関係数を計算する関数がない。 探したらあるのかもしれないがめんどうなので、作ることにした。 ちなみに GNU octave にはあるので、検証に使える。 最初に見た Wikipedia/JP の説明ではどうもよくわからないので、 Numerical Recepes in Fortran 77 を参考に Perl でスクリプトを書いた。 入力ファイルは二つ。どっちも形式は同じで、１カラム目に名前、 ２カラム目にその点数または順位が書いてあってタブで区切られているテキスト形式。 たとえば遺伝子とその発現量などだったら、一つの遺伝子が１行。こんな感じ。
dna_a 132937
clock 634
per1 3283
...
 二つのファイルは行数が同じで、同じ名前のものが含まれていなければならない。 必ずしも整列されていなくてもよい。点数や順位は同位（タイ tie）があってもよい。 順位がファイルに書いてある場合、同位は同じ順位になっていればそれでよい。 スクリプトは下。

用意するもの

やってみよう

カレントディレクトリの s1.txt と s2.txt という二つをのファイルからデータを読む。 この二つのデータファイル中には UTF-8 で日本語が書き込んであるので、 それが文字化けしたりする場合は、こっちの s1e.txt と s2e.txt がよい。 そしてダウンロードしたスクリプト中のファイル名を書き直す。

スクリプトの中身

#!/usr/bin/env perl

# spearman.perl - calculating spearman's rank correlation value
#                 from s1.txt and s2.txt
#
# September 29, 2006 TOMINAGA Daisuke
#
# Input: {s1,s2}.txt
# Output: stdout
#         tab separated text. the second column is the calculated result.
# ref. numerical receipes in fortran 77 second ed.

open(S1, "s1.txt"); @S1 = ; close(S1); $s1 = @S1;
open(S2, "s2.txt"); @S2 = ; close(S2); $s2 = @S2;
if ($s1 != $s2) { print "contens of s1.txt and s2.txt differ.\n"; exit; }

for ($i=0; $i<@S1; $i++) { @e = split(/\s+/, $S1[$i]); $ps1{$e[0]} = $e[1]; }
for ($i=0; $i<@S1; $i++) { @e = split(/\s+/, $S2[$i]); $ps2{$e[0]} = $e[1]; }

@ss1 = sort {$ps1{$a} <=> $ps1{$b}} keys %ps1;
@ss2 = sort {$ps2{$a} <=> $ps2{$b}} keys %ps2;

$prep = $m = 0;
for ($i = 0; $i < @ss1; $i++) {
	if ($ps1{$ss1[$i]} == $prep) { $m++; }
	else {
		$rs1{$ss1[$i]} = $i;
		$prep = $ps1{$ss1[$i]};
		if ($m != 0) {
			$r = $i - 1 - $m / 2;
			for ($j = $i - $m - 1; $j < $i; $j++) { $rs1{$ss1[$j]} = $r; }
			@ts1 = (@ts1, $m); $m = 0;
		}
	}
}
$i = @ss1;
if ($m != 0) {
	$r = $i - 1 - $m / 2;
	@ts1 = (@ts1, $m);
	for ($j = $i - $m - 1; $j < $i; $j++) { $rs1{$ss1[$j]} = $r; }
}
# same process as the above, s1
$prep = $m = 0;
for ($i = 0; $i < @ss2; $i++) {
	if ($ps2{$ss2[$i]} == $prep) { $m++; }
	else {
		$rs2{$ss2[$i]} = $i;
		$prep = $ps2{$ss2[$i]};
		if ($m != 0) {
			$r = $i - 1 - $m / 2;
			for ($j = $i - $m - 1; $j < $i; $j++) { $rs2{$ss2[$j]} = $r; }
			@ts2 = (@ts2, $m); $m = 0;
		}
	}
}
if ($m != 0) {
	$r = $i - 1 - $m / 2;
	@ts2 = (@ts2, $m); 
	for ($j = $i - $m - 1; $j < $i; $j++) { $rs2{$ss2[$j]} = $r; }
}

$D = 0;
foreach $key (keys %rs2) { $D += ($rs1{$key} - $rs2{$key})**2; }

$N = @ss1; $Ni = $N**3 - $N;
$tx = $ty = 0;
for ($i = 0; $i < @ts1; $i++) { $tx += ($ts1[$i]**3 - $ts1[$i]); }
for ($i = 0; $i < @ts2; $i++) { $ty += ($ts2[$i]**3 - $ts2[$i]); }
$rho = -(2*(6*$D-$Ni)+$tx+$ty)/(2*sqrt(($Ni-$tx)*($Ni-$ty)));
print "spearman's rank correlation value:\t$rho\n";