formats.tex

hmmer源程序

% --------------------------------------------------------------
% squid:formats.tex
% SRE, Wed Jul 14 17:54:59 1999
% $CVS Id$
% --------------------------------------------------------------

\chapter {Sequence file formats}

\section{Summary}

The software can handle a number of different file formats. By
default, it autodetects the file format, so you don't have to worry
about converting formats. Most common file formats are recognized,
including FASTA, Genbank, EMBL, Swissprot, PIR, and FASTA for
unaligned sequences, and GCG MSF, Clustal, Phylip, and Stockholm
format for multiple sequence alignments. Some parts of the source code
call the autodetector the ``Babelfish''. 

The Babelfish has three drawbacks. First, it takes a small amount of
time to do the autodetection. Second, the Babelfish is aggressive, and
it makes mistakes when a file isn't one of the known formats -- in
particular, it can recognize plain text files as SELEX alignments,
because the SELEX format is so free-form. Third, because the Babelfish
works by reading the first part of the file then rewinding it before
starting to process it, you can't use the Babelfish on a nonrewindable
stream: e.g. when you're taking sequence input from a UNIX pipe
instead of a file, or when the file is gzipped and has to be
decompressed before processing.  In normal use, when you're using the
software interactively from the command line on sequence files that
you're familiar with, the Babelfish is very convenient and
(relatively) safe.

However, you'll find that there are times when you want to override
the Babelfish -- particularly in high-throughput analysis, when you
know the format your files are supposed to be in, and you'd rather
increase robustness and sacrifice interactive flexibility. All the
programs have an \verb+--informat <format>+ option that lets you
specify the format and shut off the Babelfish. You \emph{must} use
\verb+--informat+ to use compressed files, or to read sequence from a
UNIX pipe... see below for more details on these tricks.

\section{Formats recognized by the Babelfish}

Recognized unaligned sequence file formats:

\begin{tabular}{ll}\hline
Format name   & Note \\ \hline
fasta         & BLAST flatfile databases, etc.\\
genbank       & NCBI Genbank flat file format.\\
embl          & Includes both EMBL (DNA) and SWISSPROT (protein) databases.\\
pir           & Protein Information Resource database (NBRF/Georgetown)\\
gcg           & Wisconsin Genetics Computer Group; only allows one sequence per file.\\
gcgdata       & I think this GCG database format is obsolete now.\\ \hline
\end{tabular}

Recognized multiple sequence alignment file formats:

\begin{tabular}{ll}\hline
Format name   & Note \\ \hline
stockholm     & Pfam format. Allows databases of more than one alignment per file\\
selex         & Old NeXagen RNA alignment format, adopted by early HMMER releases.\\
msf           & GCG's alignment format.\\
clustal       & ClustalV, ClustalW, and friends.\\
a2m           & Aligned FASTA format; see comment below.\\
phylip        & Format used by Felsenstein's PHYLIP phylogenetic inference software\\\hline
\end{tabular}

Aligned FASTA format (here called ``A2M'', though I believe that what
Haussler's group at UCSC started calling A2M is yet another variant of
aligned FASTA that's incompatible with this A2M) is only autodetected
when an alignment file is expected. Otherwise an A2M file will be
recognized as unaligned FASTA, and its gap characters (if any) will be
parsed as sequence characters -- often not what you want.

Alignment files may be used when unaligned files are expected -- the
sequences will silently be de-aligned and read sequentially. The
converse is not true; you can't give an unaligned sequence format when
an alignment is expected (makes sense, right?).

There is no provision for enforcing that single unaligned sequence
formats really do contain just a single sequence.  An attempt to
convert a multisequence file to GCG format will silently ``succeed'',
and the file may look ok to your eye, but that multisequence ``GCG''
file is illegal. The data will be corrupted if you try to read that
file back in, possibly without generating any error messages.

It turns out that other formats work too, but they're undocumented,
not subjected to any quality control testing at software release time,
and prone to change without notice at my slightest whim. (In other
words, even less supported than the software already is.) The brave,
curious, or desperate are invited to peruse
\prog{seqio.c} and \prog{squid.h}.

\section{Special tricks}

\subsection{Reading from standard input (probably UNIX-only)}

If you give ``-'' as a sequence filename, the software will read the
sequences from standard input rather than from a file. You will need
to specify the format of the incoming data using the
\verb+--informat+ option.
Any format except SELEX can be read from standard input. This lets you
use any program downstream in a standard UNIX pipe.

There is one limitation: you can't use ``-'' more than once on a
command line, for obvious reasons. (How is it supposed to read more
than one file from one standard input stream?) If you do, behavior of
the software is undefined -- in other words, the software don't check
for whether you're making this mistake, so God help you if you do.

\subsection{Reading from gzip'ed files (probably UNIX-only)}

A sequence file in any format except SELEX can be compressed by gzip,
and read in its compressed form. The software looks for the suffix
\prog{.gz} to detect gzip'ed files. This allows you to save disk space
by keeping sequence files gzip'ed, if you like. gzip is not built in;
the software needs to find a gzip executable in your current PATH.

If for some reason you name a file with a \prog{.gz} suffix and it's
\emph{not} a gzip-compressed file, the software will still try to
decompress it, and peculiar things may happen.

\section{FASTA format, the recommended unaligned format}

FASTA is probably the simplest of formats for unaligned sequences.
FASTA files are easily created in a text editor.  Each sequence is
preceded by a line starting with \verb+>+. The first word on this line
is the name of the sequence. The rest of the line is a description of
the sequence (free format). The remaining lines contain the sequence
itself. You can put as many letters on a sequence line as you want.

\textbf{Example of a simple FASTA file:}
\begin{verbatim}
>seq1 This is the description of my first sequence.
AGTACGTAGTAGCTGCTGCTACGTGCGCTAGCTAGTACGTCA CGACGTAGATGCTAGCTGACTCGATGC
>seq2 This is a description of my second sequence.
CGATCGATCGTACGTCGACTGATCGTAGCTACGTCGTACGTAG CATCGTCAGTTACTGCATGCTCG
\end{verbatim}

For better or worse, FASTA is not a documented standard. Minor (and
major) variants are in widespread use in the bioinformatics community,
all of which are called ``FASTA format''. My software attempts to
cater to all of them, and is tolerant of common deviations in FASTA
format. Certainly anything that is accepted by the database formatting
programs in NCBI BLAST or WU-BLAST (e.g. setdb, pressdb, xdformat)
will also be accepted by my software. Blank lines in a FASTA file are
ignored, and so are spaces or other gap symbols (dashes, underscores,
periods) in a sequence. Other non-amino or non-nucleic acid symbols in
the sequence are also silently ignored, mostly because some people
seem to think that ``*'' or ``.'' should be added to protein sequences
to (redundantly) indicate the end of the sequence. The parser will
also accept unlimited line lengths, which allows it to accomodate the
enormous description lines in the NCBI NR databases.

On the other hand, any FASTA files \emph{generated} by my software
adhere closely to community standards, and should be usable by other
software packages (BLAST, FASTA, etc.) that are more picky about
parsing their input files. That means you can run a sloppy FASTA file
thru \prog{sreformat} to clean it up.

Partly because of this tolerance, the software may have a difficult
time dealing with files that are \textit{not} in FASTA format,
especially if you're relying on the Babelfish to do format
autodetection.  Some (now mercifully uncommon) file formats are so
similar to FASTA format that they be erroneously called FASTA by the
Babelfish and then quietly but lethally misparsed. An example is the
old NBRF file format. If you're using \verb+--informat+, things will
be more robust, and the software should simply refuse to accept a
non-FASTA file -- but you shouldn't count on this, because files
perversely similar to FASTA will still confuse the parser.  (The gist
of these caveats applies to all formats, not just FASTA.)

\section{SELEX, the quick and dirty alignment format}

An example of a simple SELEX alignment file:

\begin{verbatim}
# Example selex file

seq1     ACGACGACGACG.
seq2     ..GGGAAAGG.GA
seq3     UUU..AAAUUU.A

seq1  ..ACG
seq2  AAGGG
seq3  AA...UUU
\end{verbatim}

SELEX is an interleaved multiple alignment format that arose as a
simple, intuitive format that was easy to write and manipulate
manually in a text editor. It is usually easy to convert other
alignment formats into SELEX format, even with a couple of lines of
Perl, but it can be harder to go the other way, since SELEX is more
free-format than other alignment formats. For instance, GCG's MSF
format and the output of the CLUSTALV multiple alignment program are
similar interleaved formats that can be converted to SELEX just by
stripping a small number of non-sequence lines out. Because SELEX
evolved to accomodate different user input styles, it is very tolerant
of various inconsistencies such as different gap symbols, varying line
lengths, etc.

Each line contains a name, followed by the aligned sequence. A space,
dash, underscore, or period denotes a gap. If the alignment is too
long to fit on one line, the alignment is split into multiple blocks,
separated by blank lines. The number of sequences, their order, and
their names must be the same in every block (even if a sequence has no
residues in a given block!) Other blank lines are ignored. You can add
comments to the file on lines starting with a \verb+#+.

SELEX stands for ``Systematic Evolution of Ligands by Exponential
Enrichment'' -- it refers to the Tuerk and Gold technology for
evolving families of small RNAs for particular functions
\cite{Tuerk90b}. SELEX files were what we used to keep track of
alignments of these small RNA families, at a company then called
NeXagen, in Boulder. It's an interesting piece of historical baggage.
With the development of HMMER and more need for annotated alignments
in Pfam, SELEX format later evolved into ``extended SELEX'', with a
reserved comment style that allowed structural markup and other
annotations, but that became unwieldy. We now use Stockholm format
(see below) for highly annotated alignments. (Extended SELEX is
deprecated and undocumented.) Still, the basic SELEX format remains a
useful ``lowest common denominator'' alignment format, and has been
retained.

\subsubsection {Detailed specification of a SELEX file}

\begin{enumerate}
\item
Any line beginning with a \verb+#=+ as the first two characters is a
parsed machine comment in extended SELEX, and is now deprecated. 

\item
All other lines beginning with a \verb+%+ or \verb+#+ as the first
character are user comments.  User comments are ignored by all
software. Anything may appear on these lines. Any number of comments
may be included in a SELEX file, and at any point.

\item
Lines of data consist of a name followed by a sequence. The total
length of the line must be smaller than 4096 characters.

\item
Names must be a single word. Any non-whitespace characters are
accepted.  No spaces are tolerated in names: names MUST be a
single word. Names must be less than 32 characters long.

\item In the sequence, any of the characters \verb+-_.+ or a space are
recognized as gaps. Any other characters are interpreted as sequence.
Sequence is case-sensitive. There is a common assumption by my
software that upper-case symbols are used for consensus (match)
positions and lower-case symbols are used for inserts. This language
of ``match'' versus ``insert'' comes from the hidden Markov model
formalism \cite{Krogh94}. To almost all of my software, this isn't
important, and it immediately converts the sequence to all upper-case
after it's read.

\item
Multiple different sequences are grouped in a block of data lines.
Blocks are separated by blank lines. No blank lines are tolerated