-
Notifications
You must be signed in to change notification settings - Fork 32
/
plot_consensus_clustering.py
97 lines (75 loc) · 4.05 KB
/
plot_consensus_clustering.py
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
# -*- coding: utf-8 -*-
"""
Consensus clustering with modularity maximization
=================================================
This example demonstrates how to generate "consensus" clustering assignments
(Bassett et al., 2013) from potentially disparate clustering solutions. This is
particularly useful when the clustering algorithm used is stochastic / greedy
and returns different results when run on the same dataset multiple times
(e.g., Louvain modularity maximization).
"""
###############################################################################
# First let's grab some data to work with. We're going to download one session
# of parcellated functional MRI data from the MyConnectome dataset (Laumann et
# al., 2015). We'll pick session 73 (though any session would do):
from urllib.request import urlopen
import numpy as np
url = 'https://s3.amazonaws.com/openneuro/ds000031/ds000031_R1.0.2/' \
'uncompressed/derivatives/sub-01/ses-073/' \
'sub-01_ses-073_task-rest_run-001_parcel-timeseries.txt'
data = np.loadtxt(urlopen(url).readlines())
print(data.shape)
###############################################################################
# The data has 518 timepoints for each of 630 regions of interest (ROIs)
# defined by the original authors. We'll use this to construct an ROI x ROI
# correlation matrix and then plot it:
import matplotlib.pyplot as plt
corr = np.corrcoef(data.T)
fig, ax = plt.subplots(1, 1)
coll = ax.imshow(corr, vmin=-1, vmax=1, cmap='viridis')
ax.set(xticklabels=[], yticklabels=[])
fig.colorbar(coll)
###############################################################################
# We can see some structure in the data, but we want to define communities or
# networks (groups of ROIs that are more correlated with one another than ROIs
# in other groups). To do that we'll use the Louvain algorithm from the `bctpy`
# toolbox.
#
# Unfortunately the defaults for the Louvain algorithm cannot handle negative
# data, so we will make a copy of our correlation matrix and zero out all the
# negative correlations:
import bct
nonegative = corr.copy()
nonegative[corr < 0] = 0
ci, Q = bct.community_louvain(nonegative, gamma=1.5)
num_ci = len(np.unique(ci))
print('{} clusters detected with a modularity of {:.2f}.'.format(num_ci, Q))
###############################################################################
# We'll take a peek at how the correlation matrix looks when sorted by these
# communities. We can use the :func:`~.plotting.plot_mod_heatmap` function,
# which is a wrapper around :func:`plt.imshow()`, to do this easily:
from netneurotools import plotting
plotting.plot_mod_heatmap(corr, ci, vmin=-1, vmax=1, cmap='viridis')
###############################################################################
# The Louvain algorithm is greedy so different instantiations will return
# different community assignments. We can run the algorithm ~100 times to see
# this discrepancy:
ci = [bct.community_louvain(nonegative, gamma=1.5)[0] for n in range(100)]
fig, ax = plt.subplots(1, 1, figsize=(6.4, 2))
ax.imshow(ci, cmap='Set1')
ax.set(ylabel='Assignments', xlabel='ROIs', xticklabels=[], yticklabels=[])
###############################################################################
# We'll provide these different assignments to our consensus-finding algorithm
# which will generate one final community assignment vector:
from netneurotools import cluster
consensus = cluster.find_consensus(np.column_stack(ci), seed=1234)
plotting.plot_mod_heatmap(corr, consensus, cmap='viridis')
###############################################################################
# The :func:`netneurotools.modularity.consensus_modularity` function provides a
# wrapper for this process of generating multiple community assignmenta via the
# Louvain algorithm and finding a consensus. It also generates and returns some
# metrics for assessing the quality of the community assignments.
#
# Nevertheless, the :func:`~.cluster.find_consensus` function is useful for
# generating a consensus clustering solution from the results of _any_
# clustering algorithm (not just Louvain).