\chapter{Graph clustering}

\section{Introduction}

Abruptely, graph clustering consists in grouping the nodes of the networks into different
classes or clusters. The groupment of the nodes can be done according to various different criteria, i.e., nodes of 
the same color, nodes of the same type, etc. Commonly, nodes are grouped according to the fact they present a relatively high number of connections between them compared to the number of connections with the other nodes 
composing the network. In the following, we will only consider clustering methods aiming at retrieving highly interconnected groups of nodes in a network.

In bioinformatics, a lot of clustering approaches have already been applied to various types of network, e.g. protein-protein interaction network (see among others \cite{Spirin2003,King2004,Krogan2006}), metabolic graphs \cite{Gagneur2003}, biological sequences (\cite{Enright2002,Leplae2004}), etc. 

Clustering of protein interaction network may be of valuable help in order to retrieve in a large graphs real biological complexes in the cell. Moreover, if in the detected complexes some of proteins are of unknown function but the rest of the proteins present all present a similar function, this may give insights in the function of the unknown protein.

In the following, we will apply different graph based clustering approaches on the yeast protein - protein interaction network published by Gavin \textit{et al} \cite{Gavin2006} and obtained by multiple co-immunoprecipitation experiments with each yeast protein used as bait followed by a mass spectrometry procedure to identify all the proteins that precipitated with the baits.

The clustering algorithms we will apply are the \program{MCL} \cite{VanDongenPHD2000,Enright2002} and \program{RNSC} \cite{King2004}. Hereafter, follows a short description of both clustering algorithms copied from \cite{Brohee2006}.
 
The Markov Cluster algorithm (\program{MCL}) simulates a flow on the graph by
calculating successive powers of the associated adjacency matrix. At
each iteration, an \textit{inflation step} is applied to enhance the
contrast between regions of strong or weak flow in the graph. The
process converges towards a partition of the graph, with a set of 
high-flow regions (the clusters) separated by boundaries with no 
flow. The value of the \textit{inflation parameter} strongly 
influences the number of clusters.

The second algorithm, Restricted Neighborhood Search Clustering (\program{RNSC}), 
is a cost-based local search
algorithm that explores the solution space to minimize a cost
function, calculated according to the numbers of intra-cluster and
inter-cluster edges. Starting from an initial random solution, \program{RNSC}
iteratively moves a vertex from one cluster to another if this
move reduces the general cost. When a (user-specified) number of moves
has been reached without decreasing the cost function, the program
ends up.

In order to dispose of a negative control, we advice the reader to read the next chapter about graph randomization and alteration. 

\section{Network clustering comparison}
\subsection{Study case}

In this demonstration, we will compare the performances of two graph based clustering algorithms
\program{MCL} and \program{RNSC}. First, we will apply them to the protein - protein interaction described by Gavin \textit{et al} \cite{Gavin2006}, secondly we will compare the resulting clusters to the complexes annotated for the yeast in the MIPS database \cite{Mewes2007}.

Note that as the 
interaction network and the MIPS complexes are different dataset (i.e. different proteins), the performances of the algorithm
will be rather low.

To run this tutorial on the command line, you need to have both \program{RNSC} and \program{MCL} installed on your computer. You can find the MCL source code on \textit{http://micans.org/mcl/} and \program{RNSC} on \textit{http://rsat.ulb.ac.be/~rsat/rnsc/rnsc\_rewritten\_compiled32.zip}. 

\subsection{Protocol for the web server}

\subsubsection{Dataset download}
Go on the NeAT demo dataset web page (\textit{http://rsat.ulb.ac.be/rsat/data/neat\_tuto\_data/}) and download the MIPS complexes (\file{mips\_complexes.tab})..

\subsubsection{Network clustering with MCL}

\begin{enumerate}

\item In the \neat left menu, select the command \program{graph-based clustering (MCL)}. 

  In the right panel, you should now see a form entitled
  ``MCL''.

\item Click on the button \option{DEMO}. 

  The form is now filled with the Gavin co-immunoprecipitation protein interaction network graph in the tab-delimited format, and the parameters have been set up to their appropriate value for the demonstration, i.e., the inflation value (the MCL main parameter) is set to 1.8, the optimal value for \program{MCL} protein interaction network clustering \cite{Brohee2006}. 
  
  The inflation acts mainly on the number of clusters resulting from the clustering, i.e., by increasing the inflation, you will obtain a larger number of smaller clusters.
  
  Note that MCL accepts weighted networks (which is not the case here), a higher weight on an edge will reinforce the strength of the link between two nodes.

\item Click on the button \option{GO}. 

  The computation should take less than one minute. 
  On one hand, the result page displays a link to the result file and on the other hand a graphic showing the size distribution of the obtained complexes is also available. These will be discussed in the \textit{Interpretation of the results} section. 
  
\item Save the resulting file under the name \file{gavin\_2006\_mcl\_inf\_1.8\_clusters.tab} by right clicking on the resulting file and choosing \textit{Save as ...}.
\end{enumerate}
\subsubsection{Network clustering with RNSC}

\begin{enumerate}

\item In the \neat left menu, select the command \program{graph-based clustering (RNSC)}. 

  In the right panel, you should now see a form entitled
  ``RNSC''.

\item Click on the button \option{DEMO}. 

  The form is now filled with the Gavin co-immunoprecipitation protein interaction network graph in the tab-delimited format, and the parameters have been set up to their appropriate value for the demonstration, i.e., the numerous RNSC parameters are set to the optimal values for \program{RNSC} protein interaction network clustering determined in \cite{Brohee2006}. However, in this study, we found that the \program{RNSC} performances were not strongly affected by the parameters values.
  
  Note that, unlike MCL, RNSC does not accept weighted networks.

\item Click on the button \option{GO}. 

  The computation should take less than one minute. 
  On one hand, the result page displays a link to the result file and on the other hand a graphic showing the size distribution of the obtained complexes is also avaible. These will be discussed in the \textit{Interpretation of the results} section. 
  
\item Save the resulting file under the name \file{gavin\_2006\_rnsc\_clusters.tab} by right clicking on the resulting file and choosing \textit{Save as ...}.

\end{enumerate}
\subsubsection{Clustering quality assessment}
In the following, we will only describe the procedure to quantify the performances of 
the clustering algorithms by comparing the \program{MCL} obtained clusters to the 
complexes annotated in the MIPS database.
You will thus have to redo this whole section with the \program{RNSC} clustering results.

\begin{enumerate}

\item In the \neat left menu, select the command \program{Compare classes/clusters}. 

  In the right panel, you should now see a form entitled
  ``compare-classes''. This program will build a contigency table, i.e., a table where each line represents
  the annotated complexes and each column the clusters of highly connected proteins. This matrix will then be used to 
  compute quality statistics.
  
\item In the ``Upload query classes from file'' menu, select the file \file{gavin\_2006\_mcl\_inf\_1.8\_clusters.tab} we just computed.
\item In the ``Upload reference classes from file'' menu, select the file \file{mips\_complexes.tab} we just downloaded.
\item Select ``Matrix file'' as output format
\item Click on the button \option{GO}. 
\item The contigency table (see the resulting links as text or HTML file) can now be used in the next process by clicking on the button
\option{contingency table statistics}.

  In the right panel, you should now see a form entitled
  ``contingency-stats''. This program will compute the statistics described in \cite{Brohee2006}, namely the $PPV$, the \textit{sensitivity} and the \textit{Separation} statistics in order 
  to estimate the quality of a clustering results to predict the complexes annotated in the MIPS. 

\item Click on the button \option{GO}.

The resulting statistics will be described in the following section \textit{Interpretation of the results}, save them under the name \file{gavin\_2006\_mcl\_inf\_1.8\_vs\_mips\_stats.tab}. 

Re-do the whole procedure with the file obtained with \program{RNSC} and save the contingency-stats output under the name \file{gavin\_2006\_rnsc\_vs\_mips\_stats.tab}. 

\end{enumerate}

\subsection{Protocol for the command-line tools}

If you have installed a stand-alone version of the NeAT distribution,
you can also use all the programs on the 
command-line. This requires to be familiar with the Unix shell
interface. If you do not have the stand-alone tools, you can skip this
section and read the next section (\textit{Interpretation of the results}).

The explanation of the parameters used for \program{RNSC} and \program{MCL} in this approach are described in the \textit{Web server} section of this chapter.

We will now describe the use of \program{RNSC}, \program{MCL}, \program{compare-classes}, \program{convert-classes}, \program{convert-graph}  and \program{contigency-stats} as command line tools. 
As a preliminary step, go on the NeAT demo dataset web page (\textit{http://rsat.ulb.ac.be/rsat/data/neat\_tuto\_data/}) and download the MIPS complexes (\file{mips\_complexes\_names.tab}) and the Gavin interaction dataset (\file{gavin\_2006\_names.tab}).



\subsubsection{Network clustering with MCL}
\begin{enumerate}
\item The first step consist in applying \program{MCL} on the co-immunoprecipitation dataset. 
To this, go into the directory where you downloaded the file \file{uetz\_2001.tab} and use this command. 
	
	{\color{Blue} \begin{footnotesize} 
		\begin{verbatim}
	mcl gavin_2006_names.tab -I 1.8 --abc -o gavin_2006_mcl_inf_1.8_clusters.mcl
		\end{verbatim} \end{footnotesize}
	}		
The file \file{gavin\_2006\_mcl\_inf\_1.8\_clusters.mcl} is created and contains the clusters of highly connected node in the interaction dataset. However, this file is formatted in the \program{MCL} format that is not usable by the NeAT / RSAT tools. We will thus use the program \program{convert-classes} to convert this file in a tab delimited format with
the following command.


	{\color{Blue} \begin{footnotesize} 
		\begin{verbatim}
	convert-classes -i gavin_2006_mcl_inf_1.8_clusters.mcl  \
			-o gavin_2006_mcl_inf_1.8_clusters.tab \
			-from mcl -to tab
		\end{verbatim} \end{footnotesize}
	}
The resulting file is a two column file containing for each node (first column) the cluster to which it belongs (second column).
\end{enumerate}
\subsubsection{Network clustering with RNSC}
\begin{enumerate}
\item The first step will consist in converting the tab delimited format in which the protein interaction dataset
is encoded into a format readable by the RNSC clustering algorithm. To this, we will use the convert-graph programwith
the following command.

	{\color{Blue} \begin{footnotesize} 
		\begin{verbatim}
	convert-graph -from tab -to rnsc \
	              -i gavin_2006_names.tab -o gavin_2006_rnsc
		\end{verbatim} \end{footnotesize}
	}

Two files are created, \file{gavin\_2006\_rnsc.rnsc} and \file{gavin\_2006\_rnsc\_node\_names.rnsc}. The first one contains the graph in itself, under the format of an adjacency list. However, each node is identified by a number. The protein names corresponding to the nodes identifiers are encoded in the second file (two column tab delimited file).

\item We can now apply \program{RNSC} on the network with the following command.

	{\color{Blue} \begin{footnotesize} 
		\begin{verbatim}
	rnsc -g gavin_2006_rnsc.rnsc -t 50 -T 1 -n 15 -N 15 -e 3 -D 50  \
	     -d 3 -o gavin_2006_rnsc_clusters.rnsc
		\end{verbatim} \end{footnotesize}
	}
		
		
The file \file{gavin\_2006\_rnsc\_clusters.rnsc} is created and contains the clusters of highly connected node in the interaction dataset. However, this file is formatted in the \program{RNSC} format that is not usable by the \textit{NeAT / RSAT} tools. We will thus use the program \program{convert-classes} to convert this file in a tab delimited format with the following command.

	{\color{Blue} \begin{footnotesize} 
		\begin{verbatim}
	convert-classes -i gavin_2006_rnsc_clusters.rnsc  \
                        -o gavin_2006_rnsc_clusters.tab  \
                        -from rnsc -to tab \
                        -names gavin_2006_rnsc_node_names.rnsc
		\end{verbatim} \end{footnotesize}
	}
The resulting file is a two column file containing for each node (first column) the cluster to which it belongs (second column).

\end{enumerate}

\subsubsection{Assessing clustering quality}

In this section, we will describe how to build a contingency table by comparing the clusters extracted from the networks by \program{MCL} and \program{RNSC} to annotated complexes and the way to compute statistics on this contingency-table.

We will only describe the procedure for the \program{MCL} results. You should redo this section for the \program{RNSC} clustering results. 


\begin{enumerate} 

\item The program \program{compare-classes} can build (among other things) a contingency table, i.e., a table where each line represents the annotated complexes and each column the clusters of highly connected proteins. This table will then be used to compute quality statistics.

	{\color{Blue} \begin{footnotesize} 
		\begin{verbatim}
	compare-classes -q gavin_2006_inf_1.8.tab \
                        -r mips_complexes_names.tab -matrix QR \
                        -o gavin_2006_inf_1.8_cc_complexes_matrix.tab
		\end{verbatim} \end{footnotesize}
	}

The file \file{gavin\_2006\_inf\_1.8\_cc\_complexes\_matrix.tab} now contains a contigency table in a tab delimited format.

\item We can now study the quality of the clustering with the \program{contingency-stats} tool that was used in \cite{Brohee2006} to computed standard evaluation statistics like the $PPV$, sensitivity and the accuracy that will be precisely described in the following.

	{\color{Blue} \begin{footnotesize} 
		\begin{verbatim}
	contingency-stats -i gavin_2006_inf_1.8_cc_complexes_matrix.tab  \
                          -o gavin\_2006\_mcl\_inf\_1.8\_vs\_mips\_stats.tab 
		\end{verbatim} \end{footnotesize}
	}

Re-do this section with the \file{gavin\_2006\_rnsc\_clusters.rnsc} to obtain a file called \file{gavin\_2006\_rnsc\_vs\_mips\_stats.tab}.
\end{enumerate} 

\subsection{Interpretation of the results}
\subsubsection{Files description}
\paragraph{Contingency table}
As already explained in a previous section, having $n$ MIPS complexes and $m$ clusters, the contingency table $T$ is
a $n \cdot m$ matrix where row $i$ corresponds to the $i^{th}$
annotated complex, and column $j$ to the $j^{th}$ cluster. The value
of a cell $T_{i,j}$ indicates the number of proteins found in common
between complex $i$ and cluster $j$.

The clustering quality will be evaluated from this table by calculating the Sensitivity ($Sn$), the Positive predictive
value ($PPV$), the row wise separation ($Sep_{r}$) and the cluster separation ($Sep_{c}$).

\paragraph{Contingency table metrics}
  A list of metrics and their value. These will be described in the next section.

\subsubsection{Metrics description}
\paragraph{Sensitivity, Positive predictive value and geometric accuracy}

For each complex, we can calculate a sensitivity value. This corresponds to the maximal fraction of protein of a complex that are attributed by a clustering algorithm to the same cluster. $Sn$ measures how well proteins belonging to the same complex are grouped within the same cluster.
\[Sn_{i.}=\frac{max_{i.}(T_{ij})}{N_i}\]
where $N_i$ corresponds to the size of the complex.

Moreover, for each cluster $j$, we calculated the Positive Predictive Value ($PPV$) which corresponds to the maximal fraction of a cluster belonging to the same complex. This reflects the ability of this cluster to detect one complex.
\[PPV_{.j}=\frac{max_{.j}(T_{ij})}{M_j}\]
where $M_j$ corresponds to the cluster size.

To summarize these values at the level of the confusion table, we calculated the average of these values. 
First, we calculated their classical mean by averaging all the $PPV_{.j}$ and $Sn_{i.}$ values. We also calculated a weighted mean where the clusters and complexes have a weight proportional to their relative size on the the calculation of the mean.

\[Sn = \frac{\sum_{i=1}^nSn_{i.}}{n}\]
\[PPV = \frac{\sum_{j=1}^mPPV_{.j}}{m}\]
\[Sn_w = \frac{\sum_{i=1}^nN_iSn_{i.}}{\sum_{i=1}^nN_i}\]
\[PPV_w = \frac{\sum_{j=1}^mM_jPPV_{.j}}{\sum_{j=1}^mM_j}\]

Sensitivity and $PPV$ reflect two contradictory tendencies of the clustering. $Sn$ increases when all the proteins of the same complex are grouped in the same cluster and $PPV$ decreases when proteins coming from different complexes are grouped in the same cluster. If all the proteins of the network are grouped in the same cluster, we maximize the $Sn$ but the $PPV$ is almost 0. On the other hand, if each protein is placed in a different cluster, the $PPV$ is maximized  but the sensitivity is very low. A compromise must be found between these two cases by using another statistics. We defined the geometric accuracy as the geometrical mean of the $PPV$ and the $Sn$.

\[Acc_g = \sqrt{PPV \cdot Sn}\]

\paragraph{Separation}

We also defined another metrics called \textit{Separation} ($Sep$). High $Sep$ values indicated a high bidirectionnal correspondance between a cluster and a complex.

The row-wise separation estimates how a complex is isolated from the others. Its maximal value is 1 if this correspondance is perfect, i.e., when all the protein of a complex are grouped in one cluster and if this cluster does not contain any other protein. This maximal value may also be reached when the complex is separated between many clusters containing only members of the complex.

\[Sep_{r_{i.}} = \sum_{j=1}^{m}{(\frac{T_{i,j}}{\sum_{j=1}^{m}{T_{i,j}}} \cdot \frac{T_{i,j}}{\sum_{i=1}^{n}{T_{i,j}}})}\]

The column-wise separation indicates how well a cluster isolates one or more complex from the other clusters. The maximal value 1 indicates that a cluster contains all the elements of one or more complexes.  

\[Sep_{c_{.j}} = \sum_{i=1}^{n}{(\frac{T_{i,j}}{\sum_{j=1}^{m}{T_{i,j}}} \cdot \frac{T_{i,j}}{\sum_{i=1}^{n}{T_{i,j}}})}\]

As for the sensitivity and the $PPV$, for each clustering result, all values of $Sep_{c_{.j}}$ and $Sep_{r_{i.}}$ are averaged over all clusters and all complexes. We then calculate a global separation value by calculating the geometrical mean of the average row wise separation and of the average column wise separation.

\[Sep = \sqrt{Sep_c \cdot Sep_r}\]


\subsubsection{Score comparaison}

In the following, we can observe the statistics described in the previous paragraph computed for the clustering results of \program{RNSC} and \program{MCL}.

We can observe that \program{MCL} seems to produce slightly more valuable results as 
\begin{enumerate}
  \item The unweighted sensitivity is a bit higher for \program{MCL} than for \program{RNSC} and the weighted sensitivity is much higher.
  \item The $PPV$ is a only bit lower for \program{MCL} than for \program{RNSC}.

These results might certainly be explained by the large number of clusters found by \program{RNSC} compared to \program{MCL}. Indeed, the $PPV$ increases and the sensitivity decreases with the number of a clusters.
We can observe the same tendencies for the other metrics.

\item Global metrics (accuracy and separation) are generally higher for MCL than for RNSC

\end{enumerate}

\begin{tabular}{|l|l|l|}
\hline
metrics & RNSC & MCL \\ \hline
ncol & 470 & 189 \\ \hline
nrow & 220 & 220 \\ \hline
min & 0 & 0 \\ \hline
max & 18 & 27 \\ \hline
mean & 0.0086 & 0.0214 \\ \hline
sum & 889 & 889 \\ \hline
Sn & 0.603 & 0.652 \\ \hline
PPV & 0.424 & 0.472 \\ \hline
acc & 0.513 & 0.562 \\ \hline
acc\_g & 0.505 & 0.555 \\ \hline
Sn\_w & 0.622 & 0.767 \\ \hline
PPV\_w & 0.642 & 0.549 \\ \hline
acc\_w & 0.632 & 0.658 \\ \hline
acc\_g\_w & 0.632 & 0.649 \\ \hline
sep & 0.303 & 0.353 \\ \hline
sep\_c & 0.207 & 0.381 \\ \hline
sep\_r & 0.443 & 0.327 \\ \hline

\end{tabular}

\textit{Remark}: The following table was generated using the RSAT program \program{compare-scores}, see the help of this command line tool for more information