// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
// ! This script is provided "as is" by Applied Maths.              !
// ! You are free to use and modify this script for your own needs. !
// ! Redistribution or reproduction of the script is prohibited.    !
// ! DISCLAIMER:                                                    !
// ! Improper use of scripts may corrupt your database.             !
// ! Running this script is entirely at your own responsibility.    !
// ! Applied Maths accepts no lialibility for any consequences      !
// ! resulting from its use.                                        !
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

//$DISPLAY_COMPWIN

integer i,j,totlen,linedrawn,seqcount;
string seq,experlist,expername,ch;
integer nA[],nC[],nG[],nT[],n[];
float conserv[],xp,yp,xsep;
DIALOG dlg;


if not(CmpIsPresent) then CmpAttach;
if not(CmpIsPresent) then {
   message("Unable to attach a comparison to this script");
   stop;
}

for i=1 to DbGetExperCount do if DbGetExperClass(i)="SEQ" then {
   if expername="" then expername=DbGetExpername(i);
   experlist=experlist+DbGetExperName(i)+"	";
}
if expername="" then { message("No sequence type experiments found!"); stop; }
DlgAddList(dlg,experlist,expername,15,15,150,150,"LIST");

if not(DlgShow(dlg,"Sequence diversity",300,200)) then stop;

seqcount=CmpGetEntryCount;
setbusy("Calculating statistics");
for i=1 to seqcount do {
   if i/5=floor(i/5) then setbusy("Calculating statistics "+str(i,0,0)+" / "+str(CmpGetEntryCount,0,0));
   seq=CmpSeqGet(expername,i);
   if length(seq)>totlen then totlen=length(seq);
   for j=1 to length(seq) do {
      ch=substring(seq,j,j);
      if ch="A" then { nA[j]=nA[j]+1; n[j]=n[j]+1; }
      if ch="C" then { nC[j]=nC[j]+1; n[j]=n[j]+1; }
      if ch="G" then { nG[j]=nG[j]+1; n[j]=n[j]+1; }
      if ch="T" then { nT[j]=nT[j]+1; n[j]=n[j]+1; }
   }
}
setbusy("");

for i=1 to totlen do {
   if n[i]>0 then {
      conserv[i]=((nA[i]*nA[i]+nC[i]*nC[i]+nG[i]*nG[i]+nT[i]*nT[i])/(n[i]*n[i])-0.25)/0.75;
   }
}



PrnStartScreen(18,24);
xsep=16*7/totlen*0.9;
if xsep>0.3 then xsep=0.3;
xp=1;yp=1;linedrawn=0;
for i=1 to totlen do {
   if not(linedrawn) then {
      PrnLinestyle(0.05,0,0,0);
      PrnMoveto(1,yp+3-2.4-0.025);PrnLineto(17+xsep,yp+3-2.4-0.025);
      PrnMoveto(1,yp+3+0.025);PrnLineto(17+xsep,yp+3+0.025);
      linedrawn=1;
   }
   PrnFillRect(xp,yp+3,xp+xsep*1.05,yp+3-conserv[i]*2.4,60,60,60);
   PrnLinestyle(0.01,0,0,0);
   PrnCircle(xp+xsep/2,yp+3-n[i]/seqcount*2.38,xsep/3,0,0,0);
   xp=xp+xsep;
   if xp>=17 then { xp=1; yp=yp+3; linedrawn=0; }
}
PrnEnd;