program seq5c;
{to show g-c content of sequences and restriction site locations.
Sequence must be in GENOME directory with '.txt' extension.
Use PARSGENE program to convert continuous string data to lined file
For restriction enzyme list, use file enzymes2.txt}

uses crt,graph,dos;
const maxWindow = 500;
var i,x,y,xMax,yMax,lineIndex,windowIndex,windowsize,searchLength,eNumber:integer;
    baseIndex,startIndex,displayWidth:longInt;
    cgContent:byte;
    dataline,base,reply,message,searchStringF,searchStringR,snert,eName:string;
    s,fileName:pathStr;
    dird:dirStr; name:nameStr; ext:extStr;
    seqdata, outFile, enzymes :text;
    fileopen,endOfDataline,baseFound,recentBitmap:boolean;
    character,ch:char;
    plotdata:text;
    windowEntry:array[1..maxWindow] of byte;

{$I Bitmap.pas}

procedure openfile(var filename:pathstr; var fileopen:boolean);
{open an existing file for input}
begin
     writeln;
     writeln('Which genome do you want to look at?');
     readln(reply);
     message:='GENOME\'+reply+'.txt';
     filename:=fExpand(message);
     fSplit(fileName,dird,name,ext);
     s:=fSearch(filename,dird);
     if s<>'' then fileopen:=true
     else
       begin
       writeln('Unfortunately, file name '+message+' does not exist');
       fileopen:=false;
       end;
end;

procedure initWindow;
var i:integer;
begin
for i:=1 to windowSize do if odd(i) then windowentry[i]:=1
  else windowEntry[i]:=0;
end;

procedure  getBaseType;
begin
if endOfDataLine then
   begin readln(seqData,dataLine); lineIndex :=0; endOfDataLine:=false; end;
inc(lineIndex);
baseFound:=false;
base:=copy(dataLine,lineIndex,1);
if (base = 'a') or (base='A') then begin cgContent:=0;baseFound:=true; end;
if (base = 't') or (base='T') then begin cgContent:=0;baseFound:=true; end;
if (base = 'c') or (base='C') then begin cgContent:=1;baseFound:=true; end;
if (base = 'g') or (base='G') then begin cgContent:=1;baseFound:=true; end;
if lineIndex >= 80 then endOfDataLine:=true;
end;

procedure placeInWindow;
begin
inc(baseIndex);
inc(windowIndex);
if windowIndex>windowsize then windowIndex:=0;
windowEntry[windowIndex]:=cgContent;
end;

procedure plotGCcontent;
var i,cgSum:integer;
begin
cgSum:=0;
for i:=1 to windowSize do cgSum:=cgSum+windowEntry[i];
y:=round(yMax *(1- cgSum/windowSize));
x:=round(xMax *((baseIndex-startIndex)/displayWidth));
setColor(red);    if baseindex<2 then setColor(black);
if (x>=0) and (x<=xMax) and (y>=0) and (y<=yMax) then lineTo(x,y);
end;

procedure searchForSequence;
begin
snert:=snert+base;
if length(snert)>searchLength then snert:=copy(snert,2,searchLength);
if snert=searchStringF then
  begin setColor(cyan);
  setTextStyle(0,1,1); {write sideways}
  line(x,yMax div 2-100,x,yMax div 2);
  str(baseIndex,message);
  outTextXY(x,50,message);
  moveTo(x,y);
  setTextStyle(0,0,1); {text back to normal}
  end;
if snert=searchStringR then
  begin setColor(yellow);
  setTextStyle(0,1,1); {write sideways}
  line(x,yMax div 2,x,yMax div 2+100);
  str(baseIndex,message);
  outTextXY(x,yMax-50,message);
  moveTo(x,y);
  setTextStyle(0,0,1); {text back to normal}
  end;
end;


procedure showMenu;
var i:integer;
begin
for i:=0 to 10 do
  begin
  if i=5 then setColor(white) else setColor(blue);
  line(0,i*yMax div 10,xMax,i*yMax div 10);
  end;
setColor(white);
outTextXY(10,yMax-10,'Z/z..zoom  >/<..right/left  W/w..window   s..search   E/e..enzyme   q..quit');
str(windowSize,message);
outTextXY(10,10,'windowSize = '+message);
setTextStyle(0,1,1); {write sideways}
str(startIndex,message);
outTextXY(10,yMax div 2,message);
str(startIndex+displayWidth,message);
outTextXY(xMax,yMax div 2,message);
setTextStyle(0,0,1); {back to normal}
end;

function reverseCompliment (inputString:string):string;
{to search for sequences that would be hit going reverse direction
on complimentary strand}
var i:integer; base,temp:string; strLength:byte;
begin
temp:='';
strLength := length(inputString);
for i:=strLength downTo 1 do
  begin
  base:=copy(searchStringF,i,1);
  if base = 'a' then temp:=temp+'t';
  if base = 'c' then temp:=temp+'g';
  if base = 'g' then temp:=temp+'c';
  if base = 't' then temp:=temp+'a';
  if base = 'A' then temp:=temp+'T';
  if base = 'C' then temp:=temp+'G';
  if base = 'G' then temp:=temp+'C';
  if base = 'T' then temp:=temp+'A';
  end;
reverseCompliment:=temp
end;

begin
  clrscr;
  repeat openfile(filename,fileopen)
  until (fileopen or (reply='q'));
  assign(seqdata,filename);
  reset(seqdata);
  readln(seqData,displayWidth);  {reads total bases and sets to display all}
  startIndex:=0; {starting at the beginning}
  assign(enzymes,'enzymes2.txt');
  x:=detect;
  initgraph(x,y,'');
  xMax:=getMaxX; yMax:=getMaxY;
  windowSize:=60;
  repeat
    initWindow;
    reset(seqdata);
    readln(seqData);  {skip base total} endOfDataline:=true;
    clearDevice;
    if eName<>'' then outTextXY(400,100,eName+'   '+searchStringF);
    if recentBitmap then outTextXY(100,100,'Image saved as '+reply+'.bmp');
    recentBitmap:=false;
    baseIndex:=0;
      repeat
      getBaseType;
      if baseFound then placeInWindow;
      plotGCcontent;
      if (searchLength>0) and baseFound then searchForSequence;
      until (eof(seqData) and endOfDataline) or keypressed;
    showMenu;
    ch:=readkey;
    case ch of
      'Z': displayWidth := round(displayWidth*0.8);
      'z': displayWidth := round(displayWidth*1.2);
      '<': startIndex:=startIndex+round(displayWidth*0.2);
      '>': startIndex:=startIndex-round(displayWidth*0.2);
      'W': begin windowSize := round(windowSize*1.5);
           if windowSize>maxWindow then windowSize:=maxWindow; end;
      'w': begin windowSize := round(windowSize/1.5);
           if windowSize<1 then windowSize:=1; end;
      's':begin outTextXY(100,100,'input search sequence [atattcg...] from 5-prime toward 3-prime');
          readln(searchStringF); searchLength:=length(searchStringF);
          searchStringR := reverseCompliment(searchStringF);
          end;
      'e','E':begin
           if ch='e' then inc(eNumber) else dec(eNumber);
           if eNumber<0 then eNumber :=0;
           reset(enzymes);
           eName:='';    searchStringF :='';
           for i:=1 to eNumber do
               begin
               readln(enzymes,eName);
               readln(enzymes,searchStringF);
               end;
           if EOF(enzymes) then dec(eNumber);
           searchLength:=length(searchStringF);
           searchStringR := reverseCompliment(searchStringF);
           end;
       'b': begin
            assign(outFile,reply+'.bmp');
            rewrite(outFile);
            writeBitmap;
            recentBitmap:=true;
            end;

      chr(0):  {also arrow keys for pan right/left}
        begin
        ch:=readkey;
        if ord(ch)=75 then startIndex:=startIndex+round(displayWidth*0.2);
        if ord(ch)=77 then startIndex:=startIndex-round(displayWidth*0.2);
        end;
      end; {of case}
  until ch='q'  ;
  close(outFile);
end.
