\documentclass{article}
\setlength{\textheight}{9in}
\setlength{\topmargin}{0in}
\begin{document}
\title{CRITICA Documentation}
\date{August 1, 2000}
\author{Jonathan Badger}
\maketitle
\tableofcontents
\newpage
\section{Introduction}

CRITICA is a series of programs designed to find protein coding
regions in intron-less genomic DNA. It doesn't require (or use) any
sequence annotation for its predictions. However, it does require that
at least some of the genes in the sequence have homologs in data-banks
for CRITICA's comparative analysis and that that the overall dicodon
bias in these genes is typical of the organism for CRITICA's
non-comparative analysis. In practice, this means that one should
either analyze complete genomes with CRITICA or a sizeable chunk
($>100,000$ nucleotides). It is not required that the sequence be in one
contig, however.

\section{Installation}

CRITICA is implemented as a series of Perl 5.0 scripts and ANSI C
code. To make the C programs simply type make (if your compiler is not
GCC, you'll have to edit the Makefile). Put the generated programs
somewhere on your path. The scripts in the scripts directory also need
to be on your path (Critica.pm should be placed somewhere on your
PERLLIB path). The location of perl in the scripts may have to be
edited if your perl is not in /usr/bin. Also, you may have to adjust
several environmental variables for the script blast-contigs:
CRITICA\_BLASTN (default value ``blastn'') is the name of the BLASTN
program used by CRITICA. CRITICA\_BLASTPARM (default value ``-warnings
-nogaps E=1e-4 E2=1e-4'' for use with WU-BLAST2; a good set of values
for NCBI BLAST2 is ``-gF -e1e-4'') is a set of parameters to be sent
to BLASTN (if you use a version of BLASTN that generates gapped
alignments, you should turn off this feature by adding a flag such as
-nogap to CRITICA\_BLASTPARM). CRITICA\_BLASTDB (default value ``gb'')
is the name of the BLASTN database to be searched. WinCritica.mak is a
makefile for Microsoft Visual C++ and Makefile is the UNIX
makefile. The two files critica.tex and critica.html (with a
stylesheet in critica.css) are simply the documentation that you are
reading now in \LaTeX~and html format. LICENSE is a copy of the terms
of the GNU Public License (GPL).

\section{Running CRITICA -- a tutorial}

Probably the best way to demonstrate how to use CRITICA is to walk
though a typical analysis. To start with, you need the DNA to be
analyzed to be in FASTA (also called Pearson) format. \begin{samepage}
Here is an example of part of a FASTA file:

\begin{verbatim}
>ECOLI
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTC
TGATAGCAGCTTCTGAACTGGTTACCTGCCGTGAGTAAATTAAAATTTTATTGACTTAGG
TCACTAAATACTTTAACCAATATAGGCATAGCGCACAGACAGATAAAAATTACAGAGTAC
ACAACATCCATGAAACGCATTAGCACCACCATTACCACCACCATCACCATTACCACAGGT
AACGGTGCGGGCTGACGCGTACAGGAAACACAGAAAAAAGCCCGCACCTGACAGTGCGGG
\end{verbatim}
\end{samepage}

Let's say you were going to analyze {\it E. coli} with CRITICA and the
genomic data was contained in a file named contigs. The first
part of a CRITICA analysis requires running blastn on small
subsequences of the genomic data: the script to do this is called
blast-contigs. Below is an example:

\begin{verbatim}
        blast-contigs contigs > EC.blast
\end{verbatim}

If you had additional data not in the main database. you could
additionally run blast-contigs as follows:

\begin{verbatim}
        blast-contigs -db=mydata contigs >> EC.blast
\end{verbatim}

BLASTing a genome against GenBank often takes overnight; be prepared.

Next the BLAST output has to be processed into a format
readable by CRITICA.

\begin{verbatim}
	make-blastpairs EC.blast > EC.blast.pairs
\end{verbatim}

Now the program scanblastpairs can be used to tabulate the
triplets aligned to {\it E. coli} in the blast.pairs
file

\begin{verbatim}
	scanblastpairs contigs EC.blast.pairs EC.triplets
\end{verbatim}

Now we can predict coding regions. Because it is best to have CRITICA
iterate several times, a script iterate-CRITICA has been
supplied. It is run as follows:

\begin{verbatim}
	iterate-CRITICA ecorfs contigs EC.triplets
\end{verbatim}

This will take twenty minutes (on a 300 Mhz Pentium II) to several
hours (on slower machines) to run and will generate numerous files,
among them ecorfs0.cds, ecorfs1.cds, ecorfs2.cds,
and ecorfs3.cds. These files are the predictions made by
CRITICA at each iteration.  The last file (ecorfs3.cds) should
be the most accurate.
\begin{samepage}
Here is part of the ecorfs3.cds file for interpretation:

\begin{verbatim}
ECOLI  35371  34781  2.86e-35   16   3399   2708  -93 GTG    48  12 AGGA    
ECOLI  36162  35377  8.28e-50   32   2171   6818   62 ATG    40   8 GGAG    
ECOLI  37824  36271  8.18e-45   64   1002   7385   63 ATG    84   3 GAGGTG  
ECOLI  39115  37898  7.12e-53  128      0  10183   63 ATG    76   7 AGGAG   
ECOLI  40386  39244  4.09e-50  128    324   9258   63 ATG    73   7 GAGGT  
ECOLI  41931  40417  1.29e-53   64    225  10286   63 ATG  -123   0 -  
  ^      ^      ^      ^        ^      ^    ^      ^^^^^^   ^^^^^^^^^^^^
contig start   end    P-value   M     comp  dicod  initiator    SD
\end{verbatim}
\end{samepage}

The P-value is the amount of statistical support for the coding
region. Like BLAST, low scores are better. The matrix (M) is the best
comparative matrix -- high matrices mean the homologs found were more
distant. The start and end are the start and end of the coding region,
including the terminator. Comp and dicod represent the scores the orf
has from comparative support and dicodon support respectively, The
initiator information consists of the score given to the initiator
codon and the codon itself. The SD information consists of the score
given to the best SD sequence for the chosen initiator, the number of
nucleotides between the end of the SD sequence and the start of the
coding region, and the SD sequence itself.


\section{Summary of programs and scripts}

\subsection{addlongorfs}

\begin{verbatim}
Usage: addlongorfs [options] cds-file seq-file 
Valid addlongorfs options:
        -orf-aa-length=length
        -genetic-code=number
\end{verbatim}

This program lists orfs over a certain length that aren't included in
the cds file given as input. The main use of this program is that it
is called by iterate-CRITICA if the -add-longorfs=
option is given.

\subsection{annotation-compare}

\begin{verbatim}
usage: annotation-compare annotation.cds prediction.cds
\end{verbatim}

This program compares the results of a CRITICA run against a file of
known coding regions of the format contig-name start end (...).
It outputs the percentage of false positives and false negatives, as
well as the percentages of too early and late start positions. It also
generates files ending in .fp, .fn, .early and
.late which have the respective CRITICA calls in them.

\subsection{blast-contigs}

\begin{verbatim}
Usage: blast-contigs [-db=dbname] fasta-file 
\end{verbatim}

This script blasts the DNA in a fasta file against a DNA
database. Either the default database (usually GenBank or a
non-redundant DNA database) is used, or if the -db= option is
chosen, a database of the user's choice. The DNA to be blasted is
split into small fragments before BLASTing.

\subsection{compare-orf-lists}

\begin{verbatim}
Usage: compare-orf-lists list1 list2
\end{verbatim}

This script compares lists of coding regions. The lists can either be
in the format of CRITICA .cds files or in a simple format: 

\begin{verbatim}  
contig start end
\end{verbatim}

The resulting list is a list of coding regions in list1 not in list2
and the vice versa. The start position of the coding region doesn't
matter in these comparisons.

\subsection{critica}

\begin{verbatim}
Usage: critica [options] seq triplet-scores
Valid critica options:
        -scoring-matrix=file
        -dicodon-scores=file
        -init-scores=file
        -sd-scores=file
        -prom-scores=file
        -no-sdscores
        -prom-find
        -genetic-code=number
        -threshold=value
        -alpha=value
        -strict-threshold
        -frameshift-threshold=value
        -quick-stats
\end{verbatim}

This is the main CRITICA program. It takes as arguments the sequence
to be analyzed and set of comparative triplets (created by
scanblastpairs --- see above). It outputs a list of predicted
coding regions with start positions within two orders of magnitude of
the best start position for each coding region.

\subsubsection*{Options}

\noindent -scoring-matrix=file --- this causes CRITICA to use
an empirically derived scoring matrix instead of the theoretically
determined scores as described in the manuscript. An advantage of the
empirical scores is that conservative changes can be
rewarded. However, in practice I've usually found that this does not
significantly change the results, suggesting that identities are
sufficient for finding coding regions.

\bigskip \noindent -dicodon-scores --- this causes CRITICA to
use a set of dicodon scores (generated by the dicodontable
program --- see below) for noncomparative support.

\bigskip \noindent -init-scores --- this causes CRITICA to use
a set of initiator scores (generated by iterate-critica ---
see below) for initiator choice.

\bigskip \noindent -sd-scores --- this causes CRITICA to use a
set of SD scores (generated by iterate-critica --- see
below) to improve initiator choice.

\bigskip \noindent -prom-scores --- this causes CRITICA to use a
set of promoter scores (generated by iterate-critica --- see
below) to improve initiator choice.

\bigskip \noindent -no-sdscores --- this turns off the search
for SD sequences. This is useful for analyzing organisms which you
know don't use SD sequences.

\bigskip \noindent -prom-find --- This turns on the currently
experimental (and disabled by default) promoter finding routines of
CRITICA.

\bigskip \noindent -genetic code=number --- this causes CRITICA
to use a different genetic code than the ``universal'' one. The
numbers used are are the standard NCBI genetic code numbers.

\bigskip \noindent -threshold=value --- this sets the threshold
for the maximum p-value allowed for an region to be considered a
coding region. The default is $1\times 10^{-4}$, which seems to work
the best.

\bigskip \noindent -alpha=value --- this sets the parameter
``alpha'' which scales the dicodon scores down to compensate for their
non-independent nature. The default value is 1.0.

\bigskip \noindent -strict-threshold --- by default CRITICA
accepts as coding regions regions that meet the demand of the
threshold even if extending them to a terminator causes the region to
drop below the amount of support required by the threshold. Using
-strict-threshold eliminates them. In practice this tends not
to be desirable, however.

\bigskip \noindent -frameshift-threshold=value --- causes
CRITICA to check for potential frame\-shift errors in the sequences it
is analyzing. The value supplied is the number of nucleotides between
two regions of coding evidence in difference frames below which the
regions are assumed to be caused by a frameshift. A good value to use
is 10. The potential frameshifts are outputed to stderr -- in the
case of iterate-critica (below) the messages are redirected to
the .mesg files.

\bigskip \noindent -quick-stats --- one of the more time
consuming parts of a CRITICA run is calculating the parameters K and
lambda. This option causes CRITICA to use conservative values (0.2
and 0.015 respectively) and thus speeds the run (at the expense of
losing some marginal cases).

\subsection{dicodontable}

\begin{verbatim}
Usage: dicodontable [options] cds-file seq-file
Valid dicodontable options:
        -fraction-coding=fraction
\end{verbatim}

This program generates a set of dicodon scores from a list of coding
regions (generated by best-inits, above) and the FASTA file of
the DNA to be analyzed. The option -fraction-coding allows
adjustment for the fact that in the first (comparative evidence only
iteration) less than a realistic percentage of the DNA is called
coding. The default behavior of iterate CRITICA is to set this
parameter to 0.8 in the first iteration, and not set it at all in
following iterations. The table generated by dicodontable has the
following format:

\begin{verbatim}
dicodon dicodon-score frequency-in-noncoding DNA
\end{verbatim}

\subsection{empirical-matrix}

\begin{verbatim}
Usage: empirical-matrix coding-blast-pairs
\end{verbatim}

This program creates a scoring matrix including conservative changes
from a blast-pairs file created from a BLASTN run of known coding
regions in the organism of interest. The scores are based on the log
odds of two triplets being aligned in the coding frame (the 1 frame)
vs. being aligned in other frames. The matrix is used via the
-scoring-matrix option of CRITICA and
iterate-CRITICA.

\subsection{extractcoding}

\begin{verbatim}
Usage: extractcoding [options] fasta-file cds-file
valid extractcoding options:
     -desc  (include gene descriptions in headers)
     -prot  (output translated sequences)
     -genetic-code=NCBI-num
     -min=num (minimum length (in amino acids) of extracted genes)
     -max=num (maximum length (in amino acids) of extracted genes)
\end{verbatim}

This program generates FASTA files of the genes predicted by CRITICA
(actually, it can read any file that starts contig-name start end).  
The default is to output DNA sequences; the program can also
translate the genes if desired (through the -prot option.  The
-desc option is useful particularly for generating FASTA files
from annotation files of the format 
contig-name start end description....

\subsection{histo-orf}

\begin{verbatim}
Usage: histo-orf bin-size cds-file
\end{verbatim}

This script produces a histogram of coding-region lengths (in amino
acids). The arguments are the size (in amino acids) of the histogram
bin size and the list of coding regions. The list can either be in the
format of CRITICA .cds files or any file that starts 
contig-name start end.

The program produces input showing the number of orfs that have a
length less than the number on the left (and whose length is greater
or equal to the number on the left on the previous line). The percent
of all coding regions in the bin is listed, as well as the culminative
percentage.

\subsection{intergenic}

\begin{verbatim}
Usage: intergenic [options] fasta-file cds-file
valid intergenic options:
     -desc  (include gene descriptions in headers)
     -min=num (minimum length (in amino acids) of extracted genes)
     -max=num (maximum length (in amino acids) of extracted genes)
\end{verbatim}

This program generates FASTA files of the regions of the sequences
provided that are not covered by the genes listed in cds-file.

\subsection{iterate-critica}

\begin{verbatim}
iterate-critica: [options] output-name seq-file triplets-scores
Valid iterate-critica options:
        -iterations=number (default: 3)
        -scoring-matrix=file
        -initial-dicod=file
        -initial-init=file
        -initial-sd=file
        -initial-prom=file
        -no-sdscores
        -prom-find
        -threshold=value
        -alpha=value
        -fraction-coding=value
        -add-longorfs=length
        -genetic-code=number
        -strict-threshold
        -frameshift-threshold=value
        -quick-stats
\end{verbatim}

This script is the normal way of automating a CRITICA run. The options
are mostly identical to those of the CRITICA program itself except for
the ability to set the number of iterations.

Some options of note --- -fraction-coding also changing the
default value is 0.8, (80\%), which is typical of the fraction of DNA
which is actually part of a gene in prokaryotes. If you know that your
organism is quite different from this you can set your own value.

-add-longorfs=length adds to the first iteration all ORFS that
would encode proteins over a certain length (given in amino
acids). This can be helpful if there is very little comparative
information available for the organism being analyzed.

\subsection{lookat}

\begin{verbatim}
Usage: lookat [options] fasta-file [contig-name] start end
valid lookat options:
     -num include num bases upstream of given region
     num include num bases downstream of given region
\end{verbatim}

This program allows the user to examine any sub-sequence of a DNA or
protein sequence. If the start is less than end, the reverse complement
of the subsequence is given (for DNA of course!)

\subsection{make-blast-pairs}

\begin{verbatim}
Usage: make-blastpairs [-exclude=string] blast-file
\end{verbatim}

This script converts BLAST output (generated by blast -contigs,
above), into a FASTA file with alternating records representing query
DNA and its corresponding BLAST match. The option -exclude=
allows the exclusion of DNA from a particular genus or species,
generally the same genus or species of the query organism. In general,
this is unneeded because make-blastpairs automatically excludes
matches with more than 97\% identity.

\subsection{map-critica-orfs}

\begin{verbatim}
Usage: map-critica-orfs cds-list
\end{verbatim}

This script lists the amount of overlap of each coding region
predicted by CRITICA with the previous coding region. The output is
the same as the input but with two additional columns -- the first is
the number of nucleotides of overlap with the previous coding region
(negative numbers mean there is no overlap but a gap between the two
genes). The second is the amount of overlap represented as the percent
of the length of the coding region. 

\subsection{motiffind}

\begin{verbatim}
Usage: motiffind [-options] pattern fasta-file
Valid motiffind options:
        -compstrand
\end{verbatim}

This program finds instances of the pattern supplied in the FASTA
file. In the FASTA file contains DNA sequences, ambiguous codes such as
R and Y can be used (the program decides if a sequence is protein if
its GC is less than 0.20). If the compstrand option is
selected, the reverse complement is searched (for DNA sequences).

\subsection{removeoverlaps}

\begin{verbatim}
Usage: removeoverlaps seq-file cds-file percent-allowed
\end{verbatim}
 
This program removes coding regions which overlap more than the
percentage given. The coding region removed is the one with the least
support.

\subsection{reportscore}

\begin{verbatim}
Usage: reportscore [options] seq triplet-scores contig start end matrix
Valid options:
        -dicodon-scores=file
        -init-scores=file
        -sd-scores=file
        -genetic-code=number
\end{verbatim}

The program allows one to determine the score and p-value of a region of
DNA. This is useful for understanding why CRITICA did not call an expected
region a coding region.

\subsection{scanblastpairs}

\begin{verbatim}
Usage: scanblastpairs seq blast-pairs triplets
\end{verbatim}

This program takes as inputs the FASTA file of the sequence to be
analyzed, the blast-pairs generated by make-blastpairs (see above),
and generates a table of triplets aligned to each triplet in the
sequence file.

\subsection{translate}

\begin{verbatim}
Usage: translate [options] fasta-file [contig-name] start end
valid translate options:
     -num include num residues upstream of given region
     num include num residues downstream of given region
     -genetic-code=NCBI-num
\end{verbatim}

This program allows the user to translate any sub-sequence of a DNA
sequence. If the start is less than end, the reverse complement of the
subsequence is used. 

\section{What's New}

\subsection*{Version 1.05}

\begin{itemize}

\item Changes to Array.c and Array.h to make them more portable.
Thanks to Peter Li and Xiaoying Lin.

\item BLAST2 puts carriage returns in sequence names to keep lines
less than 80 characters in length. This caused make-blastpairs to get
confused. It has been fixed.

\item Pseudocount method in \verb+dicodontable.c+ could cause
noncoding count to go negative. Gary Olsen discovered this and
suggested a fix.
\end{itemize}

\subsection*{Version 1.04}

\begin{itemize}

\item CRITICA now runs under Windows (when compiled with Visual
C++). Thanks to Brian Schick.

\item blast-contigs now also works with NCBI BLAST2, which has a
slightly different output format than traditional BLAST.

\item CRITICA now handles periods, spaces, tabs, and numbers in contig
names more gracefully. There {\em shouldn't} be any need for renaming 
contigs any more.

\end{itemize}

\section{License}

CRITICA (and its source code) are copyrighted by Jonathan Badger but
are licensed under the GNU Public License (GPL). What this means is
that not only are you free to use CRITICA and have access to its
source code, but are free to incorporate parts of CRITICA into your
own programs. However, if such programs are used beyond your
institution or company, you must release their source code under the
GPL as well.

Frankly, I'm surprised I have to say this explicitly. The GPL simply
embodies what I've always understood as the scientific spirit of
cooperation which distinguishes modern scientists from secret-hoarding
medieval alchemists. Just as bench biologists make the mutants and DNA
sequences mentioned in their papers freely available to their
colleagues so that additional research can be done, computational
biologists should be expected to make their source code of the
programs mentioned in their papers freely available \footnote{I'm not
against proprietary software (sometimes proprietary software has a
nicer interface and fewer quirks than the equivalent research
software). However, writing a proprietary sequence analysis package is
no more science than the writing of a proprietary word processor and
thus is not worthy of a scientific publication.}

Yet there is a disturbing trend in computational biology to only make
programs available through web servers and make access to the source
code impossible or only possible after signing legal documents which
generally forbid use of the licensed source code in the licensee's
projects. This behavior can only impede scientific progress. So don't
do it.


\section{Contact Information}
If you wish to contact me about CRITICA, here is my address:

\begin{verbatim}
Jonathan Badger
Department of Computer Science
University of Waterloo
Waterloo, Ontario N2L 3G1
Canada
email: jhbadger@monod.uwaterloo.ca 
phone: (519) 888-4567 ext. 5321
fax: (519) 885-1208
\end{verbatim}

\end{document}
