Program CnvacVRI (input,output); (*Prg reads std txt files,LINE BY LINE, and converts it into RIBvH values which are coded on 10 levels for input to Neural wks*) type filename=string[12]; prob=array[1..650] of real; dprob=array[1..650] of real; table=array[1..101,1..2] of real; var x1,ans,o :char; z,i,j,k,l,sn,error,m,filnum :integer; brg,brc,fnbrg,fnbrc,LTG,fngbr, fncbr,fnabr,fntbr, bra,brt, fnbra,fnbrt, abr,cbr,gbr,tbr,cf, a1,a2,a3,a4,a5,a6,a7,a8,a9, a10,a11,rtot :real; rseq,rseqdb :prob; lt :table; Msg :string; length,ptn ,w,x,y :integer; ifil1,ofil,ifil2 :text; namei :string; freqar :array[1..4,1..58] of real; (*SIS*) Topval :array[1..58] of real; (*SIS*) mbvh :array[1..58] of real; (*SIS*) DBA :array[1..76,1..58] of char; (*sens to input size*) SA :string[58]; (*SIS*) DBS :string[58]; (*SIS*) (*******************************************************) procedure questions; begin writeln('Enter the name of the file of test sequences'); readln(namei); Writeln('How many sequences are in this file?'); readln(filnum); writeln('Outputs to produce 0 or 1 output? Must code true and false sets separately'); readln(o); assign(ifil1,namei); reset(ifil1); writeln('Enter the name of the file containing dbase 'prototype' sequences'); readln(namei); assign(ifil2,namei); reset(ifil2); writeln('What is the sequence length?'); readln(length); writeln('How many sequences are in dbase file?'); readln(ptn); assign(ofil,'outcode.nna'); rewrite(ofil); writeln('Use default base ratios?'); readln(ans); ans:=UpCase(ans); LTG:=2; if ans='N' then begin writeln('Enter non-std local base ratios, A C G T as decimal values.'); readln(abr, cbr,gbr,tbr); LTG:=-(gbr*1.4427*ln(gbr)+cbr*1.4427*ln(cbr)+abr*1.4427*ln(abr)+ tbr*1.4427*ln(tbr)); end else begin abr:=0; cbr:=0; gbr:=0; tbr:=0; end; end; (***************************************************************) procedure initialize; begin k:=0; w:=0; x:=0; y:=0; z:=0; i:=1; j:=0; l:=0; rtot:=0; fnbrg:=0; fnbrc:=0; fnbra:=0; fnbrt:=0; brg:=0; brc:=0; bra:=0; brt:=0; error:=0; for l:=1 to 650 do begin rseq[l]:=0; rseqdb[l]:=0; end; l:=0; a1:=0.00100; a2:=0.0100; a3:=0.10; a4:=0.30; a5:=0.60; a6:=1.100; a7:=2.200; a8:=3.300; a9:=4.400; a10:=5.500; (****************************************************************) (*LOOKUP TABLE FOR SAMPLE SIZE ERROR CORRECTION*) (****************************************************************) LT[1,1]:=0.0; LT[1,2]:=0.0; LT[2,1]:=0.750; LT[2,2]:=0.187; LT[3,1]:=1.111;(*FROM SCHNEIDER ET AL JMB 188 415 86*) LT[3,2]:=0.182; LT[4,1]:=1.324; LT[4,2]:=0.152; LT[5,1]:=1.463; LT[5,2]:=0.121; LT[6,1]:=1.559; LT[6,2]:=0.096; LT[7,1]:=1.629; LT[7,2]:=0.077; LT[8,1]:=1.681; LT[8,2]:=0.061; LT[9,1]:=1.722; LT[9,2]:=0.049; LT[10,1]:=1.753; LT[10,2]:=0.040; LT[11,1]:=1.779; LT[11,2]:=0.033; LT[12,1]:=1.800; LT[12,2]:=0.028; LT[13,1]:=1.817; LT[13,2]:=0.023; LT[14,1]:=1.832; LT[14,2]:=0.020; LT[15,1]:=1.844; LT[15,2]:=0.017; LT[16,1]:=1.855; LT[16,2]:=0.015; LT[17,1]:=1.864; LT[17,2]:=0.013; LT[18,1]:=1.872; LT[18,2]:=0.011; LT[19,1]:=1.880; LT[19,2]:=0.010; LT[20,1]:=1.886; LT[20,2]:=0.009; LT[21,1]:=1.892; LT[21,2]:=0.008; LT[22,1]:=1.897; LT[22,2]:=0.007; LT[23,1]:=1.902; LT[23,2]:=0.007; LT[24,1]:=1.906; LT[24,2]:=0.006; LT[25,1]:=1.910; LT[25,2]:=0.006; LT[26,1]:=1.914; LT[26,2]:=0.006; LT[27,1]:=1.917; LT[28,1]:=1.920; LT[29,1]:=1.923; LT[30,1]:=1.926; LT[31,1]:=1.928; LT[32,1]:=1.930; LT[33,1]:=1.932; LT[34,1]:=1.9331; LT[34,2]:=0.005; LT[35,1]:=1.9342; LT[35,2]:=0.005; LT[36,1]:=1.9353; LT[36,2]:=0.005; LT[37,1]:=1.9364; LT[38,1]:=1.9375; LT[39,1]:=1.9386; LT[40,1]:=1.94; LT[41,1]:=1.941; LT[42,1]:=1.9421; LT[43,1]:=1.9432; LT[44,1]:=1.9443; LT[45,1]:=1.9454; LT[46,1]:=1.9465; LT[47,1]:=1.9476; LT[48,1]:=1.9487; LT[49,1]:=1.9498; LT[50,1]:=1.956; LT[51,1]:=1.9564; LT[52,1]:=1.9568; LT[53,1]:=1.9572; LT[54,1]:=1.9576; LT[55,1]:=1.958; LT[56,1]:=1.9584; LT[57,1]:=1.9588; LT[58,1]:=1.9592; LT[59,1]:=1.9596; LT[60,1]:=1.96; LT[61,1]:=1.961; LT[62,1]:=1.962; LT[63,1]:=1.963; LT[64,1]:=1.964; LT[65,1]:=1.965; LT[66,1]:=1.966; LT[67,1]:=1.967; LT[68,1]:=1.968; LT[69,1]:=1.969; LT[70,1]:=1.97; LT[71,1]:=1.9703; LT[72,1]:=1.9706; LT[73,1]:=1.9709; LT[74,1]:=1.9712; LT[75,1]:=1.9715; LT[76,1]:=1.9718; LT[77,1]:=1.9721; LT[78,1]:=1.9724; LT[79,1]:=1.9727; LT[80,1]:=1.973; LT[81,1]:=1.9732; LT[82,1]:=1.9734; LT[83,1]:=1.9736; LT[84,1]:=1.9738; LT[85,1]:=1.9740; LT[86,1]:=1.9742; LT[87,1]:=1.9744; LT[88,1]:=1.9746; LT[89,1]:=1.9748; LT[90,1]:=1.975; LT[91,1]:=1.9752; LT[92,1]:=1.9754; LT[93,1]:=1.9756; LT[94,1]:=1.9758; LT[95,1]:=1.976; LT[96,1]:=1.9762; LT[97,1]:=1.9764; LT[98,1]:=1.9766; LT[99,1]:=1.9768; LT[100,1]:=1.977; LT[27,2]:=0.005; LT[28,2]:=0.005; LT[29,2]:=0.005; LT[30,2]:=0.005; LT[31,2]:=0.004; LT[32,2]:=0.004; LT[33,2]:=0.004; LT[40,2]:=0.004; LT[41,2]:=0.004; LT[50,2]:=0.003; LT[51,2]:=0.003; LT[60,2]:=0.003; LT[61,2]:=0.003; LT[62,2]:=0.003; LT[70,2]:=0.002; LT[71,2]:=0.002; LT[75,2]:=0.002; LT[80,2]:=0.002; LT[81,2]:=0.002; LT[90,2]:=0.001; LT[91,2]:=0.001; LT[100,2]:=0.00; LT[101,2]:=0.00; end; (*******READ db sequences into an array***********************) procedure getdb; begin for j:=2 to (1+ptn) do begin readln(ifil2,DBS); for k:=1 to length do DBA[j,k]:=DBS[k]; end; end; (****************CALC RI for db seq*************************) (***************************************************************) procedure calcRI; (**************************************************************) (*CODE TO READ array BY COLUMNS AND AVERAGE EACH COLUMN FOR 4 BASES*) (**************************************************************) begin for j:=1 to length do begin for i:=1 to ptn do begin x1:=DBA[(i+1),j]; x1:=upcase(x1); if (x1='G') then w:=w+1; if (x1='C') then x:=x+1; if (x1='A') then y:=y+1; if (x1='T') then z:=z+1; if not ((x1='G') or (x1='C') or (x1='A') or (x1='T')or (x1='N')) then begin WRITELN('INPUT ERROR A,T,G,C EXPECTED'); error:=succ(error); readln; end; end;{for i} m:=ptn; brg:=w/(m-error); brc:=x/(m-error); bra:=y/(m-error); brt:=z/(m-error); freqar[1,j]:=brg; freqar[2,j]:=brc; freqar[3,j]:=bra; freqar[4,j]:=brt; topval[j]:=0; for i:=1 to 4 do if freqar[i,j]>topval[j] then topval[j]:=freqar[i,j]; freqar[1,j]:=ln((topval[j]+(1/ptn))/(freqar[1,j]+1/ptn)); freqar[2,j]:=ln((topval[j]+(1/ptn))/(freqar[2,j]+1/ptn)); freqar[3,j]:=ln((topval[j]+(1/ptn))/(freqar[3,j]+1/ptn)); freqar[4,j]:=ln((topval[j]+(1/ptn))/(freqar[4,j]+1/ptn)); if brg = 0 then fnbrg:=1 else fnbrg:=ln(brg); if brc = 0 then fnbrc:=1 else fnbrc:=ln(brc); if bra = 0 then fnbra:=1 else fnbra:=ln(bra); if brt = 0 then fnbrt:=1 else fnbrt:=ln(brt); if gbr = 0 then fngbr:=1 else fngbr:=ln(gbr); if cbr = 0 then fncbr:=1 else fncbr:=ln(cbr); if abr = 0 then fnabr:=1 else fnabr:=ln(abr); if tbr = 0 then fntbr:=1 else fntbr:=ln(tbr); if j<=length then begin sn:=ptn; if sn > 100 then sn:=100; { else if sn > 90 then sn:=90 else if sn>80 then sn:=80 else if sn>70 then sn:=70 else if sn > 60 then sn:=60 else if sn>50 then sn:=50 else if sn>40 then sn:=40 else if sn>33 then sn:=33; } if LTG<>2 then rseq[j]:=(brg*1.4427*(fnbrg)+brc*1.4427*(fnbrc)+ bra*1.4427*(fnbra)+brt*1.4427*(fnbrt))- (brg*1.4427*(fngbr)+brc*1.4427*(fncbr)+ bra*1.4427*(fnabr)+brt*1.4427*(fntbr)) else rseq[j]:=LT[sn,1]+(brg*1.4427*fnbrg+brc*1.4427*fnbrc+ bra*1.4427*fnbra+brt*1.4427*fnbrt); rseqdb[j]:=rseq[j]; if rseqdb[j]<0 then rseqdb[j]:=0; w:=0; x:=0; y:=0; z:=0; error:=0; end;(*if j*) end;{for j} (********rseqdb has pos RI values for the dbs**********) (*****now add tst seq 1 at a time to top db grp and recalcRI**) for k:=1 to filnum do begin readln(ifil1,SA); write(ofil,' ',' '); for j:=1 to length do begin x1:=SA[j]; x1:=upcase(x1); if not ((x1='G') or (x1='C') or (x1='A') or (x1='T')or (x1='N')) then begin WRITELN('INPUT ERROR A,T,G,C EXPECTED'); error:=succ(error); readln; end; case x1 of 'G': mbvh[j]:=freqar[1,j]; 'C': mbvh[j]:=freqar[2,j]; 'A': mbvh[j]:=freqar[3,j]; 'T': mbvh[j]:=freqar[4,j]; end;(*case*) end;(*for j*) for l:=1 to length do begin if (l mod 6 =0) then writeln(ofil); rseq[l]:=rseqdb[l]*mbvh[l]; rtot:=rtot+rseq[l]; (****write graded binary values to output file +expected output**) if rseq[l]8.5 then write(ofil,'1 1 ') else if rtot>6.5 then write(ofil,'1 0 ') else if rtot>4.5 then write(ofil,'0 1 ') else write(ofil,'0 0 '); rtot:=0; writeln(ofil); if o='0' then writeln(ofil,'0 1') else writeln(ofil,'1 0'); end;(*for k*) end;(*calcRI*) (**********************************************************) Begin(*main*) questions; initialize; getdb; calcRI; close(ifil1); close(ifil2); close(ofil); end.