home *** CD-ROM | disk | FTP | other *** search
/ Mac-Source 1994 July / Mac-Source_July_1994.iso / Pascal / source / Dotty Folder / readseq.p < prev   
Encoding:
Text File  |  1989-06-27  |  16.8 KB  |  600 lines  |  [TEXT/MPS ]
  1. {[b+,c+,n+,u+,r+,rec+,#+,j=13-/40/1o,:+,t=2,o=80] PasMat }
  2. { readseq.p -- read nucleic/protein sequence in various formats
  3.       -- file may have multiple sequences,
  4.       -- based loosely on uwgcg & zuker ftn source
  5. }
  6.  
  7. UNIT sequenceIO;
  8.   INTERFACE
  9.  
  10.     TYPE
  11.       Sequence     = PACKED ARRAY [0..32000] OF CHAR;
  12.       SeqPtr       = ^Sequence;
  13.             
  14.     CONST
  15.       {readSeq control calls}
  16.       rseqFetch1   = 0;                {fetch 1st seq, no prior read}
  17.       rseqList     = 1;                {get list of available seq}
  18.       rseqFetchN   = 2;                {fetch nth seq, with Inquire info}
  19.  
  20.     FUNCTION readSeq(fname: String;    {file name}
  21.                      seq: SeqPtr;      {sequence storage(0,2), or list(1)}
  22.                      VAR nseq: longint; {in: storage size, out: seq size}
  23.                      VAR seqId: String; {in(2): seq to fetch, out: seq info}
  24.                      control: integer  {rseq control(0,1,2)}
  25.                      ): integer;       {error code: 0=okay}
  26.  
  27.   IMPLEMENTATION
  28.  
  29.     FUNCTION readSeq(fname: String;    {file name}
  30.                      seq: SeqPtr;      {sequence storage(0,2), or list(1)}
  31.                      VAR nseq: longint; {in: storage size, out: seq size}
  32.                      VAR seqId: String; {in(2): seq to fetch, out: seq info}
  33.                      control: integer  {rseq control(0,1,2)}
  34.                      ): integer;       {error code: 0=okay}
  35.  
  36.       VAR
  37.         f           : text;
  38.         err, i      : integer;
  39.         s, mySeq    : String;
  40.         done, gotuw, isfitch, addit, doadd: boolean;
  41.         maxseq, ninfo: longint;
  42.         seqCharSet  : SET OF CHAR;
  43.  
  44.       PROCEDURE exitseq(err: integer);
  45.  
  46.         BEGIN
  47.           IF control = rseqList THEN nseq := ninfo;
  48.           close(f);
  49.           readSeq := err;
  50.           exit(readSeq);
  51.         END;
  52.  
  53.       FUNCTION isSeqChar(c: CHAR): boolean;
  54.     {---
  55.     NOTE: (Var c:char) messes up ! (MPW Pas)
  56.  
  57.     GCG programs allow all upper and lower case
  58.     letters, periods  (.),
  59.     asterisks  (*),  pluses  (+), ampersands (&), and ats
  60.     (@) as symbols in biological seqences.
  61.     Note: Digits ['0'..'9'] are NOT allowed here.
  62.     ----}
  63.  
  64.         BEGIN
  65.           IF (c <= ' ') THEN isSeqChar := false
  66.           ELSE IF (c IN seqCharSet) THEN isSeqChar := true
  67.           ELSE isSeqChar := false;
  68.         END;                           {isSeqChar}
  69.  
  70.       PROCEDURE addseq( s: String);
  71.  
  72.         VAR
  73.           i           : integer;
  74.  
  75.         BEGIN
  76.           FOR i := 1 TO length(s) DO
  77.             IF (nseq < maxseq) AND isSeqChar(s[i]) THEN BEGIN
  78.               IF addit THEN seq^[nseq] := s[i];
  79.               nseq := nseq + 1;
  80.             END;
  81.         END;                           {addseq}
  82.  
  83.       PROCEDURE addinfo(VAR s: string);
  84.  
  85.         VAR
  86.           i, l        : integer;
  87.  
  88.         BEGIN
  89.           l := length(s) + 1;
  90.           s[l] := chr(13);
  91.           FOR i := 1 TO l DO
  92.             IF (ninfo < maxseq) THEN BEGIN
  93.               seq^[ninfo] := s[i];
  94.               ninfo := ninfo + 1;
  95.             END;
  96.         END;                           {addinfo}
  97.  
  98.       PROCEDURE readLoop(margin: integer; endadd, firstadd: boolean; FUNCTION
  99.                          endTest(VAR s: String): boolean);
  100.  
  101.         BEGIN
  102.           IF control = rseqFetchN THEN doadd := seqId = mySeq
  103.           ELSE doadd := true;
  104.           nseq := 0;
  105.           IF firstadd THEN addseq(s);  {! fitch 1st string}
  106.           REPEAT
  107.             {! check eof Before read !}
  108.             done := eof(f);
  109.             readln(f, s);
  110.             done := done OR endTest(s);
  111.             IF doadd AND (endadd OR NOT done) THEN BEGIN
  112.               IF margin > 0 THEN delete(s, 1, margin);
  113.               addseq(s);
  114.             END;
  115.           UNTIL done;
  116.  
  117.           CASE control OF
  118.             rseqFetch1: exitseq(err);
  119.             rseqList: addinfo(seqId);
  120.             rseqFetchN: IF doadd THEN exitseq(err);
  121.           END;
  122.         END;                           {readLoop}
  123.  
  124.       FUNCTION endIG(VAR s: String): boolean;
  125.  
  126.         BEGIN
  127.           endIG := (pos('1', s) > 0) OR (pos('2', s) > 0);
  128.         END;
  129.  
  130.       PROCEDURE readIG;                {IG -- many seqs/file }
  131.  
  132.         BEGIN
  133.           WHILE true DO BEGIN
  134.             REPEAT
  135.               readln(f, s);
  136.             UNTIL eof(f) OR ((s <> '') AND (pos(';', s) <> 1));
  137.             IF eof(f) THEN exitseq(err);
  138.             seqId := concat(s, '  [Stanford/IG 1]');
  139.             readLoop(0, true, false, endIG);
  140.             WHILE NOT (eof(f) OR ((s <> '') AND (pos(';', s) <> 1))) DO
  141.               readln(f, s);
  142.             IF eof(f) THEN exitseq(err);
  143.           END;
  144.         END;                           {readIG}
  145.  
  146.       FUNCTION endStrider(VAR s: String): boolean;
  147.  
  148.         BEGIN
  149.           endStrider := (pos('//', s) > 0);
  150.         END;
  151.  
  152.       PROCEDURE readStrider;
  153.  
  154.         BEGIN
  155.           WHILE true DO BEGIN
  156.             {Strider -- 1 seq only? }
  157.             readln(f, s);
  158.             seqId := concat(s, '  [DNA Strider 2]');
  159.             delete(seqId, 1, 1);
  160.             REPEAT
  161.               readln(f, s);
  162.             UNTIL eof(f) OR (pos(';', s) <> 1);
  163.             readLoop(0, false, true, endStrider);
  164.             WHILE NOT (eof(f) OR ((s <> '') AND (pos(';', s) <> 1))) DO
  165.               readln(f, s);
  166.             IF eof(f) THEN exitseq(err);
  167.           END;
  168.         END;                           {readStrider}
  169.  
  170.       FUNCTION endGB(VAR s: String): boolean;
  171.  
  172.         BEGIN
  173.           endGB := (pos('//', s) > 0) OR (pos('LOCUS', s) = 1);
  174.         END;
  175.  
  176.       PROCEDURE readGenBank;
  177.  
  178.         BEGIN
  179.           WHILE true DO BEGIN
  180.             {GenBank -- many seqs/file }
  181.             seqId := concat(s, '  [GenBank 3]');
  182.             delete(seqId, 1, 12);
  183.             REPEAT
  184.               readln(f, s);
  185.             UNTIL eof(f) OR ((s <> '') AND (pos('ORIGIN', s) = 1));
  186.             readLoop(9, false, false, endGB);
  187.             WHILE NOT (eof(f) OR ((s <> '') AND (pos('LOCUS', s) = 1))) DO
  188.               readln(f, s);
  189.             IF eof(f) THEN exitseq(err);
  190.           END;
  191.         END;                           {readGenBank}
  192.  
  193.       FUNCTION endNBRF(VAR s: String): boolean;
  194.  
  195.         BEGIN
  196.           endNBRF := (pos('>', s) = 1);
  197.         END;
  198.  
  199.       PROCEDURE readNBRF;
  200.  
  201.         BEGIN
  202.           WHILE true DO BEGIN
  203.             {NBRF -- many seqs/file }
  204.             seqId := concat(s, '  [NBRF/PIR 4]');
  205.             delete(seqId, 1, 4);
  206.             readln(f, s);              {junk line}
  207.             readLoop(0, false, false, endNBRF);
  208.             WHILE NOT (eof(f) OR ((s <> '') AND (pos('>DL;', s) = 1))) DO
  209.               readln(f, s);
  210.             IF eof(f) THEN exitseq(err);
  211.           END;
  212.         END;                           {readNBRF}
  213.  
  214.       FUNCTION endPearson(VAR s: String): boolean;
  215.  
  216.         BEGIN
  217.           endPearson := (pos('>', s) = 1);
  218.         END;
  219.  
  220.       PROCEDURE readPearson;
  221.  
  222.         BEGIN
  223.           WHILE true DO BEGIN
  224.             {Pearson -- many seqs/file }
  225.             seqId := concat(s, '  [Pearson 5]');
  226.             delete(seqId, 1, 1);
  227.             readLoop(0, false, false, endPearson);
  228.             WHILE NOT (eof(f) OR ((s <> '') AND (pos('>', s) = 1))) DO
  229.               readln(f, s);
  230.             IF eof(f) THEN exitseq(err);
  231.           END;
  232.         END;                           {readPearson}
  233.  
  234.       FUNCTION endEMBL(VAR s: String): boolean;
  235.  
  236.         BEGIN
  237.           endEMBL := (pos('ID   ', s) = 1);
  238.         END;
  239.  
  240.       PROCEDURE readEMBL;
  241.  
  242.         BEGIN
  243.           WHILE true DO BEGIN
  244.             seqId := concat(s, '  [EMBL 6]');
  245.             delete(seqId, 1, 5);
  246.             REPEAT
  247.               readln(f, s);
  248.             UNTIL eof(f) OR (pos('SQ   ', s) = 1);
  249.             readLoop(0, false, false, endEMBL);
  250.             WHILE NOT (eof(f) OR ((s <> '') AND (pos('ID   ', s) = 1))) DO
  251.               readln(f, s);
  252.             IF eof(f) THEN exitseq(err);
  253.           END;
  254.         END;                           {readEMBL}
  255.  
  256.       FUNCTION endZuker(VAR s: String): boolean;
  257.  
  258.         BEGIN
  259.           endZuker := (pos('(', s) = 1);
  260.         END;
  261.  
  262.       PROCEDURE readZuker;
  263.  
  264.         BEGIN
  265.           WHILE true DO BEGIN
  266.             {Zuker -- many seq/file ?}
  267.             {! 1st string is Zuker's Fortran format }
  268.             readln(f, s);              {s == "seqLen seqId string..."}
  269.             seqId := concat(s, '  [NRC 7]');
  270.             delete(seqId, 1, 6);
  271.             readLoop(0, false, false, endZuker);
  272.             WHILE NOT (eof(f) OR ((s <> '') AND (pos('(', s) = 1))) DO
  273.               readln(f, s);
  274.             IF eof(f) THEN exitseq(err);
  275.           END;
  276.         END;                           {readZuker }
  277.  
  278.       FUNCTION endFitch(VAR s: String): boolean;
  279.  
  280.         BEGIN
  281.           endFitch := (s[1] <> ' ');
  282.         END;
  283.  
  284.       PROCEDURE readFitch;
  285.  
  286.         VAR
  287.           first       : boolean;
  288.  
  289.         BEGIN
  290.           first := true;
  291.           WHILE true DO BEGIN
  292.             {Fitch -- many seqs/file }
  293.             IF NOT first THEN seqId := s;
  294.             seqId := concat(seqId, '  [Fitch 8]');
  295.             readLoop(0, false, first, endFitch);
  296.             IF eof(f) THEN exitseq(err);
  297.             first := false;
  298.           END;
  299.         END;                           {readFitch }
  300.  
  301.       PROCEDURE readPreer;
  302.    {Used locally at InU only
  303.     control=99=rseqList else rseqFetch1 }
  304.  
  305.         BEGIN
  306.           {readPreer -- 1seq/file}
  307.                     if control=99 then control:= rSeqList
  308.                     else control:= rseqFetch1;
  309.           seqCharSet := seqCharSet + ['0'..'9']; {really, just digits}
  310.           addit := control <> rSeqList;
  311.           readln(f, s);                { skip 1st line = line count }
  312.           seqId := concat(fname, '  [Preer format Z]');
  313.           REPEAT
  314.             readln(f, s);
  315.             done := eof(f);
  316.             IF NOT done THEN addseq(s); {skip last line=index}
  317.           UNTIL done;
  318.  
  319.           IF control = rSeqList THEN addinfo(seqId);
  320.           exitseq(err);
  321.         END;                           {readPreer }
  322.  
  323.       PROCEDURE readUnknown;
  324.  
  325.         BEGIN
  326.           {Unknown -- 1seq/file}
  327.  
  328.           addit := control <> rseqList;
  329.           addseq( seqId);               {from above..}
  330.           seqId := concat(fname, '  [Unknown format, assume all sequence 0]');
  331.           REPEAT
  332.             addseq(s);
  333.             done := eof(f);
  334.             readln(f, s);
  335.           UNTIL done;
  336.  
  337.           CASE control OF
  338.             rseqFetchN, rseqFetch1: ;
  339.             rseqList: addinfo(seqId);
  340.           END;
  341.           exitseq(err);
  342.         END;                           {readUnknown }
  343.  
  344.       PROCEDURE readUWGCG;
  345.  
  346.         BEGIN
  347.           { UWGCG -- 1seq/file}
  348.           addit := control <> rseqList;
  349.           seqId := concat(s, '  [UWGCG 9]');
  350.           REPEAT
  351.             done := eof(f);
  352.             readln(f, s);
  353.             IF NOT done THEN BEGIN
  354.               delete(s, 1, 9);         {skip margin}
  355.               addseq(s);
  356.               { if nseq=0 then sysbeep(1); !true }
  357.             END;
  358.           UNTIL done;
  359.  
  360.           CASE control OF
  361.             rseqFetchN, rseqFetch1: ;
  362.             rseqList: addinfo(seqId);
  363.           END;
  364.           exitseq(err);
  365.         END;                           {readUWGCG }
  366.  
  367.       BEGIN
  368.  
  369.         maxseq := nseq;
  370.         nseq := 0;
  371.         ninfo := 0;
  372.         mySeq := seqId;                {for rseqFetchN}
  373.         seqId := '';
  374.         readSeq := 0;
  375.         addit := control <> rseqList;
  376.         seqCharSet := ['A'..'Z', 'a'..'z', '_', '@', '+', '-', '*', '.', '&'];
  377.  
  378.         reset(f, fname);
  379.         err := ioresult;
  380.         IF err <> 0 THEN exitseq(err);
  381.  
  382.         {InU only fix for Preer data: }
  383.         IF (control IN [98, 99]) THEN readPreer;
  384.  
  385.         IF control = rseqFetchN THEN BEGIN
  386.           {don't need to re-check format}
  387.           REPEAT
  388.             readln(f, s);
  389.           UNTIL (s <> '') OR eof(f);
  390.           IF (s = '') THEN exitseq( - 1);
  391.           CASE mySeq[length(mySeq) - 1] OF
  392.             '0': readUnknown;
  393.             '1': readIG;
  394.             '2': readStrider;
  395.             '3': readGenBank;
  396.             '4': readNBRF;
  397.             '5': readPearson;
  398.             '6': readEMBL;
  399.             '7': readZuker;
  400.             '8': BEGIN
  401.               seqId := s;
  402.               readln(f, s);
  403.               readFitch;
  404.             END;
  405.             '9':
  406.               REPEAT
  407.                 gotuw := pos('..', s) > 0;
  408.                 IF gotuw THEN readUWGCG;
  409.                 readln(f, s);
  410.               UNTIL eof(f);
  411.           END;
  412.         END;
  413.  
  414.      { check for ".." of uwgcg, since it can masquerade as any
  415.      other format }
  416.         i := 0;
  417.         REPEAT
  418.           readln(f, s);
  419.           i := i + 1;
  420.           gotuw := pos('..', s) > 0;
  421.           IF gotuw THEN gotuw := pos('Check:', s) > 0;
  422.           done := gotuw OR eof(f) OR (i > 500);
  423.           {! ECOLAC UW/GenBank document header is 300 lines !}
  424.           IF (i < 5) AND (pos(';', s) = 1) THEN BEGIN
  425.             gotuw := false; done := true;
  426.             {fix for ToIG of UWGCG ... also NBRF/EMBL ?}
  427.           END
  428.         UNTIL done;
  429.  
  430.         IF gotuw THEN readUWGCG
  431.         ELSE BEGIN
  432.           reset(f);
  433.           REPEAT
  434.             readln(f, s);
  435.           UNTIL (s <> '') OR eof(f);
  436.           IF (s = '') THEN exitseq( - 1);
  437.         END;
  438.  
  439.         IF pos(';', s) = 1 THEN BEGIN
  440.           IF pos('Strider', s) > 0 THEN readStrider
  441.           ELSE readIG;
  442.         END
  443.  
  444.         ELSE IF pos('LOCUS', s) = 1 THEN readGenBank
  445.  
  446.         ELSE IF pos('>DL;', s) = 1 THEN readNBRF
  447.  
  448.         ELSE IF pos('>', s) = 1 THEN readPearson
  449.  
  450.         ELSE IF pos('ID   ', s) = 1 THEN readEMBL
  451.  
  452.         ELSE IF pos('(', s) = 1 THEN readZuker
  453.  
  454.         ELSE BEGIN
  455.           seqId := s;
  456.           readln(f, s);                {test for fitch format}
  457.           i := 1;
  458.           REPEAT
  459.             isfitch := (((i - 1) MOD 4 = 0) AND (s[i] = ' ')) OR (((i - 1) MOD
  460.                        4 <> 0) AND (s[i] <> ' '));
  461.             i := i + 1;
  462.           UNTIL (i >= length(s)) OR NOT isfitch;
  463.           IF isfitch THEN readFitch
  464.           ELSE readUnknown;
  465.         END;
  466.  
  467.       END;                             {readSeq}
  468.  
  469. END.
  470.  
  471. {
  472. sequence formats....
  473. ---------------------------------------------------
  474.  
  475. stanford/IG
  476. ;comments
  477. ;...
  478. seq1 info
  479. abcd...
  480. efgh1 (or 2 = terminator)
  481. ;another seq
  482. ;....
  483. seq2 info
  484. abcd...1
  485. --- for e.g. ----
  486. ;     Dro5s-T.Seq  Length: 120  April 6, 1989  21:22  Check: 9487  ..
  487. dro5stseq
  488. GCCAACGACCAUACCACGCUGAAUACAUCGGUUCUCGUCCGAUCACCGAAAUUAAGCAGCGUCGCGGGCG
  489. GUUAGUACUUAGAUGGGGGACCGCUUGGGAACACCGCGUGUUGUUGGCCU1
  490.  
  491. ;  TOIG of: Dro5srna.Seq  check: 9487  from: 1  to: 120
  492. ---------------------------------------------------
  493.  
  494. Genbank:
  495. LOCUS    seq1 ID..
  496. ...
  497. ORIGIN ...
  498. 123456789abcdefg....(1st 9 columns are formatting)
  499.      hijkl...
  500. //         (end of sequence)
  501. LOCUS     seq2 ID ..
  502. ...
  503. ORIGIN
  504.       abcd...
  505. //
  506. ---------------------------------------------------
  507.  
  508. PIR format
  509. > seq1 id
  510. ?? junk 2nd line
  511. abcdefg...
  512. hijkl...
  513. > seq2 ID
  514. ?? junk
  515. abcd....
  516. ---------------------------------------------------
  517.  
  518. NBRF format: (from uwgcg ToNBRF) == PIR format
  519. >DL;DRO5SRNA
  520. Iubio$Dua0:[Gilbertd.Gcg]Dro5srna.Seq;2 => DRO5SRNA
  521.  
  522.        1  GCCAACGAC CAUACCACGC UGAAUACAUC GGUUCUCGUC CGAUCACCGA
  523.       51  AAUUAAGCAG CGUCGCGGGC GGUUAGUACU UAGAUGGGGG ACCGCUUGGG
  524.      101  AACACCGCGU GUUGUUGGCC U
  525.  
  526. >DL;DROEST6
  527. Iubio$Dua0:[Gilbertd.Gcg]Droest6.Seq;1 => DROEST6
  528.  
  529.        1  GAATTCGCC GGAGTGAGGA GCAACATGAA CTACGTGGGA CTGGGACTTA
  530. ---------------------------------------------------
  531.  
  532. EMBL format
  533. ID345 seq1 id   (the 345 are spaces)
  534. ... other info
  535. SQ345Sequence   (the 3,4,5 are spaces)
  536. abcd...
  537. hijk...
  538. ID    seq2 id
  539. ...
  540. SQ   Sequence
  541. abcd...
  542. ...
  543. ---------------------------------------------------
  544.  
  545. UWGCG Format:
  546. comments of any form, up to ".." signal
  547. signal line has seq id.  only 1 seq/ file
  548. Seq ID -----  ..
  549.       1  abcd ///
  550. 123456788^seq starts in 10 chars
  551. -- e.g. ---
  552. LOCUS       DROEST6      1819 bp ss-mRNA            INV       31-AUG-1987
  553. ////////////////
  554. ORIGIN      1 bp upstream of EcoRI site; chromosome BK9 region 69A1.
  555.  
  556. INVERTEBRATE:DROEST6  Length: 1819  January 9, 1989  16:48  Check: 8008  ..
  557.  
  558.        1  GAATTCGCCG GAGTGAGGAG CAACATGAAC TACGTGGGAC TGGGACTTAT
  559.  
  560.       51  CATTGTGCTG AGCTGCCTTT GGCTCGGTTC GAACGCGAGT GATACAGATG
  561.  
  562.  
  563. ---------------------------------------------------
  564.  
  565. DNAStrider (Mac) = modified Stanford:
  566. ; ### from DNA Strider  Friday, April 7, 1989   11:04:24 PM
  567. ; DNA sequence  pBR322   4363  b.p. complete sequence
  568. ;
  569. abcd...
  570. efgh
  571. //  (end of sequence)
  572. ---------------------------------------------------
  573.  
  574. Fitch format:
  575. Dro5srna.Seq
  576.  GCC AAC GAC CAU ACC ACG CUG AAU ACA UCG GUU CUC GUC CGA UCA CCG AAA UUA AGC AGC
  577.  GUC GCG GGC GGU UAG UAC UUA GAU GGG GGA CCG CUU GGG AAC ACC GCG UGU UGU UGG CCU
  578. Droest6.Seq
  579.  GAA TTC GCC GGA GTG AGG AGC AAC ATG AAC TAC GTG GGA CTG GGA CTT ATC ATT GTG CTG
  580.  AGC TGC CTT TGG CTC GGT TCG AAC GCG AGT GAT ACA GAT GAC CCT CTG TTG GTG CAG CTG
  581. ---------------------------------------------------
  582.  
  583. W.Pearson format:
  584. >BOVPRL GenBank entry BOVPRL from omam file.  907 nucleotides.
  585. TGCTTGGCTGAGGAGCCATAGGACGAGAGCTTCCTGGTGAAGTGTGTTTCTTGAAATCAT
  586.  
  587. ---------------------------------------------------
  588.  
  589. * preer (zbasic ZDAT) format:
  590.  459   == length
  591.  3    == numeric representation of GCATU...
  592.  4
  593.  3
  594. ////
  595.  4
  596.  4
  597.  4
  598.  1070  = start index (can be negative...)
  599. }
  600.