From: SMTP%"clark@mshri.utoronto.ca" 12-FEB-1990 16:47 To: GILBERTD Subj: GenBank email shells Date: Mon, 12 Feb 90 16:36:37 EST Message-Id: <9002122136.AA12391@lash.utcs.utoronto.ca> From: clark@mshri.utoronto.ca To: @mail-shells@lash.utcs.utoronto.ca Subject: GenBank email shells Fellow GCGer, Thanks for your interest in the DCL shells I wrote to make it easier to use the GenBank services that are available via Email. The procedures are in the following order, separated by a string of 20 asterisks: FAMAIL.COM FAMAILQ.COM DBMAIL.COM SEQINFO.FOR (used by FAMAIL.COM) WHAT THESE SHELLS DO -------------------- Shortly before Christmas 1989 GenBank made available search and retrieval services that can be accessed via electronic mail. There are no charges to use these services. These shells construct messages in the format expected by GenBnak and sends them by VMS mail. The users do not have to be familiar with the format of the message that GenBank expects, nor with sending email across the network; however, they do need to know how to read the mail that GenBank returns. FAMAIL.COM initiates a Fasta-type search of one of the GenBank databases. It enquires about the sequence name and a few search parameters, converts the sequence to IG format (required by the GenBank system), constructs the message and mails it off. FAMAILQ.COM sends a message to the GenBank Fasta server, inquiring about the state of the queue. This is used to determine where your job is waiting in the queue, which is important to know because GenBank has restrictions on how many jobs each user can have waiting in the queue. DAMAIL.COM sends a request to GenBank to return a database entry. The sequence can be specified by locus name or accession number. Only one entry can be retrieved per email request. SEQINFO.FOR is a GCG-style program required by FAMAIL.COM. It checks that the sequence specified by the user is really in GCG format, and returns the length and some other information to the shell via symbols. This program will have to be compiled and linked under the GCG environment. I have found the GenBank email services to be extremely useful. Not only does submitting a Fasta serach to them offload the proecessing from our VAX, but their database is updated daily (weekly at worst), and the reply usually comes in less than than an hour when searching the entire database (or a few minutes for searching a subset, eg, primates), which is much faster than our microVAX II can manage it. The retrieval requests are returned in a matter of seconds. Note: We are on Internet and I don't know what the response time will be like for you folks on Bitnet. I would appreciate it if you Bitnetters could let me know, just to satsify my curiosity. To get more information about these services, GenBank will send you a short document describing them if you send the word HELP on the first line of an email message (with no subject: line) to SEARCH@GENBANK.BIO.NET (for Fasta help) or RETRIEVE@GENBANK.BIO.NET (for retrieval help). I strongly recommend that you do this. INSTALLATION ------------ I put the command procedures in a subdirectory of my personal account, away from the rest of the GCG stuff. This is so that 1) I can find them if I need to make changes, and 2) so they don't get touched by anything on the GCG update tapes (not likely, but always possible). The appropriate initializing GCG command procedure has to be edited so that the symbols FAMAIL, FAMAILQ and DBMAIL point to the shells. You will also have to add a symbol for the program SEQINFO, such as $ SEQINFO :== $device:[directory]SEQINFO Don't forget the $ sign in front of the device name or logical. SEQINFO must be compiled and linked. This is done by initiating the GCG support environment with the command GCGSUPPORT, compiling the program (FORTRAN SEQINFO) and linking it (GENLINK SEQINFO). The genlink command knows where all the appropriate object libraries are kept. The three command procedures use a symbol called SEARCH_ADDRESS or RETRIEVAL_ADDRESS that is equated to the network address of the GenBank servers. These WILL have to be changed to accomodate your local conditions, such as gateways, or whatever. OPERATION --------- If the installation is done properly, you should be able to run the shells just by typing their names, as with other GCG programs. I see no point in going over their use point by point. The best way to see how they work is to try them. Both FAMAIL and DBMAIL allow you to specify on the command line the sequence in question. For example, to have GenBank perform a Fasta search of a sequence called nobel_sexy.seq, typing "famail nobel_sexy.seq" will result in the shell skipping the prompt asking for the sequence name. Note: The GenBank Fasta server only searches one orientation of the sequence - to search the opposite strand, it is necessary to send a second search request, specifying that the sequence should be reversed. MODIFICATIONS, UPDATES, BUGS ---------------------------- If you find that the shells don't support all the options you need, feel free to add to them. It should be easy to find where to make changes because the programs are pretty well commented, as these things go, and besides, isn't it true that DCL is the only real self-documenting language? I only request that you leave the lines at the top of the files that reference me as the author, and the institution where I work, unless you make such hash out of them that they no longer work, in which case I would be happy to relinquish credit. Alternatively, send suggestions for enhancements to me, and I'll consider them. If you find any bugs, please let me know so I can fix them and alert the other users of these shells. Stephen Clark, Ph.D. Computer Resources Manager Samuel Lunenfeld Research Institute Mount Sinai Hospital 600 University Avenue Toronto, Ontario Canada M5G 1X5 clark@mshri.utoronto.ca (Internet) sinai@utoroci (Netnorth/Bitnet) ******************** $ ! FAMAIL.COM $ $ ! December 22, 1989 $ $ ! Written by Steve Clark $ ! Mt. Sinai Hospital Research Institute, Toronto, Canada $ $ ! Command procedure to send a sequence to GenBank to have a FASTA $ ! search performed on it. THis procedure asks all the relevant $ ! questions, constructs a text file with the sequence in Intelligenetics $ ! format, and mails it to GenBank. It accepts the name of the query $ ! sequence on the command line as P1. $ ! Note: the symbol SEARCH_ADDRESS is the network address for the Genbank $ ! Search service. This will have to be changed to accomodate local $ ! gateways, etc. $ $ on control_y then goto terminate $ bell[0,7] = 7 $ ws := "write sys$output" $ iq := inquire/nopunctuation $ $ ! The Internet address for sending the search file is $ ! SEARCH@GENBANK.BIO.NET $ $ search_address := "lash::genbank.bio.net::search" $ $ ws "" $ ws "This procedure initiates a FASTA search for similarity between" $ ws "your query sequence and one of the databases maintained by GenBank." $ ws "The information required for executing the search is sent to" $ ws "GenBank via electronic mail and is executed by the GenBank people" $ ws "themselves. Their databases are much more current than our local" $ ws "ones, and their computer is very fast. You must specify which" $ ws "strand to search, since GenBank only searches ONE of the strands" $ ws "of your sequence. The results of the search will be returned to" $ ws "you via e-mail." $ ws "" $ ws "Please remember that if you submit a second search against ALL of" $ ws "GenBank or EMBL while a previous search is still waiting to be" $ ws "executed, the previous search will be aborted." $ $ seqname := "''p1'" $ if(seqname.NES."") then goto no_query $ $get_query: $ $ ws "" $ iq seqname "GenBank FASTA with what query sequence? " $ if(seqname.EQS."") then goto get_query $ $no_query: $ $ ! See if the sequence exists $ $ assign/user_mode nl: sys$output $ seqinfo/infile='seqname' $ if(seqinfotype.NES."NONE") then goto check_gcg $ ws "" $ ws "''bell'''seqname' doesn't exist. Please try again." $ goto get_query $ $check_gcg: $ $ ! Check if the sequence is in GCG format. $ $ if(seqinfotype.NES."NOGCG") then goto get_start $ ws "" $ ws "''bell'''seqname' is not a legitimate GCG sequence file!" $ ws "" $ ws "Select option by number -" $ ws "" $ ws "1) Specify another sequence" $ ws "2) Quit" $ ws "" $ iq choice "Choice (* 1 *) ? " $ if(choice.EQS."2") then exit $ goto get_query $ $get_start: $ $ ws "" $ iq begin "Begin (* 1 *) ? " $ if (begin.EQS."") then begin := 1 $ ibegin = f$integer(begin) $ if((ibegin.GE.1).AND.(ibegin.LT.f$integer(seqinfolength))) then - goto get_end $ ws "''bell'The start must be between 1 and ''seqinfolength'." $ goto get_start $ $get_end: $ $ iq end "End (* ''seqinfolength' *) ? " $ if (end.EQS."") then end := 'seqinfolength' $ iend = f$integer(end) $ if((iend.GT.ibegin).AND.(iend.LE.f$integer(seqinfolength))) then - goto get_reverse $ ws "''bell'The start must be between 1 and ''seqinfolength'." $ goto get_end $ $get_reverse: $ $ iq reverse "Reverse (* No *)? " $ if(reverse.EQS."") then reverse := NO $ if(f$extract(0,1,reverse).EQS."Y") then reverse := "YES" $ if(f$extract(0,1,reverse).EQS."N") then reverse := "NO" $ if(reverse.EQS."YES") then goto get_database $ if(reverse.EQS."NO") then goto get_database $ ws "''bell'Please answer Yes or No." $ goto get_reverse $ $get_database: $ $ ! Find out which database to search. If the sequence is DNA, the $ ! default is the Genbank database. The default for proteins is $ ! the SWISS-PROT database $ $ ws "" $ ws "Database to search:" $ ws "" $ ws " 1) ALL of GenBank" $ ws " 2) GenBank Primate sequences" $ ws " 3) GenBank Rodent sequences" $ ws " 4) Other GenBank Mammalian sequences" $ ws " 5) Other GenBank Vertebrate sequences" $ ws " 6) GenBank Invertebrate sequences" $ ws " 7) GenBank Plant sequences" $ ws " 8) GenBank Bacterial sequences" $ ws " 9) GenBank Organelle sequences" $ ws "10) GenBank Phage sequences" $ ws "11) GenBank Viral sequences" $ ws "12) GenBank Structural RNA sequences" $ ws "13) GenBank Synthetic sequences" $ ws "14) GenBank Unannotated sequences" $ ws "15) New GenBank sequences since the last quarterly release" $ ws "16) ALL of EMBL" $ ws "17) New EMBL sequences since the last release $ ws "18) ALL of SWISS-PROT $ ws "" $ default := 1 $ if(seqinfotype.EQS."PROTEIN") then default := 18 $ iq choice "Please enter choice (* ''default' *): " $ if(choice.EQS."") then choice := 'default' $ database := "" $ if(choice.EQS."1") then database := GENBANK/ALL $ if(choice.EQS."2") then database := GENBANK/PRIMATE $ if(choice.EQS."3") then database := GENBANK/RODENT $ if(choice.EQS."4") then database := GENBANK/OTHER_MAMMALIAN $ if(choice.EQS."5") then database := GENBANK/OTHER_VERTEBRATE $ if(choice.EQS."6") then database := GENBANK/INVERTEBRATE $ if(choice.EQS."7") then database := GENBANK/PLANT $ if(choice.EQS."8") then database := GENBANK/BACTERIAL $ if(choice.EQS."9") then database := GENBANK/ORGANELLE $ if(choice.EQS."10") then database := GENBANK/PHAGE $ if(choice.EQS."11") then database := GENBANK/VIRAL $ if(choice.EQS."12") then database := GENBANK/STRUCTURAL_RNA $ if(choice.EQS."13") then database := GENBANK/SYNTHETIC $ if(choice.EQS."14") then database := GENBANK/UNANNOTATED $ if(choice.EQS."15") then database := GENBANK/NEW $ if(choice.EQS."16") then database := EMBL/ALL $ if(choice.EQS."17") then database := EMBL/NEW $ if(choice.EQS."18") then database := SWISS-PROT/ALL $ if(database.NES."") then goto get_wordsize $ ws "''bell'Valid responses are 1 - 18, inclusive." $ goto get_database $ $get_wordsize: $ $ ! Find out how long the word should be $ $ if(seqinfotype.EQS."PROTEIN") then goto get_pwordsize $ $ ws "" $ iq wordsize "What word size (* 4 *)? " $ if(wordsize.EQS."") then wordsize := 4 $ if((f$integer(wordsize).GT.2).AND.(f$integer(wordsize).LT.7)) - then goto get_nscores $ ws "''bell'Word size must be in the range from 3 to 6." $ goto get_wordsize $ $get_pwordsize: $ $ ws "" $ iq wordsize "What word size (* 1 *)? " $ if(wordsize.EQS."") then wordsize := 1 $ if((f$integer(wordsize).GT.0).AND.(f$integer(wordsize).LT.3)) - then goto get_nscores $ ws "''bell'Word size must be in the range from 1 to 2." $ goto get_pwordsize $ $get_nscores: $ $ ws "" $ iq nscores "List how many best scores (* 100 *)? " $ if(nscores.EQS."") then nscores := 100 $ if(f$integer(nscores).GE.10) then goto get_nalign $ ws "''bell'At least 10 scores should be listed." $ goto get_nscores $ $get_nalign: $ $ ws "" $ iq nalign "Align how many matches (* 20 *)? " $ if(nalign.EQS."") then nalign := 20 $ if(nalign.GT.4) then goto summarize $ ws "''bell'At least 5 alignments should be shown." $ goto get_nalign $ $summarize: $ $ ws "" $ ws "" $ ws "The following FASTA search will be executed:" $ ws "" $ if(reverse.EQS."NO") then ws "Query sequence: ''seqname' from ", - "''begin' to ''end' (''seqinfotype')" $ if(reverse.EQS."YES") then ws "Query sequence: Reverse of ''seqname'", - " from ''begin' to ''end' (''seqinfotype')" $ ws "Database to be searched: ''database'" $ ws "Word size: ''wordsize'" $ ws "Scores to list: ''nscores'" $ ws "Alignments to show: ''nalign'" $ ws "" $ iq choice "Are these parameters correct (* Yes *)? " $ choice = f$extract(0, 1, choice) $ if(choice.EQS."") then goto do_it $ if(choice.EQS."Y") then goto do_it $ $ ! Something is wrong. Give the chance to correct it, or give up. $ $ask_repeat: $ $ ws "" $ ws "Do you want to" $ ws "" $ ws "1) Try again" $ ws "2) Give up" $ ws "" $ iq choice "Please enter the number of your choice (* 1 *): " $ if(choice.eqs."") then goto get_query $ if(choice.eqs."1") then goto get_query $ if(choice.eqs."2") then exit $ ws "''bell'Wasn't the question simple enough for you?" $ goto ask_repeat $ $do_it: $ $ ! Determine the root name of the sequence for comparison. This in used $ ! in specifying the output file name, which has the extension .FAM. $ $ ! Check if the sequence is a file. If not, assume it is a database $ ! entry and get the locus name to use as a root. $ $ seqroot = seqname $ root = f$parse(seqroot,,,"NAME") ! root name of sequence $ if(root.NES."") then goto set_outname $ $ ! Remove database name $ $ pos = 'f$locate(":", seqroot)' $ len = 'f$length(seqroot)' $ if(pos.NE.len) then root = f$extract(pos+1, len-pos, seqroot) $ $ set_outname: $ $ outname := "''root'.FAM" $ $ ! Convert the sequence to Intelligenetics format. Since ToIG does not $ ! allow command line input, need to construct a little command $ ! procedure to substitute in the appropriate sequence names. $ $ ws "" $ ws "Converting ''seqname' to IntelliGenetics format..." $ ws "" $ $ open/write comfile fam.cmd $ wc := "write comfile" $ wc "$ set noon" $ wc "$ set verify" $ wc "ToIG" $ wc "''seqname'" $ wc "''begin'" $ wc "''end'" $ wc "''reverse'" ! Reverse? $ wc "''root'.IG" ! Output file name $ wc "$ set noverify" $ close comfile $ @fam.cmd $ del fam.cmd;0 $ $ ! Write the text file that will be mailed to GenBank $ $ ws "Creating the file to be mailed to GenBank..." $ $ open/write comfile 'outname' $ wc "DATALIB ''database'" $ wc "KTUP ''wordsize'" $ wc "SCORES ''nscores'" $ wc "ALIGNMENTS ''nalign'" $ wc "BEGIN" $ $ ! The IG format files that GCG makes starts with a blank line, which $ ! chokes the GenBank FASTA program. Therefore it is not possible to $ ! just append the sequence file to the text file. Read it line by $ ! line, copying over only those records that aren't zero length. (At $ ! the present time there is just one blank line - the first one - but $ ! who knows what the future has in store for us?) $ $ open/read igfile 'root'.IG $ $read_loop: $ $ read/end_of_file=eof igfile record $ if(f$length(record).EQ.0) then goto read_loop $ write comfile record $ $ goto read_loop $ $eof: $ $ close igfile $ close comfile $ delete 'root'.IG;0 $ $ ! Mail the file away. $ $ ws "Mailing the file to GenBank..." $ $ mail 'outname' 'search_address' $ $ ws "The file ''outname' has been sent to GenBank." $ ws "The results will be mailed back to you shortly." $ ws "" $ $ delete 'outname';0 $ $ terminate: ! Jump here on ^y. Don't do any cleanup $ $ exit ******************** $ ! FAMAILQ.COM $ $ ! December 27, 1989 $ $ ! Written by Steve Clark $ ! Mt. Sinai Hospital Research Institute, Toronto, Canada $ $ ! Command procedure to mail to the GenBank FASTA server a query about $ ! the status of the searches waiting in the queues. $ ! It works by constructing a text file containing the single word $ ! "QUEUE" and mailing it to SEARCH@GENBANK.BIO.NET $ $ ! NOTE: The symbol SEARCH_ADDRESS is the network address for the $ ! GenBank search service. This will have to be changed to accomodate $ ! local gateways, etc. $ $ search_address := "lash::genbank.bio.net::search" $ $ ws := "write sys$output" $ ws "" $ ws "This procedure inquires about the status of the FASTA server" $ ws "at GenBank via electronic mail so you can determine where your" $ ws "searches are waiting in the queues." $ ws "" $ inquire/nopunctuation choice "Do you want to continue (* Yes *)? " $ if(choice.EQS."") then choice := YES $ if(f$extract(0,1,choice).NES."Y") then exit $ ws "" $ ws "Constructing message file..." $ open/write comfile search_queue.txt $ write comfile "QUEUE" $ close comfile $ ws "Mailing message file..." $ mail search_queue.txt 'search_address' $ delete search_queue.txt;* $ ws "" $ ws "GenBank will mail back the search queue status shortly." $ exit ******************** $ !DBMAIL.COM $ $ ! December 27, 1989 $ $ ! Written by Steve Clark $ ! Mt. Sinai Hospital Research Institute, Toronto, Canada $ $ !Command procedure to mail a request to GenBank for a database sequence. $ ! GenBank will return the sequence via email. The sequence can be $ ! specified by either its locus name or accession number, from either $ ! the GenBank or EMBL DNA databases. The sequence to be retrieved can $ ! specified on the command line as parameter 1. $ $ !NOTE: The symbol RETRIEVAL_ADDRESS is the GenBank address for the $ ! retrieval service: RETRIEVE@GENBANK.BIO.NET. This will have to be $ ! changed to accomodate local gateways, etc. $ $ retrieval_address := "lash::genbank.bio.net::retrieve" $ $ ws := "write sys$output" $ $ ws "" $ ws "This procedure retrieves from GenBank a single nucleotide sequence" $ ws "via electronic mail. The EMBL or GenBank sequence can be specified" $ ws "either by its ACCESSION NUMBER or LOCUS NAME, but not both." $ ws "" $ $ sequence := "''p1'" $ if(sequence.NES."") then goto no_query $ $ask_sequence: $ $ inquire/nopunctuation sequence "Sequence to retrieve: " $ if(sequence.EQS."") then goto ask_sequence $ $no_query: $ $ ws "Constructing message file..." $ open/write comfile retrieval_sequence.txt $ write comfile "''sequence'" $ close comfile $ ws "Mailing the request to GenBank..." $ mail retrieval_sequence.txt 'retrieval_address' $ delete retrieval_sequence.txt;0 $ ws "" $ ws "GenBank will mail the sequence back to you shortly. When it" $ ws "arrives, you will have to EXTRACT it from mail into a disk file," $ ws "then convert it to the GCG format with the program 'FromGenBank'." $ ws "" $ exit ******************** !*** SEQINFO *********************************************************** !* !* This program reads the sequence specified on the command line with the !* qualifier /INfile=... and sets two DCL symbols: SEQINFOTYPE ("PROTEIN", !* "DNA", "NONE" (if the sequence doesn't exist) or "NOGCG" (if the sequence !* is not in the GCG format)), and SEQINFOLENGTH (the length of the sequence). !* It is intended to be executed from a command procedure, so the !* command procedure can get this information about any specified sequence. !* !* To be used with the command procedure shells, the symbol !* SEQINFO :== $LOGICAL:SEQINFO !* has to be set, where LOGICAL is of course replaced with the logical name !* of the appropriate device and directory. !* !* Written by Steve Clark November 15, 1989 !* !************************************************************************* program seqinfo implicit none character seqfname(256), strand(100001), lengthstr(64) integer infile, checksum, length, inttostr, istatus logical isprotein, readsq, clgetoldfname, dclsetsymbol logical openf, logstatus ! Check for the existance of the sequence file name on the command line if( .not. clgetoldfname('INfile', 1, seqfname)) then logstatus = dclsetsymbol('seqinfotype', 'NONE') stop ' ' endif ! Open the file. if( .not. openf(infile, seqfname, 'rdb')) then logstatus = dclsetsymbol('seqinfotype', 'NONE') stop ' ' endif ! Read the sequence to determine its length. if( .not. readsq(infile, strand, length, checksum)) then logstatus = dclsetsymbol('seqinfotype', 'NOGCG') stop ' ' endif istatus = inttostr(length, lengthstr) logstatus = dclsetsymbol('seqinfolength', lengthstr) ! Find out if it is protein or DNA if(isprotein(strand)) then logstatus = dclsetsymbol('seqinfotype', 'PROTEIN') else logstatus = dclsetsymbol('seqinfotype', 'DNA') endif stop ' ' end