%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Agostino Dovier, August 30, 2012
%% TESTED with BProlog
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Input = [ a, g, c, u, a, ...]
rna_align(Input,Pairing,Energy) :-
length(Input,N),
length(Pairing,N),
rnadomain(Pairing,1,Input),
constraint(Input,Pairing),
energy(Pairing,Input,Energy, N),
statistics(runtime,_),
labeling([ff,minimize(Energy)],Pairing),
statistics(runtime,[_,T]),
format("Running time: ~d ms~n",T),
myprint(Input,Pairing).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% SET THE DOMAIN, ACCODING TO THE NUCLEOTIDES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rnadomain([],_,_).
rnadomain([A|R],I,Input) :- %% 0 stands for no pair
nth1(I,Input,Base),
(Base==a ->
findall(J,(nth1(J,Input,u),abs(I-J)>1),DOMLIST );
Base==c ->
findall(J,(nth1(J,Input,g),abs(I-J)>1),DOMLIST );
Base==g ->
findall(J,(or(J,Input,c,u),abs(I-J)>1),DOMLIST );
Base==u ->
findall(J,(or(J,Input,a,g),abs(I-J)>1),DOMLIST );
true -> (write('wrong input!'),nl,fail)),
A in [0|DOMLIST],
I1 is I+1,
rnadomain(R,I1,Input).
or(V,LIST,A,_) :- nth1(V,LIST,A).
or(V,LIST,_,B) :- nth1(V,LIST,B).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% CONSTRAINTS (main procedure)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
constraint(_,Pairing) :-
different(Pairing),
symmetric(Pairing),
%contiguity(Pairing),
no_knots(Pairing),
write('Constraints added. Starting search.'),nl.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Injectivity of contacts
%%%% (it is not all_different due to 0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
different([]).
different([_]).
different([A,B|R]) :-
different_aux(A,[B|R]),
different([B|R]).
different_aux(_A,[]).
different_aux(A,[B|R]) :-
A #> 0 #=> A #\= B,
different_aux(A,R).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Symmetry of contacts
%%%% If Pairing[i]=j > 0 Then Pairing[j]=i :
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Expressed using "element"
%%% Assume Pairing[i] = 0 (P below)
%%% Then element(1,[i,...],i) is true
%%% Assume Pairing[i] > 0 (e.g. 2 = P)
%%% Then element(3,[i,_,Pairing[2],_, ...],i) is true iff Pairing[2]=i
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
symmetric(Pairing) :-
symmetric(Pairing,Pairing,1).
symmetric([],_,_).
symmetric([P|Airing],Pairing,I) :-
P1 #= P+1, element(P1,[I|Pairing],I), %%% <-Symmetry
I1 is I + 1,
symmetric(Airing,Pairing,I1).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Contiguity (partial) of contacts
%%%% If Pairing[i]=j > 0
%%%% Then Pairing[i+1]=0 or it is j+1 or j-1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Observe that it is symmetrical...
%%% Difficult to be stated with the monotonic encoding
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
contiguity([]).
contiguity([_]).
contiguity([A,B|Pairing]) :-
A #>0 #=> (B #= A-1 #\/ B #= A+1 #\/ B #= 0),
contiguity([B|Pairing]).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% no_pseudo_knots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
no_knots(Pairings) :-
no_knots(Pairings,1,Pairings).
no_knots(L,_,_) :- length(L,NL),NL<3,!.
no_knots([End|Pairing],Init,Pairs) :-
%% I only consider increasing pairs (not the symmetrical)
truthlist(Init,End,Pairs,Truth),
(End #> Init) #=> Truth #= 1,
I1 is Init + 1,
no_knots(Pairing,I1, Pairs).
truthlist(Init,End,Pairs,Flag) :-
myselect(Init,Pairs,InterestingPart),
J is Init+2,
loopcontrol(InterestingPart,J,Init,End,Truthlist),
Flag #= min(Truthlist).
loopcontrol([],_,_,_,[]).
loopcontrol([P|Rest],J,Init,End,[T|Truth]) :-
%%% (we know that (Init,End) in P), with Init < End
%%% and J > Init.
%%% We do not consider here backwards jumps
P #< J #=> T #= 1,
(P #>= J #/\ P #< End) #=> T #= 1,
(P #>= J #/\ P #>= End) #=> T #= 0,
J1 is J+1,
loopcontrol(Rest, J1, Init, End, Truth).
myselect(1,[_,_|Pairs],Pairs).
myselect(N,[_|Pairs],InterestingPart) :-
N>1,
N1 is N-1,
myselect(N1,Pairs,InterestingPart).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% energy function (used a statistical one, as
%%% done by Dahl/Bavarian )
%%% AU~35%, CG~53%, UG~12%
%%% NIL = number of zeros
%%% Energy ~ k1 | AU / (N-NIL) - 0.35 | +
%%% k2 | CG / (N-NIL) - 0.53 | +
%%% k3 NIL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
energy(Pairing,Input,Energy, N) :-
findall(J,nth1(J,Input,a),Alist),
findall(J,nth1(J,Input,c),Clist),
findall(J,nth1(J,Input,g),Glist),
findall(J,nth1(J,Input,u),Ulist),
%%%%%%%
cons_energy(Pairing,Input,Pairing, 1, 0,AU, 0,CG, Alist,Clist,Glist,Ulist),
count(0,Pairing,#=,NonContacts),
%%%% PARAMETERS
C1=400,
C2=5,
C3=3,
%%%%
Energy #= C1*NonContacts +
C2*abs(100*AU - 70*N + 70*NonContacts) +
C3*abs(100*CG - 106*N + 106*NonContacts).
cons_energy([],_,_,_, AU,AU,CG,CG,_,_,_,_).
cons_energy([P|Airing],Input,Ass, I, AUo,AUn,CGo,CGn,
Aset,Cset,Gset,Uset) :-
% write('#######'),write(I),
nth1(I,Input,Base),
% write('#######'),write(Base),nl,
(Base==a ->
(P #= 0 #=> AUi #= AUo,
P #> 0 #=> AUi #= AUo+1,
CGo #= CGi);
Base==c ->
(P #= 0 #=> CGi #= CGo,
P #> 0 #=> CGi #= CGo+1,
AUo #= AUi);
Base==g ->
(P #= 0 #=> CGi #= CGo,
P in Cset #=> CGi #= CGo+1,
P in Uset #=> CGi #= CGo,
AUo #= AUi);
Base==u ->
(P #= 0 #=> AUi #= AUo,
P in Aset #=> AUi #= AUo+1,
P in Gset #=> AUi #= AUo,
CGo #= CGi)),
I1 is I+1,
cons_energy(Airing,Input,Ass, I1, AUi,AUn,CGi,CGn,
Aset,Cset,Gset,Uset).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Printing predicate
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
myprint(Input,Pairing) :-
myprint(1,Input,Pairing,Input).
myprint(_,[],[],_) :- nl.
myprint(N,[I|Nput],[0|Pairing], Bases) :-
!,
write(N),write(':'),write(I),nl,
N1 is N+1,
myprint(N1,Nput,Pairing,Bases).
myprint(N,[I|Nput],[J|Pairing], Bases) :-
nth1(J,Bases,B),
write(N),write(':'),write(I),write('--'),
write(B),write(':'),write(J),nl,
N1 is N+1,
myprint(N1,Nput,Pairing,Bases).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% call as
%%% :- rna(0,L),rna_align(L,Pairing,Energy).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rna(0,[g,c,g,a,a,a,c,g,c]).
rna(1,[u,c,g,a,a,a,c,g,c,c,c,a]).
rna(2,[a,c,g,a,a,a,u,c,g,a,a,a,c,g,c,c,c,a,u,u,u,u,g,u]).
%%% prion
rna(3,[c,c,a,a,g,a,u,g,u,g,g,a,g,g,c,u,g,g,g,g,u,c,a,g,c,c,c,u,u,g,g,u,g,g,a,g,g,u,u,g,c]).
%%% top riboswitch:
rna(4,[g,a,g,c,a,a,g,g,g,g,u,u,c,u,a,g,u,g,c,a,g,c,u,c,g,a,c,u,a,g,a,a,c,c,c,u,g,c,a,a,u,g,g,c,a,c,g,g,a,c,a,u,c,g,a,u,a,a,g,u,g,u,g,u,u,c,g,u,g,c,u,u,u,u,g,a,c,a,u]).