Le associazioni fisiche tra proteine sono oggetto di interesse medico principalmente per due aspetti: un interesse sulla fisiologia del corpo e un aspetto patologico. Il primo approccio cerca di individuare i pattern che si creano tra proteine interagenti e i conseguenti processi biologici che avvengono all'interno di cellule; il secondo cerca le interazioni con le proteine che sono "viziate" (proteine intatte ma la cui genesi non si è svolta correttamente): queste interagiscono con altre proteine viziate, tendono a legarsi e a creare patologie all'interno dell'organismo.
L'obiettivo è scoprire le strutture di clustering nella rete PPI, cioè determinare una struttura di cluster di proteine in cui ogni proteina è più vicina alle altre all'interno dello stesso insieme rispetto alle proteine al di fuori dell'insieme.
Non solo per l'importanza di questi studi ma anche per la loro complessità, il cercare le interazioni tra proteine (Protein-Protein Interactions, PPI) è un compito molto stimolante e impegnativo. Molti infatti sono gli algoritmi scritti, di cui uno dei più popolari è quello per l'edge-betweeness di Girvan-Newmans.
Questo progetto propone una nuova implementazione di quest'ultimo, e successivamente analizza la struttura in cluster tramite la modularità (indice Q). Ci siamo basati sull'articolo "A MapReduce-Based Parallel Clustering Algorithm for Large Protein-Protein Interaction Networks", di Li Liu e altri.
L'algoritmo è implementato in java-Spark, con il paradigma Map-Reduce. Inoltre alcuni risultati di questo algoritmo verranno visualizzati con un database NoSQL Neo4J.
I dataset utilizzati sono disponibile al link https://dip.doe-mbi.ucla.edu/dip/Download.cgi?SM=3, con data 18/05/2012.
Sono formati da 16 variabili, di cui abbiamo preso solo le prime due, rappresentanti il nodo di interazione A e il nodo di interazione B. Le relazioni descritte nel dataset sono indirette.
Il dataset Mus Musculus è composto da 728 nodi e 634 archi.
Il dataset Drosophila Melanogaster è composto da 284 nodi e 297 archi.
| Metodi | Contenuto |
|---|---|
| AccoppiaPath | Interfaccia che individua i singoli archi dei path |
| CalcolaVicini | Interfaccia che trova i vicini dei nodi |
| ComputaCluster | Interfaccia che associa a ogni nodo il cluster di appartenenza |
| DisaccoppiaNodi | Interfaccia che da ogni arco ricava i nodi dagli archi e ci associa il cluster |
| EstrapolaNodi | Interfaccia che ricava i nodi distinti dagli archi |
| MapForward | Interfaccia che implementa il MapForward |
| ReduceForward | Interfaccia che implementa il ReduceForward |
| Dataset Java | Contenuto |
|---|---|
| DatasetGrande | Mus Musculus |
| DatasetPiccolo | Drosophila Melanogaster |
| Dataset Neo4J | Contenuto |
|---|---|
| archiBest | file con lista degli archi del dataset con migliore modularità |
| archiFirst | file con lista degli archi del dataset iniziale |
| archiLast | file con lista degli archi del dataset all'ultima iterazione dell'algoritmo |
| nodiBesteFirst | file con lista dei nodi del dataset iniziale e con migliore modularità |
| nodiLast | file con lista dei nodi del dataset all'ultima iterazione dell'algoritmo |
Per poter importare i dataset sia in Java sia in Neo4J sono state apportate delle modifiche ai file di partenza:
Su Excel abbiamo eliminato tutte le colonne diverse dai due nodi di interazione. Il risultato è una lista di archi con i nodi di interazione separati da uno spazio.
L'algoritmo Clustering-MR lavora su una tabella composta da nxn righe (dove n è il numero di nodi distinti che compongono il grafo). Questa tabella viene definita Adjacency List ed ogni record è registrato mediante una coppia Chiave,Valore.
La chiave è composta dal nome del nodo e dal nome di un altro nodo definito root: le nxn righe rappresentano tutte le possibili combinazioni tra i nodi (funzione cartesian). I valori invece conservano, per ogni coppia nodo-root, i nodi vicini alla proteina in chiave, la distanza tra il nodo e la root, il colore del nodo (status) e il percorso che collega il nodo alla root (path).
Il colore è WHITE se il record non è stato considerato, GRAY se è considerato ma non ha concluso la fase di aggiornamento delle caratteristiche, BLACK se è stato completamente aggiornato.
Di default la distanza è impostata pari a MAX, il colore WHITE e il path NULL. Nei record il cui nodo è uguale alla root la distanza è pari a 0 e il colore GRAY; questi saranno i nodi da cui l'algoritmo Forward-MR comincerà ad aggiornare i valori.
Se un record presenta un nodo e una root che non sono collegati da nessun arco (si trovano in cluster diversi), l'Adjacency List conserverà i valori di default (dist = MAX, col = WHITE).
I vicini sono stati trovati con il metodo CalcolaVicini, la cui logica è quella di duplicare ogni record invertendo NodeID e root, eliminare i record uguali e ridurli con ReduceByKey. Per tutti gli altri valori abbiamo eseguito un mapToPair:
JavaRDD<String> nodi = proteine.flatMap(new EstrapolaNodi()).distinct();
JavaPairRDD<String, String> nodiroot = nodi.cartesian(nodi);
JavaPairRDD<String, String> vicini = proteine.flatMapToPair(new CalcolaVicini()).distinct().reduceByKey((x, y) -> (x + "," + y));
JavaPairRDD<String, Tuple2<String, String>> nodirootvicini = nodiroot.join(vicini);
JavaPairRDD<String,String> finale = nodirootvicini.mapToPair(x-> {
if(x._1.equalsIgnoreCase(x._2._1))
return new Tuple2<String, String>(x._1 +" "+ x._2._1, x._2._2 +" "+ 0 +" "+ "GRAY" +" "+ "null");
else return new Tuple2<String,String> (x._1 +" "+x._2._1 , x._2._2 +" "+ Integer.MAX_VALUE +" "+ "WHITE" +" "+ "null");
});Il Forward-MR è un algoritmo che partendo dalle root percorre ed aggiorna i valori (distanza, colore e path) in parallelo di ogni vicino del nodo in chiave.
I primi record toccati sono quelli dove il NodeID è uguale al root. L'algoritmo genera una coppia (K,V) per ogni vicino del nodo, in cui la distanza è la distanza del nodo madre + 1, il colore è GRAY e il path è composto dal path del nodo madre + il NodeID del nodo madre. Il colore della riga del nodo madre viene aggiornato in BLACK affinchè non venga più considerata nelle iterazioni successive.
La fase di reduce seleziona tra le diverse coppie (K,V) dei vicini che sono state generate quella dove la distanza è minima.
Successivamente l'algoritmo ricomincia dagli archi grigi, finché questi non sono più presenti, ovvero fino a quando sono stati processati tutti i vicini dei vicini possibili per i root in parallelo. L'output finale è una tabella composta da nxn righe, come l' input iniziale, dove però è registrata la distanza minima e il path tra il NodeID e la Root.
do {
// Map
JavaPairRDD<String, String> map = finale.flatMapToPair(new MapForward());
// Reduce
finale = map.reduceByKey((x,y)->x +"_"+ y).mapToPair(new ReduceForward());
}while(finale.filter(x->x._2.split(" ")[2].equalsIgnoreCase("GRAY")).isEmpty() == false);Il Backward-MR è un algoritmo che prende in input l'output del Forward-MR, e determina quanti shortest path passano per ogni arco del nostro grafo.
E' composto da una fase Map in cui si applica il metodo AccoppiaPath per ottenere il path del record frammentato in tutti gli archi che lo compongono, accompagnati dalla frequenza 1.
Nella fase Reduce l'algoritmo va a sommare le frequenze dei record con la stessa chiave e restituisce ogni arco con il numero di volte in cui è presente nei path dell'adjacency list.
// Map
JavaPairRDD<String,Integer> archetti = finale.flatMapToPair(new AccoppiaPath());
// Reduce
JavaPairRDD<String,Integer> archettiridotto = archetti.reduceByKey((x,y)-> x+y);Successivamente l'algoritmo calcola per ogni arco il valore di betweenness centrality, che ci dà una misura di centralità di ogni arco all'intero di un grafo:
Long denominatore = finale.filter(x-> !x._2.split("")[2].equalsIgnoreCase("WHITE")&& x._1.split(" ")[0]!=x._1.split(" ")[1]).count();
JavaPairRDD<Float,String> BCeOrdinato = archettiridotto.mapToPair(x->new Tuple2<Float,String>((float) x._2/denominatore,x._1.split(",")[0] + " " + x._1.split(",")[1])).sortByKey(false);Si elimina l'arco con il valore di BC(e) maggiore, si aggiorna il dataset, si calcola l'indice Q e si ripete l'algoritmo.
L'indice Q va a misurare la modularità del grafo, ovvero la qualità della suddivisione delle proteine in cluster a seconda delle loro interazioni. Un valore di Q prossimo all'1 indica una suddivisione buona, dove le proteine più connesse si trovano all'interno dello stesso cluster, ma sono debolmente connesse con i nodi delle altre comunità della rete.
JavaPairRDD<String,Float> conta = assegnoClusterB.flatMapToPair(new ComputaCluster()).reduceByKey((x,y)-> x+y).sortByKey().mapToPair(x->{
if(x._1.contains(" "))
return new Tuple2<String,Float> (x._1.split(" ")[0],x._2/den);
else return new Tuple2<String,Float> (x._1, (float) Math.pow(x._2/den,2));
}).reduceByKey((x,y)-> y-x);
Float Q = conta.values().reduce((x,y)->x+y);Per il calcolo di Q bisogna conoscere il cluster di appartenenza di ogni nodo di un arco, così da discernere archi i cui nodo fanno parte dello stesso cluster e archi i cui nodi connettono due comunità distinte. Per fare ciòabbiamo applicato agli archi il cui colore è rimasto WHITE una serie di trasformazioni (i nodi di questi archi si trovano in cluster diversi per definizione):
JavaPairRDD<String,String> passo1 = finale.filter(x->x._2.split(" ")[2].equalsIgnoreCase("WHITE")).sortByKey().mapToPair(x-> new Tuple2<String,String> (x._1.split(" ")[0],x._1.split(" ")[1])).reduceByKey((x,y)->x+" "+y);
JavaPairRDD<String,String> passo2 = passo1.mapToPair(x->new Tuple2<String, String>(x._2,x._1)).reduceByKey((x,y)->x+" "+y);
JavaPairRDD<String, Long> passo3 = passo2.zipWithIndex().flatMapToPair(new DisaccoppiaNodi());
JavaPairRDD<String, Tuple2<String,Long>> assegnoClusterA = copiaproteine.join(passo3).mapToPair(x-> new Tuple2<String,Tuple2<String,Long>> (x._2._1, new Tuple2<String,Long>( x._1 , x._2._2)));
JavaPairRDD<String, String> assegnoClusterB = assegnoClusterA.join(passo3).mapToPair(x-> new Tuple2 <String,String> (x._1 +" "+ x._2._1._1 , x._2._2 +" "+ x._2._1._2 ));L'algoritmo si interrompe quando le sue azioni diventano superflue, ovvero quando l'eliminazione dell'arco più attraversato restituisce un numero di cluster minore rispetto a quello del passaggio precedente. In una situazione di questo tipo ci troviamo di fronte ad un grafo strutturato in cluster composti da uno o due proteine l'uno.
Per l'importazione in Neo4J i file txt sono stati trasformati in csv. Attraverso il metodo EstrapolaNodi abbiamo ottenuto la lista dei nodi che è richiesta, insieme alla lista degli archi, per la creazione del grafo in Neo4J.
Sono stati caricati i file dei grafi sia del datset di partenza, sia del dataset che presenta il valore della migliore modularità, sia del dataset all'ultima iterazione.
Le query per l'importazione sono:
LOAD CSV WITH HEADERS FROM 'file:///nodi.csv' AS row
CREATE (p:Protein)
SET p = row
RETURN pLOAD CSV WITH HEADERS FROM 'file:///archi.csv' AS line
MATCH (p1:Protein {nodi:line.nodoA}), (p2:Protein {nodi:line.nodoB})
MERGE (p1)-[:INTERACTION]-(p2)L'output del Forward MR è una lista di tutte le combinazioni dei nodi con annesse distanze e shortest path. Per verificare la correttezza di questi risultati abbiamo applicato un algoritmo presente nella libreria Graph Data Science:
CALL gds.alpha.allShortestPaths.stream
({nodeProjection: 'Protein', relationshipProjection:
{INTERACTION: {type: 'INTERACTION' }}})
YIELD sourceNodeId, targetNodeId, distance
WITH sourceNodeId, targetNodeId, distance
WHERE gds.util.isFinite(distance) = true
MATCH (source:Protein) WHERE id(source) = sourceNodeId
MATCH (target:Protein) WHERE id(target) = targetNodeId
RETURN source.nodi AS source, target.nodi AS target, distanceCon questa interrogazione otteniamo uno tra i percorsi più brevi e la distanza minima tra tutte le coppie di nodi.
Per controllare la distanza e il path più breve tra due nodi dati abbiamo utilizzato la query:
MATCH (p:Protein {nodi: 'DIP-27777N'}), (q:Protein {nodi: 'DIP-29672N'}),
ShortestPath=shortestPath((p)-[*]-(q))
RETURN ShortestPath, length(ShortestPath)Per verificare la correttezza dell'assegnazione dei cluster in Java abbiamo utilizzato l'algoritmo Weakly Connected Components della libreria Graph Data Science:
CALL gds.wcc.stream({nodeProjection: 'Protein',
relationshipProjection: {INTERACTION: {type: 'INTERACTION' }}}')
YIELD nodeId, componentId
RETURN gds.util.asNode(nodeId).nodi AS Nodo, componentId
ORDER BY componentId, NodoCon questa interrogazione otteniamo il cluster di appartenenza di ogni nodo.
Per assegnare una nuova etichetta con il cluster ai nodi nel grafo utilizziamo la query:
CALL gds.wcc.write(
{nodeProjection: 'Protein', relationshipProjection: {INTERACTION: {type: 'INTERACTION' }}}', { writeProperty: 'componentId' })
YIELD nodePropertiesWritten, componentCount;Così da poter controllare i nodi appartenenti a un cluster scelto con:
MATCH n:Protein{componentID:'22'})
RETURN n