-
Notifications
You must be signed in to change notification settings - Fork 0
/
clustering.rb
115 lines (98 loc) · 2.43 KB
/
clustering.rb
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
#!/usr/bin/env ruby
require_relative 'disease.rb'
# do not merge cluster that are less similar than this value
MERGE_THRESHOLD=0.9
class Cluster < Array
def initialize(disease)
super([disease])
end
def average_similarity(other)
sum = 0.0
count = 0
self.each { |d1|
other.each { |d2|
tmp = d1.sorensen_dice(d2)
sum += tmp
count += 1
}
}
sum / count
end
end
diseases = {}
clusters = []
File.open('Disease-Gene-Net.txt') { |file|
file.readline # skip header line
file.each_line { |line|
if line =~ /^([\w\d]+)\s+(.+)/
gene = $1
disease = $2
# initialize gene list with [] if needed
diseases[disease] ||= Disease.new(disease)
# and append this gene to the list
diseases[disease] << gene
end
}
}
print "#{diseases.size} diseases\n"
# build cluster list by first merging diseases with
# similarity 1
clusters = []
diseases.values.each { |d|
new_c = Cluster.new(d)
c = clusters.find { |tmp| tmp.average_similarity(new_c) == 1.0 }
if c
c.concat(new_c)
else
clusters << new_c
end
}
print "#{clusters.size} clusters\n"
diseases.values.map { |d| Cluster.new(d) }
# then perform hierarchical clustering
# keep a cache of the pairwise similarity values
$cache = Array.new(clusters.size) { Array.new(clusters.size) { nil } }
# function to find the pair of closest clusters
def find_closest(list)
max_sim = 0.0
best_pair = nil
0.upto(list.size-2) { |i|
(i+1).upto(list.size-1) { |j|
d = if $cache[i][j].nil?
$cache[i][j] = list[i].average_similarity(list[j])
else
$cache[i][j]
end
if d > max_sim && d > MERGE_THRESHOLD
max_sim = d
best_pair = [i, j]
end
}
}
if max_sim > MERGE_THRESHOLD
print "merging #{best_pair.join('-')}, with similarity #{'%.2f' % (100*max_sim)}\n"
$stdout.flush
best_pair
else
nil
end
end
i = 0
while (pair = find_closest(clusters)) && i < 1000
a = pair.min
b = pair.max
# invalidate cache for a, and the values of other clusters towards a
$cache[a].size.times { |j| $cache[a][j] = nil }
clusters.size.times { |j| $cache[j][a] = nil }
clusters[a].concat(clusters[b])
clusters.delete_at(b)
$cache.delete_at(b)
i += 1
end
top = clusters.sort { |a, b| b.size <=> a.size }[0, 10]
top.each_with_index { |c, idx|
print "Cluster #{idx}:\n"
c.each { |disease|
print " #{disease}\n"
}
}