+1 for suggesting igraph
Hello,
Using STRING-DB, I need to build a data frame containing for each pair of proteins the number of interactors between them...
For example, if protein A and protein B have a common line in the node_node_links table, they are in direct interaction, and the number of intermediates will be 0. If protein A and protein D are not in direct interaction, but protein A interacts with protein C which in turn interacts with protein D, the number of intermediates between A and D will be 1. And so on, for the 20770 proteins...
A B 0
A D 1
...
I wonder if any of you ever faced this problem and would be able to give me a time and memory efficient solution.
Thanks !
4 answers
If I were to do it now, I would pick either of the following two methods:
1) Calculate the distance for each node pair using a graph library
This approach requires you to convert the data into a graph model for your library of choice. I previously did a similar thing with Jung, which is Java-based, and was really happy with it. It is hard to estimate the required memory, but I was able to calculate the length of shortest path for all possible node pairs on a PPI network I obtained from GeneMania (human) in a few hours on a standard laptop -- so it was not that bad.
If you want to see an example, you can check out this repository for an old project of mine and you will find my Java implementation that utilizes Jung (the documentation is a little bit scarce, but the code is simple):
https://bitbucket.org/armish/pathwayguy
and here is the code that calculates the pairwise shortest-path distances:
AbstractPairwiseDistanceMethod.java
I am guessing you can also go with the R-library igraph to do this calculation; I never tried it, but from the things I heard from my colleagues, I can tell it is a good alternative to Jung (especially if you are an R-person).
2) Work with adjancecy matrix
This requires you to convert the interaction file into an adjacency matrix (all nodes vs all nodes) and then multiplying it by itself, i.e. taking powers. I believe this is also not that memory intensive, especially if you go with a matrix multiplication method that is optimized for sparse matrices (R and Matlab both have many of those). For this alternative, the idea goes like this:
- You create the adjancency matrix and every non-zero value in a cell indicates a distance of 0 for the corresponding node pair (A^1)
- You multiply the matrix by itself and now every non-zero value in a cell indicates a distance of 1 for the corresponding node pair (A^2)
- You multiply the matrix by itself and now every non-zero value in a cell indicates a distance of 2 for the corresponding node pair (A^3)
- You multiply the matrix by itself and now every non-zero value in a cell indicates a distance of 3 for the corresponding node pair (A^4) ...
It is not a great approach, but surprisingly people who are familiar with matrix-based operations in either R or Matlab find this approach relatively easier to implement.
I don't have the node_node_links table, but let's assume that you can obtain the partners as a tab-delimited text file like this:
A B
A C
B C
C D
Then a graph can be created using R igraph:
library(igraph)
d <- read.table("myfile.txt", sep = "\t", header = F, stringsAsFactors = T)
g <- graph.data.frame(d)
A list of all igraph methods is here. Something like what you want can be achieved using shortest.paths():
shortest.paths(g)
A B C D
A 0 1 1 2
B 1 0 1 2
C 1 1 0 1
D 2 2 1 0
I don't know how long this would take using the STRING database; it could potentially be a very long time.
the input would be ftp://string-db.org/STRING/9.1/human.links.txt.gz
prior warning: pfbusson works in my lab.
I've implemented a java idea: https://github.com/lindenb/jvarkit/wiki/Biostar92368
Put all the interactions in a berkeleyDB with key=protein, value=set(interactors)
(...)
while(r.hasNext())
{
String line=r.next();
String tokens[]=tab.split(line);
for(int i=0;i< 2;++i)
{
String a=tokens[i==0?0:1];
String b=tokens[i==0?1:0];
StringBinding.stringToEntry(a, key);
OperationStatus status=c.getSearchKey(key, data, LockMode.DEFAULT);
Set<String> partners;
if(status==OperationStatus.SUCCESS)
{
partners=bindingSet.entryToObject(data);
}
else
{
partners=new HashSet<String>(1);
}
partners.add(b);
bindingSet.objectToEntry(partners,data);
if(status==OperationStatus.SUCCESS)
{
status=c.putCurrent(data);
}
else
{
status=c.put(key, data);
}
}
}
(...)
and then recursively scan all the paths to go from protein1 to protein2:
private int recursive(
Transaction txn,
final String endProt,
Stack<String> path,
int bestDepth,
final int maxDepth
)
{
DatabaseEntry key=new DatabaseEntry();
DatabaseEntry data=new DatabaseEntry();
StringBinding.stringToEntry(path.lastElement(), key);
if(this.database.get(txn,key,data,null)==OperationStatus.SUCCESS)
{
Set<String> partners=this.bindingSet.entryToObject(data);
for(String prot2:partners)
{
int depth=-1;
if(prot2.equals(endProt))
{
depth=path.size()-1;
}
else if(path.size()<maxDepth && !path.contains(prot2))
{
path.push(prot2);
depth=recursive(txn,endProt,path,bestDepth,maxDepth);
path.pop();
}
if(depth!=-1 && (bestDepth==-1 || bestDepth>depth))
{
bestDepth=depth;
}
}
}
return bestDepth;
}
Example:
$ cat input.txt
A1 A2
A2 A3
A3 A4
A4 A5
A1 A6
$ mkdir -p tmp
$ java -jar dist/biostar92368.jar -D tmp input.txt
A1 A2 0
A1 A3 1
A1 A4 2
A1 A6 0
A2 A3 0
A2 A4 1
A2 A5 2
A2 A6 1
A3 A4 0
A3 A5 1
A3 A6 2
A4 A5 0
The networkx library for python has a method for it: networkx.shortest_path_length(self.network) http://networkx.github.io/documentation/latest/reference/generated/networkx.algorithms.shortest_paths.generic.shortest_path.html
Log in to answer this question.
"data frame" suggests you want a solution in R?
yes Neil, I know him: that's what pfbusson wants.