Protein Folding (mpex28)
/***************************************************************/
/* */
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: mpex28 */
/* TITLE: Protein Folding (mpex28) */
/* PRODUCT: OR */
/* SYSTEM: ALL */
/* KEYS: OR */
/* PROCS: OPTMODEL, SGPLOT */
/* DATA: */
/* */
/* SUPPORT: UPDATE: */
/* REF: */
/* MISC: Example 28 from the Mathematical Programming */
/* Examples book. */
/* */
/***************************************************************/
data hydrophobic_data;
input position @@;
datalines;
2 4 5 6 11 12 17 20 21 25 27 28 30 31 33 37 44 46
;
%let num_acids = 50;
proc optmodel;
num n = &num_acids;
set POSITIONS = 1..n;
set HYDROPHOBIC;
read data hydrophobic_data into HYDROPHOBIC=[position];
set PAIRS
= {i in HYDROPHOBIC, j in HYDROPHOBIC: i + 1 < j and mod(j-i-1,2) = 0};
/* IsFold[i] = 1 if fold occurs between positions i and i + 1; 0 otherwise */
var IsFold {1..n-1} binary;
/* IsMatch[i,j] = 1 if hydrophobic pair at positions i and j are matched;
0 otherwise */
var IsMatch {PAIRS} binary;
/* maximize number of matches */
max NumMatches = sum {<i,j> in PAIRS} IsMatch[i,j];
/* if IsMatch[i,j] = 1 then IsFold[k] = 0 */
con DoNotFold {<i,j> in PAIRS, k in i..j-1 diff {(i+j-1)/2}}:
IsMatch[i,j] + IsFold[k] <= 1;
/* if IsMatch[i,j] = 1 then IsFold[k] = 1 */
con FoldHalfwayBetween {<i,j> in PAIRS}:
IsMatch[i,j] <= IsFold[(i+j-1)/2];
solve;
print {i in 1..n-1: IsFold[i].sol > 0.5} IsFold;
print {<i,j> in PAIRS: IsMatch[i,j].sol > 0.5} IsMatch;
num x {POSITIONS};
num y {POSITIONS};
num xx init 0;
num yy init 0;
num dir init 1;
for {i in POSITIONS} do;
xx = xx + dir;
x[i] = xx;
y[i] = yy;
if i = n or IsFold[i].sol > 0.5 then do;
xx = xx + dir;
dir = -dir;
yy = yy - 1;
end;
end;
create data plot_data from [i] x y is_hydrophobic=(i in HYDROPHOBIC);
create data edge_data from [i]=(1..n-1)
x1=x[i] y1=y[i] x2=x[i+1] y2=y[i+1] linepattern=1;
create data match_data from [i j]={<i,j> in PAIRS: IsMatch[i,j].sol > 0.5}
x1=x[i] y1=y[i] x2=x[j] y2=y[j] linepattern=2;
quit;
data sganno(drop=i j);
retain drawspace "datavalue" linethickness 1;
set edge_data match_data;
function = 'line';
run;
proc sgplot data=plot_data sganno=sganno;
scatter x=x y=y / group=is_hydrophobic datalabel=i;
xaxis display=none;
yaxis display=none;
run;