% analyze the network data %:-use_module(library(lists)). %:-['./2010_9_25/25sep2010_cancer_ver1.pl']. %:-['day14_1500_kegg_no_equ.xml-Original.pl']. %:-['../observations/day14/metabolomicInfo_lycocyc201Version.pl']. %%%%%%%%%%%%%% geneIDs_for_oneReaction(reversed-RID,GeneIDs):- !, concept_name(RID,GeneIDs,false). geneIDs_for_oneReaction(RID,GeneIDs):- concept_name(RID,GeneIDs,false). topCla0([RID],T) :- topCla(RID,T). interestMet(RID-EITIs):- member(EI-TI,EITIs), member(EI,[23,31,33,34,49,40,42]). analyzeIsolated(IsolatedEIExplanations):- maplist(analyzeIsolated_for_oneMet,IsolatedEIExplanations). analyzeIsolated_for_oneMet([EI]-[MinScore-MinCandidates|_]):- maplist(getCandidateEITI,MinCandidates,MinCandidates_EITIs), maplist(outputEachH,MinCandidates_EITIs). getCandidateEITI(MinCandidate,EI-TI):- last(MinCandidate,endNode-[EI-TI]-[EI-TI]). analyzeShared(CombineExplanation):- maplist(extractSharedReactions,CombineExplanation,Shares0), flatten(Shares0,Shares1), sort(Shares1,Shares), %remove_duplicates(Shares1,Shares), maplist(outputSharedInformation,Shares). outputSharedInformation(ClaID-EIs):- ClaID=[rs-{RID,enzymeLimiting,Hstate,Time}],!, outputRegulatedMets(EIs), outputEachCla(ClaID),nl. outputSharedInformation(ClaID-EIs). % do nothing /*outputSharedInformation(ClaID-EIs):- outputRegulatedMets(EIs), outputEachCla(ClaID),nl. */ outputRegulatedMets([EI]):- !, ex(EI,concentration(MID,Change,Time),1), concept_accession(MID,'KEGG',MName,true), write(MName),write(' ('),write(Change),write(') are regulated together by '). outputRegulatedMets([EI|EIs]):- ex(EI,concentration(MID,Change,Time),1), concept_accession(MID,'KEGG',MName,true), write(MName),write(' ('),write(Change),write(') and '), outputRegulatedMets(EIs). extractSharedReactions(Score-EIs-CombineExplanation,Shared):- maplist(extractSharedReaction,CombineExplanation,Shared0), flatten(Shared0,Shared). extractSharedReaction(startNode-EITIs,[]):-!. extractSharedReaction(endNode-EITIs-EITIs1,[]):-!. extractSharedReaction(ClaID-EITIs,ClaID-EIs):- length(EITIs,N),N>1,!,exExtraction(EITIs,EIs). extractSharedReaction(X,[]). /* outputEachH(EI-TI):- ex(EI,concentration(MID,Change,Time),1), concept_accession(MID,'KEGG',MName,true), candidateH0(EI,TI,Tr), Tr=[[rs-{RID,LimitingType,Hstate,Time}]|Rest], concept_accession(RID,'KEGG',RName,true), portray_clause(control_points(MName,RName,LimitingType,Hstate,Time)). */ outputEachH(EI-TI):- candidateH0(EI,TI,Tr), %write({EI-TI,Tr}),nl, ex(EI,concentration(MID,Change,Time),1), concept_accession(MID,'KEGG',MName,true), nl,write('#Explanation for '),write(MName),write(' :'),nl, outputEachH0(Tr), write(MName),write('|'),write(Change). %,write(')'),nl. /* write('Kegg Draw Input: '),nl, maplist(colour0,Tr), concept_accession(MID,'KEGG',MName0,true),atom_chars(MName0,[A,B,C,D|Code]),atom_chars(MCode,Code),write(MCode),write(' '),colourScheme(Change,Colour),write(Colour),nl,nl. */ outputEachH0(Tr):- Tr=[[rs-{RID,enzymeLimiting,Hstate,Time}]], geneIDs_for_oneReaction(RID,GeneIDs),write(GeneIDs),%print_list_oneLine(GeneIDs), write('|'), %concept_name(RID,RName,true),write(RName),write('('), %concept_accession(RID,'KEGG',KeggRID0,true),atom_chars(KeggRID0,[A,B,C|Code]),atom_chars(KeggRID,Code),write(KeggRID),%write(') '), concept_accession(RID,'KEGG',RName,true),write(RName), outputPathway(RID), %concept_name(RID,RName,true),write(RName), write('|'),write(Hstate),write('|'),write('EC'), (catalyzed_by_ECclass(RID,EC_Number) -> write(EC_Number),write('|'), (concentration_e(EC_Number,GEChange,Degree,Time)-> write(GEChange); write(?) ); write('?') ), (produced_by(MID,RID,_) -> write('|->|'); write('|<-|') ). outputEachH0(Tr):- Tr=[[rs-{RID,enzymeLimiting,Hstate,Time}],[met-{MID,Change,Time}]|Tr_rest], geneIDs_for_oneReaction(RID,GeneIDs),write(GeneIDs),%print_list_oneLine(GeneIDs), write('|'), %concept_name(RID,RName,true),write(RName),write('('), %concept_accession(RID,'KEGG',KeggRID0,true),atom_chars(KeggRID0,[A,B,C|Code]),atom_chars(KeggRID,Code),write(KeggRID),%write(') '), %concept_name(RID,RName,true),write(RName), concept_accession(RID,'KEGG',RName,true),write(RName),outputPathway(RID), write('|'),write(Hstate),write('|'),write('EC'), (catalyzed_by_ECclass(RID,EC_Number) -> write(EC_Number),write('|'), (concentration_e(EC_Number,GEChange,Degree,Time)-> write(GEChange); write(?) ); write('0.0.0.0') ), (produced_by(MID,RID,_) -> write('|->|'); write('|<-|'), consumed_by(MID,RID,_) ), concept_accession(MID,'KEGG',MName,true),write(MName),write('|'),write(Change),%write(')'), maplist(outputEachCla,Tr_rest). outputEachH0(Tr):- Tr=[[rs-{RID,substrateLimiting,Hstate,Time}]|Tr_rest], consumed_by(SID,RID,_),ex(EI,concentration(SID,Change,Time),1), concept_accession(SID,'KEGG',SName,true),write(SName),write('|'),write(Change),write('|->|'), concept_accession(RID,'KEGG',RName,true),write(RName),outputPathway(RID), write('|'),write(Hstate),write('|'),write('EC'), (catalyzed_by_ECclass(RID,EC_Number) -> write(EC_Number),write('|'), (concentration_e(EC_Number,GEChange,Degree,Time)-> write(GEChange); write(?) ); write('0.0.0.0') ), write('|->|'), maplist(outputEachCla,Tr_rest). outputPathway(RID):- member_of(RID,PathwayID,_),concept_accession(PathwayID,'KEGG',PathwayName,true),!, write('|'),write(PathwayName). outputPathway(RID):- write('|PATHWAY?'). outputEachCla([rs-{RID,LimitingType,Hstate,Time}]):- %concept_name(RID,RName,true),write(RName),write('('), write('|->|'), %concept_accession(RID,'KEGG',KeggRID0,true),atom_chars(KeggRID0,[A,B,C|Code]),atom_chars(KeggRID,Code),write(KeggRID),%write(')'), %concept_name(RID,RName,true), concept_accession(RID,'KEGG',RName,true), write(RName),outputPathway(RID), write('|'),write(Hstate),write('|'),write('EC'), (catalyzed_by_ECclass(RID,EC_Number) -> write(EC_Number),write('|'), (concentration_e(EC_Number,GEChange,Degree,Time)-> write(GEChange); write(?) ); write('0.0.0.0') ), write('|->|'). outputEachCla([met-{MID,Change,Time}]):- concept_accession(MID,'KEGG',MName,true),write(MName),write('|'),write(Change). %write('('),write(Change),write(')'). /********************** This part analyze their module belonging **********************/ mm1:- numbersList(1,66,PosEIs0), UbqitousCompounds=[1,2,12,25,39,47,49,60,61,62,62,63,64,65,66], % -- in fact some of them are not ubqitous ord_subtract(PosEIs0,UbqitousCompounds,PosEIs), Day1_PosEIs=PosEIs, metabolicModule(Day1_PosEIs,day1), nl. %maplist(add(66),PosEIs,Day3_PosEIs),metabolicModule(Day3_PosEIs,day3), %maplist(add(132),PosEIs,Day7_PosEIs), metabolicModule(Day7_PosEIs,day7), %maplist(add(198),PosEIs,Day14_PosEIs), metabolicModule(Day14_PosEIs,day14). add(X,Y,Z):- Z is X+Y. mm3:- numbersList(1,66,PosEIs0), UbqitousCompounds=[1,2,12,25,39,47,49,60,61,62,62,63,64,65,66], % -- in fact some of them are not ubqitous ord_subtract(PosEIs0,UbqitousCompounds,PosEIs), maplist(add(66),PosEIs,Day3_PosEIs),metabolicModule(Day3_PosEIs,day3). mm7:- numbersList(1,66,PosEIs0), UbqitousCompounds=[1,2,12,25,39,47,49,60,61,62,62,63,64,65,66], % -- in fact some of them are not ubqitous ord_subtract(PosEIs0,UbqitousCompounds,PosEIs), maplist(add(132),PosEIs,Day7_PosEIs), metabolicModule(Day7_PosEIs,day7). mm14:- numbersList(1,66,PosEIs0), UbqitousCompounds=[1,2,12,25,39,47,49,60,61,62,62,63,64,65,66], % -- in fact some of them are not ubqitous ord_subtract(PosEIs0,UbqitousCompounds,PosEIs), maplist(add(198),PosEIs,Day14_PosEIs), metabolicModule(Day14_PosEIs,day14). metabolicModule(PosEIs,Day):- % metabolite module %numbersList(1,66,PosEIs0), %UbqitousCompounds=[1,2,12,25,39,47,49,60,61,62,62,63,64,65,66], % -- in fact 47 is not ubqitous %ord_subtract(PosEIs0,UbqitousCompounds,PosEIs), prefixTreeBuilder(PosEIs), consult('candidateH.pl'), tell('pathwayModules.pl'), forAll(candidateH0(EI,TI,T),outputModules(EI,T)), told, consult('pathwayModules.pl'), forAll(member(EI,PosEIs),(setof(EIPathwayName,withinPathway(EI,PathwayName),EIPathways),length(EIPathways,EIPathwayNums),portray_clause(oneMetabolite_PathwayNums(EI,EIPathwayNums)))), setof(PathwayName,EI^withinPathway(EI,PathwayName),AllPathways), tell('metabolitePathwayModule.pl'), forAll(member(OneWithinPathway,AllPathways),(setof(EI_any,withinPathway(EI_any,OneWithinPathway),EIs_SamePathway),portray_clause(metabolitePathwayModule(OneWithinPathway,EIs_SamePathway)))), %print_list(Pathways). told, consult('metabolitePathwayModule.pl'), tell('metabolitePathwayModuleDetail.pl'), forAll(metabolitePathwayModule(OneWithinPathway,EIs_SamePathway),seeMetabolitePathwayModule(OneWithinPathway,EIs_SamePathway)), told, consult('metabolitePathwayModuleDetail.pl'), findall(OneWithinPathway, (consistentModule(OneWithinPathway,[up-UpList,no_change-NoChangeList,down-DownList]), portray_clause(consistentModule(OneWithinPathway,[up-UpList,no_change-NoChangeList,down-DownList],Day))), InterestingPathways). %forAll(candidateH0(EI_any,TI_any,T),--score and sort. -- count-> how many agree consistentModule(OneWithinPathway,[up-UpList,no_change-NoChangeList,down-DownList]):- seeMetabolitePathwayModule0(OneWithinPathway,[up-UpList,no_change-NoChangeList,down-DownList]), ((UpList==[],NoChangeList==[]);(UpList==[],DownList==[]);(NoChangeList==[],DownList==[])), metabolitePathwayModule(OneWithinPathway,EIs_SamePathway), length(EIs_SamePathway,EINum),EINum>1. % for early time point, there is more agree than non-agree seeMetabolitePathwayModule(OneWithinPathway,EIs_SamePathway):- findall(EI-M,(member(EI,EIs_SamePathway),ex0(EI,concentration(M,up,_,Time),1)),UpList), findall(EI-M,(member(EI,EIs_SamePathway),ex0(EI,concentration(M,no_change,_,Time),1)),NoChangeList), findall(EI-M,(member(EI,EIs_SamePathway),ex0(EI,concentration(M,down,_,Time),1)),DownList), portray_clause(seeMetabolitePathwayModule0(OneWithinPathway,[up-UpList,no_change-NoChangeList,down-DownList])). outputModules(EI,T):- maplist(oneReactionModule(EI),T). oneReactionModule(EI,[rs-{ReactionID,LimitingType,HypothesizedState,Time}]):- findall(PathwayID, (member_of(ReactionID,PathwayID,_), once((pathway(PathwayID,_,_,'KEGG','IMPD'),concept_accession(PathwayID,'KEGG',PathwayLycocycsID,true),concept_name(PathwayID,PathwayName,true))), %concept_accession(PathwayID,'KEGG',PathwayLycocycsID,true), %concept_name(PathwayID,PathwayName,true)), portray_clause(withinPathway(EI,PathwayLycocycsID-PathwayName))), %PathwayID))), % PathwayIDs). mEC:- findall(RID,reaction(RID2,_,_,'KEGG','IMPD'),RIDs), length(RIDs,TotalReactionNum),write('total number of reactions is '), write(TotalReactionNum),nl, findall(RID2, (reaction(RID2,_,_,'KEGG','IMPD'), %part_of_catalyzing_class(RID2,EC_Number1,'IMPD'), % catalyzed_by_ECclass(RID2,EC_Number1), %part_of_catalyzing_class(RID2,EC_Number2,'IMPD'), % catalyzed_by_ECclass(RID2,EC_Number2), EC_Number1\==EC_Number2 ), RIDs2), length(RIDs2,MultiECReactionNum),write('total number of reactions is '), write(MultiECReactionNum). /************************************* OLD **********************************************/ pathDetail([X],[X]). pathDetail([S,P|Path],[S,'<-',K,'->'|Record]):- pathDetail([P|Path],Record), pair(S,P,Reactions,K). pair(M1,M2,ReactionIDs,K):- findall(ReactionID, ( consumed_by(M1,ReactionID,'IMPD'), produced_by(M2,ReactionID,'IMPD') ), ReactionIDs ), length(ReactionIDs,K). /*% no need for now -- the data is directly about enzyme, rather than EC number % link reactions to their EC_Number catalyzed_by_ECclass(ReactionID,EC_Number):- setof(EC_Num, EnzID^(catalyzed_by(ReactionID,EnzID,'IMPD'), part_of_catalyzing_class(EnzID,EnzClassID,'IMPD'), concept_accession(EnzClassID,'EC',EC_Num,true)), EC_Nums), member(EC_Number,EC_Nums). */ catalyzed_by_ECclass(ReactionID,EC_Num):- catalyzed_by(ReactionID,EnzClassID,'IMPD'), %part_of_catalyzing_class(EnzID,EnzClassID,'IMPD'),!, concept_accession(EnzClassID,'EC',EC_Num,false). % concept_accession(EnzClassID,'KEGG',EC_Num,true). /*catalyzed_by_ECclass(ReactionID,EC_Num):- catalyzed_by(ReactionID,EnzID,'IMPD'), part_of_catalyzing_class(EnzID,EnzClassID,'IMPD'),!, concept_accession(EnzClassID,'EC',EC_Num,true). */ /* current Lycocyc file don't have this catalyzed_by_ECclass2(ReactionID,EC_Number):- part_of_catalyzing_class(ReactionID,EC_Number,'IMPD'). ec:- findall(RID, (reaction(ReactionID,'','','KEGG','IMPD'), catalyzed_by_ECclass1(ReactionID,EC_Number1),catalyzed_by_ECclass2(ReactionID,EC_Number2), EC_Number1\==EC_Number2 ), RIDs), write(RIDs). */ %enzyme_state(EC_Number,no_change):- \+enzyme_state(EC_Number,inhibited). /*writeReaction(ReactionID):- catalyzed_by_ECclass(ReactionID,EC_Number), findall(SID,consumed_by(SID,ReactionID,'IMPD'),SIDs), findall(PID,produced_by(PID,ReactionID,'IMPD'),PIDs), nl,write(EC_Number-ReactionID-{SIDs,PIDs}),nl. */ %for each example, with the other example exPath(Ex1,Ex2,Path):- ex(Ex1,concentration(M1,Change1),1), ex(Ex2,concentration(M2,Change2),1), path(M1,M2,[Ex1],Path). % detailPath(Path) -- output together. Node, enzyme on each node /* path(M,M,Path,Path). path(M1,M2,PrePath,Path):- pair(M1,M), \+member(M,PrePath), % avoid the loop append(PrePath,[M],NewPath), path(M,M2,NewPath,Path). */ path(M,M,Path,Path). path(M1,M2,PrePath,Path):- pair(M1,M), \+member(M,PrePath), % avoid the loop (ex(EI,concentration(M,Change),1)-> Middle=EI; Middle=M), path(M,M2,[Middle|PrePath],Path). % pair('41065: UDP-GLUCOSE','38003: SUCROSE-6P',RID,EC). pair(M1,M2,ReactionID,EC_Number):- consumed_by(M1,ReactionID,'IMPD'), produced_by(M2,ReactionID,'IMPD'), catalyzed_by_ECclass(ReactionID,EC_Number). pair(M1,M2):- consumed_by(M2,ReactionID,'IMPD'), produced_by(M1,ReactionID,'IMPD'). writePairDetail(Ex1,Ex2):- ex(Ex1,concentration(M1,Change1),1), ex(Ex2,concentration(M2,Change2),1), findall(ReactionID-EC_Nums, (consumed_by(M1,ReactionID,'IMPD'), produced_by(M2,ReactionID,'IMPD'), writeReaction(ReactionID), catalyzed_by_ECclass(ReactionID,EC_Nums)), AllDetail), print_list(AllDetail). %{EC_Number,ReactionID}):- print_list([]). print_list([X|List]):- portray_clause(X), %writeCla(X),nl, print_list(List). %:- depth_bound_call(user:path('27616: ILE','27565: GLT',[],Path),3),write(Path). t:- tell('linking.txt'),%,6,7,8,9,10 List=[1,2,3,4,5], findall(Ex1-Ex2-Path, (member(Ex1,List),member(Ex2,List),Ex1\==Ex2, depth_bound_call(user:exPath(Ex1,Ex2,Path),6)), Paths0 ), remove_duplicates(Paths0,Paths), sort(Paths,SortedPath), print_list(SortedPath), told. exclude(EIs):- findall(M_name, (ex0(EI,concentration(M_name,Change,_,Day),1), concept_name(MID,M_name,true), \+produced_by(MID,RID1,_), \+consumed_by(MID,RID2,_) ), EIs), write(EIs). %:- depth_bound_call(user:exPath(2,4,Path),5),write(Path).