path="c:/users/dre2/desktop/DRYGINresubmission/";

(* Create my epistasis matrix for group B: common x arrayonly *)

singlekor=Import[StringJoin[path,"FILES2012/a_common.dat"]];
singlekoc=Import[StringJoin[path,"FILES2012/a_arrayonly.dat"]];
doubleko=Import[StringJoin[path,"FILES2012/b_doublekoB.dat"]];

eps=Table[0,{Length[singlekor]},{Length[singlekoc]}];

For[n1=1,n1≤Length[singlekor],n1++,
For[n2=1,n2≤Length[singlekoc],n2++,

a=doubleko[[n1,n2]];

If[a=="NaN",
eps[[n1,n2]]=0,
If[a≤0,
Print["Non positive growthrate for the pair (",n1,",",n2,"): ",doubleko[[n1,n2]]],
eps[[n1,n2]]=Abs[Log[2,(a/(singlekor[[n1,2]]singlekoc[[n2,2]]))]];
];
];

];
]
Export[StringJoin[path,"FILES2012/y_myepsB.dat"],eps]

(* Create the literature epistasis matrix for group B: common x arrayonly *)

singlekor=Import[StringJoin[path,"FILES2012/a_common.dat"]];
singlekoc=Import[StringJoin[path,"FILES2012/a_arrayonly.dat"]];
doubleko=Import[StringJoin[path,"FILES2012/b_doublekoB.dat"]];

eps=Table[0,{Length[singlekor]},{Length[singlekoc]}];

For[n1=1,n1≤Length[singlekor],n1++,
For[n2=1,n2≤Length[singlekoc],n2++,

a=doubleko[[n1,n2]];

If[a=="NaN",
eps[[n1,n2]]=0,
If[a≤0,
Print["Non positive growthrate for the pair (",n1,",",n2,"): ",doubleko[[n1,n2]]],
eps[[n1,n2]]=Abs[a-singlekor[[n1,2]]singlekoc[[n2,2]]];
];
];

];
]
Export[StringJoin[path,"FILES2012/y_litepsB.dat"],eps]

(* Group B *)

singlekor=Import[StringJoin[path,"FILES2012/a_common.dat"]];
singlekoc=Import[StringJoin[path,"FILES2012/a_arrayonly.dat"]];
doubleko=Import[StringJoin[path,"FILES2012/b_doublekoB.dat"]];
doublekosigma=Import[StringJoin[path,"FILES2012/b_doublekosigmaB.dat"]];

(* Build the matrices of ranks *)

(* ..using our definition *)

interactionstrength=Import[StringJoin[path,"FILES2012/y_myepsB.dat"]];


rank=Table[0,{Length[singlekor] Length[singlekoc]}];
flat=Flatten[interactionstrength];
ord=Reverse[Ordering[flat]];

For[i=1,i≤Length[singlekor] Length[singlekoc],i++,
rank[[ord[[i]]]]=i;
If[Mod[i,100000]==1,Print[i]];
]
interactionrank=Partition[rank,Length[singlekoc]];

Total[Total[interactionrank]]
Sum[i,{i,1,Length[flat]}]

Export[StringJoin[path,"FILES2012/y_myinteractionrankB.dat"],interactionrank];

(* ..using the standard definition *)

interactionstrength=Import[StringJoin[path,"FILES2012/y_litepsB.dat"]];

rank=Table[0,{Length[singlekor] Length[singlekoc]}];
flat=Flatten[interactionstrength];
ord=Reverse[Ordering[flat]];

For[i=1,i≤Length[singlekor] Length[singlekoc],i++,
rank[[ord[[i]]]]=i;
If[Mod[i,100000]==1,Print[i]];
]
interactionrank=Partition[rank,Length[singlekoc]];

Total[Total[interactionrank]]
Sum[i,{i,1,Length[flat]}]

Export[StringJoin[path,"FILES2012/y_litinteractionrankB.dat"],interactionrank];

(* Compare the rankings obtained using my definition and the standard definition. *)

singlekor=Import[StringJoin[path,"FILES2012/a_common.dat"]];
singlekoc=Import[StringJoin[path,"FILES2012/a_arrayonly.dat"]];
litinteractionrank=Import[StringJoin[path,"FILES2012/y_litinteractionrankB.dat"]];
myinteractionrank=Import[StringJoin[path,"FILES2012/y_myinteractionrankB.dat"]];

nn=100000;
k=1;
list=Table[0,{2nn}];
For[i=1,i≤Length[singlekor],i++,
For[j=1,j≤Length[singlekoc],j++,
If[litinteractionrank[[i,j]]<nn||myinteractionrank[[i,j]]<nn,
list[[k]]={litinteractionrank[[i,j]],myinteractionrank[[i,j]]};
k++;
];
];
]
list=DeleteCases[list,0];

p1=ListPlot[list,AspectRatio->1,PlotStyle->{Green,PointSize[0.001]}];
p2=Plot[{x},{x,0,Length[list]},PlotStyle->Red];
Show[p1,p2]

p1=ListPlot[list,AspectRatio->1,PlotStyle->{Green,PointSize[0.001]},PlotRange->{{0,1000000},{0,1000000}}];
p2=Plot[{x},{x,0,Length[list]},PlotStyle->Red];
Show[p1,p2]

(* Compare the rankings obtained using my definition and the standard definition. *)

singlekor=Import[StringJoin[path,"FILES2012/a_common.dat"]];
singlekoc=Import[StringJoin[path,"FILES2012/a_arrayonly.dat"]];
doubleko=Import[StringJoin[path,"FILES2012/b_doublekoB.dat"]];
doublekosigma=Import[StringJoin[path,"FILES2012/b_doublekosigmaB.dat"]];
litinteractionrank=Import[StringJoin[path,"FILES2012/y_litinteractionrankB.dat"]];
myinteractionrank=Import[StringJoin[path,"FILES2012/y_myinteractionrankB.dat"]];

nn=100000;
k=1;
list=Table[0,{2nn}];
For[i=1,i≤Length[singlekor],i++,
For[j=1,j≤Length[singlekoc],j++,
If[doubleko[[i,j]]<0.75||doubleko[[i,j]]>1.00,Continue[]];
If[((doublekosigma[[i,j]])/(doubleko[[i,j]]))>0.1,Continue[]];
If[((singlekor[[i,3]])/(singlekor[[i,2]]))>0.1,Continue[]];
If[((singlekoc[[j,3]])/(singlekoc[[j,2]]))>0.1,Continue[]];
If[litinteractionrank[[i,j]]<nn||myinteractionrank[[i,j]]<nn,
list[[k]]={litinteractionrank[[i,j]],myinteractionrank[[i,j]]};
k++;
];
];
]
list=DeleteCases[list,0];

p1=ListPlot[list,AspectRatio->1,PlotStyle->{Blue,PointSize[0.001]},PlotRange->{{0,100000},{0,100000}}]

Show[,,,,Plot[x,{x,0,100000},PlotStyle->{Thick,Red}]]
Show[Plot[x,{x,0,1000},AspectRatio->1,PlotStyle->{Thick,Red}],,,,,Plot[x,{x,0,1000},PlotStyle->{Thick,Red}]]

litinteractionstrength=Import[StringJoin[path,"FILES2012/y_litepsB.dat"]];
myinteractionstrength=Import[StringJoin[path,"FILES2012/y_myepsB.dat"]];
litinteractionrank=Import[StringJoin[path,"FILES2012/y_litinteractionrankB.dat"]];
myinteractionrank=Import[StringJoin[path,"FILES2012/y_myinteractionrankB.dat"]];
litinteractionsign=Import[StringJoin[path,"FILES2012/y_litinteractionsignB.dat"]];
myinteractionsign=Import[StringJoin[path,"FILES2012/y_myinteractionsignB.dat"]];


kor=Import[StringJoin[path,"FILES2012/a_common.dat"]];
koc=Import[StringJoin[path,"FILES2012/a_arrayonly.dat"]];
kor=Transpose[kor][[1]];
koc=Transpose[koc][[1]];
Print["Number of rows in group B: ",Length[kor]]
Print["Number of columns in group B: ",Length[koc]]

zz=Import["c:/users/dre2/desktop/interactome2.dat"];
Print["Number of entries in the interactome database: ",Length[zz]]
genes=Union[Transpose[zz][[1]],Transpose[zz][[2]]];
Print["Number of unique genes in the interactome database: ",Length[genes]]

counter=0;
For[i=1,i≤Length[genes],i++,
If[Count[kor,genes[[i]]]==1,counter++];
]
Print["Number of genes in the interactome database indexing rows in group B: ",counter]

counter=0;
For[i=1,i≤Length[genes],i++,
If[Count[koc,genes[[i]]]==1,counter++];
]
Print["Number of genes in the interactome database indexing columns in group B: ",counter]

lista={};
For[i=1,i≤Length[zz],i++,
If[Count[kor,zz[[i,1]]]==1&&Count[koc,zz[[i,2]]]==1,
AppendTo[lista,{Position[kor,zz[[i,1]]][[1,1]],Position[koc,zz[[i,2]]][[1,1]]}];
];
If[Count[kor,zz[[i,2]]]==1&&Count[koc,zz[[i,1]]]==1,
AppendTo[lista,{Position[kor,zz[[i,2]]][[1,1]],Position[koc,zz[[i,1]]][[1,1]]}];
];
]
Print["Number of pairs in group B: ",Length[lista]]

we=0;
they=0;
none=0;
abovethreshold=0;
For[i=1,i≤Length[lista],i++,
If[litinteractionstrength[[lista[[i,1]],lista[[i,2]]]]==0&&myinteractionstrength[[lista[[i,1]],lista[[i,2]]]]==0,
none++;
Continue[]
];
If[litinteractionrank[[lista[[i,1]],lista[[i,2]]]]>1000000&&myinteractionrank[[lista[[i,1]],lista[[i,2]]]]>1000000,abovethreshold++;Continue[]];
If[litinteractionrank[[lista[[i,1]],lista[[i,2]]]]<myinteractionrank[[lista[[i,1]],lista[[i,2]]]],
they++,
we++;
];
Print[litinteractionrank[[lista[[i,1]],lista[[i,2]]]],"\t\t\t",myinteractionrank[[lista[[i,1]],lista[[i,2]]]],"\t\t\t",myinteractionsign[[lista[[i,1]],lista[[i,2]]]]];
]
Print["They: ",they,"    We: ",we,"    Above threshold: ",abovethreshold,"    None: ",none];

Length[Flatten[litinteractionstrength]]-Count[Flatten[litinteractionstrength],0]
Length[Flatten[myinteractionstrength]]-Count[Flatten[myinteractionstrength],0]

singlekor=Import[StringJoin[path,"FILES2012/a_common.dat"]];
singlekoc=Import[StringJoin[path,"FILES2012/a_arrayonly.dat"]];
doubleko=Import[StringJoin[path,"FILES2012/b_doublekoB.dat"]];
doublekosigma=Import[StringJoin[path,"FILES2012/b_doublekosigmaB.dat"]];

all=Table[0,{13},{13}];
mytop=Table[0,{13},{13}];
littop=Table[0,{13},{13}];

For[i=1,i≤Length[singlekor],i++,
For[j=1,j≤Length[singlekoc],j++,

If[doubleko[[i,j]]=="NaN",Continue[]];
If[doubleko[[i,j]]<0,Continue[]];

x=Ceiling[10singlekor[[i,2]]];
y=Ceiling[10singlekoc[[i,2]]];

all[[x,y]]++;
If[myinteractionrank[[i,j]]≤100000,mytop[[x,y]]++];
If[litinteractionrank[[i,j]]≤100000,littop[[x,y]]++];

];
]

myfraction=Table[0,{13},{13}];
litfraction=Table[0,{13},{13}];

all=((all+Transpose[all])/2.);
mytop=((mytop+Transpose[mytop])/2.);
littop=((littop+Transpose[littop])/2.);

For[x=1,x≤13,x++,
For[y=1,y≤13,y++,

If[all[[x,y]]>0,
myfraction[[x,y]]=((mytop[[x,y]])/(all[[x,y]]));
litfraction[[x,y]]=((littop[[x,y]])/(all[[x,y]]));
];

];
]

Total[Total[all]]

p0=Plot3D[0.5,{x,0,13},{y,0,13},PlotStyle->Yellow];
p1=ListPlot3D[((423all)/(Total[Total[all]])),InterpolationOrder->1,PlotRange->All];
Show[p1,p0]
p2=ListPlot3D[myfraction,InterpolationOrder->1,PlotStyle->Green]
p3=ListPlot3D[litfraction,InterpolationOrder->1,PlotStyle->Orange]
Show[p2,p3]

logall=Log[2.,all];
logmytop=Log[2.,mytop];
loglittop=Log[2.,littop];
logmyfraction=Log[2.,myfraction];
loglitfraction=Log[2.,litfraction];

p1=ListPlot3D[logall,InterpolationOrder->1]
p2=ListPlot3D[logmyfraction,InterpolationOrder->1,PlotStyle->Green]
p3=ListPlot3D[loglitfraction,InterpolationOrder->1,PlotStyle->Orange]
Show[p2,p3]

(* 1000 *)

(* 10000 *)

(* 100000 *)

(* Study network properties *)

kor=Import[StringJoin[path,"FILES2012/a_common.dat"]];
koc=Import[StringJoin[path,"FILES2012/a_arrayonly.dat"]];

litinteractionrank=Import[StringJoin[path,"FILES2012/y_litinteractionrankB.dat"]];
myinteractionrank=Import[StringJoin[path,"FILES2012/y_myinteractionrankB.dat"]];

interactionrank=myinteractionrank;
adj=Table[0,{Length[kor]},{Length[koc]}];   (* adjacency matrix *)

For[i=1,i≤Length[kor],i++,
For[j=1,j≤Length[koc],j++,
If[interactionrank[[i,j]]<10000,adj[[i,j]]=1];
];
]

v=Join[((Total[adj])/(Length[koc])),((Total[Transpose[adj]])/(Length[kor]))];
v=Reverse[Sort[v]];
lit=Transpose[{Range[Length[v]],v}];

ListLogLogPlot[Take[lit,1000],PlotRange->All,PlotStyle->Blue]

ListLogLogPlot[Take[my,1000],PlotRange->All,PlotStyle->Red]

Show[,]

v=Join[((Total[adj])/(Length[koc])),((Total[Transpose[adj]])/(Length[kor]))];
lit=Transpose[{Join[Transpose[koc][[2]],Transpose[kor][[2]]],v}];

ListPlot[lit,PlotRange->All,PlotStyle->Blue]

v=Join[((Total[adj])/(Length[koc])),((Total[Transpose[adj]])/(Length[kor]))];
my=Transpose[{Join[Transpose[koc][[2]],Transpose[kor][[2]]],v}];

ListPlot[my,PlotRange->All,PlotStyle->Red]

lm=LinearModelFit[my,x,x]
lm=LinearModelFit[lit,x,x]

pmy=Show[ListPlot[my,PlotRange->All,PlotStyle->Red],Plot[LinearModelFit[my,x,x][y],{y,0,1.3},PlotStyle->Red]]

plit=Show[ListPlot[lit,PlotRange->All,PlotStyle->Blue],Plot[LinearModelFit[lit,x,x][y],{y,0,1.3},PlotStyle->Blue]]

Show[pmy,plit]

list=Table[0,{13}];
tot=Table[0,{13}];
For[i=1,i≤Length[my],i++,
list[[Ceiling[10 my[[i,1]]]]]+=my[[i,2]];
tot[[Ceiling[10 my[[i,1]]]]]++;
]

For[i=1,i≤13,i++,
If[tot[[i]]≠0,list[[i]]=((list[[i]])/(tot[[i]]))];
]

ListPlot[list,PlotRange->All,PlotStyle->Red,Filling->Axis]

list=Table[0,{13}];
tot=Table[0,{13}];
For[i=1,i≤Length[lit],i++,
list[[Ceiling[10 lit[[i,1]]]]]+=lit[[i,2]];
tot[[Ceiling[10 lit[[i,1]]]]]++;
]

For[i=1,i≤13,i++,
If[tot[[i]]≠0,list[[i]]=((list[[i]])/(tot[[i]]))];
]

ListPlot[list,PlotRange->All,PlotStyle->Blue,Filling->Axis]

Total[list]//N

list//N

Total[my]
Total[lit]