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

(* Clean up the interactions using p-values *)
(* Group C *)

SeedRandom[2];

singlekor=Import[StringJoin[path,"FILES2012/a_queryonly.dat"]];
singlekoc=Import[StringJoin[path,"FILES2012/a_genuinearray.dat"]];
doubleko=Import[StringJoin[path,"FILES2012/b_doublekoC.dat"]];
doublekosigma=Import[StringJoin[path,"FILES2012/b_doublekosigmaC.dat"]];
eps=Import[StringJoin[path,"FILES2012/c_epsC.dat"]];

nn=10000;

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

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[]];
If[doublekosigma[[i,j]]<(10^(-6)),doublekosigma[[i,j]]=(10^(-6))];

gi=RandomReal[StudentTDistribution[singlekor[[i,2]],((singlekor[[i,3]])/2),3],{nn}];
gj=RandomReal[StudentTDistribution[singlekoc[[j,2]],((singlekoc[[j,3]])/2),3],{nn}];
gij=RandomReal[StudentTDistribution[doubleko[[i,j]],((doublekosigma[[i,j]])/2),3],{nn}];

gi=DeleteCases[gi,x_/;x≤0];
gj=DeleteCases[gj,x_/;x≤0];
gij=DeleteCases[gij,x_/;x≤0];
mm=Min[{Length[gi],Length[gj],Length[gij]}];
gi=Take[gi,mm];
gj=Take[gj,mm];
gij=Take[gij,mm];


(*interaction=gij-gi gj;*)
interaction=Log[2,(gij/(gi gj))];

If[eps[[i,j]]>0,   (* regardless of the def of "interaction" the sign of eps is the same *)
q=Quantile[interaction,0.05];
If[q>0, 
interactionstrength[[i,j]]=q;
interactionsign[[i,j]]=1;
];
];

If[eps[[i,j]]<0,
q=Quantile[interaction,0.95];
If[q<0, 
interactionstrength[[i,j]]=-q;
interactionsign[[i,j]]=-1;
];
];

];
If[Mod[i,10]==1,Print[i]];
]
(*Export[StringJoin[path,"FILES2012/xC_litinteractionstrengthC10000.dat"],interactionstrength];
Export[StringJoin[path,"FILES2012/xC_litinteractionsignC10000.dat"],interactionsign];*)
(*Export[StringJoin[path,"FILES2012/xC_myinteractionstrengthC10000.dat"],interactionstrength];
Export[StringJoin[path,"FILES2012/xC_myinteractionsignC10000.dat"],interactionsign];*)
Export[StringJoin[path,"FILES2012/xC_myinteractionstrengthC10000_2.dat"],interactionstrength];
Export[StringJoin[path,"FILES2012/xC_myinteractionsignC10000_2.dat"],interactionsign];

(* Build the matrix of ranks *)

singlekor=Import[StringJoin[path,"FILES2012/a_queryonly.dat"]];
singlekoc=Import[StringJoin[path,"FILES2012/a_genuinearray.dat"]];
(*interactionstrength=Import[StringJoin[path,"FILES2012/xC_litinteractionstrengthC10000.dat"]];*)
(*interactionstrength=Import[StringJoin[path,"FILES2012/xC_myinteractionstrengthC10000.dat"]];*)
interactionstrength=Import[StringJoin[path,"FILES2012/xC_myinteractionstrengthC10000_2.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/xC_litinteractionrankC10000.dat"],interactionrank]*)
(*Export[StringJoin[path,"FILES2012/xC_myinteractionrankC10000.dat"],interactionrank]*)
Export[StringJoin[path,"FILES2012/xC_myinteractionrankC10000_2.dat"],interactionrank]

(* Look at the noise interoduced by the cleaning procedure in determining the p-values *)

litinteractionrank=Import[StringJoin[path,"FILES2012/xC_myinteractionrankC10000_2.dat"]];
myinteractionrank=Import[StringJoin[path,"FILES2012/xC_myinteractionrankC10000.dat"]];

nn=10000;
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,PlotRange->{{0,20000},{0,20000}},AxesStyle->18];
p2=Plot[{x},{x,0,Length[list]},PlotStyle->Red];
Show[p1,p2]

(* Compare the rankings obtained with the two definitions *)
(* Prepare a "list" comparing all the interactions according to the two definitions *)
(* In the following I will compile "list2" with high-quality GO interactions only.  *)

litinteractionrank=Import[StringJoin[path,"FILES2012/xC_litinteractionrankC10000.dat"]];
myinteractionrank=Import[StringJoin[path,"FILES2012/xC_myinteractionrankC10000.dat"]];

nn=10000;
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];

p3=ListPlot[list,AspectRatio->1,PlotRange->{{0,20000},{0,20000}},AxesStyle->18];
p4=Plot[{x},{x,0,Length[list]},PlotStyle->Red];
Show[p3,p4]

(* GO analysis *)

(* Selection of high-quality GO interactions *)

(* Import files *)

a1=Import[StringJoin[path,"GO_queryonly.dat"]];
a2=Import[StringJoin[path,"GO_genuinearray.dat"]];
singlekor=Import[StringJoin[path,"FILES2012/a_queryonly.dat"]];
singlekoc=Import[StringJoin[path,"FILES2012/a_genuinearray.dat"]];
doubleko=Import[StringJoin[path,"FILES2012/b_doublekoC.dat"]];
doublekosigma=Import[StringJoin[path,"FILES2012/b_doublekosigmaC.dat"]];
litinteractionrank=Import[StringJoin[path,"FILES2012/xC_litinteractionrankC10000.dat"]];
myinteractionrank=Import[StringJoin[path,"FILES2012/xC_myinteractionrankC10000.dat"]];
interactionsign=Import[StringJoin[path,"FILES2012/xC_myinteractionsignC10000.dat"]];


(* "tot" is a table in which the matrix element (i,j) is a list of genetic pairs in which gene 1 has i GO terms and gene 2 has j GO terms  *)

tot=Table[{},{35},{35}];
num=Table[0,{35},{35}];
For[i=1,i≤Length[a1],i++,
For[j=1,j≤Length[a2],j++,
c1=Length[a1[[i]]]-1;   (* I don't want to count the gene name *)
c2=Length[a2[[j]]]-1;   (* I don't want to count the gene name *)
ci=Length[Intersection[a1[[i]],a2[[j]]]];
If[c1<c2,
AppendTo[tot[[c1,c2]],ci];
num[[c1,c2]]++,
AppendTo[tot[[c2,c1]],ci];
num[[c2,c1]]++;
];
];
If[Mod[i,100]==1,Print[i]];
]

(* in "thresh" I indicate how many common GO terms a pair of genes should have to be in the top 10% of pairs with common GO terms. *)
(* in "num" I indicate how many pairs with i and j GO terms there are in the db. *)

thresh=Table[0,{35},{35}];
For[i=1,i≤35,i++,
For[j=i,j≤35,j++,
If[num[[i,j]]==1,
thresh[[i,j]]=tot[[i,j]][[1]];
thresh[[j,i]]=thresh[[i,j]];   (* to symmetrize thresh *)
];
If[num[[i,j]]≥2,
thresh[[i,j]]=Quantile[tot[[i,j]],0.90];
thresh[[j,i]]=thresh[[i,j]];   (* to symmetrize thresh *)
];
];
];
MatrixForm[N[thresh]]
MatrixForm[num]

(* Group B only: in "numhighq" I indicate how many high-quality GO pairs there are in the db. *)

numhighq=Table[0,{35},{35}];
nn=10000000;
k=1;
list2=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,

c1=Length[a1[[i]]]-1;   (* I don't want to count the gene name *)
c2=Length[a2[[j]]]-1;   (* I don't want to count the gene name *)
ci=Length[Intersection[a1[[i]],a2[[j]]]];

If[ci>thresh[[c1,c2]],
list2[[k]]={litinteractionrank[[i,j]],myinteractionrank[[i,j]]};
If[c1<c2,
numhighq[[c1,c2]]++;,
numhighq[[c2,c1]]++;;
];

k++;
];
];
];
If[Mod[i,100]==1,Print[i]];
]
list2=DeleteCases[list2,0];

MatrixForm[numhighq]


(* Fraction of selected high-quality GO interactions (to be compared with the 5% objective). *)

ratio=Table[0,{35},{35}];
For[i=1,i≤35,i++,
For[j=i,j≤35,j++,
If[num[[i,j]]≥1,
ratio[[i,j]]=100N[((numhighq[[i,j]])/(num[[i,j]]))];
ratio[[j,i]]=100N[((numhighq[[i,j]])/(num[[i,j]]))];   (* to symmetrize thresh *)
];
];
];
N[ratio]//MatrixForm

<<PlotLegends`
ShowLegend[MatrixPlot[ratio-6,ColorFunction->"Rainbow",FrameStyle->18],{ColorData["Rainbow"][#1]&,10,"0%","10%",LegendPosition->{1.1,-0.4}}]

list3={};
For[i=1,i≤Length[list2],i++,
If[list2[[i,1]]>10000&&list2[[i,2]]>10000,Continue[]];
AppendTo[list3,list2[[i]]];
]

(* Plot high-quality interactions stored in "list2", computed right after thresh above. *)

p5=ListPlot[list,AspectRatio->1,PlotStyle->{Green,PointSize[0.005]},PlotRange->{{0,20000},{0,20000}},AxesStyle->18];
p6=ListPlot[list3,AspectRatio->1,PlotStyle->{Red,PointSize[0.01]},PlotRange->{{0,1000},{0,1000}}];
Show[p5,p6]

(* For the top-ranking pairs decide if they are GO interacting or not.  *)

bingo={};
interactingpairs={};
For[k=1,k≤10000,k++,
i=Position[myinteractionrank,k][[1,1]];
j=Position[myinteractionrank,k][[1,2]];

If[Length[Intersection[a1[[i]],a2[[j]]]]>  thresh[[Length[a1[[i]]]-1,Length[a2[[j]]]-1]],
AppendTo[bingo,1];
AppendTo[interactingpairs,{k,singlekor[[i,1]],singlekoc[[j,1]],singlekor[[i,2]],singlekoc[[j,2]]}],
AppendTo[bingo,0];
];

If[Mod[k,100]==1,Print[k]];
]

Export[StringJoin[path,"FILES2012/xC_bingomy_p010.dat"],bingo];
Export[StringJoin[path,"FILES2012/xC_interactingpairsmy_p010.dat"],interactingpairs];


bingo={};
interactingpairs={};
For[k=1,k≤10000,k++,
i=Position[litinteractionrank,k][[1,1]];
j=Position[litinteractionrank,k][[1,2]];

If[Length[Intersection[a1[[i]],a2[[j]]]]>  thresh[[Length[a1[[i]]]-1,Length[a2[[j]]]-1]],
AppendTo[bingo,1];
AppendTo[interactingpairs,{k,singlekor[[i,1]],singlekoc[[j,1]],singlekor[[i,2]],singlekoc[[j,2]]}],
AppendTo[bingo,0];
];

If[Mod[k,100]==1,Print[k]];
]

Export[StringJoin[path,"FILES2012/xC_bingolit_p010.dat"],bingo];
Export[StringJoin[path,"FILES2012/xC_interactingpairslit_p010.dat"],interactingpairs];

myinteractingpairs=Import[StringJoin[path,"FILES2012/xC_interactingpairsmy_p010.dat"]];
litinteractingpairs=Import[StringJoin[path,"FILES2012/xC_interactingpairslit_p010.dat"]];

myinteractingpairs2=Drop[myinteractingpairs,None,1];
litinteractingpairs2=Drop[litinteractingpairs,None,1];

Length[myinteractingpairs2]
Length[litinteractingpairs2]

inters=Intersection[litinteractingpairs2,myinteractingpairs2];

litonly=Complement[litinteractingpairs2,inters];
myonly=Complement[myinteractingpairs2,inters];

Length[inters]
Length[litonly]
Length[myonly]

inters=Drop[inters,None,2];
myonly=Drop[myonly,None,2];
litonly=Drop[litonly,None,2];

Show[ListPlot[{inters,myonly,litonly},PlotStyle->{Directive[Black,PointSize->0.012],Directive[Red,PointSize->0.014],Directive[Blue,PointSize->0.012]},AspectRatio->1,AxesOrigin->{0.2,0.2}],Plot[1,{x,0,1.2},PlotStyle->Black],Plot[1000(x-1),{x,0,1.2},PlotStyle->Black],AxesStyle->18]

(* GO enrichment by picking random pairs *)

Total[Total[numhighq]]
Total[Total[num]]
slope=((Total[Total[numhighq]])/(Total[Total[num]]))//N

(* GO interactions with p < 0.05 ; all pairs considered, regardless of the single knockout effects *)

bingolit=Import[StringJoin[path,"FILES2012/xC_bingolit_p010.dat"]];
bingomy=Import[StringJoin[path,"FILES2012/xC_bingomy_p010.dat"]];
bingolit=Flatten[bingolit];
bingomy=Flatten[bingomy];
bingolit=Accumulate[bingolit];
bingomy=Accumulate[bingomy];
plit=ListPlot[bingolit,PlotStyle->Blue,AxesStyle->18];
pmy=ListPlot[bingomy,PlotStyle->Red];
prandom=Plot[slope x,{x,0,10000},PlotStyle->Directive[Green,Dashed,Thick]];
scale=ListPlot[{0,0},PlotRange->{{0,1000},{0,Max[{bingolit[[1000]],bingomy[[1000]]}]+1}},AxesStyle->18];
Print["Number of 'GO interactions' detected by the traditional definition: ",bingolit[[1000]],"/1000"]
Print["Number of 'GO interactions' detected by the our definition: ",bingomy[[1000]],"/1000"]
Show[scale,plit,pmy,prandom]
Print["Number of 'GO interactions' detected by the traditional definition: ",bingolit[[10000]],"/10000"]
Print["Number of 'GO interactions' detected by the our definition: ",bingomy[[10000]],"/10000"]
Show[plit,pmy,prandom]