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

(* GO analysis *)

(* Selection of high-quality GO interactions *)

(* Import files *)

a1=Import[StringJoin[path,"GO_common.dat"]];
a2=Import[StringJoin[path,"GO_arrayonly.dat"]];
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/x_litinteractionrankB10000.dat"]];
myinteractionrank=Import[StringJoin[path,"FILES2012/x_myinteractionrankB10000.dat"]];
interactionsign=Import[StringJoin[path,"FILES2012/x_myinteractionsignB10000.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 x% 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.80];
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]

(* See how many high-quality GO interactions are selected if I pick interactions with n. of common terms > thresh.. *)

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

(* ..or ≥ thresh. *)
Total[Total[numhighq]]
Total[Total[num]]
((Total[Total[numhighq]])/(Total[Total[num]]))//N

(* Fraction of selected high-quality GO interactions with n. of common terms > thresh.. *)

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-11,ColorFunction->"Rainbow"],{ColorData["Rainbow"][#1]&,10,"0%","20%",LegendPosition->{1.1,-0.4}}]

coarsegrainedhighq=Table[0,{5},{5}];
coarsegrainednum=Table[0,{5},{5}];
For[i=1,i≤35,i++,
For[j=1,j≤35,j++,
coarsegrainedhighq[[Floor[(i-1)/7]+1,Floor[(j-1)/7]+1]]+=numhighq[[i,j]];
coarsegrainednum[[Floor[(i-1)/7]+1,Floor[(j-1)/7]+1]]+=num[[i,j]];
];
]
MatrixForm[coarsegrainedhighq]
MatrixForm[coarsegrainednum]

coarsegrainedratio=Table[0,{5},{5}];
For[i=1,i≤5,i++,
For[j=i,j≤5,j++,
If[coarsegrainednum[[i,j]]≥2,
coarsegrainedratio[[i,j]]=100N[((coarsegrainedhighq[[i,j]])/(coarsegrainednum[[i,j]]))];
coarsegrainedratio[[j,i]]=100N[((coarsegrainedhighq[[i,j]])/(coarsegrainednum[[i,j]]))];
];
];
]
MatrixForm[coarsegrainedratio]
MatrixPlot[coarsegrainedratio]

(* Growth rate and number of high-quality GO interactions *)

(* study the number of GO interactions as a function of the growth rate effect of a gene *)

GOint=Table[0,{Length[singlekor]},{Length[singlekoc]}];
For[i=1,i≤Length[singlekor],i++,
For[j=1,j≤Length[singlekoc],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[ci>thresh[[c1,c2]],
GOint[[i,j]]=1;
];

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

Needs["ErrorBarPlots`"]

Length[Total[GOint]]
Length[singlekoc]

growthrates=Transpose[singlekoc][[2]];
interactions=N[((Total[GOint])/1138)];     (* each column has 1138 elements, so i'm counting the number of interactions _per gene_ *)
listino=Transpose[{growthrates,interactions}];

ListPlot[listino,PlotRange->All]

inter=Table[0,{13}];
counter=Table[0,{13}];
For[i=1,i≤Length[listino],i++,
inter[[Round[10listino[[i,1]]+1]]]+=listino[[i,2]];
counter[[Round[10listino[[i,1]]+1]]]++;
]
inter
counter
meaninter={};
For[i=1,i≤13,i++,
If[counter[[i]]>1,
AppendTo[meaninter,{{0.1i-0.05,N[((inter[[i]])/(counter[[i]]))]},ErrorBar[(1/(Sqrt[counter[[i]]]))Sqrt[(((inter[[i]])^2)/(counter[[i]]-1))-((counter[[i]])/(counter[[i]]-1))(((((inter[[i]])/(counter[[i]]))))^2)]]}];
];
If[counter[[i]]==1,
AppendTo[meaninter,{{0.1i-0.05,N[((inter[[i]])/(counter[[i]]))]},ErrorBar[0]}];
];
]
meaninter
ErrorListPlot[meaninter,Filling->Axis,AxesOrigin->{0,0},PlotRange->All,PlotStyle->Purple]

Length[Total[Transpose[GOint]]]
Length[singlekor]

growthrates=Transpose[singlekor][[2]];
interactions=N[((Total[Transpose[GOint]])/2742)];     (* each row has 2742 elements, so i'm counting the number of interactions _per gene_ *)
listino=Transpose[{growthrates,interactions}];

ListPlot[listino,PlotRange->All]

inter=Table[0,{13}];
counter=Table[0,{13}];
For[i=1,i≤Length[listino],i++,
inter[[Round[10listino[[i,1]]+1]]]+=listino[[i,2]];
counter[[Round[10listino[[i,1]]+1]]]++;
]
inter
counter
meaninter={};
For[i=1,i≤13,i++,
If[counter[[i]]>1,
AppendTo[meaninter,{{0.1i-0.05,N[((inter[[i]])/(counter[[i]]))]},ErrorBar[(1/(Sqrt[counter[[i]]]))Sqrt[(((inter[[i]])^2)/(counter[[i]]-1))-((counter[[i]])/(counter[[i]]-1))(((((inter[[i]])/(counter[[i]]))))^2)]]}];
];
If[counter[[i]]==1,
AppendTo[meaninter,{{0.1i-0.05,N[((inter[[i]])/(counter[[i]]))]},ErrorBar[0]}];
];
]
meaninter
ErrorListPlot[meaninter,Filling->Axis,AxesOrigin->{0,0},PlotRange->All]

Show[,]

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

p0=ListPlot[list,AspectRatio->1,PlotStyle->{Green,PointSize[0.005]},PlotRange->{{0,200000},{0,200000}}];
p1=ListPlot[list2,AspectRatio->1,PlotStyle->{Red,PointSize[0.005]},PlotRange->{{0,1000},{0,1000}}];
Show[p0,p1]

(* p<0.05 *)

(* Need to chekc how many total pairs and how many high quality pairs I identify *)

Length[list]
Length[list2]

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

If[singlekor[[i,2]]>2.0&&singlekoc[[j,2]]>2.0,   (* only if I want to see how things change with growth rate *)
AppendTo[bingo,0];
If[Mod[k,100]==1,Print[k]];
Continue[];
];

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/x20_bingomy_p005_all.dat"],bingo];

litinteractingpairs=interactingpairs

myinteractingpairs=interactingpairs

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

inters=Intersection[litinteractingpairs2,myinteractingpairs2];
Length[inters]

traditiononly=Complement[litinteractingpairs2,inters];
myonly=Complement[myinteractingpairs2,inters];
Length[traditiononly]
Length[myonly]
Length[inters]

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

Show[ListPlot[{inters,myonly,traditiononly},PlotStyle->{Directive[Green,PointSize->0.01],Directive[Red,PointSize->0.01],Directive[Blue,PointSize->0.01]},AspectRatio->1,AxesOrigin->{0,0}],Plot[1,{x,0,1.2}],Plot[1000(x-1),{x,0,1.2}]]

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

bingolit=Import[StringJoin[path,"FILES2012/x20_bingolit_p005_all.dat"]];
bingomy=Import[StringJoin[path,"FILES2012/x20_bingomy_p005_all.dat"]];
bingolit=Flatten[bingolit];
bingomy=Flatten[bingomy];
bingolit=Accumulate[bingolit];
bingomy=Accumulate[bingomy];
plit=ListPlot[bingolit,PlotStyle->Blue];
pmy=ListPlot[bingomy,PlotStyle->Red];
prandom=Plot[0.10952554537748022`x,{x,0,10000},PlotStyle->Directive[Green,Dashed,Thick]];
scale=ListPlot[{0,0},PlotRange->{{0,1000},{0,Max[{bingolit[[1000]],bingomy[[1000]]}]+1}}];
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,pmy,plit,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[pmy,plit,prandom]

all=0;
highq=0;
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[]];

all++;
If[singlekor[[i,2]]>2.0&&singlekoc[[j,2]]>2.0,
Continue[];
];
If[Length[Intersection[a1[[i]],a2[[j]]]]>thresh[[Length[a1[[i]]]-1,Length[a2[[j]]]-1]],
highq++;
];
];
If[Mod[i,100]==1,Print[i]];
]
all

(* all growth rates *)
highq
all
N[(highq/all)]

selected={};
For[i=1,i≤Length[singlekor],i++,
For[j=1,j≤Length[singlekoc],j++,
If[((litinteractionrank[[i,j]])/(myinteractionrank[[i,j]]))>30,AppendTo[selected,{Round[((litinteractionrank[[i,j]])/(myinteractionrank[[i,j]]))],singlekor[[i,1]],singlekoc[[j,1]],myinteractionrank[[i,j]],litinteractionrank[[i,j]]}]];
];
If[Mod[i,100]==1,Print[i]];
]

selected//MatrixForm

Reverse[Sort[Flatten[Drop[selected,None,-4]]]]

rankratio=Table[0,{Length[singlekor]},{Length[singlekoc]}];
For[i=1,i≤Length[singlekor],i++,
For[j=1,j≤Length[singlekoc],j++,
rankratio[[i,j]]=((litinteractionrank[[i,j]])/(myinteractionrank[[i,j]]));
];
If[Mod[i,100]==1,Print[i]];
]

sorted=Sort[Flatten[rankratio]];

Round[Reverse[Take[sorted,-20]]]
Round[(1/(Take[sorted,20]))]

otherranks={};
For[k=1,k≤100,k++,
i=Position[myinteractionrank,k][[1,1]];
j=Position[myinteractionrank,k][[1,2]];

AppendTo[otherranks,litinteractionrank[[i,j]]];

If[litinteractionrank[[i,j]]>1000,
Print[singlekor[[i,1]]," ",singlekoc[[j,1]],"  ",litinteractionrank[[i,j]],"  ",myinteractionrank[[i,j]],"  ",Round[((litinteractionrank[[i,j]])/(myinteractionrank[[i,j]]))]];
];

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

Sort[otherranks]

Sort[otherranks]

Sort[otherranks]

Sort[otherranks]