From:	SMTP%"clark@mshri.utoronto.ca"  4-MAY-1990 14:44
To:	ARCHIVE
Subj:	Update GCG shells for searching databases in batch

Date: Fri, 4 May 90 15:38:21 EDT
Message-Id: <9005041938.AA06948@lash.utcs.utoronto.ca>
From: clark@mshri.utoronto.ca
To: @gcg-servers@lash.utcs.utoronto.ca
Subject: Update GCG shells for searching databases in batch

Fellow GCGer,                                               May 4, 1990

   Thanks for your interest in the DCL shells I wrote to make it easier
to run the Wisconsin GCG database searching programs in batch mode. The
procedures are in the following order, separated by a string of 20
asterisks: 

WORDSEARCH.COM
FIND.COM
FASTA.COM
TFASTA.COM
STRINGS.COM
MULTIDOT.COM
CHK_PROCESS.COM		(Called by WORDSEARCH.COM)
SEQINFO.FOR


WHAT THESE SHELLS DO
--------------------

   These are DCL programs that interrogate the user as to the kind of 
database search he wants to do: what sequence, which database, what 
options, etc. The shells are completely menu-driven so the user doesn't 
have to know before-hand what commands to use. Intelligent default answers 
are provided most of the time, and all the answers are checked as much as 
possible for errors. Simple error-checking is to make sure the specified 
answer was in the allowable range. Also, all sequences and data files are 
checked for their existance and, whenever possible, conformity to what the 
GCG programs expect. My feeling is that there is no point in submitting 
something to batch and have to search through the .LOG file to discover 
that it gave up after 5 minutes because the name of the sequence was typed 
incorrectly.

   After the user has answered all the questions, a summary of the search 
to be executed is displayed, and the user can then proceed or change the 
search or quit. To start the search, the shell creates a command procedure 
(with the extension .CMD) which it submits to the batch queue SYS$BATCH. 
The submitted batch job then invokes the real GCG program with the 
appropriate command line switches.

   WORDSEARCH.COM and FIND.COM allow the searches to be done interactively 
if the search database is a list of user files, such as "*.seq" or 
"@seqnames.list".

   STRINGS.COM executes the searches interactively if the short definitions 
are being searched, but in batch if the complete records are being 
searched.

   I have avoided the use of IF...THEN...ELSE...ENDIF constructs so they 
should run properly under VMS version 4.x.

   The shells STRINGS.COM, WORDSEARCH.COM, FIND.COM and (T)FASTA.COM 
control the respective GCG program. (WORDSEARCH.COM also runs SEGMENTS to 
align the best matching sequences.) MULTIDOT.COM is a shell that runs 
multiple sets of COMPARE/DOTPLOTs under different conditions without having 
to re-enter all the sequence information for each analysis. The option is 
given first for the type of search - word-style or window/stringency. If 
word-stype is chosen, the user then specifies a starting word size, final 
word size, and increment. For example, if the starting word size is 2, 
final is 8 and increment is 3, three runs will be executed with word sizes 
of 2, 5 and 8.

   Since word-style searches are so fast, the analysis can either be 
executed interactively or in batch.

   Similarly, for window/stringency searches, the user specifies starting, 
final, and increment values for window size and percent similarity. This 
can generate a lot of different analyses to best find or show the homology 
of interest (or perhaps hide the uninteresting coincidnetal sequence 
similarity ;-)). Note that, unlike the original GCG program, the shell asks 
for PERCENT MATCHES in the window, instead of stingency. I think PERCENT 
MATCHES is an easier value to deal with, especially with a window of 
varying size. Also, I have replaced the COMPARE similarity matrix with one 
that might be better called an identity matrix - all 1s on the diagonal and 
zeros everywhere else. This makes it much easier to interpret the resulting 
dot plot. Furthermore, I seriously question the biological relevance of 
some of the similarity that is allowed in the standard GCG matrices. (You 
will have to make your own identity matrix and add the switch /DATA=... to 
the symbol COMPARE.)

   The window/stringency analyses will only run in batch mode.

   At our site, the graphics output is spooled to SYS$PRINT, a DEC LN03 PLUS 
laser printer that emulates a Tektronix device. If you cannot spool your 
graphics output (for example, it expects to go to the terminal), I don't 
know what changes will have to be made to MULTIDOT to make it useful. 
Perhaps the easiest would be to send the graphics to files which could be 
displayed or plotted later.

INSTALLATION
------------

   I put the command procedures in a subdirectory of my personal account, 
away from the rest of the GCG stuff. This is to 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 WORDSEARCH, STRINGS, FIND, FASTA and TFASTA point to the 
shells, rather than to the executable programs. The shells expect the 
executables to live in $GENUTIL:. 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.

   The shell WORDSEARCH.COM has a reference to the command procedure 
CLARK$COMMANDS:CHK_PROCESS. You will have to change the logical 
CLARK$COMMANDS to point to the directory where CHK_PROCESS.COM lives. (This 
is required to make sure that two wordsearches with the same name aren't 
running at the same time. Since the shell also runs SEGMENTS, and the 
wordsearch output file name is based on the name of the sequence, a problem 
could arise if a short search is started after a long search with the same 
query sequence. Both searches will run SEGMENTS on the second WORDSEARCH 
output file, since it has the same name but a higher version number. 
CHK_PROCESS isn't infallible, since it looks at the process name of the 
running batch jobs. If the system is busy, a second job could be submitted 
before an earlier one has progressed far enough to set its process name.)

   (T)FASTA.COM submits the batch jobs it creates to the batch queue 
FASTA$BATCH. This is set up with a priority of 2 and a job_limit of 1 to 
avoid multiple concurrent FASTA searches, which freeze our machine when 
SWAPPER takes over. You will either have to change the queue to which the 
batch job is submitted, or else create one like ours.

   Finally, you will have to compile and link the Fortran program 
SEQINFO.FOR, and a symbol for it must be created as described above. This 
program sets some symbols that the shells use to tell if the specified 
sequence exists, if it is in the GCG format, if it is DNA or protein, and 
how long it is. To link the program, use the GENLINK command under the 
GCG support environment.

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. If you want 
to make their use optional, give the symbol that executes them a different 
name. Operation is pretty straightforward, since everything is menu-driven 
and checked for idiot input. I see no advantage in going over their use point
by point. The best way to see how they work is to try them. I might point out 
that not all options provided by the GCG programs are supported. For 
example, the /PLOT option for WORDSEARCH is ommitted. I used to have the 
plot sent to the laser printer by default, but got tired of throwing out 
the histograms that sat in the output tray for weeks, so I removed it.

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.

   Finally, if I make any major updates (for example, if the GCG programs 
change substantially in future releases) I'll send you some e-mail to 
inquire whether you would like them sent to you.

---------------------------------------------------------------------------

For my next act...

   You may have heard that GenBank now provides a service to do a FASTA 
search of a sequence that you e-mail to them. Not only is the service very 
fast (it takes less than a half hour to get the results), but it offloads 
processing from your VAX and uses the most up-to-date database. The down
side is that it expects a message of a certain strict format, and the
sequence must be in Intelligenetics format. To make it easy for the people
here, I have written a menu-driven shell that asks all the apropriate
questions, creates an ASCII file of the necessary format, and mails it to
GenBank. If you would like a copy of this shell, send me a note. 


Stephen Clark

Computer Resources Manager
Samuel Lunenfeld Research Institute
Mt. Sinai Hospital 
Toronto, Ontario, Canada

clark@mshri.utoronto.ca  (Internet)
sinai@utoroci            (bitnet)


********************
$	! WORDSEARCH.COM
$
$	! May 1, 1990
$
$	! Written by Steve Clark
$	! Samuel Lunenfeld Research Institute
$	! Mt. Sinai Hospital, Toronto, Canada
$
$	! Command procedure to execute tthe GCG program WORDSEARCH in the batch 
$	! stream. THis procedure is interactive and asks the user for the
$	! various parameters.
$
$	on control_y then goto terminate
$	bell[0,7] = 7
$	ws := "write sys$output"
$	iq := inquire/nopunctuation
$
$	ws ""
$	ws "WORDSEARCH does a Wilbur and Lipman search for similarity between a"
$	ws "new sequence and any group of sequences. WORDSEARCH answers the"
$	ws "question 'What sequences in the database are similar to my sequence?'"
$
$get_query:
$
$	mindiags = 10 ! minimum number of matches to list
$	mode := "BATCH"
$	ws ""
$	iq seqname "WORDSEARCH with what query sequence? "
$	if(seqname.EQS."") then goto get_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_type
$	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_type:
$
$	! Determine if this is a DNA or protein sequence.
$
$	type := NUCLEOTIDE
$	if(seqinfotype.EQS."PROTEIN") then type := PEPTIDE
$
$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 NBRF database
$
$	userfile := ""
$	ws ""
$	ws "Database to search:"
$	ws ""
$	ws "1) GenEMBL"
$	ws "2) GenBank"
$	ws "3) EMBL"
$	ws "4) Swiss-Prot"
$	ws "5) NBRF (protein)"
$	ws "6) NBRF (nucleic)"
$	ws "7) VecBase"
$	ws "8) Your own sequence(s)"
$	ws ""
$	if(type.EQS."NUCLEOTIDE") then iq choice "Please enter choice (* 1 *): "
$	if(type.EQS."PEPTIDE") then iq choice "Please enter choice (* 4 *): "
$	if((type.EQS."NUCLEOTIDE").AND.(choice.EQS."")) then choice := 1
$	if((type.EQS."PEPTIDE").AND.(choice.EQS."")) then choice := 4
$	if(choice.EQS."1") then database := genEMBL
$	if(choice.EQS."2") then database := genbank
$	if(choice.EQS."3") then database := EMBL
$	if(choice.EQS."4") then database := swissprot
$	if(choice.EQS."5") then database := NBRF
$	if(choice.EQS."6") then database := Nucleic
$	if(choice.EQS."7") then database := VecBase
$	if(choice.EQS."8") then goto get_userfile
$	if("''database'".EQS."GENBANK") then goto get_genset
$	if("''database'".NES."") then goto ask_subset
$	ws "''bell'Valid choices are 1,2,3,4,5,6,7 or 8."
$	goto get_database
$
$get_userfile:
$
$	ws ""
$	iq userfile "Sequence name(s): "
$	if(userfile.EQS."") then goto get_userfile
$
$	! The filename can be a list of filenames, in which case it will
$	! start with @. If so, strip itoff,then check for the existence of
$	! the file.
$
$	list := FALSE
$	if(f$extract(0,1,userfile).EQS."@") then list := TRUE
$	if(list.EQS."TRUE") then userfile = userfile - "@"
$	if(f$search(userfile).NES."") then goto fix_filename
$	ws "''bell'''userfile' doesn't exist. Please try again."
$	goto get_userfile
$
$fix_filename:
$
$	if(list.EQS."TRUE") then userfile := @'userfile'
$	if(list.EQS."FALSE") then mindiags =1
$
$get_mode:
$
$	! For user-defined databases, the search can be done interactively
$
$	ws ""
$	ws "How do you want the search done?"
$	ws ""
$	ws "1) Interactively (faster but more expensive)"
$	ws "2) In Batch (might be slower, but cheaper)"
$	ws ""
$	iq choice "Please enter your choice (* 2 *): "
$	if(choice.EQS."") then choice := 2
$	if(choice.EQS."1") then mode := INTERACTIVE
$	if(choice.EQS."2") then mode := BATCH
$	if(mode.NES."") then goto ask_mask
$	ws "''bell'Please enter 1 or 2."
$	goto get_mode
$
$get_genset:
$
$	! Find out which of the GenBank databases should be searched, or
$	! if all of them should be.
$
$	ws ""
$	ws "GenBank sequences to search:"
$	ws ""
$	ws " 1) All of them"
$	ws " 2) Primate sequences"
$	ws " 3) Rodent sequences"
$	ws " 4) Other Mammalian sequences"
$	ws " 5) Other Vertebrate sequences"
$	ws " 6) Invertebrate sequences"
$	ws " 7) Plant sequences"
$	ws " 8) Bacterial sequences"
$	ws " 9) Organelle sequences"
$	ws "10) Phage sequences"
$	ws "11) Viral sequences"
$	ws "12) Structural RNA sequences"
$	ws "13) Synthetic sequences"
$	ws "14) Unannotated sequences"
$	ws ""
$	iq choice "Please enter choice (* 1 *): "
$	if(choice.EQS."") then choice := 1
$	database := ""
$	if(choice.EQS."1") then database := GENBANK
$	if(choice.EQS."2") then database := PRIMATE
$	if(choice.EQS."3") then database := RODENT
$	if(choice.EQS."4") then database := OTHERMAMMAL
$	if(choice.EQS."5") then database := OTHERVERTEBRATE
$	if(choice.EQS."6") then database := INVERTEBRATE
$	if(choice.EQS."7") then database := PLANT
$	if(choice.EQS."8") then database := BACTERIAL
$	if(choice.EQS."9") then database := ORGANELLE
$	if(choice.EQS."10") then database := PHAGE
$	if(choice.EQS."11") then database := VIRAL
$	if(choice.EQS."12") then database := STRUCTURAL
$	if(choice.EQS."13") then database := SYNTHETIC
$	if(choice.EQS."14") then database := UNANNOTATED
$	if(database.NES."") then goto ask_subset
$	ws "''bell'Valid responses are 1 - 14, inclusive."
$	goto get_genset
$
$ask_subset:
$
$	! Find out if a subset of the database should be searched
$
$	subset := "*"
$	ws ""
$	ws "How much of ''database' do you want to search?"
$	ws ""
$	ws "1) All of it."
$	ws "2) Only part of it."
$	ws ""
$	iq choice "Please enter your choice (* 1 *): "
$	if(choice.EQS."") then goto ask_mask
$	if(choice.EQS."1") then goto ask_mask
$	if(choice.EQS."2") then goto get_subset
$	ws "''bell'Please enter 1 or 2."
$	goto ask_subset
$
$get_subset:
$
$	! Find out what subset of the database to search
$
$	ws ""
$	iq subset "Specify the subset of the database to search (eg HUM*): "
$	if(subset.NES."") then goto ask_mask
$	ws "''bell'Invalid answer."
$	ws "If you want to search the whole database, type *"
$	goto get_subset
$
$ask_mask:
$
$	! Find out if a mask should be used
$
$	mask := ""
$	ws ""
$	iq choice "Do you want to use a mask (* No *)? "
$	if(choice.EQS."") then goto get_wordsize
$	choice = f$extract(0, 1, choice)
$	if(choice.EQS."Y") then goto get_mask
$	if(choice.EQS."N") then goto get_wordsize
$	ws "''bell'Please answer Yes or No."
$	goto ask_mask
$
$get_mask:
$
$	! Get the mask to be used
$
$	ws ""
$	if(type.EQS."NUCLEOTIDE") then iq mask "What word mask (* ++-++-++ *)? "
$	if(type.EQS."PEPTIDE") then iq mask "What word mask? "
$	if((type.EQS."NUCLEOTIDE").AND.(mask.EQS."")) then mask := "++-++-++"
$	if(mask.EQS."") then goto ask_mask ! Give a chance to back out
$	goto get_ndiag
$
$get_wordsize:
$
$	! Find out how long the word should be
$
$	wordsize := ""
$	ws ""
$	if(type.EQS."NUCLEOTIDE") then iq wordsize "What word size (* 6 *)? "
$	if(type.EQS."PEPTIDE") then iq wordsize "What word size (* 2 *)? "
$	if((type.EQS."NUCLEOTIDE").AND.(wordsize.EQS."")) then wordsize := 6
$	if((type.EQS."PEPTIDE").AND.(wordsize.EQS."")) then wordsize := 2
$	if((f$integer(wordsize).GT.0).AND.(f$integer(wordsize).LT.31)) -
			then goto get_ndiag
$	ws "''bell'Word size must be in the range from 1 to 30."
$	goto get_wordsize
$
$get_ndiag:
$
$	ws ""
$	iq ndiag "List how many best diagonals (* 20 *)? "
$	if(ndiag.EQS."") then ndiag := 20
$	if(f$integer(ndiag).GE.mindiags) then goto get_nadjacent
$	ws "''bell'At least ''mindiags' diagonals should be listed."
$	goto get_ndiag
$
$get_nadjacent:
$
$	ws ""
$	iq nadjacent "Integrate how many adjacent diagonals (* 3 *)? "
$	if(nadjacent.EQS."") then nadjacent := 3
$	if(f$integer(nadjacent).GT.0) then goto summarize
$	ws "''bell'The minimum value is 1."
$	goto get_nadjacent
$
$summarize:
$
$	if(userfile.EQS."") then subset := ":''subset'"
$	if(userfile.NES."") then subset := ""
$	if(userfile.NES."") then database := 'userfile'
$	database := "''database'''subset'"
$	ws ""
$	ws ""
$	ws "The following WORDSEARCH will be executed:"
$	ws ""
$	ws "Query sequence: ''seqname' (''type')"
$	ws "Database to be searched: ''database'"
$	if("''mask'".EQS."") then ws "No mask"
$	if("''mask'".NES."") then ws "Mask: ''mask'"
$	if("''mask'".EQS."") then ws "Word size: ''wordsize'"
$	ws "Diagonals to list: ''ndiag'"
$	ws "Integrating ''nadjacent' adjacent diagonals"
$	if(mode.EQS."INTERACTIVE") then -
   		ws "The search will be done INTERACTIVELY"
$	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:
$
$	c_dir = f$environment("default") ! Current directory
$
$	! Determine the root name of the sequence for comparison. This in used
$	! in specifying the process name of the batch job, and the output
$	! file name. Also check that another process with the same name isn't
$	! running. If it is, the WORDSEARCH and SEGMENTS files could get 
$	! mixed up, so change the process name and the output file name for
$	! WORDSEARCH, and the corresponding input file name for SEGMENTS.
$
$	! 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.
$
$	version = 0
$	seqroot = seqname
$	root = f$parse(seqroot,,,"NAME") ! root name of sequence
$	if(root.NES."") then goto check_process
$
$	! Remove database name
$
$	pos = 'f$locate(":", seqroot)'
$	len = 'f$length(seqroot)'
$	if(pos.NE.len) then root = f$extract(pos+1, len-pos, seqroot)
$
$ check_process:
$
$	version = version + 1
$	procname := "WRD''f$string(version)'_''root'"
$
$	! The process name must be no more than 15 characters long for this
$	! routine to work properly.
$
$	procname = f$extract(0, 15, procname)
$	@clark$commands:chk_process "''procname'"
$	if(exist_process.eq.1) then goto check_process
$
$	! Now set up the switches before executing the command
$
$	if(mask.EQS."") then wordsize := "/WOR="'wordsize'
$	if(mask.NES."") then mask := "/MAS="'mask'
$	if(mask.NES."") then wordsize := ""
$	seqname := "/INFILE1="'seqname'
$	if(type.EQS."PEPTIDE") then seqname := 'seqname'"/PRO"
$	if(type.EQS."NUCLEOTIDE") then seqname := 'seqname'"/NOPRO"
$	database := "/INFILE2="'database'
$	ndiag := "/LIS="'ndiag'
$	nadjacent := "/INT="'nadjacent'
$	if(version.eq.1) then version := ""
$	newroot := 'root''version'
$
$	! Create the command file that will be submitted to the batch
$	! processor. This procedure also runs SEGMENTS and PLOT to show
$	! the histograms.
$
$	open/write comfile wordsearch.cmd
$	wc := "write comfile"
$	wc "$ set noon"
$	if(mode.EQS."BATCH") then wc "$ set process/name=''procname'"
$	wc "$ set on"
$	wc "$ set default ''c_dir'"
$	wc "$ wsearch := $genutil:wordsearch" ! Won't work directly
$	wc "$ wsearch ''seqname'''database'''mask' -"
$	wc "          ''wordsize'''ndiag'''nadjacent' -"
$	wc "          /OUT=''newroot'.WRD/nomonitor"
$	wc "$ if($severity.EQS."0") then exit ! Program crapped out
$
$	! Run SEGMENTS. The output filename has the extension .SEG
$	! If the sequence is a protein, use the local datafile
$	! DNA$CMP:SEGPEP.HI, which requires strict matches.
$
$	datafile := ""
$	if(type.EQS."PEPTIDE") then datafile :="/DATA=DNA$CMP:SEGPEP.HI"
$	wc "$ sg := $genutil:segments"
$	wc "$ sg /WIDTH=80/INFILE=''newroot'.WRD''DATAFILE' -"
$	wc "          /OUTFILE=''newroot'.SEG/NOMONITOR"
$
$	close comfile
$
$	if(mode.EQS."BATCH") then goto batch_mode
$
$	@wordsearch.cmd
$	ws ""
$	delete wordsearch.cmd;0
$	exit
$
$batch_mode:
$
$	! A friendly message for the user
$
$	ws ""
$	ws "The name of the file containing the aligned sequences will be"
$	ws "''newroot'.SEG. You will be notified when this is ready."
$	ws "Thank-you for your patronage."
$	ws ""
$
$	submit/notify/noprinter/name='procname' wordsearch.cmd/delete
$
$	exit
$
$ terminate:  ! Jump here on ^y
$
$	if(f$search("wordsearch.cmd").EQS."") then exit
$	close comfile
$	delete wordsearch.cmd;*
$
$	exit
********************
$	! FIND.COM
$
$	! May 1, 1990
$
$	! Written by Steve Clark
$	! Samuel Lunenfeld Research Institute
$	! Mt. Sinai Hospital, Toronto, Canada
$
$	! Command procedure to execute the GCG program FIND in the batch 
$	! stream. This procedure is interactive and asks the user for the
$	! various parameters.
$
$	on control_y then goto terminate
$	bell[0,7] = 7
$	ws := "write sys$output"
$	iq := inquire/nopunctuation
$
$	ws ""
$	ws "FIND looks for patterns like GAATTC or YRYRYRYR in sequences. You can"
$	ws "define the patterns ambiguously and allow mismatches. You can provide"
$	ws "the patterns in a file or simply type them in from the terminal."
$	ws ""
$
$ ask_findseq:
$
$	mode := "BATCH"
$	ws ""
$	ws "Get the sequence to FIND from"
$	ws ""
$	ws "1) The terminal"
$	ws "2) A file"
$	ws ""
$	iq choice "Please enter choice (* 1 *): "
$	if(choice.EQS."") then choice := 1
$	if(choice.EQS."1") then goto get_findseq
$	if(choice.EQS."2") then goto get_findfile
$	ws "''bell'Valid choices are 1 and 2. Please try again."
$	goto ask_findseq
 $
$ get_findseq:
$
$	! Ask the user to type in the sequence to be found, one at a time.
$	! The individual patterns are concatenated into one string called
$	! FINDSEQ.
$
$	ws ""
$	ws "Enter the search patterns individually, one per line."
$	ws "End the list with a blank line."
$	ws ""
$	ipattern = 1
$	source := "/PATTERN="
$
$startloop:
$
$		temp := pattern'ipattern'
$		iq 'temp' "Pattern ''ipattern': "
$		if('temp'.EQS."") then goto endloop
$		if(ipattern.EQ.1) then findseq := p,'temp'
$		if(ipattern.GT.1) then findseq := 'findseq',q,'temp'
$		ipattern = ipattern + 1
$		if(ipattern.GT.5) then goto endloop ! Can do 5 patterns max
$		goto startloop
$
$endloop:
$
$	findseq:='findseq',p
$	ipattern = ipattern - 1
$	if(ipattern.GT.0) then goto get_database
$	ws "''bell'It's too easy to find nothing!"
$	goto get_findseq
$
$ get_findfile:
$
$	! Ask for the name of the of the file that contains sequences to search for
$
$	findseq := temp
$	source := "/DATA="
$	ws ""
$	ws "Name of the file containing the"
$	iq temp "sequence(s) to FIND (* FIND.DAT *): "
$	if(temp.EQS."") then temp := FIND.DAT
$	if(f$search(temp).NES."") then goto get_database
$	ws "''bell'''temp' doesn't exist. Please try again."
$	goto get_findfile
$
$get_database:
$
$	! Find out which database to search.
$
$	userfile := ""
$	type := NUCLEOTIDE
$	ws ""
$	ws "Database to search:"
$	ws ""
$	ws "1) GenEMBL"
$	ws "2) GenBank"
$	ws "3) EMBL"
$	ws "4) Swiss-Prot"
$	ws "5) NBRF (Protein)"
$	ws "6) NBRF (Nucleic)"
$	ws "7) VecBase"
$	ws "8) Your own sequence(s)"
$	ws ""
$	iq choice "Please enter choice (* 1 *): "
$	if(choice.EQS."") then database := GenEMBL
$	if(choice.EQS."1") then database := GenEMBL
$	if(choice.EQS."2") then database := genbank
$	if(choice.EQS."3") then database := EMBL
$	if(choice.EQS."4") then database := SwissProt
$	if(choice.EQS."5") then database := NBRF
$	if(choice.EQS."6") then database := Nucleic
$	if(choice.EQS."7") then database := VecBase
$	if(choice.EQS."8") then goto get_userfile
$	if(choice.EQS."4") then type := PEPTIDE
$	if(choice.EQS."5") then type := PEPTIDE
$	if("''database'".EQS."GENBANK") then goto get_genset
$	if("''database'".NES."") then goto ask_subset
$	ws "''bell'Valid choices are 1,2,3,4,5,6,7 or 8."
$	goto get_database
$
$get_userfile:
$
$	ws ""
$	iq userfile "Sequence name(s): "
$	if(userfile.EQS."") then goto get_userfile
$
$	! The filename can be a list of filenames, in which case it will
$	! start with @. If so, strip itoff,then check for the existence of
$	! the file.
$
$	list := FALSE
$	if(f$extract(0,1,userfile).EQS."@") then list := TRUE
$	if(list.EQS."TRUE") then userfile = userfile - "@"
$	if(f$search(userfile).NES."") then goto fix_filename
$	ws "''bell'''userfile' doesn't exist. Please try again."
$	goto get_userfile
$
$fix_filename:
$
$	if(list.EQS."TRUE") then userfile := @'userfile'
$
$get_type:
$
$	! Ask if the sequences are DNA or protein
$	ws ""
$	ws "Are these sequences"
$	ws ""
$	ws "1) DNA / RNA"
$	ws "2) protein"
$	ws ""
$	iq choice "Please enter choice: "
$	if(choice.EQS."2") then type := PEPTIDE
$	if((choice.EQS."1") .OR. (choice.EQS."2")) then goto get_mode
$	ws "''bell'Enter either 1 or 2."
$	goto get_type
$
$ get_mode:
$
$	! For user-defined databases, the search can be done interactively
$
$	ws ""
$	ws "How do you want the search done?"
$	ws ""
$	ws "1) Interactively (faster but more expensive)"
$	ws "2) In Batch (might be slower, but cheaper)"
$	ws ""
$	iq choice "Please enter your choice (* 2 *): "
$	if(choice.EQS."") then choice := 2
$	if(choice.EQS."1") then mode := INTERACTIVE
$	if(choice.EQS."2") then mode := BATCH
$	if(mode.NES."") then goto get_strand
$	ws "''bell'Please enter 1 or 2."
$	goto get_mode
$
$get_genset:
$
$	! Find out which of the GenBank databases should be searched, or
$	! if all of them should be.
$
$	ws ""
$	ws "GenBank sequences to search:"
$	ws ""
$	ws " 1) All of them"
$	ws " 2) Primate sequences"
$	ws " 3) Rodent sequences"
$	ws " 4) Other Mammalian sequences"
$	ws " 5) Other Vertebrate sequences"
$	ws " 6) Invertebrate sequences"
$	ws " 7) Plant sequences"
$	ws " 8) Bacterial sequences"
$	ws " 9) Organelle sequences"
$	ws "10) Phage sequences"
$	ws "11) Viral sequences"
$	ws "12) Structural sequences"
$	ws "13) Synthetic sequences"
$	ws "14) Unannotated sequences"
$	ws ""
$	iq choice "Please enter choice (* 1 *): "
$	if(choice.EQS."") then choice := 1
$	database := ""
$	if(choice.EQS."1") then database := GENBANK
$	if(choice.EQS."2") then database := PRIMATE
$	if(choice.EQS."3") then database := RODENT
$	if(choice.EQS."4") then database := OTHERMAMMAL
$	if(choice.EQS."5") then database := OTHERVERTEBRATE
$	if(choice.EQS."6") then database := INVERTEBRATE
$	if(choice.EQS."7") then database := PLANT
$	if(choice.EQS."8") then database := BACTERIAL
$	if(choice.EQS."9") then database := ORGANELLE
$	if(choice.EQS."10") then database := PHAGE
$	if(choice.EQS."11") then database := VIRAL
$	if(choice.EQS."12") then database := STRUCTURAL
$	if(choice.EQS."13") then database := SYNTHETIC
$	if(choice.EQS."14") then database := UNANNOTATED
$	if(database.NES."") then goto ask_subset
$	ws "''bell'Valid responses are 1 - 14, inclusive."
$	goto get_genset
$
$ask_subset:
$
$	! Find out if a subset of the database should be searched
$
$	subset := "*"
$	ws ""
$	ws "How much of ''database' do you want to search?"
$	ws ""
$	ws "1) All of it."
$	ws "2) Only part of it."
$	ws ""
$	iq choice "Please enter your choice (* 1 *): "
$	if(choice.EQS."") then goto get_strand
$	if(choice.EQS."1") then goto get_strand
$	if(choice.EQS."2") then goto get_subset
$	ws "''bell'Please enter 1 or 2."
$	goto ask_subset
$
$get_subset:
$
$	! Find out what subset of the database to search
$
$	ws ""
$	iq subset "Specify the subset of the database to search (eg HUM*): "
$	if(subset.NES."") then goto get_strand
$	ws "''bell'Invalid answer."
$	ws "If you want to search the whole database, type *"
$	goto get_subset
$
$get_strand:
$
$	! Search one or both strands?
$
$	strand = 1
$	if(type.EQS."PEPTIDE") then goto get_mismatch
$
$	ws ""
$	ws "Which strands should be searched?"
$	ws ""
$	ws "1) Top strand only"
$	ws "2) Both strands"
$	ws ""
$	iq strand "Please enter your choice (* 2 *): "
$	if(strand.EQS."") then strand := 2
$	if(strand.EQS."1") then goto get_mismatch
$	if(strand.EQS."2") then goto get_mismatch
$	ws "''bell'Please enter 1 or 2."
$	goto get_strand
$
$get_mismatch:
$
$	! Find out if any mismatches should be allowed
$
$	ws ""
$	iq nmismatch "Number of allowable mismatches (* 0 *): "
$	nmismatch = f$integer(nmismatch) ! Convert to 0 or 1 if string
$
$ get_outfile:
$
$	! Determine the name of the output file
$
$	ws ""
$	iq outfile "Name of the output file (* ???.FND *): "
$	if(outfile.NES."") then goto fix_outname
$	ws "''bell'It will be easier for you if the name is longer."
$	goto get_outfile
$
$fix_outname:
$
$	! Make sure the extension is .FND. Strip out the root and add extension.
$
$	pos = f$locate(".", outfile)
$	outfile := 'f$extract(0, pos, outfile)'.FND
$
$ summarize:
$
$	if(userfile.EQS."") then subset := ":''subset'"
$	if(userfile.NES."") then subset := ""
$	if(userfile.NES."") then database := "''userfile'"
$	database := "''database'''subset'"
$	p := """
$	q := ""","""
$	ws ""
$	ws ""
$	ws "The following FIND will be executed:"
$	ws ""
$	if(source.EQS."/PATTERN=") then ws -
		"Sequence to FIND: ",'findseq'," (''type')"
$	if(source.EQS."/DATA=") then ws -
		"File of sequences to FIND: ''temp' (''type')"
$	ws "Database to be searched: ''database'"
$	ws "Number of strands to search: ''strand'"
$	ws "Number of allowable mismatches: ''nmismatch'"
$	ws "Output file name: ''outfile'"
$	if(mode.EQS."INTERACTIVE") then -
   		ws "The search will be done INTERACTIVELY."
$	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 ask_findseq
$	if(choice.eqs."1") then goto ask_findseq
$	if(choice.eqs."2") then exit
$	ws "''bell'Wasn't the question simple enough for you?"
$	goto ask_repeat
$
$do_it:
$
$	c_dir = f$environment("default") ! Current directory
$
$	! Determine the root name of the output file. This in used
$	! in specifying the process name of the batch job.
$
$	root = f$parse(outfile,,,"NAME") ! root name of sequence
$	procname := "FIND_''root'"
$
$	! Now set up the switches that will be appended to the command
$	! The individual patterns have to be enclosed in double quotes
$	! in case there are any commas in the search patterns.
$
$do_file:
$
$	if(strand.EQS."1") then strand := "/ONE"
$	if(strand.EQS."2") then strand := ""
$	database := "/INFILE=''database'"
$	outfile := "/OUTFILE=''outfile'"
$	mismatch := "/MISMATCH=''nmismatch'"
$	if(type.EQS."PEPTIDE") then type := "/PEPTIDE"
$	if(type.NES."/PEPTIDE") then type := ""
$
$	! Create the command file that will be submitted to the batch queue
$
$	open/write comfile find.cmd
$	wc := "write comfile"
$	wc "$ set noon"
$	if(mode.EQS."BATCH") then wc "$ set process/name=''procname'"
$	wc "$ set on"
$	wc "$ set default ''c_dir'"
$	wc "$ f := $genutil:find"	! won't work directly
$	wc "$ f",source,'findseq'," -"
$	wc "   ",database,outfile,mismatch,type,strand,"/NOMONITOR"
$	close comfile
$
$	if(mode.EQS."BATCH") then goto batch_mode
$
$	@find.cmd
$	delete find.cmd;0
$	exit
$
$batch_mode:
$
$	submit/notify/noprinter/name='procname' find.cmd/delete
$	exit
$
$ terminate:  ! Jump here on ^y
$
$	if("''f$search("find.cmd")'".EQS."") then exit
$	close comfile
$	delete find.cmd;*
$
$	exit
********************
$	! FASTA.COM
$
$	! May 1, 1990
$
$	! Written by Steve Clark
$	! Mt. Sinai Hospital Research Institute, Toronto, Canada
$
$	! Command procedure to execute tthe GCG program FASTA in the batch 
$	! stream. This procedure is interactive and asks the user for the
$	! various parameters.
$
$	on control_y then goto terminate
$	bell[0,7] = 7
$	ws := "write sys$output"
$	iq := inquire/nopunctuation
$
$	ws ""
$	ws "FASTA does a Pearson and Lipman search for similarity between a"
$	ws "query sequence and any group of sequences. FASTA answers the"
$	ws "question, What sequences in the database are similar to my sequence?"
$	ws "The relationship between FASTA and WORDSEARCH has not been"
$	ws "characterized, but FASTA is faster, and, for some searches, more"
$	ws "sensitive."
$
$get_query:
$
$	ws ""
$	iq seqname "FASTA with what query sequence? "
$	if(seqname.EQS."") then goto get_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_type
$	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_type:
$
$	! Determine if this is a DNA or protein sequence.
$
$	type := NUCLEOTIDE
$	if(seqinfotype.EQS."PROTEIN") then type := PEPTIDE
$
$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 NBRF database
$
$	userfile := ""
$	ws ""
$	ws "Database to search:"
$	ws ""
$	ws "1) GenEMBL"
$	ws "2) GenBank"
$	ws "3) EMBL"
$	ws "4) Swiss-Prot"
$	ws "5) NBRF (protein)"
$	ws "6) NBRF (nucleic)"
$	ws "7) VecBase"
$	ws "8) Your own sequence(s)"
$	ws ""
$	if(type.EQS."NUCLEOTIDE") then iq choice "Please enter choice (* 1 *): "
$	if(type.EQS."PEPTIDE") then iq choice "Please enter choice (* 4 *): "
$	if((type.EQS."NUCLEOTIDE").AND.(choice.EQS."")) then choice := 1
$	if((type.EQS."PEPTIDE").AND.(choice.EQS."")) then choice := 4
$	if(choice.EQS."1") then database := genEMBL
$	if(choice.EQS."2") then database := genbank
$	if(choice.EQS."3") then database := EMBL
$	if(choice.EQS."4") then database := SwissProt
$	if(choice.EQS."5") then database := NBRF
$	if(choice.EQS."6") then database :=Nucleic
$	if(choice.EQS."7") then database := VecBase
$	if(choice.EQS."8") then goto get_userfile
$	if("''database'".EQS."GENBANK") then goto get_genset
$	if("''database'".NES."") then goto ask_subset
$	ws "''bell'Valid choices are 1,2,3,4,5,6,7 or 8."
$	goto get_database
$
$get_userfile:
$
$	ws ""
$	iq userfile "Sequence name(s): "
$	if(userfile.EQS."") then goto get_userfile
$
$	! The filename can be a list of filenames, in which case it will
$	! start with @. If so, strip itoff,then check for the existence of
$	! the file.
$
$	list := FALSE
$	if(f$extract(0,1,userfile).EQS."@") then list := TRUE
$	if(list.EQS."TRUE") then userfile = userfile - "@"
$	if(f$search(userfile).NES."") then goto fix_filename
$	ws "''bell'''userfile' doesn't exist. Please try again."
$	goto get_userfile
$
$fix_filename:
$
$	if(list.EQS."TRUE") then userfile := @'userfile'
$	goto get_wordsize  !  FASTA only runs in BATCH mode
$
$get_genset:
$
$	! Find out which of the GenBank databases should be searched, or
$	! if all of them should be.
$
$	ws ""
$	ws "GenBank sequences to search:"
$	ws ""
$	ws " 1) All of them"
$	ws " 2) Primate sequences"
$	ws " 3) Rodent sequences"
$	ws " 4) Other Mammalian sequences"
$	ws " 5) Other Vertebrate sequences"
$	ws " 6) Invertebrate sequences"
$	ws " 7) Plant sequences"
$	ws " 8) Bacterial sequences"
$	ws " 9) Organelle sequences"
$	ws "10) Phage sequences"
$	ws "11) Viral sequences"
$	ws "12) Structural RNA sequences"
$	ws "13) Synthetic sequences"
$	ws "14) Unannotated sequences"
$	ws ""
$	iq choice "Please enter choice (* 1 *): "
$	if(choice.EQS."") then choice := 1
$	database := ""
$	if(choice.EQS."1") then database := GENBANK
$	if(choice.EQS."2") then database := PRIMATE
$	if(choice.EQS."3") then database := RODENT
$	if(choice.EQS."4") then database := OTHERMAMMAL
$	if(choice.EQS."5") then database := OTHERVERTEBRATE
$	if(choice.EQS."6") then database := INVERTEBRATE
$	if(choice.EQS."7") then database := PLANT
$	if(choice.EQS."8") then database := BACTERIAL
$	if(choice.EQS."9") then database := ORGANELLE
$	if(choice.EQS."10") then database := PHAGE
$	if(choice.EQS."11") then database := VIRAL
$	if(choice.EQS."12") then database := STRUCTURAL
$	if(choice.EQS."13") then database := SYNTHETIC
$	if(choice.EQS."14") then database := UNANNOTATED
$	if(database.NES."") then goto ask_subset
$	ws "''bell'Valid responses are 1 - 14, inclusive."
$	goto get_genset
$
$ask_subset:
$
$	! Find out if a subset of the database should be searched
$
$	subset := "*"
$	ws ""
$	ws "How much of ''database' do you want to search?"
$	ws ""
$	ws "1) All of it."
$	ws "2) Only part of it."
$	ws ""
$	iq choice "Please enter your choice (* 1 *): "
$	if(choice.EQS."") then goto get_wordsize
$	if(choice.EQS."1") then goto get_wordsize
$	if(choice.EQS."2") then goto get_subset
$	ws "''bell'Please enter 1 or 2."
$	goto ask_subset
$
$get_subset:
$
$	! Find out what subset of the database to search
$
$	ws ""
$	iq subset "Specify the subset of the database to search (eg HUM*): "
$	if(subset.NES."") then goto get_wordsize
$	ws "''bell'Invalid answer."
$	ws "If you want to search the whole database, type *"
$	goto get_subset
$
$get_wordsize:
$
$	! Find out how long the word should be
$
$	wordsize := ""
$	ws ""
$	if(type.EQS."NUCLEOTIDE") then iq wordsize "What word size (* 6 *)? "
$	if(type.EQS."PEPTIDE") then iq wordsize "What word size (* 2 *)? "
$	if((type.EQS."NUCLEOTIDE").AND.(wordsize.EQS."")) then wordsize := 6
$	if((type.EQS."PEPTIDE").AND.(wordsize.EQS."")) then wordsize := 2
$	if((f$integer(wordsize).GT.0).AND.(f$integer(wordsize).LT.7)) -
			then goto get_nscores
$	ws "''bell'Word size must be in the range from 1 to 6."
$	goto get_wordsize
$
$get_nscores:
$
$	ws ""
$	iq nscores "List how many best scores and alignments (* 20 *)? "
$	if(nscores.EQS."") then nscores := 20
$	if(f$integer(nscores).GE.10) then goto summarize
$	ws "''bell'At least 10 scores should be listed."
$	goto get_nscores
$
$summarize:
$
$	if(userfile.EQS."") then subset := ":''subset'"
$	if(userfile.NES."") then subset := ""
$	if(userfile.NES."") then database := 'userfile'
$	database := "''database'''subset'"
$	ws ""
$	ws ""
$	ws "The following FASTA search will be executed:"
$	ws ""
$	ws "Query sequence: ''seqname' (''type')"
$	ws "Database to be searched: ''database'"
$	ws "Word size: ''wordsize'"
$	ws "Alignments to list: ''nscores'"
$	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:
$
$	c_dir = f$environment("default") ! Current directory
$
$	! Determine the root name of the sequence for comparison. This in used
$	! in specifying the process name of the batch job, and the output
$	! file name, which has the extension .FA.
$
$	! 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_procname
$
$	! 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_procname:
$
$	outname := "''root'.FA"
$
$	procname := "FA_''root'"
$
$	! The process name must be no more than 15 characters long for this
$	! routine to work properly.
$
$	procname = f$extract(0, 15, procname)
$
$	! Now set up the switches before executing the command
$
$	wordsize := "/WORDSIZE="'wordsize'
$	seqname := "/INFILE1="'seqname'
$	if(type.EQS."PEPTIDE") then seqname := 'seqname'"/PRO"
$	if(type.EQS."NUCLEOTIDE") then seqname := 'seqname'"/NOPRO"
$	database := "/INFILE2="'database'
$	nscores := "/LIS="'nscores'
$
$	! Create the command file that will be submitted to the batch
$	! processor. 
$
$	open/write comfile fasta.cmd
$	wc := "write comfile"
$	wc "$ set noon"
$	wc "$ set process/name=''procname'"
$	wc "$ set on"
$	wc "$ set default ''c_dir'"
$	wc "$ fast := $genutil:fasta" ! Won't work directly
$	wc "$ fast ''seqname'''database' -"
$	wc "          ''wordsize'''nscores' -"
$	wc "          /OUTFILE=''outname'/NODOCL/NOMONITOR/LINESIZE=80"
$	wc "" ! Response to "How many scores would you like to see?"
$
$	close comfile
$
$	! A friendly message for the user
$
$	ws ""
$	ws "The name of the file containing the aligned sequences will be"
$	ws "''outname' You will be notified when this is ready."
$	ws "Thank-you for your patronage."
$	ws ""
$
$	! FASTA really hogs system resources, so it has a queue of its own,
$	! job_limit=1
$
$	submit/que=fasta$batch/notify/noprinter/name='procname' fasta.cmd/delete
$
$	exit
$
$ terminate:  ! Jump here on ^y
$
$	if(f$search("fasta.cmd").EQS."") then exit
$	close comfile
$	delete fasta.cmd;*
$
$	exit
********************
$	! TFASTA.COM
$
$	! May 1, 1990
$
$	! Written by Steve Clark
$	! Samuel Lunenfeld Research Institute
$	! Mt. Sinai Hospital, Toronto, Canada
$
$	! Command procedure to execute tthe GCG program FASTA in the batch 
$	! stream. This procedure is interactive and asks the user for the
$	! various parameters.
$
$	on control_y then goto terminate
$	bell[0,7] = 7
$	ws := "write sys$output"
$	iq := inquire/nopunctuation
$
$	ws ""
$	ws "TFASTA does a Pearson and Lipman search for similarity between a"
$	ws "query peptide sequence and any group of nucleotide sequences. TFASTA"
$	ws "translates the nucleotide sequence in all six reading frames before"
$	ws "performing the comparison. It is designed to answer the question,"
$	ws "What implied peptide sequences in a nucleotide sequence database"
$	ws "are similar to my peptide sequence?"
$
$get_query:
$
$	ws ""
$	iq seqname "TFASTA with what query sequence? "
$	if(seqname.EQS."") then goto get_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_type
$	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_type:
$
$	! If this is a DNA sequence, quit the program.
$
$	if(seqinfotype.EQS."PROTEIN") then goto get_database
$	ws ""
$	ws "''bell'''seqname' is a nucleotide sequence."
$	ws "Tfasta requires a PROTEIN query sequence."
$	exit
$
$get_database:
$
$	! Find out which database to search.
$
$	userfile := ""
$	ws ""
$	ws "Database to search:"
$	ws ""
$	ws "1) GenEMBL"
$	ws "2) GenBank"
$	ws "3) EMBL"
$	ws "4) NBRF (nucleic)"
$	ws "5) VecBase"
$	ws "6) Your own sequence(s)"
$	ws ""
$	iq choice "Please enter choice (* 1 *): "
$	if(choice.EQS."") then choice := 1
$	if(choice.EQS."1") then database := genEMBL
$	if(choice.EQS."2") then database := genbank
$	if(choice.EQS."3") then database := EMBL
$	if(choice.EQS."4") then database := Nucleic
$	if(choice.EQS."5") then database := VecBase
$	if(choice.EQS."6") then goto get_userfile
$	if("''database'".EQS."GENBANK") then goto get_genset
$	if("''database'".NES."") then goto ask_subset
$	ws "''bell'Valid choices are 1,2,3,4,5 or 6."
$	goto get_database
$
$get_userfile:
$
$	ws ""
$	iq userfile "Sequence name(s): "
$	if(userfile.EQS."") then goto get_userfile
$
$	! The filename can be a list of filenames, in which case it will
$	! start with an @. If so, strip it off, then check for the existence
$	! of the file.
$
$	list := FALSE
$	if(f$extract(0,1,userfile).EQS."@") then list := TRUE
$	if(list.EQS."TRUE") then userfile = userfile - "@"
$	if(f$search(userfile).NES."") then goto fix_filename
$	ws "''bell'''userfile' doesn't exist. Please try again."
$	goto get_userfile
$
$fix_filename:
$
$	if(list.EQS."TRUE") then userfile := @'userfile'
$	goto get_wordsize	! TFASTA only runs in batch mode
$
$get_genset:
$
$	! Find out which of the GenBank databases should be searched, or
$	! if all of them should be.
$
$	ws ""
$	ws "GenBank sequences to search:"
$	ws ""
$	ws " 1) All of them"
$	ws " 2) Primate sequences"
$	ws " 3) Rodent sequences"
$	ws " 4) Other Mammalian sequences"
$	ws " 5) Other Vertebrate sequences"
$	ws " 6) Invertebrate sequences"
$	ws " 7) Plant sequences"
$	ws " 8) Bacterial sequences"
$	ws " 9) Organelle sequences"
$	ws "10) Phage sequences"
$	ws "11) Viral sequences"
$	ws "12) Structural RNA sequences"
$	ws "13) Synthetic sequences"
$	ws "14) Unannotated sequences"
$	ws ""
$	iq choice "Please enter choice (* 1 *): "
$	if(choice.EQS."") then choice := 1
$	database := ""
$	if(choice.EQS."1") then database := GENBANK
$	if(choice.EQS."2") then database := PRIMATE
$	if(choice.EQS."3") then database := RODENT
$	if(choice.EQS."4") then database := OTHERMAMMAL
$	if(choice.EQS."5") then database := OTHERVERTEBRATE
$	if(choice.EQS."6") then database := INVERTEBRATE
$	if(choice.EQS."7") then database := PLANT
$	if(choice.EQS."8") then database := BACTERIAL
$	if(choice.EQS."9") then database := ORGANELLE
$	if(choice.EQS."10") then database := PHAGE
$	if(choice.EQS."11") then database := VIRAL
$	if(choice.EQS."12") then database := STRUCTURAL
$	if(choice.EQS."13") then database := SYNTHETIC
$	if(choice.EQS."14") then database := UNANNOTATED
$	if(database.NES."") then goto ask_subset
$	ws "''bell'Valid responses are 1 - 14, inclusive."
$	goto get_genset
$
$ask_subset:
$
$	! Find out if a subset of the database should be searched
$
$	subset := "*"
$	ws ""
$	ws "How much of ''database' do you want to search?"
$	ws ""
$	ws "1) All of it."
$	ws "2) Only part of it."
$	ws ""
$	iq choice "Please enter your choice (* 1 *): "
$	if(choice.EQS."") then goto get_wordsize
$	if(choice.EQS."1") then goto get_wordsize
$	if(choice.EQS."2") then goto get_subset
$	ws "''bell'Please enter 1 or 2."
$	goto ask_subset
$
$get_subset:
$
$	! Find out what subset of the database to search
$
$	ws ""
$	iq subset "Specify the subset of the database to search (eg HUM*): "
$	if(subset.NES."") then goto get_wordsize
$	ws "''bell'Invalid answer."
$	ws "If you want to search the whole database, type *"
$	goto get_subset
$
$get_wordsize:
$
$	! Find out how long the word should be
$
$	wordsize := ""
$	ws ""
$	iq wordsize "What word size (* 2 *) ? "
$	if(wordsize.EQS."") then wordsize := 2
$	if((f$integer(wordsize).GT.0).AND.(f$integer(wordsize).LT.7)) -
			then goto get_nscores
$	ws "''bell'Word size must be in the range from 1 to 6."
$	goto get_wordsize
$
$get_nscores:
$
$	ws ""
$	iq nscores "List how many best scores and alignments (* 20 *) ? "
$	if(nscores.EQS."") then nscores := 20
$	if(f$integer(nscores).GE.10) then goto summarize
$	ws "''bell'At least 10 scores should be listed."
$	goto get_nscores
$
$summarize:
$
$	if(userfile.EQS."") then subset := ":''subset'"
$	if(userfile.NES."") then subset := ""
$	if(userfile.NES."") then database := 'userfile'
$	database := "''database'''subset'"
$	ws ""
$	ws ""
$	ws "The following TFASTA search will be executed:"
$	ws ""
$	ws "Query sequence: ''seqname'"
$	ws "Database to be searched: ''database'"
$	ws "Word size: ''wordsize'"
$	ws "Alignments to list: ''nscores'"
$	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:
$
$	c_dir = f$environment("default") ! Current directory
$
$	! Determine the root name of the sequence for comparison. This in used
$	! in specifying the process name of the batch job, and the output
$	! file name, which has the extension .TFA.
$
$	! 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_procname
$
$	! Remove the database name
$
$	pos = f$locate(":", seqroot)
$	len = f$length(seqroot)
$	if(pos.NE.len) then root = f$extract(pos+1, len-pos, seqroot)
$
$set_procname:
$
$	outname := "''root'.TFA"
$
$	procname := "TFA_''root'"
$
$	! The process name must be no more than 15 characters long for this
$	! routine to work properly.
$
$	procname = f$extract(0, 15, procname)
$
$	! Now set up the switches before executing the command
$
$	wordsize := "/WORDSIZE="'wordsize'
$	seqname := "/INFILE1="'seqname'
$	database := "/INFILE2="'database'
$	nscores := "/LIS="'nscores'
$
$	! Create the command file that will be submitted to the batch
$	! processor. 
$
$	open/write comfile tfasta.cmd
$	wc := "write comfile"
$	wc "$ set noon"
$	wc "$ set process/name=''procname'"
$	wc "$ set on"
$	wc "$ set default ''c_dir'"
$	wc "$ tfast := $genutil:fasta" ! Won't work directly
$	wc "$ tfast ''seqname'''database' -"
$	wc "          ''wordsize'''nscores'/OUTFILE=''outname' -"
$	wc "          /TRANSLATE/NODOCL/NOMONITOR/LINESIZE=80"
$	wc "" ! Response to "How many scores would you like to see?"
$
$	close comfile
$
$	! A friendly message for the user
$
$	ws ""
$	ws "The name of the file containing the aligned sequences will be"
$	ws "''outname' You will be notified when this is ready."
$	ws "Thank-you for your patronage."
$	ws ""
$
$	! FASTA really hogs system resources, so it has a queue of its own,
$	! with job_limit=1.
$
$	submit/queue=fasta$batch/notify/noprinter/name='procname' -
		tfasta.cmd/delete
$
$	exit
$
$ terminate:  ! Jump here on ^y
$
$	if(f$search("tfasta.cmd").EQS."") then exit
$	close comfile
$	delete tfasta.cmd;*
$
$	exit
********************
$ 	! STRINGS.COM
$
$	! May 1, 1990
$
$	! Written by Steve Clark
$	! Samuel Lunenfeld Research Institute
$	! Mt. Sinai Hospital, Toronto, Canada
$
$	! This procedure helps the user execute the GCG program STRINGS.
$	! The program is run interactively, except a search of the complete
$	! sequence records is done in batch mode.
$
$	on control_y then goto terminate
$	bell[0,7] = 7
$	ws := "write sys$output"
$	iq := inquire/nopunctuation
$
$	ws ""
$	ws "STRINGS finds sequences by searching sequence library documentation"
$	ws "for character patterns such as 'globin' or 'Human'."
$	ws ""
$	ws "STRINGS has been changed to search either definitions or complete"
$	ws "sequence records. The records take 100 times as long to search as"
$	ws "the catalogue, and is done in batch mode. If you specify more than"
$	ws "one word to search for, STRINGS will only find entries that contain"
$	ws "all the specified words."
$	ws ""
$	ws "STRINGS search through which database?"
$	ws ""
$	ws "1) GenEMBL"
$	ws "2) GenBank"
$	ws "3) EMBL"
$	ws "4) Swiss-Prot
$	ws "5) NBRF (protein)"
$	ws "6) NBRF (nucleic)"
$	ws "7) VecBase"
$	ws ""
$
$get_database:
$
$	iq choice "Please choose one (* 1 *): "
$	if(choice.EQS."") then choice := 1
$	if(choice.EQS."1") then database := genEMBL
$	if(choice.EQS."2") then database := genbank
$	if(choice.EQS."3") then database := EMBL
$	if(choice.EQS."4") then database := swissprot
$	if(choice.EQS."5") then database := NBRF
$	if(choice.EQS."6") then database :=Nucleic
$	if(choice.EQS."7") then database := VecBase
$	if("''database'".EQS."GENBANK") then goto get_genset
$	if("''database'".NES."") then goto ask_subset
$	ws "''bell'Valid choices are 1,2,3,4,5,6, or 7."
$	goto get_database
$
$get_genset:
$
$	! Find out which of the GenBank databases should be searched, or
$	! if all of them should be.
$
$	ws ""
$	ws "GenBank sequences to search:"
$	ws " 1) All of them"
$	ws " 2) Primate sequences"
$	ws " 3) Rodent sequences"
$	ws " 4) Other Mammalian sequences"
$	ws " 5) Other Vertebrate sequences"
$	ws " 6) Invertebrate sequences"
$	ws " 7) Plant sequences"
$	ws " 8) Bacterial sequences"
$	ws " 9) Organelle sequences"
$	ws "10) Phage sequences"
$	ws "11) Viral sequences"
$	ws "12) Structural RNA sequences"
$	ws "13) Synthetic sequences"
$	ws "14) Unannotated sequences"
$	iq choice "Please enter choice (* 1 *): "
$	if(choice.EQS."") then choice := 1
$	database := ""
$	if(choice.EQS."1") then database := GENBANK
$	if(choice.EQS."2") then database := PRIMATE
$	if(choice.EQS."3") then database := RODENT
$	if(choice.EQS."4") then database := OTHERMAMMAL
$	if(choice.EQS."5") then database := OTHERVERTEBRATE
$	if(choice.EQS."6") then database := INVERTEBRATE
$	if(choice.EQS."7") then database := PLANT
$	if(choice.EQS."8") then database := BACTERIAL
$	if(choice.EQS."9") then database := ORGANELLE
$	if(choice.EQS."10") then database := PHAGE
$	if(choice.EQS."11") then database := VIRAL
$	if(choice.EQS."12") then database := STRUCTURAL
$	if(choice.EQS."13") then database := SYNTHETIC
$	if(choice.EQS."14") then database := UNANNOTATED
$	if(database.NES."") then goto ask_subset
$	ws "''bell'Valid responses are 1 - 14, inclusive."
$	goto get_genset
$
$ask_subset:
$
$	! Find out if a subset of the database should be searched
$
$	subset := "*"
$	ws ""
$	ws "How much of ''database' do you want to search?"
$	ws ""
$	ws "1) All of it."
$	ws "2) Only part of it."
$	ws ""
$	iq choice "Please enter your choice (* 1 *): "
$	if(choice.EQS."") then goto get_searcharea
$	if(choice.EQS."1") then goto get_searcharea
$	if(choice.EQS."2") then goto get_subset
$	ws "''bell'Please enter 1 or 2."
$	goto ask_subset
$
$get_subset:
$
$	! Find out what subset of the database to search
$
$	ws ""
$	iq subset "Specify the subset of the database to search (eg HUM*): "
$	if(subset.NES."") then goto get_searcharea
$	ws "''bell'Invalid answer."
$	ws "If you want to search the whole database, type *"
$	goto get_subset
$
$get_searcharea:
$
$	! See which part of the sequence headers to peruse.
$
$	ws ""
$	ws "Do you want to search through:"
$	ws ""
$	ws "1) Definitions (takes a few seconds while you watch)"
$	ws "2) Complete sequence records (done in batch mode)"
$	ws ""
$
$ask_searcharea:
$
$	iq choice "Please choose one (* 1 *): "
$	if(choice.EQS."") then choice := 1
$	if(choice.EQS."1") then mode := INTERACTIVE
$	if(choice.EQS."2") then mode := BATCH
$	if("''mode'".EQS."") then goto ask_searcharea ! No response
$
$get_text:
$
$	! Find out what strings to search for.
$
$	ws ""
$	iq text "Search for what text patterns? "
$	if(text.EQS."") then goto get_text
$
$	! There can't be any spaces in the search text
$
$	if(f$locate(" ",text).EQ.f$length(text)) then goto ask_outfile
$	ws "''bell'Spaces are not allowed. Use commas instead."
$	goto get_text
$
$ask_outfile:
$
$	! See if the output should be to the terminal or to a file
$
$	if(mode.EQS."BATCH") then goto get_outfile
$	outfile := "TERM:"
$	ws ""
$	ws "The outfile should go to..."
$	ws ""
$	ws "1) The terminal screen alone."
$	ws "2) A disk file as well as the terminal."
$	ws ""
$	iq choice "Please enter your choice (* 1 *): "
$	if(choice.EQS."") then goto do_it
$	if(choice.EQS."1") then goto do_it
$	if(choice.EQS."2") then goto get_outfile
$	ws "''bell'Please enter 1 or 2."
$	goto ask_outfile
$
$get_outfile:
$
$	! The default name of the output file is the name of the database with
$ 	! the extension .STR.
$
$	outfile := "''database'.STR"
$	ws ""
$	iq choice "Name of output file (* ''outfile' *): "
$	if(choice.EQS."") then choice := 'outfile'
$	outfile := 'choice'
$
$do_it:
$
$	if(mode.EQS."BATCH") then goto batch_mode
$
$	str := $genutil:strings
$	str/INFILE='database':'subset'/STRING='text'/MENU=a/OUTFILE='outfile'
$	exit
$
$batch_mode:
$
$	! Write a command procedure that can be sent to the batch queue.
$
$	c_dir = f$environment("default")	! Current directory
$	open/write comfile strings.cmd
$	wc := "write comfile"
$	wc "$ set noon"
$	wc "$ set process/name=STRINGS"
$	wc "$ set on"
$	wc "$ set default ''c_dir'"
$	wc "$ str := $genutil:strings"
$	wc "$ str/INFILE=''database':''subset'/STRING=''text'/MENU=b-"
$	wc "         /OUTFILE=''outfile'/NOSCREEN/NOMONITOR"
$	close comfile
$
$	! A friendly message.
$
$	ws ""
$	ws "You will be notified when your search is complete."
$	ws "The names will be written to the file ''outfile'"
$	ws ""
$
$	submit/notify/noprinter/name=STRINGS strings.cmd/delete
$	exit
$
$terminate:	! Jump to here on ^Y
$
$	if(f$search("strings.cmd").EQS."") then exit
$	close comfile
$	delete strings.cmd;*
$	exit
$
********************
$    ! MULTIDOT.COM

$    ! May 4, 1990

$    ! Written by Steve Clark
$    ! Samuel Lunenfeld Research Institute
$    ! Mt. Sinai Hospital, Toronto, Canada

$    ! Command procedure to run multiple sets of COMPARE/DOTPLOT with a
$    ! range of stringencies/windows or wordsizes. If a window/stingency
$    ! analysis is performed, it is executed in batch mode. A word analysis
$    ! can be executed either by batch or interactively, depending on how
$    ! the user feels about it.
$    
$    on control_y then goto terminate
$    bell[0,7] = 7
$    ws := "write sys$output"
$    iq := inquire/nopunctuation
$    
$    ws ""
$    ws "MULTIDOT runs a series of diagonal dot-matirx analyses at varying degrees"
$    ws "of similarity."
$    ws ""
$    
$get_seq1:
$    
$    ! Get the two sequences to compare. 
$    
$    ws ""
$    iq seq1 "Compare what horizontal sequence? "
$    if (seq1.EQS."") then goto get_seq1
$
$	! See if the sequence exists
$
$	assign/user_mode nl: sys$output
$	seqinfo/infile='seq1'
$	if(seqinfotype.NES."NONE") then goto check_gcg1
$	ws ""
$	ws "''bell'''seq1' doesn't exist. Please try again."
$	goto get_seq1
$
$check_gcg1:
$
$	! Check if the sequence is in GCG format.
$
$	if(seqinfotype.NES."NOGCG") then goto get_reg1
$	ws ""
$	ws "''bell'''seq1' 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_seq1
$
$get_reg1:
$
$	! Save the values of SEQINFOTYPE and SEQINFOLENGTH in case an error
$	! is made specifying the second sequence, and then the default
$	! is chosen the next time through.
$
$	seqtype := 'seqinfotype'
$	seqlength := 'seqinfolength'
$
$	! Get the regions to compare
$
$    ws ""
$    iq begin1 "Begin (* 1 *) ? "
$    if (begin1.EQS."") then begin1 := 1
$    iq end1 "End (* ''seqinfolength' *) ? "
$    if (end1.EQS."") then end1 := 'seqinfolength'
$    
$get_r1:
$    
$    iq reverse1 "Reverse (* No *) ? "
$    if (reverse1.EQS."") then reverse1 := NO
$    if (f$extract(0,1,reverse1).EQS."Y") then reverse1 := YES
$    if (f$extract(0,1,reverse1).EQS."N") then reverse1 := NO
$    if (reverse1.EQS."YES") then goto get_seq2
$    if (reverse1.EQS."NO") then goto get_seq2
$    ws bell
$    goto get_r1
$    
$get_seq2:
$    
$    ws ""
$    iq seq2 "To what vertical sequence (* ''seq1' *) ? "
$    if (seq2.EQS."") then goto same_seq
$
$	! See if the sequence exists
$
$	assign/user_mode nl: sys$output
$	seqinfo/infile='seq2'
$	if(seqinfotype.NES."NONE") then goto check_gcg2
$	ws ""
$	ws "''bell'''seq2' doesn't exist. Please try again."
$	goto get_seq2
$
$check_gcg2:
$
$	! Check if the sequence is in GCG format.
$
$	if(seqinfotype.NES."NOGCG") then goto get_reg2
$	ws ""
$	ws "''bell'''seq2' 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_seq2
$
$same_seq:
$
$	! It's not necessary to check for the validity of the second sequence
$	! if it is the same as the first.
$
$	seq2 := 'seq1'
$	seqinfotype := 'seqtype'
$	seqinfolength := 'seqlength'
$
$get_reg2:
$
$	! Get the regions to compare
$
$    ws ""
$    iq begin2 "Begin (* 1 *) ? "
$    if (begin2.EQS."") then begin2 := 1
$    iq end2 "End (* ''seqinfolength' *) ? "
$    if (end2.EQS."") then end2 := 'seqinfolength
$    
$get_r2:
$    
$    iq reverse2 "Reverse (* No *) ? "
$    if (reverse2.EQS."") then reverse2 := NO
$    if (f$extract(0,1,reverse2).EQS."Y") then reverse2 := YES
$    if (f$extract(0,1,reverse2).EQS."N") then reverse2 := NO
$    if (reverse2.EQS."YES") then goto get_root
$    if (reverse2.EQS."NO") then goto get_root
$    ws bell
$    goto get_r2
$
$get_root:
$    
$    ! Get the root for the output file name. Try to take it from the name of 
$    ! the sequence. If it is a database, it will have a colon. If not, assume
$    ! it is a file and let the lexicals get the root.
$    
$    len = f$length(seq1)
$    pos = f$locate(":", seq1)
$    if (pos.EQ.len) then goto parse_filename
$    root = f$extract(pos+1, len, seq1)
$    goto get_outword
$    
$parse_filename:
$    
$    root = f$parse(seq1,,,"NAME")
$    
$get_outword:
$    
$    ! Determine the output file name, current directory and process name for 
$    ! the batch process.
$    
$    outname := 'root'.PNT
$    c_dir = f$environment("default")
$    procname = f$extract(0,15,"MD_''root'")
$    
$    ws ""
$    iq choice "What should I call the output file (* ''outname' *) ? "
$    if (choice.NES."") then outname = choice
$    
$get_type:
$    
$    ws ""
$    ws "COMPARE should use which type of analysis?"
$    ws ""
$    ws "1) Word (Very fast - can be done in interactive mode)"
$    ws "2) Window/Stringency (Much slower but more sensitive - batch mode only)
$    ws ""
$    iq choice "Please enter you preference (* 1 *): "
$    if (choice.EQS."") then choice := 1
$    if (choice.EQS."1") then goto use_word
$    if (choice.EQS."2") then goto use_window
$    ws "''bell'Valid choices are 1 or 2. Please try again."
$    goto get_type
$    
$use_word:
$    
$    ! Find out the start and finishing word sizes, and the increment size. 
$    ! Then build a command procedure which can be submitted to the batch
$    ! queue or run interactively. It is necessary that the plotting goes to
$    ! a spooled printer or a file.
$	! The default initial word sizes are 2 for protein and 3 for DNA.
$	! The final size is three times the initial sizse, and the default
$	! increment is 1/4 of the difference.
$
$	dsize = 2
$	if(seqinfotype.EQS."DNA") then dsize = 3
$    
$get_wrd1:
$    
$	ws ""
$    iq wrd1 "Initial word size (* ''dsize' *) ? "
$	if(wrd1.EQS."") then wrd1 := 'dsize'
$    wrd1 = f$integer(wrd1)
$    if (wrd1.GT.0) then goto get_wrd2
$    ws "''bell'You must enter a value greater than zero!"
$    goto get_wrd1
$    
$get_wrd2:
$    
$	dsize = 3 * wrd1
$    iq wrd2 "Final word size (* ''dsize' *) ? "
$	if(wrd2.EQS."") then wrd2 := 'dsize'
$    wrd2 = f$integer(wrd2)
$    if (wrd2.EQ.wrd1) then goto just_one ! Analysis should be done only once
$    if (wrd2.GT.wrd1) then goto get_wrd3
$    ws "''bell'You must enter a value greater than or equal to the initial size."
$    goto get_wrd2
$    
$get_wrd3:
$    
$	dsize = (wrd2 - wrd1) / 4
$	if(dsize.EQ.0) then dsize = 1
$    iq wrd3 "Increment word size (* ''dsize' *) ? "
$	if(wrd3.EQS."") then wrd3 := 'dsize'
$    wrd3 = f$integer(wrd3)
$    if (wrd3.GT.0) then goto get_mode
$    ws "''bell'The increment must be a positive integer!"
$    goto get_wrd3
$    
$just_one:
$    
$    wrd3 = 0
$    
$get_mode:
$    
$    ! Find out if it should be done interactively or in batch.
$    
$    ws ""
$    ws "Do this analysis in"
$    ws ""
$    ws "1) Interactive mode (Faster but more costly)"
$    ws "2) Batch mode (Slower but less costly)"
$    ws ""
$    iq choice "Please select your choice now (* 1 *): "
$    if (choice.EQS."") then choice :=1
$    if (choice.EQS."1") then mode := INTERACTIVE
$    if (choice.EQS."2") then mode := BATCH
$    if ("''mode'".NES."") then goto summarize_word
$    ws "''bell'Please enter 1 or 2."
$    goto get_mode
$    
$summarize_word:
$    
$    ws ""
$    ws "The following MULTIDOT analysis (word-type) will be executed:"
$    ws ""
$    if (reverse1.EQS."NO") then -
	ws "Horizontal sequence: ''seq1' from ''begin1' to ''end1'."
$    if (reverse1.EQS."YES") then -
	ws "Horizontal sequence: REVERSE of ''seq1' from ''begin1' to ''end1'."
$    if (reverse2.EQS."NO") then -
	ws "Vertical sequence: ''seq2' from ''begin2' to ''end2'."
$    if (reverse2.EQS."YES") then -
	ws "Vertical sequence: REVERSE of ''seq2' from ''begin2' to ''end2'."
$    if (wrd3.EQ.0) then ws "Word size is ''wrd1'."
$    if (wrd3.GT.0) then ws "Word size: ''wrd1' to ''wrd2' by ''wrd3'."
$    ws "Intermediate point file: ''outname'"
$    ws ""
$    iq choice "Is this okay (* Yes *) ? "
$    choice = f$extract(0, 1, choice)
$    if (choice.EQS."") then goto do_word
$    if (choice.EQS."Y") then goto do_word
$    
$    ! Something is wrong. Give a chance to correct it, or to give up.
$    
$repeat_word:
$    
$    ws ""
$    ws "Do you want to -"
$    ws ""
$    ws "1) Try again"
$    ws "2) Give up"
$    ws ""
$    iq choice "Choice (* 1 *): "
$    if (choice.EQS."") then goto get_seq1
$    if (choice.EQS."1") then goto get_seq1
$    if (choice.EQS."2") then exit
$    ws "''bell'Either 1 or 2."
$    goto repeat_word
$    
$do_word:
$    
$    ! Open the command procedure file that is to be executed.
$    
$    open/write comfile multidot.cmd
$    wc := "write comfile"
$    if (mode.EQS."BATCH") then wc "$ set noon"
$    if (mode.EQS."BATCH") then wc "$ set process/name=''procname'"
$    if (mode.EQS."BATCH") then wc "$ set default ''c_dir'"
$    if (mode.EQS."BATCH") then wc "$ set on"
$    if (mode.EQS."INTERACTIVE") then wc "$ set verify"	! Echo it to the screen
$    
$    wrd1 = wrd1 - wrd3
$    
$loop_word:
$    
$    		wrd1 = wrd1 + wrd3
$    		if (wrd1.GT.wrd2) then wrd1 = wrd2
$    		wc "$ compare/WORD/OUTFILE=''outname'"
$    		wc "''seq1'"
$    		wc "''begin1'"
$    		wc "''end1'"
$    		wc "''reverse1'"
$    		wc "''seq2'"
$    		wc "''begin2'"
$    		wc "''end2'"
$    		wc "''reverse2'"
$    		wc "''wrd1'"
$    		wc "$ dotplot/INFILE=''outname'"
$    		wc ""	! What point density?
$    		wc "P"	! Plot
$    		wc "Q"	! Quit
$    if (wrd1.LT.wrd2) then goto loop_word
$    
$    if (mode.EQS."INTERACTIVE") then wc "$ set noverify"
$    close comfile
$    
$    ! Now execute the command procedure we just wrote.
$    
$    if (mode.EQS."BATCH") then goto batch_word
$    @multidot.cmd
$    delete multidot.cmd;0	! Clean up the droppings
$    exit
$    
$batch_word:
$    
$    submit/noprint/notify/name='procname' multidot.cmd/delete
$    ws ""
$    ws "Your diagonal dot matrix analyses have been submitted to the batch queue"
$    ws "for processing. You will be notified when it is finished."
$    exit
$    
$use_window:
$    
$    ! There are two parameters that have to be varied here: window size,
$    ! and stringency. Unlike the native GCG program, this procedure expects
$    ! the stringency to be expressed in percent. Thus, to see 10 or more 
$    ! matches when the window size is 20, the stringency should be set to
$    ! 50%. TO get a meaningful value for stringency, the global scoring
$    ! matrix is replaced by a local data file that consists entirely of
$    ! 1's on the diagonal and 0 everywhere else. Therefore a perfect match
$    ! scores 1 and a mismatch scores zero. This helps with analysis of
$    ! the output, too.
$    
$get_win1:
$
$	! The default initial window size is 5 for protein and 10 for DNA.
$	! The default final window size is five times the initial window,
$	! and the default increment is 1/4 ofthe difference.
$
$	dsize = 5
$	if(seqinfotype.EQS."DNA") then dsize = 10
$    
$    ws ""
$    iq win1 "Initial window size (* ''dsize' *) ? "
$	if(win1.EQS."") then win1 := 'dsize'
$    win1 = f$integer(win1)
$    if (win1.GT.0) then goto get_win2
$    ws "''bell'You must enter a value greater than zero!"
$    goto get_win1
$    
$get_win2:
$    
$	dsize = 5 * win1
$    iq win2 "Final window size (* ''dsize' *) ? "
$	if(win2.EQS."") then win2 := 'dsize'
$    win2 = f$integer(win2)
$    if (win2.EQ.win1) then goto one_window
$    if (win2.GT win1) then goto get_win3
$    ws "''bell'You must enter a value greater than or equal to ''win1'."
$    goto get_win2
$    
$one_window:
$    
$    win3 = 0
$    goto get_per1
$    
$get_win3:
$    
$	dsize = (win2 - win1) / 4
$    iq win3 "Increment window by how much (* ''dsize' *) ? "
$	if(win3.EQS."") then win3 := 'dsize'
$    win3 = f$integer(win3)
$    if (win3.GT.0) then goto get_per1
$    ws "''bell'The increment must be greater than zero!"
$    goto get_win3
$    
$get_per1:
$    
$	! The default initial percent is 20% for protein and 40% for DNA. The
$	! final default is double this. The increment is 10% for proteins and
$	! 20% for DNA
$
$	dsize = 20
$	if(seqinfotype.EQS."DNA") then dsize = 40
$    ws ""
$    iq per1 "What initial PERCENT matches (* ''dsize' *) ? "
$	if(per1.EQS."") then per1 := 'dsize'
$    per1 = f$integer(per1)
$    if ((per1.GE.10) .AND. (per1.LT.100)) then goto get_per2
$    ws "''bell'The percentage should be between 10 and 100."
$    goto get_per1
$    
$get_per2:
$    
$	dsize = 2 * per1
$	if(dsize.GT.100) then dsize = 100
$    iq per2 "What final PERCENT matches (* ''dsize' *) ? "
$	if(per2.EQS."") then per2 := 'dsize'
$    per2 = f$integer(per2)
$    if (per2.EQ.per1) then goto one_percent
$    if ((per2.GT.per1) .AND. (per2.LE.100)) then goto get_per3
$    ws "''bell'The final percentage must be between ''per1' and 100, inclusive."
$    goto get_per2
$    
$one_percent:
$    
$    per3 = 0
$    goto summarize_win
$    
$get_per3:
$    
$	dsize = 10
$	if(seqinfotype.EQS."DNA") then dsize = 20
$    iq per3 "What PERCENT increment in matches (* ''dsize' *) ? "
$	if(per3.EQS."") then per3 = 'dsize'
$    per3 = f$integer(per3)
$    if (per3.GT.5) then goto summarize_win
$    ws "''bell'The smallest increment is 5."
$    goto get_per3
$    
$summarize_win:
$    
$    ws ""
$    ws "The following MULTIDOT analysis (window/stringency) will be executed:"
$    ws ""
$    if (reverse1.EQS."NO") then -
	ws "Horizontal sequence: ''seq1' from ''begin1' to ''end1'."
$    if (reverse1.EQS."YES") then -
	ws "Horizontal sequence: REVERSE of ''seq1' from ''begin1' to ''end1'."
$    if (reverse2.EQS."NO") then -
	ws "Vertical sequence: ''seq2' from ''begin2' to ''end2'."
$    if (reverse2.EQS."YES") then -
	ws "Vertical sequence: REVERSE of ''seq2' from ''begin2' to ''end2'."
$    if (win3.EQ.0) then ws "Window size: ''win1'."
$    if (win3.GT.0) then ws "Window size: ''win1' to ''win2' by ''win3'."
$    if (per3.EQ.0) then ws "Percent matches: ''per1'."
$    if (per3.GT.0) then ws "Percent matches: ''per1' to ''per2' by ''per3'."
$    ws "Intermediate point file: ''outname'"
$    ws ""
$    iq choice "Is this okay (* Yes *) ? "
$    choice = f$extract(0, 1, choice)
$    if (choice.EQS."") then goto do_window
$    if (choice.EQS."Y") then goto do_window
$    
$    ! Something is wrong. Give a chance to correct it, or to give up.
$    
$repeat_window:
$    
$    ws ""
$    ws "Do you want to -"
$    ws ""
$    ws "1) Try again"
$    ws "2) Give up"
$    ws ""
$    iq choice "Choice (* 1 *): "
$    if (choice.EQS."") then goto get_seq1
$    if (choice.EQS."1") then goto get_seq1
$    if (choice.EQS."2") then exit
$    ws "''bell'Either 1 or 2."
$    goto repeat_window
$    
$do_window:
$    
$    ! Open and write the command procedure file that is to be executed.
$    
$    open/write comfile multidot.cmd
$    wc := write comfile
$    
$    wc "$ set noon"
$    wc "$ set process/name=''procname'"
$    wc "$ set on"
$    wc "$ set default ''c_dir'"
$    
$    win1 = win1 - win3
$    
$win_loop:
$    
$    	win1 = win1 + win3
$    	if (win1.GT.win2) then win1 = win2
$    	per = per1 - per3
$    
$per_loop:
$    
$    		per = per + per3
$    
$    		! Calculate the actual stingency that the program COMPARE expects.
$    
$    		matches = win1 * per / 100
$    
$    		wc "$ compare/OUTFILE=''outname'"
$    		wc "''seq1'"
$    		wc "''begin1'"
$    		wc "''end1'"
$    		wc "''reverse1'"
$    		wc "''seq2'"
$    		wc "''begin2'"
$    		wc "''end2'"
$    		wc "''reverse2'"
$		wc "''win1'"
$		wc "''matches'"
$    		wc "$ dotplot/INFILE=''outname'"
$    		wc ""	! What point density?
$    		wc "P"	! Plot
$    		wc "Q"	! Quit
$	if (per.LT.per2) then goto per_loop
$    if (win1.LT.win2) then goto win_loop
$
$    close comfile
$    
$    submit/noprint/notify/name='procname' multidot.cmd/delete
$    ws ""
$    ws "Your diagonal dot matrix analyses have been submitted to the batch queue"
$    ws "for processing. You will be notified when it is finished."
$    exit
$
$terminate: ! Jump here on ^Y
$    
$    if(f$search("multidot.cmd").EQS."") then exit
$    close comfile
$    delete multidot.cmd;*
$    
$    exit
********************

$! CHK_PROCESS.COM
$
$! October 11, 1988
$
$! Written by Steve Clark
$! Mt. Sinai Hospital Research Institute, Toronto, Canada.
$
$! This command procedure dtermines whether there is a process running on
$! the system with the same process name as that specified by P1. It sets
$! the global integer symbol EXIST_PROCESS to 1 if the process exists,
$! to 0 if the process doesn't exist, and to -1 if no parameter P1 is
$! specified.
$
$	if(p1.nes."") then goto get_first_pid
$	exist_process == -1
$	exit
$
$ get_first_pid:
$
$!	To get the process name, you must get each PID number in turn, determine
$!	the name associated with it, then compare it to P1.
$
$	next = 0
$
$ get_next_pid:
$
$	pid = f$pid(next)
$	if(pid.eqs."") then goto done
$	pname = f$getjpi(pid, "PRCNAM")
$	if("''p1'".nes."''pname'") then goto get_next_pid
$
$	exist_process == 1
$	exit
$
$ done:
$
$	exist_process == 0
$	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