[biopieces.git] / bp_doc / biotools_cookbook.tex

%% LyX 1.5.1 created this file.  For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass[english]{scrartcl}
\usepackage[T1]{fontenc}
\usepackage[latin9]{inputenc}
\setlength{\parskip}{\medskipamount}
\setlength{\parindent}{0pt}
\usepackage{amsmath}
\usepackage{graphicx}
\IfFileExists{url.sty}{\usepackage{url}}
                      {\newcommand{\url}{\texttt}}

\makeatletter

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LyX specific LaTeX commands.
%% Because html converters don't know tabularnewline
\providecommand{\tabularnewline}{\\}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Textclass specific LaTeX commands.
\newenvironment{lyxcode}
{\begin{list}{}{
\setlength{\rightmargin}{\leftmargin}
\setlength{\listparindent}{0pt}% needed for AMS classes
\raggedright
\setlength{\itemsep}{0pt}
\setlength{\parsep}{0pt}
\normalfont\ttfamily}%
 \item[]}
{\end{list}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
\usepackage[colorlinks=true, urlcolor=blue, linkcolor=black]{hyperref}

\usepackage{babel}
\makeatother

\begin{document}

\title{Biotools Cookbook}


\author{Martin Asser Hansen}


\publishers{John Mattick Group\\
Institute for Molecular Bioscience\\
University of Queensland\\
Australia\\
E-mail: mail@maasha.dk}

\maketitle
\thispagestyle{empty}

\newpage{}

\tableofcontents{}

\listoffigures


\newpage{}


\section{Introduction}

Biotools is a selection of simple tools that can be linked together
(piped as we shall call it) in a very flexible manner to perform both
simple and complex tasks. The fundamental idea is that biotools work
on a data stream that will only terminate at the end of an analysis
and that this data stream can be passed through several different
biotools, each performing one specific task. The advantage of this
approach is that a user can perform simple and complex tasks without
having to write advanced code. Moreover, since the data format used
to pass data between biotools is text based, biotools can be written
by different developers in their favorite programming language ---
and still the biotools will be able to work together.

In the most simple form bioools can be piped together on the command
line like this (using the pipe character '|'):

\begin{lyxcode}
read\_data~|~calculate\_something~|~write\_result
\end{lyxcode}
However, a more comprehensive analysis could be composed:

\begin{lyxcode}
read\_data~|~select\_entries~|~convert\_entries~|~search\_database~~

evaluate\_results~|~plot\_diagram~|~plot\_another\_diagram~|

load\_to\_database
\end{lyxcode}
The data stream that is piped through the biotools consists of records
of key/value pairs in the same way a hash does in order to keep as
simple a structure as possible. An example record can be seen below:

\begin{lyxcode}


REC\_TYPE:~PATSCAN

MATCH:~AGATCAAGTG

S\_BEG:~7

S\_END:~16

ALIGN\_LEN:~9

S\_ID:~piR-t6

STRAND:~+

PATTERN:~AGATCAAGTG

-{}-{}-
\end{lyxcode}
The '-\/-\/-' denotes the delimiter of the records, and each key
is a word followed by a ':' and a white-space and then the value.
By convention the biotools only uses upper case keys (a list of used
keys can be seen in Appendix~\ref{sec:Keys}). Since the records
basically are hash structures this mean that the order of the keys
in the stream is unordered, and in the above example it is pure coincidence
that HIT\_BEG is displayed before HIT\_END, however, when the order
of the keys is importent, the biotools will automagically see to that.

All of the biotools are able to read and write a data stream to and
from file as long as the records are in the biotools format. This
means that if you are undertaking a lengthy analysis where one of
the steps is time consuming, you may save the stream after this step,
and subsequently start one or more analysis from that last step%
\footnote{It is a goal that the biotools at some point will be able to dump
the data stream to file in case one of the tools fail critically.%
}. If you are running a lengthy analysis it is highly recommended that
you create a small test sample of the data and run that through the
pipeline --- and once you are satisfied with the result proceed with
the full data set (see~\ref{sub:How-to-select-a-few-records}).


\section{The Data Stream}


\subsection{How to read the data stream from file?\label{sub:How-to-read-stream}}

You want to read a data stream that you previously have saved to file
in biotools format. This can be done implicetly or explicitly. The
implicit way uses the 'stdout' stream of the Unix terminal:

\begin{lyxcode}
cat~|~<biotool>
\end{lyxcode}
cat is the Unix command that reads a file and output the result to
'stdout' --- which in this case is piped to any biotool represented
by the <biotool>. It is also possible to read the data stream using
'<' to direct the 'stdout' stream into the biotool like this:

\begin{lyxcode}
<biotool>~<~<file>
\end{lyxcode}
However, that will not work if you pipe more biotools together. Then
it is much safer to read the stream from a file explicitly like this:

\begin{lyxcode}
<biotool>~-{}-stream\_in=<file>
\end{lyxcode}
Here the filename <file> is explicetly given to the biotool <biotool>
with the switch -\/-stream\_in. This switch works with all biotools.
It is also possible to read in data from multiple sources by repeating
the explicit read step:

\begin{lyxcode}
<biotool>~-{}-stream\_in=<file1>~|~<biotool>~-{}-stream\_in=<file2>
\end{lyxcode}

\subsection{How to write the data stream to file?\label{sub:How-to-write-stream}}

In order to save the output stream from a biotool to file, so you
can read in the stream again at a later time, you can do one of two
things:

\begin{lyxcode}
<biotool>~>~<file>
\end{lyxcode}
All, the biotools write the data stream to 'stdout' by default which
can be written to a file by redirecting 'stdout' to file using '>'
, however, if one of the biotools for writing other formats is used
then the both the biotools records as well as the result output will
go to 'stdout' in a mixture causing havock! To avoid this you must
use the switch -\/-stream\_out that explictly tells the biotool to
write the output stream to file:

\begin{lyxcode}
<biotool>~-{}-stream\_out=<file>
\end{lyxcode}
The -\/-stream\_out switch works with all biotools.


\subsection{How to terminate the data stream?}

The data stream is never stops unless the user want to save the stream
or by supplying the -\/-no\_stream switch that will terminate the
stream:

\begin{lyxcode}
<biotool>~-{}-no\_stream
\end{lyxcode}
The -\/-no\_stream switch only works with those biotools where it
makes sense that the user might want to terminale the data stream,
\emph{i.e}. after an analysis step where the user wants to output
the result, but not the data stream.


\subsection{How to write my results to file?\label{sub:How-to-write-result}}

Saving the result of an analysis to file can be done implicitly or
explicitly. The implicit way:

\begin{lyxcode}
<biotool>~-{}-no\_stream~>~<file>
\end{lyxcode}
If you use '>' to redirect 'stdout' to file then it is important to
use the -\/-no\_stream switch to avoid writing a mix of biotools
records and result to the same file causing havock. The safe way is
to use the -\/-result\_out switch which explicetly tells the biotool
to write the result to a given file:

\begin{lyxcode}
<biotool>~-{}-result\_out=<file>
\end{lyxcode}
Using the above method will not terminate the stream, so it is possible
to pipe that into another biotool generating different results:

\begin{lyxcode}
<biotool1>~-{}-result\_out=<file1>~|~<biotool2>~-{}-result\_out=<file2>
\end{lyxcode}
And still the data stream will continue unless terminated with -\/-no\_stream:

\begin{lyxcode}
<biotool>~-{}-result\_out=<file>~-{}-no\_stream
\end{lyxcode}
Or written to file using implicitly or explicity \eqref{sub:How-to-write-result}.
The explicit way:

\begin{lyxcode}
<biotool>~-{}-result\_out=<file1>~-{}-stream\_out=<file2>
\end{lyxcode}

\subsection{How to read data from multiple sources?}

To read multiple data sources, with the same type or different type
of data do:

\begin{lyxcode}
<biotool1>~-{}-data\_in=<file1>~|~<biotool2>~-{}-data\_in=<file2>
\end{lyxcode}
where type is the data type a specific biotool reads.


\section{Reading input}


\subsection{How to read biotools input?}

See \eqref{sub:How-to-read-stream}.


\subsection{How to read in data?}

Data in different formats can be read with the appropriate biotool
for that format. The biotools are typicalled named 'read\_<data type>'
such as \textbf{read\_fasta}, \textbf{read\_bed}, \textbf{read\_tab},
etc., and all behave in a similar manner. Data can be read by supplying
the -\/-data\_in switch and a file name to the file containing the
data:

\begin{lyxcode}
<biotool>~-{}-data\_in=<file>
\end{lyxcode}
It is also possible to read in a saved biotools stream (see \ref{sub:How-to-read-stream})
as well as reading data in one go:

\begin{lyxcode}
<biotool>~-{}-stream\_in=<file1>~-{}-data\_in=<file2>
\end{lyxcode}
If you want to read data from several files you can do this:

\begin{lyxcode}
<biotool>~-{}-data\_in=<file1>~|~<biotool>~-{}-data\_in=<file2>
\end{lyxcode}
If you have several data files you can read in all explicitly with
a comma separated list:

\begin{lyxcode}
<biotool>~-{}-data\_in=file1,file2,file3
\end{lyxcode}
And it is also possible to use file globbing:

\begin{lyxcode}
<biotool>~-{}-data\_in={*}.fna
\end{lyxcode}
Or in a combination:

\begin{lyxcode}
<biotool>~-{}-data\_in=file1,/dir/{*}.fna
\end{lyxcode}
Finally, it is possible to read in data in different formats using
the appropriate biotool for each format:

\begin{lyxcode}
<biotool1>~-{}-data\_in=<file1>~|~<biotool2>~-{}-data\_in=<file2>~...
\end{lyxcode}

\subsection{How to read FASTA input?}

Sequences in FASTA format can be read explicitly using \textbf{read\_fasta}:

\begin{lyxcode}
read\_fasta~-{}-data\_in=<file>
\end{lyxcode}

\subsection{How to read alignment input?}

If your alignment if FASTA formatted then you can \textbf{read\_align}.
It is also possible to use \textbf{read\_fasta} since the data is
FASTA formatted, however, with \textbf{read\_fasta} the key ALIGN
will be omitted. The ALIGN key is used to determine which sequences
belong to what alignment which is required for \textbf{write\_align}.

\begin{lyxcode}
read\_align~-{}-data\_in=<file>
\end{lyxcode}

\subsection{How to read tabular input?\label{sub:How-to-read-table}}

Tabular input can be read with \textbf{read\_tab} which will read
in all rows and chosen columns (separated by a given delimter) from
a table in text format.

The table below:

\noindent \begin{center}
\begin{tabular}{lll}
Human & ATACGTCAG & 23524\tabularnewline
Dog & AGCATGAC & 2442\tabularnewline
Mouse & GACTG & 234\tabularnewline
Cat & AAATGCA & 2342\tabularnewline
\end{tabular}
\par\end{center}

Can be read using the command:

\begin{lyxcode}
read\_tab~-{}-data\_in=<file>
\end{lyxcode}
Which will result in four records, one for each row, where the keys
V0, V1, V2 are the default keys for the organism, sequence, and count,
respectively. It is possible to select a subset of colums to read
by using the -\/-cols switch which takes a comma separated list of
columns numbers (first column is designated 0) as argument. So to
read in only the sequence and the count so that the count comes before
the sequence do:

\begin{lyxcode}
read\_tab~-{}-data\_in=<file>~-{}-cols=2,1
\end{lyxcode}
It is also possible to name the columns with the -\/-keys switch:

\begin{lyxcode}
read\_tab~-{}-data\_in=<file>~-{}-cols=2,1~-{}-keys=COUNT,SEQ
\end{lyxcode}

\subsection{How to read BED input?}

The BED (Browser Extensible Data%
\footnote{\url{http://genome.ucsc.edu/FAQ/FAQformat}%
}) format is a tabular format for data pertaining to one of the Eukaryotic
genomes in the UCSC genome brower%
\footnote{\url{http://genome.ucsc.edu/}%
}. The BED format consists of up to 12 columns, where the first three
are mandatory CHR, CHR\_BEG, and CHR\_END. The mandatory columns and
any of the optional columns can all be read in easily with the \textbf{read\_bed}
biotool.

\begin{lyxcode}
read\_bed~-{}-data\_in=<file>
\end{lyxcode}
It is also possible to read the BED file with \textbf{read\_tab} (see~\ref{sub:How-to-read-table}),
however, that will be more cumbersome because you need to specify
the keys:

\begin{lyxcode}
read\_tab~-{}-data\_in=<file>~-{}-keys=CHR,CHR\_BEG,CHR\_END~...
\end{lyxcode}

\subsection{How to read PSL input?}

The PSL format is the output from BLAT and contains 21 mandatory fields
that can be read with \textbf{read\_psl}:

\begin{lyxcode}
read\_psl~-{}-data\_in=<file>
\end{lyxcode}

\section{Writing output}

All result output can be written explicitly to file using the -\/-result\_out
switch which all result generating biotools have. It is also possible
to write the result to file implicetly by directing 'stdout' to file
using '>', however, that requires the -\/-no\_stream swich to prevent
a mixture of data stream and results in the file. The explicit (and
safe) way:

\begin{lyxcode}
...~|~<biotool>~-{}-result\_out=<file>
\end{lyxcode}
The implicit way:

\begin{lyxcode}
...~|~<biotool>~-{}-no\_stream~>~<file>
\end{lyxcode}

\subsection{How to write biotools output?}

See \eqref{sub:How-to-write-stream}.


\subsection{How to write FASTA output?\label{sub:How-to-write-fasta}}

FASTA output can be written with \textbf{write\_fasta}.

\begin{lyxcode}
...~|~write\_fasta~-{}-result\_out=<file>
\end{lyxcode}
It is also possible to wrap the sequences to a given width using the
-\/-wrap switch allthough wrapping of sequence is generally an evil
thing:

\begin{lyxcode}
...~|~write\_fasta~-{}-no\_stream~-{}-wrap=80
\end{lyxcode}

\subsection{How to write alignment output?\label{sub:How-to-write-alignment}}

Pretty alignments with ruler%
\footnote{'.' for every 10 residues, ':' for every 50, and '|' for every 100%
} and consensus sequence can be created with \textbf{write\_align},
what also have the optional -\/-wrap switch to break the alignment
into blocks of a given width: 

\begin{lyxcode}
...~|~write\_align~-{}-result\_out=<file>~-{}-wrap=80
\end{lyxcode}
If the number of sequnces in the alignment is 2 then a pairwise alignment
will be output otherwise a multiple alignment. And if the sequence
type, determined automagically, is protein, then residues and symbols
(+,~:,~.) will be used to show consensus according to the Blosum62
matrix.


\subsection{How to write tabular output?\label{sub:How-to-write-tab}}

Outputting the data stream as a table can be done with \textbf{write\_tab},
which will write generate one row per record with the values as columns.
If you supply the optional -\/-comment switch, when the first row
in the table will be a 'comment' line prefixed with a '\#':

\begin{lyxcode}
...~|~write\_tab~-{}-result\_out=<file>~-{}-comment
\end{lyxcode}
You can also change the delimiter from the default (tab) to \emph{e.g.}
',':

\begin{lyxcode}
...~|~write\_tab~-{}-result\_out=<file>~-{}-delimit=','
\end{lyxcode}
If you want the values output in a specific order you have to supply
a comma separated list using the -\/-keys switch that will print
only those keys in that order:

\begin{lyxcode}
...~|~write\_tab~-{}-result\_out=<file>~-{}-keys=SEQ\_NAME,COUNT
\end{lyxcode}
Alternatively, if you have some keys that you don't want in the tabular
output, use the -\/-no\_keys switch. So to print all keys except
SEQ and SEQ\_TYPE do:

\begin{lyxcode}
...~|~write\_tab~-{}-result\_out=<file>~-{}-no\_keys=SEQ,SEQ\_TYPE
\end{lyxcode}
Finally, if you have a stream containing a mix of different records
types, \emph{e.g.} records with sequences and records with matches,
then you can use \textbf{write\_tab} to output all the records in
tabluar format, however, the -\/-comment, -\/-keys, and -\/-no\_keys
switches will only respond to records of the first type encountered.
The reason is that outputting mixed records is probably not what you
want anyway, and you should remove all the unwanted records from the
stream before outputting the table: \textbf{grab} is your friend (see~\ref{sub:How-to-grab}).


\subsection{How to write a BED output?\label{sub:How-to-write-BED}}

Data in BED format can be output if the records contain the mandatory
keys CHR, CHR\_BEG, and CHR\_END using \textbf{write\_bed}. If the
optional keys are also present, they will be output as well:

\begin{lyxcode}
write\_bed~-{}-result\_out=<file>
\end{lyxcode}

\subsection{How to write PSL output?\label{sub:How-to-write-PSL}}

Data in PSL format can be output using \textbf{write\_psl:}

\begin{lyxcode}
write\_psl~-{}-result\_out=<file>
\end{lyxcode}

\section{Manipulating Records}


\subsection{How to select a few records?\label{sub:How-to-select-a-few-records}}

To quickly get an overview of your data you can limit the data stream
to show a few records. This also very useful to test the pipeline
with a few records if you are setting up a complex analysis using
several biotools. That way you can inspect that all goes well before
analyzing and waiting for the full data set. All of the read\_<type>
biotools have the -\/-num switch which will take a number as argument
and only that number of records will be read. So to read in the first
10 FASTA entries from a file:

\begin{lyxcode}
read\_fasta~-{}-data\_in=test.fna~-{}-num=10
\end{lyxcode}
Another way of doing this is to use \textbf{head\_records} will limit
the stream to show the first 10 records (default):

\begin{lyxcode}
...~|~head\_records
\end{lyxcode}
Using \textbf{head\_records} directly after one of the read\_<type>
biotools will be a lot slower than using the -\/-num switch with
the read\_<type> biotools, however, \textbf{head\_records} can also
be used to limit the output from all the other biotools. It is also
possible to give \textbf{head\_records} a number of records to show
using the -\/-num switch. So to display the first 100 records do:

\begin{lyxcode}
...~|~head\_records~-{}-num=100
\end{lyxcode}

\subsection{How to count all records in the data stream?}

To count all the records in the data stream use \textbf{count\_records},
which adds one record (which is not included in the count) to the
data stream. So to count the number of sequences in a FASTA file you
can do this:

\begin{lyxcode}
cat~test.fna~|~read\_fasta~|~count\_records~-{}-no\_stream
\end{lyxcode}
Which will write the last record containing the count to 'stdout':

\begin{lyxcode}
-{}-{}-

count\_records:~630
\end{lyxcode}
It is also possible to write the count to file using the -\/-result\_out
switch.


\subsection{How to grab specific records?\label{sub:How-to-grab}}

The biotool \textbf{grab} is related to the Unix grep and locates
records based on matching keys and/or values using either a pattern,
a Perl regex, or a numerical evaluation. To easily \textbf{grab} all
records in the stream that has any mentioning of the pattern 'human'
just pipe the data stream through \textbf{grab} like this:

\begin{lyxcode}
...~|~grab~-{}-pattern=human
\end{lyxcode}
This will search for the pattern 'human' in all keys and all values.
The -\/-pattern switch takes a comma separated list of patterns,
so in order to match multiple patterns do:

\begin{lyxcode}
...~|~grab~-{}-pattern=human,mouse
\end{lyxcode}
It is also possible to use the -\/-pattern\_in switch instead of
-\/-pattern. -\/-pattern\_in is used to read a file with one pattern
per line:

\begin{lyxcode}
...~|~grab~-{}-pattern\_in=patterns.txt
\end{lyxcode}
If you want the opposite result --- to find all records that does
not match the patterns, add the -\/-invert switch, which not only
works with the -\/-pattern switch, but also with -\/-regex and -\/-eval:

\begin{lyxcode}
...~|~grab~-{}-pattern=human~-{}-invert
\end{lyxcode}
If you want to search the record keys only, \emph{e.g.} to find all
records containing the key SEQ you can add the -\/-keys\_only switch.
This will prevent matching of SEQ in any record value, and in fact
SEQ is a not uncommon peptide sequence you could get an unwanted record.
Also, this will give an increase in speed since only the keys are
searched:

\begin{lyxcode}
...~|~grab~-{}-pattern=SEQ~-{}-keys\_only
\end{lyxcode}
However, if you are interested in finding the peptide sequence SEQ
and not the SEQ key, just add the -\/-vals\_only switch instead:

\begin{lyxcode}
...~|~grab~-{}-pattern=SEQ~-{}-vals\_only
\end{lyxcode}
Also, if you want to grab for certain key/value pairs you can supply
a comma separated list of keys whos values will then be searched using
the -\/-keys switch. This is handy if your records contain large
genomic sequences and you dont want to search the entire sequence
for \emph{e.g.} the organism name --- it is much faster to tell \textbf{grab}
which keys to search the value for:

\begin{lyxcode}
...~|~grab~-{}-pattern=human~-{}-keys=SEQ\_NAME


\end{lyxcode}
It is also possible to invoke flexible matching using regex (regular
expressions) instead of simple pattern matching. In \textbf{grab}
the regex engine is Perl based and allows use of different type of
wild cards, alternatives, \emph{etc}%
\footnote{\url{http://perldoc.perl.org/perlreref.html}%
}. If you want to \textbf{grab} records withs the sequence ATCG or
GCTA you can do this:

\begin{lyxcode}
...~|~grab~-{}-regex='ATCG|GCTA'
\end{lyxcode}
Or if you want to find sequences beginning with ATCG:

\begin{lyxcode}
...~|~grab~-{}-regex='\textasciicircum{}ATCG'
\end{lyxcode}
You can also use \textbf{grab} to locate records that fulfill a numerical
property using the -\/-eval switch witch takes an expression in three
parts. The first part is the key that holds the number we want to
evaluate, the second part holds one the six operators:

\begin{enumerate}
\item Greater than: >
\item Greater than or equal to: >=
\item Less than: <
\item Less than or equal to: <=
\item Equal to: =
\item Not equal to: !=
\end{enumerate}
And finally comes the number used in the evaluation. So to \textbf{grab}
all records with a sequence length greater than 30:

\begin{lyxcode}
...~length\_seq~|~grab~-{}-eval='SEQ\_LEN~>~30'
\end{lyxcode}
If you want to locate all records containing the pattern 'human' and
where the sequence length is greater that 30, you do this by running
the stream through \textbf{grab} twice:

\begin{lyxcode}
...~|~grab~-{}-pattern='human'~|~length\_seq~|~grab~-{}-eval='SEQ\_LEN~>~30'
\end{lyxcode}
To get the best speed performance, use the most restrictive \textbf{grab}
first.


\subsection{How to remove keys from records?}

To remove one or more specific keys from all records in the data stream
use \textbf{remove\_keys} like this:

\begin{lyxcode}
...~|~remove\_keys~-{}-keys=SEQ,SEQ\_NAME
\end{lyxcode}
In the above example SEQ and SEQ\_NAME will be removed from all records
if they exists in these. If all keys are removed from a record, then
the record will be removed.


\subsection{How to rename keys in records?}

Sometimes you want to rename a record key, \emph{e.g.} if you have
read in a two column table with sequence name and sequence in each
column (see \ref{sub:How-to-read-table}) without specifying the key
names, then the sequence name will be called V0 and the sequence V1
as default in the \textbf{read\_tab} biotool. To rename the V0 and
V1 keys we need to run the stream through \textbf{rename\_keys} twice
(one for each key to rename):

\begin{lyxcode}
...~|~rename\_keys~-{}-keys=V0,SEQ\_NAME~|~rename\_keys~-{}-keys=V1,SEQ
\end{lyxcode}
The first instance of \textbf{rename\_keys} replaces all the V0 keys
with SEQ\_NAME, and the second instance of \textbf{rename\_keys} replaces
all the V1 keys with SEQ. \emph{Et viola} the data can now be used
in the biotools that requires these keys.


\section{Manipulating Sequences}


\subsection{How to get sequence lengths?}

The length for sequences in records can be determined with \textbf{length\_seq},
which adds the key SEQ\_LEN to each record with the sequence length
as the value. It also generates an extra record that is emitted last
with the key TOTAL\_SEQ\_LEN showing the total length of all the sequences.

\begin{lyxcode}
read\_fasta~-{}-data\_in=<file>~|~length\_seq
\end{lyxcode}
It is also possible to determine the sequence length using the generic
tool \textbf{length\_vals} (see \#\#\#), which determines the length
of the values for a given list of keys:

\begin{lyxcode}
read\_fasta~-{}-data\_in=<file>~|~length\_vals~-{}-keys=SEQ
\end{lyxcode}
To obtain the total length of all sequences use \textbf{sum\_vals}
like this:

\begin{lyxcode}
read\_fasta~-{}-data\_in=<file>~|~length\_vals~-{}-keys=SEQ

|~sum\_vals~-{}-keys=SEQ\_LEN
\end{lyxcode}
The biotool \textbf{analyze\_seq} will also determine the length of
each sequence (see~\ref{sub:How-to-analyze}).


\subsection{How to analyze sequence composition?\label{sub:How-to-analyze}}

If you want to find out the sequence type, composition, length, as
well as GC content, indel content and proportions of soft and hard
masked sequence, then use \textbf{analyze\_seq}. This handy biotool
will determine all these things per sequence from which it is easy
to get an overview using the \textbf{write\_tab} biotool to output
a table (see~\ref{sub:How-to-write-tab}). So in order to determine
the sequence composition of a FASTA file with just one entry containing
the sequence 'ATCG' we just read the data with \textbf{read\_fasta}
and run the output through \textbf{analyze\_seq} which will add the
analysis to the record like this:

\begin{lyxcode}
read\_fasta~-{}-data\_in=test.fna~|~analyze\_seq~...



-{}-{}-

GC\%:~50.00

HARD\_MASK\%:~0.00

RES:-:~0

RES:.:~0

RES:A:~1

RES:B:~0

RES:C:~1

RES:D:~0

RES:G:~1

RES:H:~0

RES:K:~0

RES:M:~0

RES:N:~0

RES:R:~0

RES:S:~0

RES:T:~1

RES:U:~0

RES:V:~0

RES:W:~0

RES:Y:~0

RES:\textasciitilde{}:~0

SEQ:~ATCG

SEQ\_LEN:~4

SEQ\_NAME:~test

SEQ\_TYPE:~DNA

SOFT\_MASK\%:~0.00
\end{lyxcode}
Now to make a table of how may As, Ts, Cs, and Gs you can add the
following:

\begin{lyxcode}
...~|~analyze\_seq~|~write\_tab~-{}-keys=RES:A,RES:T,RES:C,RES:G
\end{lyxcode}
Or if you want to see the proportions of hard and soft masked sequence:

\begin{lyxcode}
...~|~analyse\_seq~|~write\_tab~-{}-keys=HARD\_MASK\%,SOFT\_MASK\%
\end{lyxcode}
If you have a stack of sequences in one file and you want to determine
the mean GC content you can do it using the \textbf{mean\_vals} biotool:

\begin{lyxcode}
read\_fasta~-{}-data\_in=test.fna~|~analyze\_seq~|~mean\_vals~-{}-keys=GC\%
\end{lyxcode}
Or if you want the total count of Ns you can use \textbf{sum\_vals}
like this:

\begin{lyxcode}
read\_fasta~-{}-data\_in=test.fna~|~analyze\_seq~|~sum\_vals~-{}-keys=RES:N
\end{lyxcode}

\subsection{How to extract subsequences?\label{sub:How-to-extract}}

In order to extract a subsequence from a longer sequence use the biotool
extract\_seq, which will replace the sequence in the record with the
subsequence (this behaviour should probably be modified to be dependant
of a -\/-replace or a -\/-no\_replace switch). So to extract the
first 20 residues from all sequences do (first residue is designated
1): 

\begin{lyxcode}
...~|~extract\_seq~-{}-beg=1~-{}-len=20
\end{lyxcode}
You can also specify a begin and end coordinate set:

\begin{lyxcode}
...~|~extract\_seq~-{}-beg=20~-{}-end=40
\end{lyxcode}
If you want the subsequences from position 20 to the sequence end
do:

\begin{lyxcode}
...~|~extract\_seq~-{}-beg=20
\end{lyxcode}
If you want to extract subsequences a given distance from the sequence
end you can do this by reversing the sequence with the biotool \textbf{reverse\_seq}
\eqref{sub:How-to-reverse-seq}, followed by \textbf{extract\_seq}
to get the subsequence, and then \textbf{reverse\_seq} again to get
the subsequence back in the original orientation:

\begin{lyxcode}
read\_fasta~-{}-data\_in=test.fna~|~reverse\_seq

|~extract\_seq~-{}-beg=10~-{}-len=10~|~reverse\_seq
\end{lyxcode}

\subsection{How to get genomic sequence?\label{sub:How-to-get-genomic-sequence}}

The biotool \textbf{get\_genomic\_seq} can extract subsequences for
a given genome specified with the -\/-genome switch explicitly using
the -\/-beg and -\/-end/-\/-len switches:

\begin{lyxcode}
get\_genome\_seq~-{}-genome=<genome>~-{}-beg=1~-{}-len=100
\end{lyxcode}
Alternatively, \textbf{get\_genome\_seq} can be used to append the
corresponding sequence to BED, PSL, and BLAST records:

\begin{lyxcode}
read\_bed~-{}-data\_in=<BED~file>~|~get\_genome\_seq~-{}-genome=<genome>
\end{lyxcode}

\subsection{How to upper-case sequences?}

Sequences can be shifted from lower case to upper case using \textbf{uppercase\_seq}:

\begin{lyxcode}
...~|~uppercase\_seq
\end{lyxcode}

\subsection{How to reverse sequences?\label{sub:How-to-reverse-seq}}

The order of residues in a sequence can be reversed using reverse\_seq:

\begin{lyxcode}
...~|~reverse\_seq
\end{lyxcode}
Note that in order to reverse/complement a sequence you also need
the \textbf{complement\_seq} biotool (see~\ref{sub:How-to-complement}).


\subsection{How to complement sequences?\label{sub:How-to-complement}}

DNA and RNA sequences can be complemented with \textbf{complement\_seq},
which automagically determines the sequence type:

\begin{lyxcode}
...~|~complement\_seq
\end{lyxcode}
Note that in order to reverse/complement a sequence you also need
the \textbf{reverse\_seq} biotool (see~\ref{sub:How-to-reverse-seq}).


\subsection{How to remove indels from sequnces?}

Indels can be removed from sequences with the \textbf{remove\_indels}
biotool. This is useful if you have aligned some sequences (see~\ref{sub:How-to-align})
and extracted (see~\ref{sub:How-to-extract}) a block of subsequences
from the alignment and you want to use these sequence in a search
where you need to remove the indels first. '-', '\textasciitilde{}',
and '.' are considered indels:

\begin{lyxcode}
...~|~remove\_indels
\end{lyxcode}

\subsection{How to split sequences into overlapping subsequences?}

Sequences can be slit into overlapping subsequences with the \textbf{split\_seq}
biotool.

\begin{lyxcode}
...~|~split\_seq~-{}-word\_size=20~-{}-uniq
\end{lyxcode}

\subsection{How to determine the oligo frequency?}

In order to determine if any oligo usage is over represented in one
or more sequences you can determine the frequency of oligos of a given
size with \textbf{oligo\_freq}:

\begin{lyxcode}
...~|~oligo\_freq~-{}-word\_size=4
\end{lyxcode}
And if you have more than one sequence and want to accumulate the
frequences you need the -\/-all switch:

\begin{lyxcode}
...~|~oligo\_freq~-{}-word\_size=4~-{}-all
\end{lyxcode}
To get a meaningful result you need to write the resulting frequencies
as a table with \textbf{write\_tab} (see~\ref{sub:How-to-write-tab}),
but first it is important to \textbf{grab} (see~\ref{sub:How-to-grab})
the records with the frequencies to avoid full length sequences in
the table:

\begin{lyxcode}
...~|~oligo\_freq~-{}-word\_size=4~-{}-all~|~grab~-{}-pattern=OLIGO~-{}-keys\_only

|~write\_tab~-{}-no\_stream
\end{lyxcode}
And the resulting frequency table can be sorted with Unix sort (man
sort).


\subsection{How to search for sequences in genomes?}

See the following biotool:

\begin{itemize}
\item \textbf{patscan\_seq} \eqref{sub:How-to-use-patscan}
\item \textbf{blat\_seq} \eqref{sub:How-to-use-BLAT}
\item \textbf{blast\_seq} \eqref{sub:How-to-use-BLAST}
\item \textbf{vmatch\_seq} \eqref{sub:How-to-use-Vmatch}
\end{itemize}

\subsection{How to search sequences for a pattern?\label{sub:How-to-use-patscan}}

It is possible to search sequences in the data stream for patterns
using the \textbf{patscan\_seq} biotool which utilizes the powerful
scan\_for\_matches engine. Consult the documentation for scan\_for\_matches
in order to learn how to define patterns (the documentation is included
in Appendix~\ref{sec:scan_for_matches-README}).

To search all sequences for a simple pattern consisting of the sequence
ATCGATCG allowing for 3 mismatches, 2 insertions and 1 deletion:

\begin{lyxcode}
read\_fasta~-{}-data\_in=<file>~|~patscan\_seq~-{}-pattern='ATCGATCG{[}3,2,1]'
\end{lyxcode}
The -\/-pattern switch takes a comma seperated list of patterns,
so if you want to search for more that one pattern do:

\begin{lyxcode}
...~|~patscan\_seq~-{}-pattern='ATCGATCG{[}3,2,1],GCTAGCTA{[}3,2,1]'
\end{lyxcode}
It is also possible to have a list of different patterns to search
for in a file with one pattern per line. In order to get \textbf{patscan\_seq}
to read these patterns use the -\/-pattern\_in switch:

\begin{lyxcode}
...~|~patscan\_seq~-{}-pattern\_in=<file>
\end{lyxcode}
To also scan the complementary strand in nucleotide sequences (\textbf{patscan\_seq}
automagically determines the sequence type) you need to add the -\/-comp
switch:

\begin{lyxcode}
...~|~patscan\_seq~-{}-pattern=<pattern>~-{}-comp
\end{lyxcode}
It is also possible to use \textbf{patscan\_seq} to output those records
that does not contain a certain pattern by using the -\/-invert switch:

\begin{lyxcode}
...~|~patscan\_seq~-{}-pattern=<pattern>~-{}-invert
\end{lyxcode}
Finally, \textbf{patscan\_seq} can also scan for patterns in a given
genome sequence, instead of sequences in the stream, using the -\/-genome
switch:

\begin{lyxcode}
patscan~-{}-pattern=<pattern>~-{}-genome=<genome>
\end{lyxcode}

\subsection{How to use BLAT for sequence search?\label{sub:How-to-use-BLAT}}

Sequences in the data stream can be matched against supported genomes
using \textbf{blat\_seq} which is a biotool using BLAT as the name
might suggest. Currently only Mouse and Human genomes are available
and it is not possible to use OOC files since there is still a need
for a local repository for genome files. Otherwise it is just:

\begin{lyxcode}
read\_fasta~-{}-data\_in=<file>~|~blat\_seq~-{}-genome=<genome>
\end{lyxcode}
The search results can then be written to file with \textbf{write\_psl}
(see~\ref{sub:How-to-write-PSL}) or \textbf{write\_bed} (see~\ref{sub:How-to-write-BED})
allthough with \textbf{write\_bed} some information will be lost).
It is also possible to plot chromosome distribution of the search
results using \textbf{plot\_chrdist} (see~\ref{sub:How-to-plot-chrdist})
or the distribution of the match lengths using \textbf{plot\_lendist}
(see~\ref{sub:How-to-plot-lendist}) or a karyogram with the hits
using \textbf{plot\_karyogram} (see~\ref{sub:How-to-plot-karyogram}).


\subsection{How to use BLAST for sequence search?\label{sub:How-to-use-BLAST}}

Two biotools exist for blasting sequences: \textbf{create\_blast\_db}
is used to create the BLAST database required for BLAST which is queried
using the biotool \textbf{blast\_seq}. So in order to create a BLAST
database from sequences in the data stream you simple run:

\begin{lyxcode}
...~|~create\_blast\_db~-{}-database=my\_database~-{}-no\_stream
\end{lyxcode}
The type of sequence to use for the database is automagically determined
by \textbf{create\_blast\_db}, but don't have a mixture of peptide
and nucleic acids sequences in the stream. The -\/-database switch
takes a path as argument, but will default to 'blastdb\_<time\_stamp>
if not set.

The resulting database can now be queried with sequences in another
data stream using \textbf{blast\_seq}:

\begin{lyxcode}
...~|~blast\_seq~-{}-database=my\_database
\end{lyxcode}
Again, the sequence type is determined automagically and the appropriate
BLAST program is guessed (see below table), however, the program name
can be overruled with the -\/-program switch.

\noindent \begin{center}
\begin{tabular}{ccc}
Subject sequence & Query sequence & Program guess\tabularnewline
\hline
Nucleotide & Nucleotide & blastn\tabularnewline
Protein & Protein & blastp\tabularnewline
Protein & Nucleotide & blastx\tabularnewline
Nucleotide & Protein & tblastn\tabularnewline
\end{tabular}
\par\end{center}

Finally, it is also possible to use \textbf{blast\_seq} for blasting
sequences agains a preformatted genome using the -\/-genome switch
instead of the -\/-database switch:

\begin{lyxcode}
...~|~blast\_seq~-{}-genome=<genome>
\end{lyxcode}

\subsection{How to use Vmatch for sequence search?\label{sub:How-to-use-Vmatch}}

The powerful suffix array software package Vmatch%
\footnote{\url{http://www.vmatch.de/}%
} can be used for exact mapping of sequences against indexed genomes
using the biotool \textbf{vmatch\_seq}, which will e.g. map 700000
ESTs to the human genome locating all 160 mio hits in less than an
hour.

\begin{lyxcode}
...~|~vmatch\_seq~-{}-genome=<genome>
\end{lyxcode}
Only nucleotide sequences and sequences longer than 11 nucleotides
will be mapped. The resulting SCORE key will hold the number of genome
matches of a given sequence (multi-mappers).


\subsection{How to find all matches between sequences?\label{sub:How-to-find-matches}}

All matches between two sequences can be determined with the biotool
\textbf{match\_seq}. The match finding engine underneath the hood
of \textbf{match\_seq} is the super fast suffix tree program MUMmer%
\footnote{\url{http://mummer.sourceforge.net/}%
}, which will locate all forward and reverse matches between huge sequences
in a matter of minutes (if the repeat count is not too high and if
the word size used is appropriate). Matching two \emph{Helicobacter
pylori} genomes (1.7Mbp) takes around 10 seconds:

\begin{lyxcode}
...~|~match\_seq~-{}-word\_size=20~-{}-direction=both
\end{lyxcode}
The output from \textbf{match\_seq} can be used to generate a dot
plot with \textbf{plot\_matches} (see~\ref{sub:How-to-generate-dotplot}).


\subsection{How to align sequences?\label{sub:How-to-align}}

Sequences in the stream can be aligned with the \textbf{align\_seq}
biotool that uses Muscle%
\footnote{\url{http://www.drive5.com/muscle/muscle.html}%
} as aligment engine. Currently you cannot change any of the Muscle
alignment parameters and \textbf{align\_seq} will create an alignment
based on the defaults (which are really good!):

\begin{lyxcode}
...~|~align\_seq
\end{lyxcode}
The aligned output can be written to file in FASTA format using \textbf{write\_fasta}
(see~\ref{sub:How-to-write-fasta}) or in pretty text using \textbf{write\_align}
(see~\ref{sub:How-to-write-alignment}).


\subsection{How to create a weight matrix?}

If you want a weight matrix to show the sequence composition of a
stack of sequences you can use the biotool create\_weight\_matrix:

\begin{lyxcode}
...~|~create\_weight\_matrix
\end{lyxcode}
The result can be output in percent using the -\/-percent switch:

\begin{lyxcode}
...~|~create\_weight\_matrix~-{}-percent
\end{lyxcode}
The weight matrix can be written as tabular output with \textbf{write\_tab}
(see~\ref{sub:How-to-write-tab}) after removeing the records containing
SEQ with \textbf{grab} (see~\ref{sub:How-to-grab}):

\begin{lyxcode}
...~|~create\_weight\_matrix~|~grab~-{}-invert~-{}-keys=SEQ~-{}-keys\_only

|~write\_tab~-{}-no\_stream
\end{lyxcode}
The V0 column will hold the residue, while the rest of the columns
will hold the frequencies for each sequence position.


\section{Plotting}

There exists several biotools for plotting. Some of these are based
on GNUplot%
\footnote{\url{http://www.gnuplot.info/}%
}, which is an extremely powerful platform to generate all sorts of
plots and even though GNUplot has quite a steep learning curve, the
biotools utilizing GNUplot are simple to use. GNUplot is able to output
a lot of different formats (called terminals in GNUplot), but the
biotools focusses on three formats only:

\begin{enumerate}
\item The 'dumb' terminal is default to the GNUplot based biotools and will
output a plot in crude ASCII text (Fig.~\ref{fig:Dumb-terminal}).
This is quite nice for a quick and dirty plot to get an overview of
your data .
\item The 'post' or 'postscript' terminal output postscript code which is
publication grade graphics that can be viewed with applications such
as Ghostview, Photoshop, and Preview.
\item The 'svg' terminal output's scalable vector graphics (SVG) which is
a vector based format. SVG is great because you can edit the resulting
plot using Photoshop or Inkscape%
\footnote{Inkscape is a really handy drawing program that is free and open source.
Availble at \url{http://www.inkscape.org}%
} if you want to add additional labels, captions, arrows, and so on
and then save the result in different formats, such as postscript
without loosing resolution.
\end{enumerate}
The biotools for plotting that are not based on GNUplot only output
SVG (that may change in the future).

%
\begin{figure}
\noindent \begin{centering}
\includegraphics[width=12cm]{lendist_ascii}
\par\end{centering}

\caption{\label{fig:Dumb-terminal}Dumb terminal}


\begin{quote}
The output of a length distribution plot in the default 'dumb terminal'
to the terminal window. 
\end{quote}

\end{figure}



\subsection{How to plot a histogram?\label{How-to-plot-histogram}}

A generic histogram for a given value can be plotted with the biotool
\textbf{plot\_histogram} (Fig.~\ref{fig:Histogram}):

\begin{lyxcode}
...~|~plot\_histogram~-{}-key=TISSUE~-{}-no\_stream
\end{lyxcode}
(Figure missing)

\noindent \begin{flushleft}
%
\begin{figure}
\noindent \begin{centering}
\includegraphics[width=12cm]{histogram}
\par\end{centering}

\caption{\label{fig:Histogram}Histogram}

\end{figure}

\par\end{flushleft}


\subsection{How to plot a length distribution?\label{sub:How-to-plot-lendist}}

Plotting of length distributions, weather sequence lengths, patterns
lengths, hit lengths, \emph{etc.} is a really handy thing and can
be done with the the biotool \textbf{plot\_lendist}. If you have a
file with FASTA entries and want to plot the length distribution you
do it like this:

\begin{lyxcode}
read\_fasta~-{}-data\_in=<file>~|~length\_seq

|~plot\_lendist~-{}-key=SEQ\_LEN~-{}-no\_stream
\end{lyxcode}
The result will be written to the default dumb terminal and will look
like Fig.~\ref{fig:Dumb-terminal}.

If you instead want the result in postscript format you can do:

\begin{lyxcode}
...~|~plot\_lendist~-{}-key=SEQ\_LEN~-{}-terminal=post~-{}-result\_out=file.ps
\end{lyxcode}
That will generate the plot and save it to file, but not interrupt
the data stream which can then be used in further analysis. You can
also save the plot implicetly using '>', however, it is then important
to terminate the stream with the -\/-no\_stream switch:

\begin{lyxcode}
...~|~plot\_lendist~-{}-key=SEQ\_LEN~-{}-terminal=post~-{}-no\_stream~>~file.ps
\end{lyxcode}
The resulting plot can be seen in Fig.~\ref{fig:Length-distribution}.

%
\begin{figure}


\noindent \begin{centering}
\includegraphics[width=12cm]{lendist}
\par\end{centering}

\caption{\label{fig:Length-distribution}Length distribution}


\begin{quote}
Length distribution of 630 piRNA like RNAs.
\end{quote}

\end{figure}



\subsection{How to plot a chromosome distribution?\label{sub:How-to-plot-chrdist}}

If you have the result of a sequence search against a multi chromosome
genome, it is very practical to be able to plot the distribution of
search hits on the different chromosomes. This can be done with \textbf{plot\_chrdist}:

\begin{lyxcode}
read\_fasta~-{}-data\_in=<file>~|~blat\_genome~|~plot\_chrdist~-{}-no\_stream
\end{lyxcode}
The above example will result in a crude plot using the 'dumb' terminal,
and if you want to mess around with the results from the BLAT search
you probably want to save the result to file first (see~\ref{sub:How-to-write-PSL}).
To plot the chromosome distribution from the saved search result you
can do:

\begin{lyxcode}
read\_bed~-{}-data\_in=file.bed~|~plot\_chrdist~-{}-terminal=post~-{}-result\_out=plot.ps
\end{lyxcode}
That will result in the output show in Fig.~\ref{fig:Chromosome-distribution}.

%
\begin{figure}


\noindent \begin{centering}
\includegraphics[angle=90,width=12cm]{chrdist}
\par\end{centering}

\caption{\label{fig:Chromosome-distribution}Chromosome distribution}

\end{figure}



\subsection{How to generate a dotplot?\label{sub:How-to-generate-dotplot}}

A dotplot is a powerful way to get an overview of the size and location
of sequence insertions, deletions, and duplications between two sequences.
Generating a dotplot with biotools is a two step process where you
initially find all matches between two sequences using the tool \textbf{match\_seq}
(see~\ref{sub:How-to-find-matches}) and plot the resulting matches
with \textbf{plot\_matches}. Matching and plotting two \emph{Helicobacter
pylori} genomes (1.7Mbp) takes around 10 seconds:

\begin{lyxcode}
...~|~match\_seq~|~plot\_matches~-{}-terminal=post~-{}-result\_out=plot.ps
\end{lyxcode}
The resulting dotplot is in Fig.~\ref{fig:Dotplot}.

%
\begin{figure}
\noindent \begin{centering}
\includegraphics[width=12cm]{dotplot}
\par\end{centering}

\caption{\label{fig:Dotplot}Dotplot}


\begin{quote}
Forward matches are displayed in green while reverse matches are displayed
in red.
\end{quote}

\end{figure}



\subsection{How to plot a sequence logo?}

Sequence logos can be generate with \textbf{plot\_seqlogo}. The sequnce
type is determined automagically and an entropy scale of 2 bits and
4 bits is used for nucleotide and peptide sequences, respectively%
\footnote{\url{http://www.ccrnp.ncifcrf.gov/~toms/paper/hawaii/latex/node5.html}%
}.

\begin{lyxcode}
...~|~plot\_seqlogo~-{}-no\_stream~-{}-result\_out=seqlogo.svg
\end{lyxcode}
An example of a sequence logo can be seen in Fig.~\ref{fig:Sequence-logo}.

%
\begin{figure}
\noindent \begin{centering}
\includegraphics[width=12cm]{seqlogo}
\par\end{centering}

\caption{\label{fig:Sequence-logo}Sequence logo}

\end{figure}



\subsection{How to plot a karyogram?\label{sub:How-to-plot-karyogram}}

To plot search hits on genomes use \textbf{plot\_karyogram}, which
will output a nice karyogram in SVG graphics:

\begin{lyxcode}
...~|~plot\_karyogram~-{}-result\_out=karyogram.svg
\end{lyxcode}
The banding data is taken from the UCSC genome browser database and
currently only Human and Mouse is supported. Fig.~\ref{fig:Karyogram}
shows the distribution of piRNA like RNAs matched to the Human genome.

%
\begin{figure}
\noindent \begin{centering}
\includegraphics[width=12cm]{karyogram}
\par\end{centering}

\caption{\label{fig:Karyogram}Karyogram}


\begin{quote}
Hits from a search of piRNA like RNAs in the Human genome is displayed
as short horizontal bars.
\end{quote}

\end{figure}



\section{Uploading Results}


\subsection{How do I display my results in the UCSC Genome Browser?}

Results from the list of biotools below can be uploaded directly to
a local mirror of the UCSC Genome Browser using the biotool \textbf{upload\_to\_ucsc}:

\begin{itemize}
\item patscan\_seq \eqref{sub:How-to-use-patscan}
\item blat\_seq \eqref{sub:How-to-use-BLAT}
\item blast\_seq \eqref{sub:How-to-use-BLAST}
\item vmatch\_seq \eqref{sub:How-to-use-Vmatch}
\end{itemize}
The syntax for uploading data the most simple way requires two mandatory
switches: -\/-database, which is the UCSC database name (such as
hg18, mm9, etc.) and-\/-table which should be the users initials
followed by an underscore and a short description of the data:

\begin{lyxcode}
...~|~upload\_to\_ucsc~-{}-database=hg18~-{}-table=mah\_snoRNAs
\end{lyxcode}
The \textbf{upload\_to\_ucsc} biotool modifies the users \textasciitilde{}/ucsc/my\_tracks.ra
file automagically (a backup is created with the name \textasciitilde{}/ucsc/my\_tracks.ra\textasciitilde{})
with default values that can be overridden using the following switches:

\begin{itemize}
\item -\/-short\_label - Short label for track - Default=database->table
\item -\/-long\_label - Long label for track - Default=database->table
\item -\/-group - Track group name - Default=<user name as defined in env>
\item -\/-priority - Track display priority - Default=1
\item -\/-color - Track color - Default=147,73,42
\item -\/-chunk\_size - Chunks for loading - Default=10000000
\item -\/-visibility - Track visibility - Default=pack
\end{itemize}
Also, data in BED or PSL format can be uploaded with \textbf{upload\_to\_ucsc}
as long as these reference to genomes and chromosomes existing in
the UCSC Genome Browser:

\begin{lyxcode}
read\_bed~-{}-data\_in=<bed~file>~|~upload\_to\_ucsc~...



read\_psl~-{}-data\_in=<psl~file>~|~upload\_to\_ucsc~...
\end{lyxcode}

\section{Trouble shooting}

Shoot the messenger!

\appendix

\section{Keys\label{sec:Keys}}

HIT

HIT\_BEG

HIT\_END

HIT\_LEN

HIT\_NAME

PATTERN


\section{Switches\label{sec:Switches}}

-\/-stream\_in

-\/-stream\_out

-\/-no\_stream

-\/-data\_in

-\/-result\_out

-\/-num


\section{scan\_for\_matches README\label{sec:scan_for_matches-README}}

\begin{lyxcode}
~~~~~~~~~~~~~~~~~~~~~~~~~~scan\_for\_matches:

~~~~A~Program~to~Scan~Nucleotide~or~Protein~Sequences~for~Matching~Patterns

~~~~~~~~~~~~~~~~~~~~~~~~Ross~Overbeek

~~~~~~~~~~~~~~~~~~~~~~~~MCS

~~~~~~~~~~~~~~~~~~~~~~~~Argonne~National~Laboratory

~~~~~~~~~~~~~~~~~~~~~~~~Argonne,~IL~60439

~~~~~~~~~~~~~~~~~~~~~~~~USA

Scan\_for\_matches~is~a~utility~that~we~have~written~to~search~for

patterns~in~DNA~and~protein~sequences.~~I~wrote~most~of~the~code,

although~David~Joerg~and~Morgan~Price~wrote~sections~of~an

earlier~version.~~The~whole~notion~of~pattern~matching~has~a~rich

history,~and~we~borrowed~liberally~from~many~sources.~~However,~it~is

worth~noting~that~we~were~strongly~influenced~by~the~elegant~tools

developed~and~distributed~by~David~Searls.~~My~intent~is~to~make~the

existing~tool~available~to~anyone~in~the~research~community~that~might

find~it~useful.~~I~will~continue~to~try~to~fix~bugs~and~make~suggested

enhancements,~at~least~until~I~feel~that~a~superior~tool~exists.

Hence,~I~would~appreciate~it~if~all~bug~reports~and~suggestions~are

directed~to~me~at~Overbeek@mcs.anl.gov.~~

I~will~try~to~log~all~bug~fixes~and~report~them~to~users~that~send~me

their~email~addresses.~~I~do~not~require~that~you~give~me~your~name

and~address.~~However,~if~you~do~give~it~to~me,~I~will~try~to~notify

you~of~serious~problems~as~they~are~discovered.

Getting~Started:

~~~~The~distribution~should~contain~at~least~the~following~programs:

~~~~~~~~~~~~~~~~README~~~~~~~~~~~~~~~~~~-~~~~~This~document

~~~~~~~~~~~~~~~~ggpunit.c~~~~~~~~~~~~~~~-~~~~~One~of~the~two~source~files

~~~~~~~~~~~~~~~~scan\_for\_matches.c~~~~~~-~~~~~The~second~source~file

~~~~~~~~~~~~~~~~

~~~~~~~~~~~~~~~~run\_tests~~~~~~~~~~~~~~~-~~~~~A~perl~script~to~test~things

~~~~~~~~~~~~~~~~show\_hits~~~~~~~~~~~~~~~-~~~~~A~handy~perl~script

~~~~~~~~~~~~~~~~test\_dna\_input~~~~~~~~~~-~~~~~Test~sequences~for~DNA

~~~~~~~~~~~~~~~~test\_dna\_patterns~~~~~~~-~~~~~Test~patterns~for~DNA~scan

~~~~~~~~~~~~~~~~test\_output~~~~~~~~~~~~~-~~~~~Desired~output~from~test

~~~~~~~~~~~~~~~~test\_prot\_input~~~~~~~~~-~~~~~Test~protein~sequences

~~~~~~~~~~~~~~~~test\_prot\_patterns~~~~~~-~~~~~Test~patterns~for~proteins

~~~~~~~~~~~~~~~~testit~~~~~~~~~~~~~~~~~~-~~~~~a~perl~script~used~for~test

~~~~Only~the~first~three~files~are~required.~~The~others~are~useful,

~~~~but~only~if~you~have~Perl~installed~on~your~system.~~If~you~do

~~~~have~Perl,~I~suggest~that~you~type

~~~~~~~~

~~~~~~~~~~~~~~~~which~perl

~~~~to~find~out~where~it~installed.~~On~my~system,~I~get~the~following

~~~~response:

~~~~~~~~

~~~~~~~~~~~~~~~~clone\%~which~perl

~~~~~~~~~~~~~~~~/usr/local/bin/perl

~~~~indicating~that~Perl~is~installed~in~/usr/local/bin.~~Anyway,~once

~~~~you~know~where~it~is~installed,~edit~the~first~line~of~files~

~~~~~~~~testit

~~~~~~~~show\_hits

~~~~replacing~/usr/local/bin/perl~with~the~appropriate~location.~~I

~~~~will~assume~that~you~can~do~this,~although~it~is~not~critical~(it

~~~~is~needed~only~to~test~the~installation~and~to~use~the~\char`\"{}show\_hits\char`\"{}

~~~~utility).~~Perl~is~not~required~to~actually~install~and~run

~~~~scan\_for\_matches.~

~~~~If~you~do~not~have~Perl,~I~suggest~you~get~it~and~install~it~(it

~~~~is~a~wonderful~utility).~~Information~about~Perl~and~how~to~get~it

~~~~can~be~found~in~the~book~\char`\"{}Programming~Perl\char`\"{}~by~Larry~Wall~and

~~~~Randall~L.~Schwartz,~published~by~O'Reilly~\&~Associates,~Inc.

~~~~To~get~started,~you~will~need~to~compile~the~program.~~~I~do~this

~~~~using~

~~~~~~~~gcc~-O~-o~scan\_for\_matches~~ggpunit.c~scan\_for\_matches.c

~~~~If~you~do~not~use~GNU~C,~use~

~~~~~~~~cc~-O~-DCC~-o~scan\_for\_matches~~ggpunit.c~scan\_for\_matches.c

~~~~which~works~on~my~Sun.~~

~~~~Once~you~have~compiled~scan\_for\_matches,~you~can~verify~that~it

~~~~works~with

~~~~~~~~clone\%~run\_tests~tmp

~~~~~~~~clone\%~diff~tmp~test\_output

~~~~You~may~get~a~few~strange~lines~of~the~sort

~~~~~~~~clone\%~run\_tests~tmp

~~~~~~~~rm:~tmp:~No~such~file~or~directory

~~~~~~~~clone\%~diff~tmp~test\_output

~~~~These~should~cause~no~concern.~~However,~if~the~\char`\"{}diff\char`\"{}~shows~that

~~~~tmp~and~test\_output~are~different,~contact~me~(you~have~a

~~~~problem).~

~~~~You~should~now~be~able~to~use~scan\_for\_matches~by~following~the

~~~~instructions~given~below~(which~is~all~the~normal~user~should~have

~~~~to~understand,~once~things~are~installed~properly).

~==============================================================

How~to~run~scan\_for\_matches:

~~~~To~run~the~program,~you~type~need~to~create~two~files

~~~~1.~~the~first~file~contains~the~pattern~you~wish~to~scan~for;~I'll

~~~~~~~~call~this~file~pat\_file~in~what~follows~(but~any~name~is~ok)

~~~~2.~~the~second~file~contains~a~set~of~sequences~to~scan.~~These

~~~~~~~~should~be~in~\char`\"{}fasta~format\char`\"{}.~~Just~look~at~the~contents~of

~~~~~~~~test\_dna\_input~to~see~examples~of~this~format.~~Basically,

~~~~~~~~each~sequence~begins~with~a~line~of~the~form

~~~~~~~~~~~>sequence\_id

~~~~~~~~and~is~followed~by~one~or~more~lines~containing~the~sequence.

~~~~Once~these~files~have~been~created,~you~just~use

~~~~~~~~scan\_for\_matches~pat\_file~<~input\_file

~~~~to~scan~all~of~the~input~sequences~for~the~given~pattern.~~As~an

~~~~example,~suppose~that~pat\_file~contains~a~single~line~of~the~form

~~~~~~~~~~~~~~~~p1=4...7~3...8~\textasciitilde{}p1

~~~~Then,

~~~~~~~~~~~~~~~~scan\_for\_matches~pat\_file~<~test\_dna\_input

~~~~should~produce~two~\char`\"{}hits\char`\"{}.~~When~I~run~this~on~my~machine,~I~get

~~~~~~~~clone\%~scan\_for\_matches~pat\_file~<~test\_dna\_input

~~~~~~~~>tst1:{[}6,27]

~~~~~~~~cguaacc~ggttaacc~gguuacg~

~~~~~~~~>tst2:{[}6,27]

~~~~~~~~CGUAACC~GGTTAACC~GGUUACG~

~~~~~~~~clone\%~

Simple~Patterns~Built~by~Matching~Ranges~and~Reverse~Complements

~~~~Let~me~first~explain~this~simple~pattern:

~~~~~~~~~~~~~~~~

~~~~~~~~~~~~~~~~p1=4...7~3...8~\textasciitilde{}p1

~~~~The~pattern~consists~of~three~\char`\"{}pattern~units\char`\"{}~separated~by~spaces.

~~~~The~first~pattern~unit~is

~~~~~~~~~~~~~~~~p1=4...7

~~~~which~means~\char`\"{}match~4~to~7~characters~and~call~them~p1\char`\"{}.~~The

~~~~second~pattern~unit~is

~~~~~~~~~~~~~~~~3...8

~~~~which~means~\char`\"{}then~match~3~to~8~characters\char`\"{}.~~The~last~pattern~unit

~~~~is~

~~~~~~~~~~~~~~~~\textasciitilde{}p1

~~~~which~means~\char`\"{}match~the~reverse~complement~of~p1\char`\"{}.~~The~first

~~~~reported~hit~is~shown~as

~~~~~~~~>tst1:{[}6,27]

~~~~~~~~cguaacc~ggttaacc~gguuacg~

~~~~which~states~that~characters~6~through~27~of~sequence~tst1~were

~~~~matched.~~\char`\"{}cguaac\char`\"{}~matched~the~first~pattern~unit,~\char`\"{}ggttaacc\char`\"{}~the

~~~~second,~and~\char`\"{}gguuacg\char`\"{}~the~third.~~This~is~an~example~of~a~common

~~~~type~of~pattern~used~to~search~for~sections~of~DNA~or~RNA~that

~~~~would~fold~into~a~hairpin~loop.

Searching~Both~Strands

~~~~Now~for~a~short~aside:~scan\_for\_matches~only~searched~the

~~~~sequences~in~the~input~file;~it~did~not~search~the~opposite

~~~~strand.~~With~a~pattern~of~the~sort~we~just~used,~there~is~not

~~~~need~o~search~the~opposite~strand.~~However,~it~is~normally~the

~~~~case~that~you~will~wish~to~search~both~the~sequence~and~the

~~~~opposite~strand~(i.e.,~the~reverse~complement~of~the~sequence).

~~~~To~do~that,~you~would~just~use~the~\char`\"{}-c\char`\"{}~command~line.~~For~example,

~~~~~~~~scan\_for\_matches~-c~pat\_file~<~test\_dna\_input

~~~~Hits~on~the~opposite~strand~will~show~a~beginning~location~greater

~~~~than~te~end~location~of~the~match.

Defining~Pairing~Rules~and~Allowing~Mismatches,~Insertions,~and~Deletions

~~~~Let~us~stop~now~and~ask~\char`\"{}What~additional~features~would~one~need~to

~~~~really~find~the~kinds~of~loop~structures~that~characterize~tRNAs,

~~~~rRNAs,~and~so~forth?\char`\"{}~~I~can~immediately~think~of~two:

~~~~~~~~a)~you~will~need~to~be~able~to~allow~non-standard~pairings

~~~~~~~~~~~(those~other~than~G-C~and~A-U),~and

~~~~~~~~b)~you~will~need~to~be~able~to~tolerate~some~number~of

~~~~~~~~~~~mismatches~and~bulges.

~~~~~~~~

~~~~Let~me~first~show~you~how~to~handle~non-standard~\char`\"{}rules~for

~~~~pairing~in~reverse~complements\char`\"{}.~~Consider~the~following~pattern,

~~~~which~I~show~as~two~line~(you~may~use~as~many~lines~as~you~like~in

~~~~forming~a~pattern,~although~you~can~only~break~a~pattern~at~points

~~~~where~space~would~be~legal):

~~~~~~~~~~~~r1=\{au,ua,gc,cg,gu,ug,ga,ag\}~

~~~~~~~~~~~~p1=2...3~0...4~p2=2...5~1...5~r1\textasciitilde{}p2~0...4~\textasciitilde{}p1~~~~~~~~

~~~~The~first~\char`\"{}pattern~unit\char`\"{}~does~not~actually~match~anything;~rather,

~~~~it~defines~a~\char`\"{}pairing~rule\char`\"{}~in~which~standard~pairings~are

~~~~allowed,~as~well~as~G-A~and~A-G~(in~case~you~wondered,~Us~and~Ts

~~~~and~upper~and~lower~case~can~be~used~interchangably;~for~example

~~~~r1=\{AT,UA,gc,cg\}~could~be~used~to~define~the~\char`\"{}standard~rule\char`\"{}~for

~~~~pairings).~~The~second~line~consists~of~six~pattern~units~which

~~~~may~be~interpreted~as~follows:

~~~~~~~~~~~~p1=2...3~~~~~match~2~or~3~characters~(call~it~p1)

~~~~~~~~~~~~0...4~~~~~~~~match~0~to~4~characters

~~~~~~~~~~~~p2=2...5~~~~~match~2~to~5~characters~(call~it~p2)

~~~~~~~~~~~~1...5~~~~~~~~match~1~to~5~characters

~~~~~~~~~~~~r1\textasciitilde{}p2~~~~~~~~match~the~reverse~complement~of~p2,

~~~~~~~~~~~~~~~~~~~~~~~~~~~~allowing~G-A~and~A-G~pairs

~~~~~~~~~~~~0...4~~~~~~~~match~0~to~4~characters~~~~~~~~

~~~~~~~~~~~~\textasciitilde{}p1~~~~~~~~~~match~the~reverse~complement~of~p1

~~~~~~~~~~~~~~~~~~~~~~~~~~~~allowing~only~G-C,~C-G,~A-T,~and~T-A~pairs

~~~~Thus,~r1\textasciitilde{}p2~means~\char`\"{}match~the~reverse~complement~of~p2~using~rule~r1\char`\"{}.

~~~~Now~let~us~consider~the~issue~of~tolerating~mismatches~and~bulges.

~~~~You~may~add~a~\char`\"{}qualifier\char`\"{}~to~the~pattern~unit~that~gives~the

~~~~tolerable~number~of~\char`\"{}mismatches,~deletions,~and~insertions\char`\"{}.

~~~~Thus,

~~~~~~~~~~~~~~~~p1=10...10~3...8~\textasciitilde{}p1{[}1,2,1]

~~~~means~that~the~third~pattern~unit~must~match~10~characters,

~~~~allowing~one~\char`\"{}mismatch\char`\"{}~(a~pairing~other~than~G-C,~C-G,~A-T,~or

~~~~T-A),~two~deletions~(a~deletion~is~a~character~that~occurs~in~p1,

~~~~but~has~been~\char`\"{}deleted\char`\"{}~from~the~string~matched~by~\textasciitilde{}p1),~and~one

~~~~insertion~(an~\char`\"{}insertion\char`\"{}~is~a~character~that~occurs~in~the~string

~~~~matched~by~\textasciitilde{}p1,~but~not~for~which~no~corresponding~character

~~~~occurs~in~p1).~~In~this~case,~the~pattern~would~match

~~~~~~~~~~~~~~ACGTACGTAC~GGGGGGGG~GCGTTACCT

~~~~which~is,~you~must~admit,~a~fairly~weak~loop.~~It~is~common~to

~~~~allow~mismatches,~but~you~will~find~yourself~using~insertions~and

~~~~deletions~much~more~rarely.~~In~any~event,~you~should~note~that

~~~~allowing~mismatches,~insertions,~and~deletions~does~force~the

~~~~program~to~try~many~additional~possible~pairings,~so~it~does~slow

~~~~things~down~a~bit.

How~Patterns~Are~Matched

~~~~Now~is~as~good~a~time~as~any~to~discuss~the~basic~flow~of~control

~~~~when~matching~patterns.~~Recall~that~a~\char`\"{}pattern\char`\"{}~is~a~sequence~of

~~~~\char`\"{}pattern~units\char`\"{}.~~Suppose~that~the~pattern~units~were

~~~~~~~~u1~u2~u3~u4~...~un

~~~~The~scan~of~a~sequence~S~begins~by~setting~the~current~position

~~~~to~1.~~Then,~an~attempt~is~made~to~match~u1~starting~at~the

~~~~current~position.~~Each~attempt~to~match~a~pattern~unit~can

~~~~succeed~or~fail.~~If~it~succeeds,~then~an~attempt~is~made~to~match

~~~~the~next~unit.~~If~it~fails,~then~an~attempt~is~made~to~find~an

~~~~alternative~match~for~the~immediately~preceding~pattern~unit.~~If

~~~~this~succeeds,~then~we~proceed~forward~again~to~the~next~unit.~~If

~~~~it~fails~we~go~back~to~the~preceding~unit.~~This~process~is~called

~~~~\char`\"{}backtracking\char`\"{}.~~If~there~are~no~previous~units,~then~the~current

~~~~position~is~incremented~by~one,~and~everything~starts~again.~~This

~~~~proceeds~until~either~the~current~position~goes~past~the~end~of

~~~~the~sequence~or~all~of~the~pattern~units~succeed.~~On~success,

~~~~scan\_for\_matches~reports~the~\char`\"{}hit\char`\"{},~the~current~position~is~set

~~~~just~past~the~hit,~and~an~attempt~is~made~to~find~another~hit.

~~~~If~you~wish~to~limit~the~scan~to~simply~finding~a~maximum~of,~say,

~~~~10~hits,~you~can~use~the~-n~option~(-n~10~would~set~the~limit~to

~~~~10~reported~hits).~~For~example,

~~~~~~~~scan\_for\_matches~-c~-n~1~pat\_file~<~test\_dna\_input

~~~~would~search~for~just~the~first~hit~(and~would~stop~searching~the

~~~~current~sequences~or~any~that~follow~in~the~input~file).

Searching~for~repeats:

~~~~In~the~last~section,~I~discussed~almost~all~of~the~details

~~~~required~to~allow~you~to~look~for~repeats.~~Consider~the~following

~~~~set~of~patterns:

~~~~~~~~p1=6...6~3...8~p1~~~(find~exact~6~character~repeat~separated

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~by~to~8~characters)

~~~~~~~~p1=6...6~3..8~p1{[}1,0,0]~~~(allow~one~mismatch)

~~~~~~~~p1=3...3~p1{[}1,0,0]~p1{[}1,0,0]~p1{[}1,0,0]~~

~~~~~~~~~~~~~~~~~~~~~~~~~~~~(match~12~characters~that~are~the~remains

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~of~a~3-character~sequence~occurring~4~times)

~~~~~~~~~~~~~~~~

~~~~~~~~p1=4...8~0...3~p2=6...8~p1~0...3~p2

~~~~~~~~~~~~~~~~~~~~~~~~~~~~(This~would~match~things~like

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ATCT~G~TCTTT~ATCT~TG~TCTTT

~~~~~~~~~~~~~~~~~~~~~~~~~~~~)

Searching~for~particular~sequences:

~~~~Occasionally,~one~wishes~to~match~a~specific,~known~sequence.

~~~~In~such~a~case,~you~can~just~give~the~sequence~(along~with~an

~~~~optional~statement~of~the~allowable~mismatches,~insertions,~and

~~~~deletions).~~Thus,

~~~~~~~~p1=6...8~GAGA~\textasciitilde{}p1~~~~(match~a~hairpin~with~GAGA~as~the~loop)

~~~~~~~~RRRRYYYY~~~~~~~~~~~~~(match~4~purines~followed~by~4~pyrimidines)

~~~~~~~~TATAA{[}1,0,0]~~~~~~~~~(match~TATAA,~allowing~1~mismatch)

~~~~~~~~

Matches~against~a~\char`\"{}weight~matrix\char`\"{}:

~~~~I~will~conclude~my~examples~of~the~types~of~pattern~units

~~~~available~for~matching~against~nucleotide~sequences~by~discussing~a

~~~~crude~implemetation~of~matching~using~a~\char`\"{}weight~matrix\char`\"{}.~~While~I

~~~~am~less~than~overwhelmed~with~the~syntax~that~I~chose,~I~think~that

~~~~the~reader~should~be~aware~that~I~was~thinking~of~generating

~~~~patterns~containing~such~pattern~units~automatically~from

~~~~alignments~(and~did~not~really~plan~on~typing~such~things~in~by

~~~~hand~very~often).~~Anyway,~suppose~that~you~wanted~to~match~a

~~~~sequence~of~eight~characters.~~The~\char`\"{}consensus\char`\"{}~of~these~eight

~~~~characters~is~GRCACCGS,~but~the~actual~\char`\"{}frequencies~of~occurrence\char`\"{}

~~~~are~given~in~the~matrix~below.~~Thus,~the~first~character~is~an~A

~~~~16\%~the~time~and~a~G~84\%~of~the~time.~~The~second~is~an~A~57\%~of

~~~~the~time,~a~C~10\%~of~the~time,~a~G~29\%~of~the~time,~and~a~T~4\%~of

~~~~the~time.~~

~~~~~~~~~~~~~C1~~~~~C2~~~~C3~~~~C4~~~C5~~~~C6~~~~C7~~~~C8

~~~~

~~~~~~~A~~~~~16~~~~~57~~~~~0~~~~95~~~~0~~~~18~~~~~0~~~~~0

~~~~~~~C~~~~~~0~~~~~10~~~~80~~~~~0~~100~~~~60~~~~~0~~~~50

~~~~~~~G~~~~~84~~~~~29~~~~~0~~~~~0~~~~0~~~~20~~~100~~~~50

~~~~~~~T~~~~~~0~~~~~~4~~~~20~~~~~5~~~~0~~~~~2~~~~~0~~~~~0~~~

~~~~

~~~~One~could~use~the~following~pattern~unit~to~search~for~inexact

~~~~matches~related~to~such~a~\char`\"{}weight~matrix\char`\"{}:

~~~~~~~~\{(16,0,84,0),(57,10,29,4),(0,80,0,20),(95,0,0,5),

~~~~~~~~~(0,100,0,0),(18,60,20,2),(0,0,100,0),(0,50,50,0)\}~>~450

~~~~This~pattern~unit~will~attempt~to~match~exactly~eight~characters.

~~~~For~each~character~in~the~sequence,~the~entry~in~the~corresponding

~~~~tuple~is~added~to~an~accumulated~sum.~~If~the~sum~is~greater~than

~~~~450,~the~match~succeeds;~else~it~fails.

~~~~Recently,~this~feature~was~upgraded~to~allow~ranges.~~Thus,

~~600~>~~\{(16,0,84,0),(57,10,29,4),(0,80,0,20),(95,0,0,5),

~~~~~~~~~(0,100,0,0),(18,60,20,2),(0,0,100,0),(0,50,50,0)\}~>~450

~~~~will~work,~as~well.

Allowing~Alternatives:

~~~~Very~occasionally,~you~may~wish~to~allow~alternative~pattern~units

~~~~(i.e.,~\char`\"{}match~either~A~or~B\char`\"{}).~~You~can~do~this~using~something

~~~~like

~~~~~~~~~~~~~~~~(~GAGA~|~GCGCA)

~~~~which~says~\char`\"{}match~either~GAGA~or~GCGCA\char`\"{}.~~You~may~take

~~~~alternatives~of~a~list~of~pattern~units,~for~example

~~~~~~~~(p1=3...3~3...8~\textasciitilde{}p1~|~p1=5...5~4...4~\textasciitilde{}p1~GGG)

~~~~would~match~one~of~two~sequences~of~pattern~units.~~There~is~one

~~~~clumsy~aspect~of~the~syntax:~to~match~a~list~of~alternatives,~you

~~~~need~to~fully~the~request.~~Thus,

~~~~~~~~(GAGA~|~(GCGCA~|~TTCGA))

~~~~would~be~needed~to~try~the~three~alternatives.

One~Minor~Extension

~~~~Sometimes~a~pattern~will~contain~a~sequence~of~distinct~ranges,

~~~~and~you~might~wish~to~limit~the~sum~of~the~lengths~of~the~matched

~~~~subsequences.~~~For~example,~suppose~that~you~basically~wanted~to

~~~~match~something~like

~~~~ARRYYTT~p1=0...5~GCA{[}1,0,0]~p2=1...6~\textasciitilde{}p1~4...8~\textasciitilde{}p2~p3=4...10~CCT

~~~~but~that~the~sum~of~the~lengths~of~p1,~p2,~and~p3~must~not~exceed

~~~~eight~characters.~~To~do~this,~you~could~add~

~~~~~~~~length(p1+p2+p3)~<~9

~~~~as~the~last~pattern~unit.~~It~will~just~succeed~or~fail~(but~does

~~~~not~actually~match~any~characters~in~the~sequence).

~~~~

Matching~Protein~Sequences

~~~~Suppose~that~the~input~file~contains~protein~sequences.~~In~this

~~~~case,~you~must~invoke~scan\_for\_matches~with~the~\char`\"{}-p\char`\"{}~option.~~You

~~~~cannot~use~aspects~of~the~language~that~relate~directly~to

~~~~nucleotide~sequences~(e.g.,~the~-c~command~line~option~or~pattern

~~~~constructs~referring~to~the~reverse~complement~of~a~previously

~~~~matched~unit).~~

~~~~You~also~have~two~additional~constructs~that~allow~you~to~match

~~~~either~\char`\"{}one~of~a~set~of~amino~acids\char`\"{}~or~\char`\"{}any~amino~acid~other~than

~~~~those~a~given~set\char`\"{}.~~For~example,

~~~~~~~~p1=0...4~any(HQD)~1...3~notany(HK)~p1

~~~~would~successfully~match~a~string~like

~~~~~~~~~~~YWV~D~AA~C~YWV

Using~the~show\_hits~Utility

~~~~When~viewing~a~large~set~of~complex~matches,~you~might~find~it

~~~~convenient~to~post-process~the~scan\_for\_matches~output~to~get~a

~~~~more~readable~version.~~We~provide~a~simple~post-processor~called

~~~~\char`\"{}show\_hits\char`\"{}.~~To~see~its~effect,~just~pipe~the~output~of~a

~~~~scan\_for\_matches~into~show\_hits:

~~~~~Normal~Output:

~~~~~~~~clone\%~scan\_for\_matches~-c~pat\_file~<~tmp

~~~~~~~~>tst1:{[}1,28]

~~~~~~~~gtacguaacc~~ggttaac~cgguuacgtac~

~~~~~~~~>tst1:{[}28,1]

~~~~~~~~gtacgtaacc~~ggttaac~cggttacgtac~

~~~~~~~~>tst2:{[}2,31]

~~~~~~~~CGTACGUAAC~C~GGTTAACC~GGUUACGTACG~

~~~~~~~~>tst2:{[}31,2]

~~~~~~~~CGTACGTAAC~C~GGTTAACC~GGTTACGTACG~

~~~~~~~~>tst3:{[}3,32]

~~~~~~~~gtacguaacc~g~gttaactt~cgguuacgtac~

~~~~~~~~>tst3:{[}32,3]

~~~~~~~~gtacgtaacc~g~aagttaac~cggttacgtac~

~~~~~Piped~Through~show\_hits:

~~~~

~~~~~~~~clone\%~scan\_for\_matches~-c~pat\_file~<~tmp~|~show\_hits

~~~~~~~~tst1:{[}1,28]:~~gtacguaacc~~~ggttaac~~cgguuacgtac

~~~~~~~~tst1:{[}28,1]:~~gtacgtaacc~~~ggttaac~~cggttacgtac

~~~~~~~~tst2:{[}2,31]:~~CGTACGUAAC~C~GGTTAACC~GGUUACGTACG

~~~~~~~~tst2:{[}31,2]:~~CGTACGTAAC~C~GGTTAACC~GGTTACGTACG

~~~~~~~~tst3:{[}3,32]:~~gtacguaacc~g~gttaactt~cgguuacgtac

~~~~~~~~tst3:{[}32,3]:~~gtacgtaacc~g~aagttaac~cggttacgtac

~~~~~~~~clone\%~

~~~~Optionally,~you~can~specify~which~of~the~\char`\"{}fields\char`\"{}~in~the~matches

~~~~you~wish~to~sort~on,~and~show\_hits~will~sort~them.~~The~field

~~~~numbers~start~with~0.~~So,~you~might~get~something~like

~~~~~~~~clone\%~scan\_for\_matches~-c~pat\_file~<~tmp~|~show\_hits~2~1

~~~~~~~~tst2:{[}2,31]:~~CGTACGUAAC~C~GGTTAACC~GGUUACGTACG

~~~~~~~~tst2:{[}31,2]:~~CGTACGTAAC~C~GGTTAACC~GGTTACGTACG

~~~~~~~~tst3:{[}32,3]:~~gtacgtaacc~g~aagttaac~cggttacgtac

~~~~~~~~tst1:{[}1,28]:~~gtacguaacc~~~ggttaac~~cgguuacgtac

~~~~~~~~tst1:{[}28,1]:~~gtacgtaacc~~~ggttaac~~cggttacgtac

~~~~~~~~tst3:{[}3,32]:~~gtacguaacc~g~gttaactt~cgguuacgtac

~~~~~~~~clone\%~

~~~~In~this~case,~the~hits~have~been~sorted~on~fields~2~and~1~(that~is,

~~~~the~third~and~second~matched~subfields).

~~~~show\_hits~is~just~one~possible~little~post-processor,~and~you

~~~~might~well~wish~to~write~a~customized~one~for~yourself.

Reducing~the~Cost~of~a~Search

~~~~The~scan\_for\_matches~utility~uses~a~fairly~simple~search,~and~may

~~~~consume~large~amounts~of~CPU~time~for~complex~patterns.~~Someday,

~~~~I~may~decide~to~optimize~the~code.~~However,~until~then,~let~me

~~~~mention~one~useful~technique.~~

~~~~When~you~have~a~complex~pattern~that~includes~a~number~of~varying

~~~~ranges,~imprecise~matches,~and~so~forth,~it~is~useful~to

~~~~\char`\"{}pipeline\char`\"{}~matches.~~That~is,~form~a~simpler~pattern~that~can~be

~~~~used~to~scan~through~a~large~database~extracting~sections~that

~~~~might~be~matched~by~the~more~complex~pattern.~~Let~me~illustrate

~~~~with~a~short~example.~~Suppose~that~you~really~wished~to~match~the

~~~~pattern~

~~~~p1=3...5~0...8~\textasciitilde{}p1{[}1,1,0]~p2=6...7~3...6~AGC~3...5~RYGC~\textasciitilde{}p2{[}1,0,0]

~~~~In~this~case,~the~pattern~units~AGC~3...5~RYGC~can~be~used~to~rapidly

~~~~constrain~the~overall~search.~~You~can~preprocess~the~overall

~~~~database~using~the~pattern:

~~~~~~~~~~31...31~AGC~3...5~RYGC~7...7

~~~~Put~the~complex~pattern~in~pat\_file1~and~the~simpler~pattern~in

~~~~pat\_file2.~~Then~use,

~~~~~~~~scan\_for\_matches~-c~pat\_file2~<~nucleotide\_database~|

~~~~~~~~scan\_for\_matches~pat\_file1

~~~~The~output~will~show~things~like

~~~~>seqid:{[}232,280]{[}2,47]

~~~~matches~pieces

~~~~Then,~the~actual~section~of~the~sequence~that~was~matched~can~be

~~~~easily~computed~as~{[}233,278]~(remember,~the~positions~start~from

~~~~1,~not~0).

~~~~Let~me~finally~add,~you~should~do~a~few~short~experiments~to~see

~~~~whether~or~not~such~pipelining~actually~improves~performance~-{}-~it

~~~~is~not~always~obvious~where~the~time~is~going,~and~I~have

~~~~sometimes~found~that~the~added~complexity~of~pipelining~actually

~~~~slowed~things~up.~~It~gets~its~best~improvements~when~there~are

~~~~exact~matches~of~more~than~just~a~few~characters~that~can~be

~~~~rapidly~used~to~eliminate~large~sections~of~the~database.

=============

Additions:

Feb~9,~1995:~~~the~pattern~units~\textasciicircum{}~and~\$~now~work~as~in~normal~regular

~~~~~~~~~~~~~~~expressions.~~That~is

~~~~~~~~~~~~~~~~~~~~~~~~TTF~\$

~~~~~~~~~~~~~~~matches~only~TTF~at~the~end~of~the~string~and~

~~~~~~~~~~~~~~~~~~~~~~~~\textasciicircum{}~TTF~

~~~~~~~~~~~~~~~matches~only~an~initial~TTF

~~~~~~~~~~~~~~~The~pattern~unit~

~~~~~~~~~~~~~~~~~~~~~~~~<p1

~~~~~~~~~~~~~~~matches~the~reverse~of~the~string~named~p1.~~That~is,

~~~~~~~~~~~~~~~if~p1~matched~GCAT,~then~<p1~would~match~TACG.~~Thus,

~~~~~~~~~~~~~~~~~~~p1=6...6~<p1

~~~~~~~~~~~~~~~matches~a~real~palindrome~(not~the~biologically~common

~~~~~~~~~~~~~~~meaning~of~\char`\"{}reverse~complement\char`\"{})


\end{lyxcode}

\end{document}