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

common=Import[StringJoin[path,"FILES2012/a_common.dat"]];
queryonly=Import[StringJoin[path,"FILES2012/a_queryonly.dat"]];
arrayonly=Import[StringJoin[path,"FILES2012/a_arrayonly.dat"]];

singlekor=Join[common,queryonly];
singlekoc=Join[common,arrayonly];

strengthsA=Import[StringJoin[path,"FILES2012/xA_litinteractionstrengthA10000.dat"]];
strengthsB=Import[StringJoin[path,"FILES2012/xB_litinteractionstrengthB10000.dat"]];
strengthsC=Import[StringJoin[path,"FILES2012/xC_litinteractionstrengthC10000.dat"]];

signsA=Import[StringJoin[path,"FILES2012/xA_litinteractionsignA10000.dat"]];
signsB=Import[StringJoin[path,"FILES2012/xB_litinteractionsignB10000.dat"]];
signsC=Import[StringJoin[path,"FILES2012/xC_litinteractionsignC10000.dat"]];

(*strengthsA=Import[StringJoin[path,"FILES2012/xA_myinteractionstrengthA10000.dat"]];
strengthsB=Import[StringJoin[path,"FILES2012/xB_myinteractionstrengthB10000.dat"]];
strengthsC=Import[StringJoin[path,"FILES2012/xC_myinteractionstrengthC10000.dat"]];

signsA=Import[StringJoin[path,"FILES2012/xA_myinteractionsignA10000.dat"]];
signsB=Import[StringJoin[path,"FILES2012/xB_myinteractionsignB10000.dat"]];
signsC=Import[StringJoin[path,"FILES2012/xC_myinteractionsignC10000.dat"]];*)

Length[singlekor]
Length[singlekoc]

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

For[i=1,i≤Length[common],i++,
For[j=1,j≤Length[common],j++,
interactionstrengths[[i,j]]=strengthsA[[i,j]];
interactionsigns[[i,j]]=signsA[[i,j]];
];
]

For[i=1,i≤Length[common],i++,
For[j=1,j≤Length[arrayonly],j++,
interactionstrengths[[i,1138+j]]=strengthsB[[i,j]];
interactionsigns[[i,1138+j]]=signsB[[i,j]];
];
]

For[i=1,i≤Length[queryonly],i++,
For[j=1,j≤Length[singlekoc],j++,
interactionstrengths[[1138+i,j]]=strengthsC[[i,j]];
interactionsigns[[1138+i,j]]=signsC[[i,j]];
];
]

Export[StringJoin[path,"FILES2012/xABC_litinteractionstrength10000.dat"],interactionstrengths];
Export[StringJoin[path,"FILES2012/xABC_litinteractionssign10000.dat"],interactionsigns];
(*Export[StringJoin[path,"FILES2012/xABC_myinteractionstrength10000.dat"],interactionstrengths];
Export[StringJoin[path,"FILES2012/xABC_myinteractionssign10000.dat"],interactionsigns];*)

(* Build the matrix of ranks *)

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

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

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

Export[StringJoin[path,"FILES2012/xABC_litinteractionrank10000.dat"],interactionrank]
(*Export[StringJoin[path,"FILES2012/xABC_myinteractionrank10000.dat"],interactionrank]*)

(* 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/xABC_litinteractionrank10000.dat"]];
myinteractionrank=Import[StringJoin[path,"FILES2012/xABC_myinteractionrank10000.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]

For[i=1,i≤Length[list],i++,
If[((list[[i,1]])/(list[[i,2]]))>80,Print[list[[i]]]];
]

For[i=1,i≤Length[list],i++,
If[((list[[i,2]])/(list[[i,1]]))>8,Print[list[[i]]]];
]

myrank=2;
{singlekor[[Position[myinteractionrank,myrank][[1,1]]]],singlekoc[[Position[myinteractionrank,myrank][[1,2]]]]}
myrank=4;
{singlekor[[Position[myinteractionrank,myrank][[1,1]]]],singlekoc[[Position[myinteractionrank,myrank][[1,2]]]]}
myrank=32;
{singlekor[[Position[myinteractionrank,myrank][[1,1]]]],singlekoc[[Position[myinteractionrank,myrank][[1,2]]]]}
myrank=52;
{singlekor[[Position[myinteractionrank,myrank][[1,1]]]],singlekoc[[Position[myinteractionrank,myrank][[1,2]]]]}
myrank=226;
{singlekor[[Position[myinteractionrank,myrank][[1,1]]]],singlekoc[[Position[myinteractionrank,myrank][[1,2]]]]}
myrank=470;
{singlekor[[Position[myinteractionrank,myrank][[1,1]]]],singlekoc[[Position[myinteractionrank,myrank][[1,2]]]]}

YHR038W Mitochondrial ribosome recycling factor,essential for mitochondrial protein synthesis and for the maintenance of the respiratory function of mitochondria
YNL082W ATP-binding protein required for mismatch repair in mitosis and meiosis;functions as a heterodimer with Mlh1p,binds double-and single-stranded DNA via its N-terminal domain,similar to E.coli MutL (1,2,3 and see Summary Paragraph)

YER086W Threonine deaminase,catalyzes the first step in isoleucine biosynthesis;expression is under general amino acid control;ILV1 locus exhibits highly positioned nucleosomes whose organization is independent of known ILV1 regulation (2,3,4)
YML131W Putative protein of unknown function with similarity to medium chain dehydrogenase/reductases;expression induced by stresses including osmotic shock,DNA damaging agents,and other chemicals;GFP-fusion protein localizes to the cytoplasm (1,2,3,4,5,6,7,8,9 and see Summary Paragraph)

YDR159W Nuclear pore-associated protein;required for biogenesis of the small ribosomal subunit;component of TREX-2 complex (Sac3p-Thp1p-Sus1p-Cdc31p) involved in transcription elongation and mRNA export from the nucleus;involved in post-transcriptional tethering of active genes to the nuclear periphery and to non-nascent mRNP;similar to the human germinal center-associated nuclear protein (GANP) (2,3,4,5,6,7,8,9)
YMR202W C-8 sterol isomerase,catalyzes the isomerization of the delta-8 double bond to the delta-7 position at an intermediate step in ergosterol biosynthesis (1,2,3 and see Summary Paragraph)

YLR396C ATP-binding protein that is a subunit of the HOPS complex and of the CORVET tethering complex;essential for protein sorting,vesicle docking and fusion at the vacuole (9,10,11,12)
YPL184C RNA-binding protein proposed to be involved in translational regulation;binds specific categories of mRNAs,including those that contain upstream open reading frames (uORFs) and internal ribosome entry sites (IRES) (2)

YAL021C Component of the CCR4-NOT transcriptional complex,which is involved in regulation of gene expression;component of the major cytoplasmic deadenylase,which is involved in mRNA poly(A) tail shortening (4,5,6)
YJR092W Anillin-like protein involved in bud-site selection;required for the axial budding pattern;localizes with septins to the bud neck in mitosis and may constitute an axial landmark for the next round of budding;required for the formation and disassembly of the double septin ring structure,and generally for septin organization;in vivo substrate of Cdc28p/Clb2p (1,2,3,4)

YLR417W Component of the ESCRT-II complex;contains the GLUE (GRAM Like Ubiquitin binding in EAP45) domain which is involved in interactions with ESCRT-I and ubiquitin-dependent sorting of proteins into the endosome (5,6 and see Summary Paragraph)
YBR043C Multidrug transporter of the major facilitator superfamily,has a role in polyamine homeostasis;expression is upregulated under polyamine stress;required for resistance to quinidine,barban,cisplatin,and bleomycin (1,3,4 and see Summary Paragraph)

(* GO analysis *)

(* Selection of high-quality GO interactions *)

litinteractionrank=Import[StringJoin[path,"FILES2012/xABC_litinteractionrank10000.dat"]];
myinteractionrank=Import[StringJoin[path,"FILES2012/xABC_myinteractionrank10000.dat"]];

(* Import files *)

GOcommon=Import[StringJoin[path,"GO_common.dat"]];
GOarrayonly=Import[StringJoin[path,"GO_arrayonly.dat"]];
GOqueryonly=Import[StringJoin[path,"GO_queryonly.dat"]];

a1=Join[GOcommon,GOqueryonly];
a2=Join[GOcommon,GOarrayonly];

interactionsign=Import[StringJoin[path,"FILES2012/xABC_myinteractionssign10000.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]

(* 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]]];
]

list4={};
For[i=1,i≤Length[list3],i++,
If[list3[[i,1]]<10000&&list3[[i,2]]<10000,
AppendTo[list4,list3[[i]]];
];
]

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

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

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

p5=ListPlot[list,AspectRatio->1,PlotStyle->{Gray,PointSize[0.003]},PlotRange->{{0,30000},{0,30000}},Frame->True,FrameStyle->Directive[18,Thick],AspectRatio->1,FrameTicks->{{0,10000,20000,30000},{0,10000,20000,30000},None,None}];
(*p6=ListPlot[list3,AspectRatio->1,PlotStyle->{Red,PointSize[0.01]},PlotRange->{{0,1000},{0,1000}}];*)
p7=ListPlot[list4,AspectRatio->1,PlotStyle->{Purple,PointSize[0.01]},PlotRange->{{0,1000},{0,1000}}];
p8=ListPlot[list5,AspectRatio->1,PlotStyle->{Blue,PointSize[0.012]},PlotRange->{{0,1000},{0,1000}}];
p9=ListPlot[list6,AspectRatio->1,PlotStyle->{Red,PointSize[0.01]},PlotRange->{{0,1000},{0,1000}}];
Show[p5,(*p6,*)p7,p8,p9,ImageSize->300]

Export[StringJoin[path,"FinalFIGS/6d.pdf"],,ImageResolution->600]

Export["GOrankings.png",]

(* 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/xABC_bingomy_p010.dat"],bingo];
Export[StringJoin[path,"FILES2012/xABC_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/xABC_bingolit_p010.dat"],bingo];
Export[StringJoin[path,"FILES2012/xABC_interactingpairslit_p010.dat"],interactingpairs];

myinteractingpairs=Import[StringJoin[path,"FILES2012/xABC_interactingpairsmy_p010.dat"]];
litinteractingpairs=Import[StringJoin[path,"FILES2012/xABC_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[Purple,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]

Export["GOwhere.png",]

(* 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/xABC_bingolit_p010.dat"]];
bingomy=Import[StringJoin[path,"FILES2012/xABC_bingomy_p010.dat"]];
bingolit=Flatten[bingolit];
bingomy=Flatten[bingomy];
bingolit=Accumulate[bingolit];
bingomy=Accumulate[bingomy];
plit=ListPlot[bingolit,PlotStyle->Red,Frame->True,FrameStyle->Directive[18,Thick],AspectRatio->1,FrameTicks->{{0,2000,4000,6000,8000,10000},{0,400,800,1200},None,None}];
pmy=ListPlot[bingomy,PlotStyle->Blue];
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}},Frame->True,FrameStyle->Directive[18,Thick],AspectRatio->1,FrameTicks->{{0,500,1000},{0,100,200,300},None,None}];
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,ImageSize->300]

Export[StringJoin[path,"FinalFIGS/6c.pdf"],,ImageResolution->600]

Export["GOfinders.png",]

common=Import[StringJoin[path,"FILES2012/a_common.dat"]];
queryonly=Import[StringJoin[path,"FILES2012/a_queryonly.dat"]];
arrayonly=Import[StringJoin[path,"FILES2012/a_arrayonly.dat"]];

singlekor=Join[common,queryonly];
singlekoc=Join[common,arrayonly];

litinteractionrank=Import[StringJoin[path,"FILES2012/xABC_litinteractionrank10000.dat"]];
myinteractionrank=Import[StringJoin[path,"FILES2012/xABC_myinteractionrank10000.dat"]];

thresholdrank=10000;
rank=litinteractionrank;

(* common *)

zcommon={};

For[i=1,i≤Length[common],i++,
counter=0;

For[j=1,j≤Length[singlekoc],j++,
If[rank[[i,j]]≤thresholdrank,counter++];
];
For[j=1,j≤Length[queryonly],j++,
If[rank[[1138+j,i]]≤thresholdrank,counter++];
];

AppendTo[zcommon,{common[[i,2]],N[(counter/(Length[singlekoc]+Length[queryonly]))]}];
]

(* arrayonly *)

zarrayonly={};

For[j=1,j≤Length[arrayonly],j++,
counter=0;

For[i=1,i≤Length[singlekor],i++,
If[rank[[i,1138+j]]≤thresholdrank,counter++];
];
AppendTo[zarrayonly,{arrayonly[[j,2]],N[(counter/(Length[singlekor]))]}];
]

(* queryonly *)

zqueryonly={};

For[i=1,i≤Length[queryonly],i++,
counter=0;

For[j=1,j≤Length[singlekoc],j++,
If[rank[[1138+i,j]]≤thresholdrank,counter++];
];

AppendTo[zqueryonly,{queryonly[[i,2]],N[(counter/(Length[singlekoc]))]}];
]

combined=Join[zcommon,zarrayonly,zqueryonly];

lowersummary={};
summary={};
uppersummary={};
For[k=1,k≤140,k++,
total={};

For[i=1,i≤Length[combined],i++,
If[0.01 (k-5.5)<combined[[i,1]]≤0.01 (k+5.5),
AppendTo[total,combined[[i,2]]];
];
];

If[Length[total]>4,
AppendTo[lowersummary,{0.01 k,Quantile[total,0.25]}];
AppendTo[summary,{0.01 k,Quantile[total,0.5]}];
AppendTo[uppersummary,{0.01 k,Quantile[total,0.75]}];
];
]

ListLinePlot[{lowersummary,summary,uppersummary},PlotStyle->{Blue,Directive[Blue,Thick],Blue}(*PlotStyle->{Red,Directive[Red,Thick],Red}*),PlotRange->All,AxesStyle->18,AxesOrigin->{0,0},Ticks->{{0.4,0.6,0.8,1.0,1.2},{0,0.002,0.004,0.006,0.008,0.010}}]
ListLinePlot[summary,PlotRange->{{0.32,1.2},{0,0.009}},PlotStyle->Directive[Red,Thick],Frame->True,FrameStyle->Directive[18,Thick],AspectRatio->0.75,FrameTicks->{{0.2,0.4,0.6,0.8,1.0,1.2},{0,0.002,0.004,0.006,0.008,0.010},None,None},ImageSize->300]

Show[,]

Export[StringJoin[path,"FinalFIGS/6b.pdf"],,ImageResolution->600]

ListLinePlot[{lowersummary,summary,uppersummary},PlotStyle->{Blue,Directive[Blue,Thick],Blue}(*PlotStyle->{Red,Directive[Red,Thick],Red}*),PlotRange->All,AxesStyle->18,AxesOrigin->{0,0},Ticks->{{0,0.2,0.4,0.6,0.8,1.0,1.2},{0,0.003,0.006,0.009,0.012,0.015}}]
ListLinePlot[summary,PlotStyle->Directive[Red,Thick],PlotRange->All,AxesStyle->18,AxesOrigin->{0,0},Ticks->{{0,0.2,0.4,0.6,0.8,1.0,1.2},{0,0.003,0.006,0.009,0.012,0.015}}]

Show[,]

Export["GOnetwork.png",]

Show[,]

(* myrank, common, arrayonly, queryonly, combined *)

(* Import files *)

GOcommon=Import[StringJoin[path,"GO_common.dat"]];
GOarrayonly=Import[StringJoin[path,"GO_arrayonly.dat"]];
GOqueryonly=Import[StringJoin[path,"GO_queryonly.dat"]];

a1=Join[GOcommon,GOqueryonly];
a2=Join[GOcommon,GOarrayonly];

interactionsign=Import[StringJoin[path,"FILES2012/xABC_myinteractionssign10000.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.95];
thresh[[j,i]]=thresh[[i,j]];   (* to symmetrize thresh *)
];
];
];
MatrixForm[N[thresh]]
MatrixForm[num]

(* 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-2.75,ColorFunction->"Rainbow",FrameStyle->18],{ColorData["Rainbow"][#1]&,10,"0%","5%",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/xABC_bingomy_p005.dat"],bingo];
Export[StringJoin[path,"FILES2012/xABC_interactingpairsmy_p005.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/xABC_bingolit_p005.dat"],bingo];
Export[StringJoin[path,"FILES2012/xABC_interactingpairslit_p005.dat"],interactingpairs];

myinteractingpairs=Import[StringJoin[path,"FILES2012/xABC_interactingpairsmy_p005.dat"]];
litinteractingpairs=Import[StringJoin[path,"FILES2012/xABC_interactingpairslit_p005.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[Purple,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/xABC_bingolit_p005.dat"]];
bingomy=Import[StringJoin[path,"FILES2012/xABC_bingomy_p005.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]

(* Import files *)

GOcommon=Import[StringJoin[path,"GO_common.dat"]];
GOarrayonly=Import[StringJoin[path,"GO_arrayonly.dat"]];
GOqueryonly=Import[StringJoin[path,"GO_queryonly.dat"]];

a1=Join[GOcommon,GOqueryonly];
a2=Join[GOcommon,GOarrayonly];

interactionsign=Import[StringJoin[path,"FILES2012/xABC_myinteractionssign10000.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.75];
thresh[[j,i]]=thresh[[i,j]];   (* to symmetrize thresh *)
];
];
];
MatrixForm[N[thresh]]
MatrixForm[num]

(* 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-15,ColorFunction->"Rainbow",FrameStyle->18],{ColorData["Rainbow"][#1]&,10,"0%","25%",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/xABC_bingomy_p025.dat"],bingo];
Export[StringJoin[path,"FILES2012/xABC_interactingpairsmy_p025.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/xABC_bingolit_p025.dat"],bingo];
Export[StringJoin[path,"FILES2012/xABC_interactingpairslit_p025.dat"],interactingpairs];

myinteractingpairs=Import[StringJoin[path,"FILES2012/xABC_interactingpairsmy_p025.dat"]];
litinteractingpairs=Import[StringJoin[path,"FILES2012/xABC_interactingpairslit_p025.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[Purple,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/xABC_bingolit_p025.dat"]];
bingomy=Import[StringJoin[path,"FILES2012/xABC_bingomy_p025.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]