% EPN, Mon Oct 21 12:57:38 2013
% EPN, Mon Jun 27 12:32:24 2016 [1.1.2 release]
% EPN, Sat Nov 16 19:29:20 2019 [1.1.3 release]
% Actual commands run on: 
% cbbdev13
% $ uname -a
% Linux cbbdev13 2.6.32-573.18.1.el6.x86_64 #1 SMP Tue Feb 9 22:46:17 UTC 2016 x86_64 x86_64 x86_64 GNU/Linux

\section{Tutorial}
\label{section:tutorial}
\setcounter{footnote}{0}

Here's a tutorial walk-through of some small projects with
Infernal. This should suffice to get you started on work of your own,
and you can (at least temporarily) skip the rest of the Guide,
such as all the nitty-gritty details of available command line
options.

\subsection {The programs in Infernal}


\begin{tabular}{ll}
\multicolumn{2}{c}{\textbf{Core programs}}\\
 & \\ 
\textbf{cmbuild}     & Build a covariance model from an input multiple alignment.\\
\textbf{cmcalibrate} & Calibrate E-value parameters for a covariance model.\\
\textbf{cmsearch}    & Search a covariance model against a sequence database.\\
\textbf{cmscan}      & Search a sequence against a covariance model database.\\
\textbf{cmalign}     & Make a multiple alignment of many sequences to a common covariance model.\\
 & \\ 
\multicolumn{2}{c}{\textbf{Other utilities}}\\ 
 & \\ 
\textbf{cmconvert} & Convert CM formats to/from Infernal v1.1 format.\\ 
\textbf{cmemit}    & Generate (sample) sequences from a covariance model.\\
\textbf{cmfetch}   & Get a covariance model by name or accession from a CM database.\\
\textbf{cmpress}   & Format a CM database into a binary format for \prog{cmscan}.\\
\textbf{cmstat}    & Show summary statistics for each model in a CM database.\\ 
\end{tabular} \\
\\

In this section, we'll show examples of running each of these
programs, using examples in the \otext{tutorial/} subdirectory of the
distribution.

\subsection{Files used in the tutorial}

The subdirectory \otext{/tutorial} in the Infernal distribution contains the
files used in the tutorial, as well as a number of examples of various
file formats that Infernal reads. The important files for the tutorial
are:

\begin{sreitems}{\emprog{minifam.i1\{m,i,f,p\}}}
\item[\otext{tRNA5.sto}] A multiple alignment of five tRNA
  sequences. This file is a simple example of \emph{Stockholm
  format} that Infernal uses for structurally-annotated alignments.
%
\item[\otext{tRNA5.c.cm}] An example CM file. Built with
  \prog{cmbuild} from \otext{tRNA5.sto} and calibrated using
  \prog{cmcalibrate}. Included so you don't need to calibrate your own
  model file, which takes about 20 minutes. 
%
\item[\otext{mrum-genome.fa}] The 3 Mb genome of the methanogenic archeaon 
  \emph{Methanobrevibacter ruminantium}, in
  FASTA format, downloaded from the NCBI Nucleotide database
  (accession: NC\_13790.1). 
%
\item[\otext{tRNA5-mrum.out}] An example \prog{cmsearch} output file,
  obtained by searching \otext{tRNA5.c.cm} against \otext{mrum-genome.fa}.
%
\item[\otext{5S\_rRNA.sto}] The Rfam 10.1 5S ribosomal RNA (RF00001) 
  ``seed'' alignment. 
%
\item[\otext{5S\_rRNA.c.cm}] A CM file built from
  \otext{5S\_rRNA.sto} using \prog{cmbuild}
  and calibrated using \prog{cmcalibrate}.
%
\item[\otext{Cobalamin.sto}] The Rfam 10.1 Cobalamin riboswitch (RF00174) 
  ``seed'' alignment. 
%
\item[\otext{Cobalamin.c.cm}] A CM file built from
  \otext{Cobalamin.sto} using \prog{cmbuild}
  and calibrated using \prog{cmcalibrate}.
%
\item[\otext{minifam.cm}] A CM file including three calibrated CMs.
  This is actually just a concatenation 
  of the files \otext{tRNA5.c.cm}, \otext{5S\_rRNA.c.cm} and
  \otext{Cobalamin.c.cm}.
%
\item[\otext{minifam.i1\{m,i,f,p\}}] Binary compressed files
  corresponding to \otext{minifam.cm}, produced by \prog{cmpress}.
%
\item[\otext{metag-example.fa}] A FASTA sequence file containing 3
  sequences hand-selected from a metagenomics dataset \citep{Tringe05},
  used for demonstrating \prog{cmscan}. 
%
\item[\otext{minifam-metag.out}] An example \prog{cmscan} output file,
  obtained by searching \otext{minifam.cm} against \otext{metag-example.fa}.
%
\item[\otext{mrum-tRNAs10.fa}] A FASTA sequence file containing 10
  tRNAs predicted using \prog{cmsearch} in the \emph{M. ruminantium} genome.
%
\item[\otext{mrum-tRNAs10.out}] An example \prog{cmalign} output
  alignment, obtained by aligning the sequences in
  \otext{mrum-tRNAs10.fa} to the model in \otext{tRNA5.c.cm} with \prog{cmalign}.
%
\item[\otext{Cobalamin.fa}] A FASTA sequence file containing 1
  \prog{cmscan}-predicted Cobalamin riboswitch extracted from \otext{metag-example.fa}.
%
\item[\otext{tRNA5-noss.sto}] A Stockholm alignment file identical
  to \otext{tRNA5.sto} except without secondary structure annotation.
  Used to demonstrate HMM searches for models without secondary
  structure.
%
\item[\otext{tRNA5-hand.sto}] A Stockholm alignment file identical
  to \otext{tRNA5.sto} except it includes column reference annotation.
  Used to demonstrate expert annotation of model positions with
  \otext{cmbuild --hand}.
%
\item[\otext{tRNA5-hand.frag.sto}] A Stockholm alignment file similar
  to \otext{tRNA5-hand.sto} but with some sequence data at the
  beginning and end of sequences removed and replaced with gaps.  Used
  to demonstrate model building with fragmentary sequence information.
%
\item[\otext{tRNA5-hand.frag2.sto}] Same as
  \otext{tRNA5-hand.frag2.sto} except with the partial sequence
  annotated as fragments. Used to demonstrate model building with
  fragmentary sequence information.
%
\item[\otext{tRNA5-hand.fraggiven.sto}] Same as
  \otext{tRNA5-hand.frag2.sto} except with the partial sequence
  annotated as fragments. Used to demonstrate model building with
  fragmentary sequence information.
%
\item[\otext{tRNA5-hand.c.cm}] A CM file built from
  \otext{tRNA5-hand.sto} with \prog{cmbuild} and calibrated with
  \prog{cmcalibrate}. 
\end{sreitems}

\subsection{Searching a sequence database with a single covariance model}

\subsubsection{Step 1: build a covariance model with cmbuild}

Infernal starts with a multiple sequence alignment file that you
provide. It must be in Stockholm format and must include consensus
secondary structure annotation. The file \otext{tutorial/tRNA5.sto} is
an example of a simple Stockholm file. It is shown below, with a
secondary structure of the first sequence shown to the right for
reference (yeast Phe tRNA, labeled as ``tRNA1'' in the file):

\vspace{1em}
\begin{minipage}{4.0in}
\begin{sreoutput}[xleftmargin=0em]
# STOCKHOLM 1.0

tRNA1             GCGGAUUUAGCUCAGUUGGG.AGAGCGCCAGACUGAAGAUCUGGAGGUCC
tRNA2             UCCGAUAUAGUGUAAC.GGCUAUCACAUCACGCUUUCACCGUGGAGA.CC
tRNA3             UCCGUGAUAGUUUAAU.GGUCAGAAUGGGCGCUUGUCGCGUGCCAGA.UC
tRNA4             GCUCGUAUGGCGCAGU.GGU.AGCGCAGCAGAUUGCAAAUCUGUUGGUCC
tRNA5             GGGCACAUGGCGCAGUUGGU.AGCGCGCUUCCCUUGCAAGGAAGAGGUCA
#=GC SS_cons      <<<<<<<..<<<<.........>>>>.<<<<<.......>>>>>.....<

tRNA1             UGUGUUCGAUCCACAGAAUUCGCA
tRNA2             GGGGUUCGACUCCCCGUAUCGGAG
tRNA3             GGGGUUCAAUUCCCCGUCGCGGAG
tRNA4             UUAGUUCGAUCCUGAGUGCGAGCU
tRNA5             UCGGUUCGAUUCCGGUUGCGUCCA
#=GC SS_cons      <<<<.......>>>>>>>>>>>>.
//
\end{sreoutput}
\end{minipage}
\begin{minipage}{1.5in}
\includegraphics[scale=0.4]{Figures/trna1-DF6280}
\end{minipage}
\vspace{1em}

This is a simple example of a multiple RNA sequence alignment with
secondary structure annotation, in \emph{Stockholm format}.  Stockholm
format, the native alignment format used by HMMER and Infernal and the
Pfam and Rfam databases, is documented in detail later in the guide in
section~\ref{section:formats}.

For now, what you need to know about the key features of the input file is:
\begin{itemize}
\item The alignment is in an interleaved format.
Lines consist of a name, followed by an aligned sequence;
long alignments are split into blocks separated by blank lines.
\item Each sequence must have a unique name that has zero spaces in it. (This is important!)
\item For residues, any one-letter IUPAC nucleotide code is accepted,
      including ambiguous nucleotides. Case is ignored; residues
      may be either upper or lower case.
\item Gaps are indicated by the characters \otext{.}, \otext{\_}, -, or \verb+~+.
      (Blank space is not allowed.)
\item A special line starting with \otext{\#=GC SS\_cons} indicates
      the secondary structure consensus. Gap characters annotate
      unpaired (single-stranded) columns. Base pairs are indicated
      by any of the following pairs: \otext{<>}, \otext{()}, \otext{[]},
      or \otext{\{\}}. No pseudoknots are allowed; the
      open/close-brackets notation is only unambiguous for strictly
      nested base-pairing interactions.
      For more on secondary structure annotation see the WUSS format
      description in section~\ref{section:formats}.
\item The alignment begins with the special tag line
      \otext{\# STOCKHOLM 1.0}, and ends with \otext{//}.
      Stockholm alignments
      can be concatenated to create an alignment database flatfile
      containing many alignments.
\end{itemize}

The \prog{cmbuild} command builds a covariance model from an alignment (or
CMs for each of many alignments in a Stockholm file), and saves the
CM(s) in a file. For example, type:

\user{cmbuild tRNA5.cm tutorial/tRNA5.sto}

and you'll see some output that looks like:

\begin{sreoutput}
# cmbuild :: covariance model construction from multiple sequence alignments
# INFERNAL 1.1.5 (Sep 2023)
# Copyright (C) 2023 Howard Hughes Medical Institute.
# Freely distributed under the BSD open source license.
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# CM file:                                            tRNA5.cm
# alignment file:                                     tutorial/tRNA5.sto
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#                                                                      rel entropy
#                                                                      -----------
# idx    name                     nseq eff_nseq   alen  clen  bps bifs    CM   HMM description
# ------ -------------------- -------- -------- ------ ----- ---- ---- ----- ----- -----------
       1 tRNA5                       5     3.73     74    72   21    2 0.783 0.489 
#
# CPU time: 0.09u 0.00s 00:00:00.09 Elapsed: 00:00:00.10
\end{sreoutput}

If your input file had contained more than one alignment, you'd get
one line of output for each model. The information on these lines is
almost self-explanatory. The \otext{tRNA5} alignment consisted of 5
%  VERIFY WHEN UPDATING                                           ^
sequences with 74 aligned columns. Infernal turned it into a model of
%              ^^
72 consensus positions, which means it defined 2 gap-containing
%^                                             ^
alignment columns to be insertions relative to consensus. The 5
%                                                             ^
sequences were only counted as an ``effective'' total sequence number
(\otext{eff\_nseq}) of 3.73. The model includes 21 basepairs and 2
%                      ^^^^                     ^^               ^
bifurcations. The model ended up with a relative entropy per position
(\otext{rel entropy, CM}; information content) of 0.783 bits. If the
%                                                 ^^^^^
secondary structure information of the model were ignored the relative
entropy per position (\otext{rel entropy, HMM}) would be 0.489 bits.
%                                                        ^^^^^         
This output format is rudimentary. Infernal knows quite a bit more
information about what it's done to build this CM, but it's not
displaying it. You don't need to know more to be able to use the
model, so we can press on here. Model construction is described in
more detail in section~\ref{section:cmbuild}.

The new CM was saved to \otext{tRNA5.cm}. You can look at it
(e.g. \otext{> more tRNA5.cm}) if you like, but it isn't really designed
to be human-interpretable. You can treat \otext{.cm} files as compiled
models of your RNA alignment. The Infernal ASCII save file format is
defined in Section~\ref{section:formats}.

\subsubsection{Step 2: calibrate the model with cmcalibrate}

The next step is to ``calibrate'' the model. This step must be
performed prior to using your model for database searches with
\prog{cmsearch} or \prog{cmscan}. In this step, statistical parameters
necessary for reporting E-values (expectation values) are estimated
and stored in the CM file. When \prog{cmsearch} or \prog{cmscan} is
later used for a database search and a hit with score $x$ is found,
the E-value of that hit is the number of hits expected to score $x$ or
more just by chance (given the size of the search you're performing).

\emph{Importantly, if you're not going to use a model for database
search, there is no need to calibrate it.} For example, if you are
only going to use a model to create structurally annotated multiple
alignments of a large family like small subunit ribosomal RNA, don't
waste time calibrating it. \prog{cmsearch} and \prog{cmscan} are the
only Infernal programs that use E-values, so if you're not going to
use them then don't calibrate your model.

Unfortunately, CM calibration takes a long time because fairly long
random sequences must be searched to determine the expected
distribution of hit scores against nonhomologous sequences, and none
of the search acceleration heuristics described in
section~\ref{section:pipeline} can be used because they rely on
primary sequence similarity which is absent in random sequence.

The amount of time required for calibration varies widely, but
depends mainly on the size of the RNA family being modeled.
So you can know what kind of a wait you're in for, the
\prog{cmcalibrate} has a \otext{--forecast} option which reports an
estimate of the running time. To get an estimate for the tRNA model, do:

\user{cmcalibrate --forecast tRNA5.cm}\\

You should see something like this:

\begin{sreoutput}
# cmcalibrate :: fit exponential tails for CM E-values
# INFERNAL 1.1.5 (Sep 2023)
# Copyright (C) 2023 Howard Hughes Medical Institute.
# Freely distributed under the BSD open source license.
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# CM file:                                     tRNA5.cm
# forecast mode (no calibration):              on
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#
# Forecasting running time for CM calibration(s) on 4 cpus:
#
#                          predicted
#                       running time
# model name            (hr:min:sec)
# --------------------  ------------
  tRNA5                     00:03:26
#
# CPU time: 0.09u 0.00s 00:00:00.09 Elapsed: 00:00:00.09
[ok]
\end{sreoutput}

The header comes first, telling you what program you ran, on what file
and with what options. This calibration will use 4 CPUs which is the
default number for Infernal's multithreaded programs, but you can
modify this to \otext{<n>} CPUs with the \otext{--nforecast <n>}
option if you want to use a different number of CPUs or are planning
on using MPI to parallelize the calibrate (see the Installation
section). Using 4 CPUs, \prog{cmcalibrate}
estimates the time required for calibration on the machine I'm using
at about three and a half minutes.

Feel free to perform the calibration yourself if you'd like (with the
command \otext{cmcalibrate tRNA5.cm}). However, we've included the file
\otext{tRNA5.c.cm}, an already calibrated version of \otext{tRNA5.cm},
for you to use if you don't want to wait. To use this calibrated
model, copy over the \otext{tRNA5.cm} file you just made with the
calibrated version:
 
\user{cp tutorial/tRNA5.c.cm tRNA5.cm}\\
 
\subsubsection{Step 3: search a sequence database with cmsearch}

You can now use your tRNA model to search for tRNA homologs in a
database. The file \otext{mrum-genome.fa} is the genome sequence of the
archaeon \emph{Methanobrevibacter rumanitium} (accession:
\otext{NC\_13790.1}). We'll use this file as our database. To perform
the search:

\user{cmsearch tRNA5.cm tutorial/mrum-genome.fa}\\

As before, the first section is the header, telling you what program
your ran, on what, and with what options:

\begin{sreoutput}
# cmsearch :: search CM(s) against a sequence database
# INFERNAL 1.1.5 (Sep 2023)
# Copyright (C) 2023 Howard Hughes Medical Institute.
# Freely distributed under the BSD open source license.
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# query CM file:                         tRNA5.cm
# target sequence database:              tutorial/mrum-genome.fa
# number of worker threads:              4
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
\end{sreoutput}

The second section is a list of ranked top hits (sorted by E-value,
most significant hit first):

\begin{sreoutput}
Query:       tRNA5  [CLEN=72]
Hit scores:
 rank     E-value  score  bias  sequence      start     end   mdl trunc   gc  description
 ----   --------- ------ -----  ----------- ------- -------   --- ----- ----  -----------
  (1) !   1.4e-18   71.4   0.0  NC_013790.1  362026  361955 -  cm    no 0.50  Methanobrevibacter rumi ...
  (2) !   3.3e-18   70.2   0.0  NC_013790.1 2585265 2585193 -  cm    no 0.60  Methanobrevibacter rumi ...
  (3) !   9.5e-18   68.7   0.0  NC_013790.1  762490  762562 +  cm    no 0.67  Methanobrevibacter rumi ...
  (4) !   9.5e-18   68.7   0.0  NC_013790.1 2041704 2041632 -  cm    no 0.67  Methanobrevibacter rumi ...
  (5) !   2.4e-17   67.5   0.0  NC_013790.1 2351254 2351181 -  cm    no 0.62  Methanobrevibacter rumi ...
  (6) !     3e-17   67.2   0.0  NC_013790.1  735136  735208 +  cm    no 0.59  Methanobrevibacter rumi ...
  (7) !   5.3e-17   66.4   0.0  NC_013790.1 2186013 2185941 -  cm    no 0.53  Methanobrevibacter rumi ...
  (8) !   1.7e-16   64.8   0.0  NC_013790.1 2350593 2350520 -  cm    no 0.66  Methanobrevibacter rumi ...
  (9) !   2.9e-16   64.1   0.0  NC_013790.1 2585187 2585114 -  cm    no 0.59  Methanobrevibacter rumi ...
 (10) !   9.3e-16   62.5   0.0  NC_013790.1  662185  662259 +  cm    no 0.61  Methanobrevibacter rumi ...
 (11) !   1.3e-15   62.0   0.0  NC_013790.1  360887  360815 -  cm    no 0.55  Methanobrevibacter rumi ...
 (12) !   1.7e-15   61.7   0.0  NC_013790.1 2350984 2350911 -  cm    no 0.53  Methanobrevibacter rumi ...
 (13) !   3.3e-15   60.7   0.0  NC_013790.1 2186090 2186019 -  cm    no 0.54  Methanobrevibacter rumi ...
 (14) !   4.3e-15   60.4   0.0  NC_013790.1 2680159 2680233 +  cm    no 0.67  Methanobrevibacter rumi ...
 (15) !   8.1e-15   59.5   0.0  NC_013790.1 2749839 2749768 -  cm    no 0.53  Methanobrevibacter rumi ...
 (16) !   8.1e-15   59.5   0.0  NC_013790.1 2749945 2749874 -  cm    no 0.53  Methanobrevibacter rumi ...
 (17) !     1e-14   59.2   0.0  NC_013790.1  361676  361604 -  cm    no 0.51  Methanobrevibacter rumi ...
 (18) !     1e-14   59.2   0.0  NC_013790.1 2585073 2584999 -  cm    no 0.60  Methanobrevibacter rumi ...
 (19) !   1.2e-14   59.0   0.0  NC_013790.1 2130422 2130349 -  cm    no 0.59  Methanobrevibacter rumi ...
 (20) !   1.3e-14   58.9   0.0  NC_013790.1  546056  545947 -  cm    no 0.61  Methanobrevibacter rumi ...
 (21) !   4.1e-14   57.3   0.0  NC_013790.1  361915  361844 -  cm    no 0.42  Methanobrevibacter rumi ...
 (22) !   5.2e-14   56.9   0.0  NC_013790.1   97724   97795 +  cm    no 0.49  Methanobrevibacter rumi ...
 (23) !   6.3e-14   56.7   0.0  NC_013790.1 2350717 2350646 -  cm    no 0.68  Methanobrevibacter rumi ...
 (24) !   8.3e-14   56.3   0.0  NC_013790.1 1873887 1873815 -  cm    no 0.64  Methanobrevibacter rumi ...
 (25) !   1.5e-13   55.5   0.0  NC_013790.1  360730  360659 -  cm    no 0.40  Methanobrevibacter rumi ...
 (26) !   3.6e-13   54.3   0.0  NC_013790.1 2680310 2680384 +  cm    no 0.52  Methanobrevibacter rumi ...
 (27) !   3.6e-13   54.3   0.0  NC_013790.1 2664806 2664732 -  cm    no 0.60  Methanobrevibacter rumi ...
 (28) !   3.8e-13   54.2   0.0  NC_013790.1  361061  360989 -  cm    no 0.41  Methanobrevibacter rumi ...
 (29) !   7.7e-13   53.3   0.0  NC_013790.1 2130335 2130262 -  cm    no 0.55  Methanobrevibacter rumi ...
 (30) !   7.7e-13   53.3   0.0  NC_013790.1 2151672 2151745 +  cm    no 0.65  Methanobrevibacter rumi ...
 (31) !   2.9e-12   51.4   0.0  NC_013790.1  319297  319370 +  cm    no 0.62  Methanobrevibacter rumi ...
 (32) !   3.9e-12   51.1   0.0  NC_013790.1  361753  361679 -  cm    no 0.55  Methanobrevibacter rumi ...
 (33) !     4e-12   51.0   0.0  NC_013790.1  360983  360912 -  cm    no 0.50  Methanobrevibacter rumi ...
 (34) !   6.1e-12   50.4   0.0  NC_013790.1  361456  361383 -  cm    no 0.50  Methanobrevibacter rumi ...
 (35) !   7.6e-12   50.1   0.0  NC_013790.1  362798  362727 -  cm    no 0.51  Methanobrevibacter rumi ...
 (36) !     9e-12   49.9   0.0  NC_013790.1  917722  917793 +  cm    no 0.61  Methanobrevibacter rumi ...
 (37) !     1e-11   49.7   0.0  NC_013790.1 2583869 2583798 -  cm    no 0.51  Methanobrevibacter rumi ...
 (38) !   1.4e-11   49.3   0.0  NC_013790.1  360811  360740 -  cm    no 0.42  Methanobrevibacter rumi ...
 (39) !   1.4e-11   49.3   0.0  NC_013790.1  362324  362252 -  cm    no 0.51  Methanobrevibacter rumi ...
 (40) !   4.3e-11   47.8   0.0  NC_013790.1 1160526 1160609 +  cm    no 0.60  Methanobrevibacter rumi ...
 (41) !     1e-10   46.6   0.0  NC_013790.1  362403  362331 -  cm    no 0.49  Methanobrevibacter rumi ...
 (42) !   1.1e-10   46.5   0.0  NC_013790.1 2327124 2327042 -  cm    no 0.63  Methanobrevibacter rumi ...
 (43) !   1.2e-10   46.4   0.0  NC_013790.1  995344  995263 -  cm    no 0.49  Methanobrevibacter rumi ...
 (44) !   2.3e-10   45.4   0.0  NC_013790.1  256772  256696 -  cm    no 0.57  Methanobrevibacter rumi ...
 (45) !   2.5e-10   45.3   0.0  NC_013790.1 2584830 2584758 -  cm    no 0.64  Methanobrevibacter rumi ...
 (46) !   6.3e-10   44.1   0.0  NC_013790.1 2351071 2350997 -  cm    no 0.59  Methanobrevibacter rumi ...
 (47) !   6.4e-10   44.1   0.0  NC_013790.1  362552  362482 -  cm    no 0.55  Methanobrevibacter rumi ...
 (48) !   5.1e-09   41.2   0.0  NC_013790.1 1064775 1064858 +  cm    no 0.63  Methanobrevibacter rumi ...
 (49) !   1.2e-08   40.0   0.0  NC_013790.1  361222  361150 -  cm    no 0.45  Methanobrevibacter rumi ...
 (50) !   1.3e-08   40.0   0.0  NC_013790.1  361369  361297 -  cm    no 0.60  Methanobrevibacter rumi ...
 (51) !   4.8e-08   38.1   0.0  NC_013790.1  361596  361513 -  cm    no 0.61  Methanobrevibacter rumi ...
 (52) !   3.1e-07   35.6   0.0  NC_013790.1 1913310 1913227 -  cm    no 0.64  Methanobrevibacter rumi ...
 (53) !   2.6e-06   32.7   0.0  NC_013790.1  363464  363381 -  cm    no 0.51  Methanobrevibacter rumi ...
 (54) !     3e-06   32.5   0.0  NC_013790.1 2584954 2584872 -  cm    no 0.58  Methanobrevibacter rumi ...
 ------ inclusion threshold ------
 (55) ?     0.026   20.1   0.0  NC_013790.1  363803  363716 -  cm    no 0.50  Methanobrevibacter rumi ...
 (56) ?       3.4   13.4   0.0  NC_013790.1  984373  984304 -  cm    no 0.53  Methanobrevibacter rumi ...
\end{sreoutput}

The first number is the rank of each hit\footnote{Ranks of hits are in
  parantheses to make it easy to jump to/from an entry in the hit list
  and the hit alignment section, described later.}. Next comes either a
\otext{!} or \otext{?} symbol and then the \emph{E-value} of the hit.
The E-value is the statistical significance of the hit: the number of
hits we'd expect to score this highly in a database of this size
(measured by the total number of nucleotides) if the database
contained only nonhomologous random sequences. The lower the E-value,
the more significant the hit. The \otext{!} or \otext{?} that
precedes the E-value indicates whether the hit does (\otext{!}) or
does not (\otext{?}) satisfy the inclusion threshold for the search.
Inclusion thresholds are used to determine what matches should be
considered to be ``true'', as opposed to reporting thresholds that
determine what matches will be reported (often including the top of
the noise, so you can see what interesting sequences might be getting
tickled by your search). By default, inclusion thresholds usually
require an E-value of 0.01 or less, and reporting E-value thresholds
are set to 10.0, but these can be changed (see the manual page for
\prog{cmsearch} toward the end of guide).

The E-value is based on the \emph{bit score}, which is in the next
column. This is the log-odds score for the hit. Some people like to
see a bit score instead of an E-value, because the bit score doesn't
depend on the size of the sequence database, only on the covariance
model and the target sequence.

The next number, the \emph{bias}, is a correction term for biased
sequence composition that has been applied to the sequence bit
score. Infernal uses an alternative null model we call \emph{null3},
described more in section~\ref{section:pipeline}, to determine the bias
bit score correction. The bias correction is often very small and is
only reported to one decimal place, after rounding. For all hits in
this example search the bias column reads 0.0 bits, indicating that
the correction is less than 0.05 bits. On very biased sequences this
correction can become significant and is helpful for lowering the
score of high-scoring false positives that achieve high scores solely
due to their biased composition. 

Next comes the target sequence name the hit is in, and the start and
end positions of the hit within the sequence. Hits can occur on either
the top (Watson) or bottom (Crick) strand of the target
sequence\footnote{You can search either only the top strand with the
\otext{--toponly} or bottom strand with the \otext{--bottomonly} options
to \prog{cmsearch} and \prog{cmscan}.}, so the start position may be
less than (if hit is on the top strand) or greater than (if hit is on
the bottom strand) the end position. After the end position, comes a single
\otext{+} or \otext{-} symbol, indicating whether the hit is on the
top (\otext{+}) or bottom (\otext{-}) strand (solely for convenience -
so you don't have to look at the start and end positions to determine
the strand the hit is on).

After the strand symbol comes the model field, which indicates whether
the hit was found using either the CM (\otext{cm}) or a profile HMM
built from the CM (\otext{hmm}). This field is necessary because for
models with zero basepairs, \prog{cmsearch} (and \prog{cmscan}) use a
profile HMM instead of a CM for final hit scoring. This is done for
reasons of speed and efficiency, because profile HMM algorithms are
more efficient than CM ones and a CM with zero basepairs is
essentially equivalent to a profile HMM. In this example, since our
tRNA model does include basepairs, \prog{cmsearch} used a CM to score
all hits and so all hits have \prog{cm} for this column. There's an
example later in the tutorial of hits found with a profile HMM.

The next column indicates whether the hit is \emph{truncated} or
not. Infernal uses special versions of its CM dynamic programming
algorithms to allow detection of structural RNAs that have been
truncated due to missing data at the beginning and/or end of a target
sequence. Truncated hits are most common in databases that include
single reads from shotgun sequencing projects. Since our database is a
complete genome, we don't expect any hits to be truncated due to
missing data. For all hits the ``trunc'' column reads
``no'' indicating that, as expected, none of the hits are
truncated. There are examples of truncated hits in the next exercise
which uses \prog{cmscan}. Section~\ref{section:pipeline} describes how
Infernal detects and aligns truncated hits in more detail.

The next column reports the GC fraction of the hit. This is the
fraction of residues in the target sequence hit that are either G or C
residues. The GC fraction is included as an additional indication of
the level of sequence bias of the hit. Some expert users may be aided
by this number when deciding if they believe a hit is a real homolog
or a false positive.

Finally comes the description of the sequence, if any. This
description is propogated from the input target sequence file.

After the hit list comes the hit alignments section. Each hit in 
the hit list will have a corresponding entry in this section, in the
same order. As an illustrative example, let's take a look at
hit number 43. First, take a look at the first four lines for this
hit: 

\begin{sreoutput}
>> NC_013790.1  Methanobrevibacter ruminantium M1 chromosome, complete genome
 rank     E-value  score  bias mdl mdl from   mdl to       seq from      seq to       acc trunc   gc
 ----   --------- ------ ----- --- -------- --------    ----------- -----------      ---- ----- ----
 (43) !   1.2e-10   46.4   0.0  cm        1       72 []      995344      995263 - .. 0.93    no 0.49
\end{sreoutput}

The first line of each hit alignment begins with \otext{>>} followed
by a single space, the name of the target sequence, then two spaces
and the description of the sequence, if any. Next comes a set of
tabular fields that is partially redundant with the information in the
hit list. The first five columns are the same as in the hit list. The
next column reports the type of model used for the alignment, as
described above for the hit list. The next four columns report the
boundaries of the alignment with respect to the query model (``mdl
from'' and ``mdl to'') and the target sequence (``seq from'' and ``seq
to''). Following the ``seq to'' column is a \otext{+} or \otext{-}
symbol indicating whether the hit is on the top (\otext{+}) or bottom
(\otext{-}) strand.

It's not immediately easy to tell from the ``to'' coordinate whether
the alignment ended internally in the query model or target sequence,
versus ran all the way (as in a full-length global alignment) to the
end(s). To make this more readily apparent, with each pair of query
and target endpoint coordinates, there's also a little symbology. For
the normal case of a non-truncated hit: \otext{..} means both ends of
the alignment ended internally, and \otext{[]} means both ends of the
alignment were full-length flush to the ends of the query or target,
and \otext{[.}  and \otext{.]} mean only the left (5') or right (3') end
was flush/full length. For truncated hits, the symbols are the same
except that either the first and/or the second symbol will be a
\verb+~+ for the query and target. If the first symbol is \verb+~+
then the left end of the alignment is truncated because the 5' end of
the hit is predicted to be missing (extend beyond the beginning of the
target sequence). Similarly, if the second symbol is \verb+~+ then the
right end of the alignment is truncated because the 3' end of the hit
is predicted to be missing (extend beyond the end of the target
sequence). These two symbols occur just after the ``mdl to'' column for
the query, and after the strand \otext{+} or \otext{-} symbol for the
target sequence.

The next column is labeled ``acc'' and is the average posterior
probability of the aligned target sequence residues; effectively, the
expected accuracy per residue of the alignment. (This field may be
absent and replaced by a column labeled ``cyksc'' in certain
situations, as explained more below.)

The final two columns indicate whether the hit is truncated or not and
the GC fraction of the hit. These are redundant with the columns of
the same name in the hit list, described above.

Next comes the alignment display. This is an ``optimal
posterior accuracy'' alignment~\citep{Holmes98}, which means it is the
alignment with the maximal summed posterior probability of all
aligned residues. Take a look at the alignment for hit number 43:

\begin{sreoutput}
                                 v         v                                                            NC
                     (((((((,,<<<<______.._>>>>,<<<<<_______>>>>>,,,........,,<<<<<_______>>>>>))))))): CS
        tRNA5      1 gCcggcAUAGcgcAgUGGu..AgcgCgccagccUgucAagcuggAGg........UCCgggGUUCGAUUCcccGUgccgGca 72    
                     :::G:CAUAGCG AG GGU  A CGCG:CAG:CU +++A:CUG: G+        UC:GGGGUUCGA UCCCC:UG:C:::A
  NC_013790.1 995344 AGAGACAUAGCGAAGCGGUcaAACGCGGCAGACUCAAGAUCUGUUGAuuaguucuUCAGGGGUUCGAAUCCCCUUGUCUCUA 995263
                     **********************************************9444444445************************** PP
\end{sreoutput}

The alignment contains six lines. Start by looking at the second line
which ends with \otext{CS}.  The line shows the predicted secondary
structure of the target sequence. The format is a little fancier and
more informative than the simple least-common-denominator format we
used in the input alignment file. It's designed to make it easier to
see the secondary structure by eye. The format is described in detail
later (see WUSS format in section~\ref{section:formats}); for now,
here's all you need to know. Basepairs in simple stem loops are
annotated with \otext{<>} characters. Basepairs enclosing
multifurcations (multiple stem loops) are annotated with \otext{()},
such as the tRNA acceptor stem in this example. In more complicated
structures, \otext{[]} and \otext{\{\}} annotations also show up, to
reflect deeper nestings of multifurcations. For single stranded
residues, \otext{\_} characters mark hairpin loops; \otext{-}
characters mark interior loops and bulges; \otext{,} characters mark
single-stranded residues in multifurcation loops; and \otext{:}
characters mark single stranded residues external to any secondary
structure. Insertions relative to this consensus are annotated by a
\otext{.} character.

The line above the \otext{CS} line ends with \otext{NC} and marks
negative scoring non-canonical basepairs in the alignment with a
\otext{v} character. All other positions of the alignment will be
blank\footnote{For anyone trying to parse this output, this means it
is possible for this line to be completely blank except for the
\otext{NC} trailer.} More specifically, the following ten types of
basepairs which are assigned a negative score by the model at their
alignment positions will be marked with a \otext{v}: \otext{A:A},
\otext{A:C}, \otext{A:G}, \otext{C:A}, \otext{C:C}, \otext{C:U},
\otext{G:A}, \otext{G:G}, \otext{U:U}, and \otext{U:C}. The \otext{NC}
annotation makes it easy to quickly identify suspicious basepairs in
a hit. Importantly, the \otext{NC} annotation will only be present in
CM hit alignments (``mdl'' column reads ``cm'') and will be absent in
HMM hit alignments (``mdl'' column reads ``hmm'') because basepairs
are not scored by an HMM.

The third line shows that consensus of the query model. The highest
scoring residue sequence is shown. Upper case residues are highly
conserved. Lower case residues are weakly conserved or unconserved.
Dots (\otext{.}) in this line indicate insertions in the target
sequence with respect to the model.

The fourth line shows where the alignment score is coming from. For a
consensus basepair, if the observed pair is the highest-scoring
possible pair according to the consensus, both residues are shown in
upper case; if a pair has a score of $\geq 0$, both residues are
annotated by : characters (indicating an acceptable compensatory
basepair); else, there is a space, indicating that a negative
contribution of this pair to the alignment score. Note that the 
\otext{NC} line will only mark a subset of these negative scoring
pairs with a \otext{v}, as discussed above.
For a single-stranded consensus residue, if the observed residue is
the highest scoring possibility, the residue is shown in upper case;
if the observed residue has a score of $\geq 0$, a \otext{+} character
is shown; else there is a space, indicating a negative contribution to
the alignment score. Importantly, for HMM hits (``mdl'' column reads
``hmm''), \emph{all} positions are considered single stranded, since
an HMM scores each half of a basepair independently.

The fifth line, beginning with \otext{NC\_013790.1} is
the target sequence. Dashes (\otext{-}) in this line indicate deletions
in the target sequence with respect to the model.

The bottom line represents the posterior probability (essentially the
expected accuracy) of each aligned residue. A 0 means 0-5\%, 1 means
5-15\%, and so on; 9 means 85-95\%, and a \otext{*} means 95-100\%
posterior probability. You can use these posterior probabilities to
decide which parts of the alignment are well-determined or not. You'll
often observe, for example, that expected alignment accuracy degrades
around locations of insertion and deletion, which you'd intuitively
expect.

Alignments for some searches may be formatted slightly differently
than this example. Longer alignments to longer models will be broken
up into blocks of six lines each - this alignment was short enough to
be entirely contained within a single block.  If your model was built
with the \otext{--hand} option in \prog{cmbuild}, then an additional
line will be included in each block, with RF annotation.  If the model
used for the alignment was an HMM (the ``mdl'' column reads ``hmm'')
then the \otext{NC} line will be absent from each alignment
block. We'll see example of all three of these cases later in the
tutorial. Finally, the \otext{PP} line will be absent if
\prog{cmsearch} has determined that calculating posterior
probabilities for the optimal posterior accuracy alignment would have
required more than a preset amount of memory (that amount is
autodetermined based on model size, but is capped at a maximum of
1024Mb for the largest models, unless the \emprog{--mxsize} option is
used to increase it). If the \otext{PP} line is absent, the tabular
``acc'' field described above will be absent and replaced with a
``cyksc'' field which reports the bit score of the highest scoring CYK
alignment, and the displayed alignment will be the alignment computed
by the CYK algorithm which requires less memory than the optimal
posterior accuracy alignment.

After hit number 43, there's 13 more hit alignments for hits number 44
through 56. 

Finally, at the bottom of the file, you'll see some summary
statistics. For example, at the bottom of the tRNA search output,
you'll find something like:

\begin{sreoutput}
Internal CM pipeline statistics summary:
----------------------------------------
Query model(s):                                                  1  (72 consensus positions)
Target sequences:                                                1  (5874406 residues searched)
Target sequences re-searched for truncated hits:                 1  (360 residues re-searched)
Windows   passing  local HMM SSV           filter:           10932  (0.2063); expected (0.35)
Windows   passing  local HMM Viterbi       filter:                  (off)
Windows   passing  local HMM Viterbi  bias filter:                  (off)
Windows   passing  local HMM Forward       filter:             136  (0.002704); expected (0.005)
Windows   passing  local HMM Forward  bias filter:             133  (0.002634); expected (0.005)
Windows   passing glocal HMM Forward       filter:              84  (0.001948); expected (0.005)
Windows   passing glocal HMM Forward  bias filter:              84  (0.001948); expected (0.005)
Envelopes passing glocal HMM envelope defn filter:              99  (0.001331); expected (0.005)
Envelopes passing  local CM  CYK           filter:              60  (0.0007629); expected (0.0001)
Total CM hits reported:                                         56  (0.0007205); includes 0 truncated hit(s)

# CPU time: 1.44u 0.03s 00:00:01.47 Elapsed: 00:00:02.91
//
[ok]
\end{sreoutput}

This gives you some idea of what's going on in Infernal's acceleration
pipeline. You've got one query CM, and the database has one target
%                    ^                                  ^
sequence. The search examined 5,874,406 residues, even though the
%                             ^^^^^^^^^
actual target sequence length is only half that, because both the top
and bottom strand of each target is searched. 360 of those residues
%                                             ^^^
were searched more than once in an effort to find truncated
hits. Ignore this for the moment, we'll revisit this later after
discussing the filter pipeline (in the subsection entitled ``Truncated
RNA detection'' below).

Each sequence goes through a multi-stage filter pipeline of four
scoring algorithms called SSV, Viterbi, Forward\footnote{Actually two
  separate Forward based filters are used, the first with the profile
  HMM in local mode and the next with the profile HMM in global
  mode. There's more detail on this in
  section~\ref{section:pipeline}.}, and CYK in order of increasing
sensitivity and increasing computational requirement. The filter
pipeline is the topic of section~\ref{section:pipeline} of this guide
but briefly, SSV, Viterbi and Forward are profile HMM algorithms which
are more efficient than CM algorithms. These three algorithms are the
same ones used by HMMER3 and are the main reason that version 1.1 of
Infernal is so much faster than previous versions. For these HMM
stages, Infernal uses a filter profile HMM that was constructed
simultaneously with the CM, from the same training alignment in
\prog{cmbuild}, and stored in the CM file. CYK is a CM scoring
algorithm, so it's slow, but it is accelerated using banded dynamic
programming with bands derived from an HMM alignment. Subsequences
that survive all filters are finally scored with the CM Inside
algorithm, agains using HMM bands. Subsequences that score
sufficiently high with Inside are then aligned using the optimal
posterior accuracy algorithm and displayed.

The score thresholds for a subsequence surviving each HMM filter stage
are dependent on the search space size (sequence database size for
\prog{cmsearch}). This differs from HMMER3 which always uses the same
filter thresholds. In general, the larger the search
space the more strict the thresholds are, because a hit must have a
higher bit score to have a significant E-value.  In this case, the
database is relatively small so the filter thresholds have been set
relatively loosely. The SSV filter has been configured to allow
subsequences with a P-value of $\leq 0.35$ through
%                                   ^^^^^^
the SSV score filter (thus, if the database contained no homologs and
P-values were accurately calculated, the highest scoring 35\% of the
%                                           ^^^^
residues will pass the filter). Here, about 21\% of the database in
%                                           ^^^^
10,932 separate windows got through the SSV filter. For a database of
%^^^^^
this size, the local Viterbi filter is turned off.  The local Forward filter
is set to allow an expected 0.5\% of the database survive. Here about
%                           ^^^^
0.3\% survives in 136 windows. Next, each surviving window is checked
%^^
to see if the target sequence is ``obviously'' so biased in its
composition that it's unlikely to be a true homolog. This is called
the ``bias filter''\footnote{There's also a bias filter step used in
  the local Viterbi filter stage, when it is used.} and applying a bit
score correction to previous filter's score for each window and
recomputing the P-value. Three of the 136 windows fail to pass
%                        ^^^^^        ^^^
the local Forward bias filter stage. Next, the Forward algorithm is
used to score each window again, but this time with the HMM configured
in glocal mode requiring a full length alignment to the
model\footnote{The use of glocal Forward is another important
  difference between Infernal and HMMER3's (v3.0) pipeline. HMMER v3.0
  only uses local HMM algorithms.}  As with the local stage, an
expected 0.5\% of the database is expected to survive. In this case,
84 of the 133 windows, comprising about 0.2\% of the database,
%^        ^^^
survive. The bias filter is run again, this time applying a correction
to the glocal Forward scores. For this search, 0 windows are removed at
this stage. The envelope definition stage is next. This stage is very
similar to the HMMER3 domain definition stage, with the difference
that the HMM is configured in glocal rather than local mode. In this
stage, the Forward and Backward algorithms are used to identify zero
or more hit envelopes in each window, where each envelope contains one
putative hit.  Often residues at the beginning and ends of windows are
determined to be nonhomologous and are not included in the
envelope. In this search,  99 envelopes are defined within the 84
%                         ^^^                                  ^^
windows. Note that the envelopes comprise only about 70\% of the
%                                                    ^^^
residues from the 84 windows, indicated by the drop of 0.1948\% to
%                                                      ^^^^^^^^
0.1331\%.
%%^^^^

After hit envelopes have been defined with the filter HMM, the two
remaining stages of the pipeline use the CM to score both the
conserved sequence and structure of each possible hit. In both of
these stages, constraints are derived from an HMM alignment of the
envelope and enforced as \emph{bands} on the CM dynamic programming
matrices (more on this in section~\ref{section:pipeline}). In the
first CM stage, the CYK algorithm (which is the SCFG analog of the
Viterbi HMM algorithm) is used to determine the best scoring maximum
likelihood alignment of any subsequence in each envelope to the CM. If
this alignment has a P-value of $\leq 0.0001$ then the envelope
survives to the final round. The envelopes passed to the final stage
may be shorter than those examined during the CYK stage. Specifically,
envelopes are redefined as starting and ending at the first and final
residues for which at least one alignment exists with a P-value $\leq
0.001$.

In the final round, the Inside algorithm (the SCFG analog of the HMM
Forward algorithm) is used to define final hit boundaries and
scores. Hits with scores above the reporting threshold were output, as
described above. In this search there were 56 such hits.
%                                          ^^

Finally, the running time of the search is reported, in CPU time and
elapsed time. This search took about 3 seconds (wall
clock time) running on 4 cores. 
%                      ^^^^^^^^^^^
\subsubsection{Truncated RNA detection}

Now, we come back to the topic of truncated hit detection.  As briefly
mentioned above, the pipeline statistics summary from our search above
reported that 360 residues were re-searched for truncated
%             ^^^
hits. Infernal explicitly looks at the 5' and 3' ends of target
sequences using specialized algorithms for detection of truncated
hits, in which part of the 5' and/or 3' end of the actual full length
sequence from the source organism's full genomic context is missing in
the input target sequence. Truncated hits will be most common in
sequence files consisting of unassembled sequencing reads. In our
search of the full archaeal genome above, no truncated hits were
found. However, there is an example of a truncated hit in the
\prog{cmscan} tutorial section below in a sequence from a metagenomics
sequencing survey.

Special dynamic programming algorithms are required for truncated hit
detection \citep{KolbeEddy09}, so sequences must be re-searched for
truncated hits after they are initially searched for standard
(non-truncated) hits. However, only the sequence ends must be
re-searched because Infernal assumes that only hits at sequence
terminii might be truncated. This is why our search above reported
that only 360 residues were re-searched for truncated hits. For each
%         ^^^
model, an expected maximum length of a hit\footnote{The expected
  maximum hit length is defined as the maximum of 1.25 times the
  consensus length of the model and the $W$ parameter computed with
  the QDB algorithm \citep{NawrockiEddy07}. $W$ is computed based on
  the transition probabilities of the model by \prog{cmbuild} and
  stored in the CM file.} 
is used to define the window length at the
beginning and end of each sequence which must be re-searched for
truncated hits. For our tRNA model the maximum expected length is
90, so exactly 360 residues were re-searched: the first and final 90
%^             ^^^                                                ^^
residues on the top strand, and the first and final 90 residues on
%                                                   ^^
the bottom strand.

There is one more aspect of truncated hit detection that is important
to mention here. Because Infernal expects that hit truncation only
occurs at the ends of sequences, 5' truncated hits are forced to start
at the first nucleotide of a sequence, 3' truncated hits are forced to
end at the final nucleotide of the sequence, and 5' \emph{and} 3'
truncated hits are forced to include the full sequence.

The annotation of truncated hits in \prog{cmsearch} output is slightly
different than for standard (non-truncated) hits. An example is
included in the next section below.

\subsection{Searching a CM database with a query sequence}

The \prog{cmscan} program is for annotating all the different
known/detectable RNAs in a given sequence. It takes a single query
sequence and a CM database as input. The CM database might be Rfam,
for example, or another collection of your choice. \prog{cmscan} is
new to version 1.1 of Infernal. It is designed to be useful for
sequence datasets from RNA-Seq experiments or metagenomics surveys,
for which one wants to know what families of RNAs are represented in
the sequence dataset.

\subsubsection{Step 1: create an CM database flatfile}

A CM ``database'' flatfile is simply a concatenation of individual CM
files. To create a database flatfile, you can either build and
calibrate individual CM files and concatenate them, or you can
concatenate Stockholm alignments and use \prog{cmbuild} to build a CM
database of all of them in one command, and then calibrate that
database with \prog{cmcalibrate}. Importantly, \prog{cmscan} can only
be used with calibrated CMs.

Let's create a tiny database called \otext{minifam.cm} containing the
tRNA model we've been working with, a 5S ribosomal RNA model, and a
Cobalamin riboswitch model. To save you time, calibrated versions of
the 5S and Cobalamin models are included in the \otext{tutorial/}
directory in the files \otext{5S\_rRNA.c.cm}, and
\otext{Cobalamin.c.cm}. These files were created using \prog{cmbuild}
and \prog{cmcalibrate} from the Rfam 10.1 seed alignments for 5S\_rRNA
(RF00001) and Cobalamin (RF00174), provided in
\otext{tutorial/5S\_rRNA.sto} and \otext{tutorial/Cobalamin.sto}. The
third model is the tRNA model from earlier in the tutorial 
(\otext{tRNA5.c.cm}). Feel free to build and calibrate these
models yourself if you'd like, but if you'd like to keep moving on
with the tutorial, use the pre-calibrated ones. To create the database,
simply concatenate the three provided files:

\user{cat tutorial/tRNA5.c.cm tutorial/5S\_rRNA.c.cm tutorial/Cobalamin.c.cm > minifam.cm}

\subsubsection{Step 2: compress and index the flatfile with cmpress}

The \prog{cmscan} program has to read a lot of CMs and their filter
HMMs in a hurry, and Infernal's ASCII flatfiles are bulky. To
accelerate this, \prog{cmscan} uses binary compression and indexing of
the flatfiles. To use \prog{cmscan}, you must first compress and
index your CM database with the \prog{cmpress} program:

\user{cmpress minifam.cm}

This will quickly produce:

\begin{sreoutput}
Working...    done.
Pressed and indexed 3 CMs and p7 HMM filters (3 names and 2 accessions).
Covariance models and p7 filters pressed into binary file:  minifam.cm.i1m
SSI index for binary covariance model file:                 minifam.cm.i1i
Optimized p7 filter profiles (MSV part)  pressed into:      minifam.cm.i1f
Optimized p7 filter profiles (remainder) pressed into:      minifam.cm.i1p
\end{sreoutput}

and you'll see these four new binary files in the directory. 

The \otext{tutorial} directory includes a copy of the
\otext{minifam.cm} file, which has already been pressed, so there
are example binary files \otext{tutorial/minifam.cm.i1\{m,i,f,p\}}
included in the tutorial.

Their format is ``proprietary'', which is an open source term of art
that means both ``We haven't found time to document them yet'' and ``We
still might decide to change them arbitrarily without telling you''.

\subsubsection{Step 3: search the CM database with cmscan}

Now we can analyze sequences using our CM database and
\prog{cmscan}. 

For example, the file \otext{tutorial/metag-example.fa} includes 3
sequences (whole genome shotgun sequencing reads) derived from samples
of ``whale fall'' carcasses in a metagenomics study
\citep{Tringe05}. These three sequences were chosen from the roughly
200,000 in the complete dataset because they include statistically
significant hits to the three models in our toy CM database.

To scan the sequences with our database: 

\user{cmscan minifam.cm tutorial/metag-example.fa}

The header and the first section of the output will look like:

\begin{sreoutput}
# cmscan :: search sequence(s) against a CM database
# INFERNAL 1.1.5 (Sep 2023)
# Copyright (C) 2023 Howard Hughes Medical Institute.
# Freely distributed under the BSD open source license.
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# query sequence file:                   tutorial/metag-example.fa
# target CM database:                    minifam.cm
# number of worker threads:              4
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Query:       AAGA01015927.1  [L=943]
Description: Metagenome sequence AHAI1002.g1, whole genome shotgun sequence
Hit scores:
 rank     E-value  score  bias  modelname  start    end   mdl trunc   gc  description
 ----   --------- ------ -----  --------- ------ ------   --- ----- ----  -----------
  (1) !   1.5e-23   92.7   0.0  5S_rRNA       59    174 +  cm    no 0.66  5S ribosomal RNA
  (2) !   9.6e-19   62.4   0.0  tRNA5        229    302 +  cm    no 0.62  -
  (3) !   6.1e-16   53.5   0.0  tRNA5        314    386 +  cm    no 0.59  -
\end{sreoutput}

\prog{cmscan} has identified three putative RNAs in the first query
sequence, one 5S rRNA and two tRNAs. The output fields are in the
same order and have the same meaning as in \prog{cmsearch}'s output.

Before we move on, this is a good opportunity to point out an
important difference between \prog{cmsearch} and \prog{cmscan} related
to the size of the search space (often referred to in the code and in
this guide as the $Z$ parameter).  The size of the search space for
\prog{cmscan} is double the length of the current (single) query
sequence (doubled because we're searching both strands) multiplied by
the number of models in the CM database (here, 3; for a Rfam search,
on the order of 1000). Because each query sequence is probably a
different size this means $Z$ changes for each query sequence. In
\prog{cmsearch}, the size of the search space is double the summed
length of all sequences in the database (again, doubled because both
strands are searched). This means that E-values may differ even for
the same individual CM vs. sequence comparison, depending on how you
do the search. The search space size also affects what filter
thresholds \prog{cmsearch} or \prog{cmscan} will use, which is
discussed more in section~\ref{section:pipeline}.

Now back to the \prog{cmscan} results. What follows the ranked list of
three hits are the hit alignments. These are constructed and annotated
the same as in \prog{cmsearch}. The 5S alignment is:

\begin{widesreoutput}
                             v               v       v         v         v             v            v                NC
                     (((((((((,,,,<<-<<<<<---<<--<<<<<<_______>>-->>>>-->>---->>>>>-->><<<-<<----<-<<-----<<____>>-- CS
         5S_rRNA   1 ccuggcggcCAUAgcggggcggAaacACccGauCCCAUCCCGaACuCggAAguUAAGcgcccuagcgccgggguAGuAcuggGgUGgGuGAcCaC 95 
                     CCUGGC:G C UAGCG:GG:GG   CACC:GA CCCAU CCG ACUC:GAAG  AA+C:CC:UAGC CCG::G  G:A: G+G  GGGU  CC C
  AAGA01015927.1  59 CCUGGCGGCCGUAGCGCGGUGGUCCCACCUGACCCCAUGCCGAACUCAGAAGUGAAACGCCGUAGCGCCGAUG--GUAGUGUG--GGGUCUCCCC 149
                     *************************************************************************..********..********** PP

                        v         v v          NC
                     --->>->->>->>>,))))))))): CS
         5S_rRNA  96 cUGggAaaccagguCgccgccaggc 120
                      UG :A:A::AGG   C:GCCAGGC
  AAGA01015927.1 150 AUGCGAGAGUAGGGAACUGCCAGGC 174
                     ************************* PP
\end{widesreoutput}

After the three sequences, you'll see a pipeline statistics summary
report:

\newpage

\begin{widesreoutput}
Internal CM pipeline statistics summary:
----------------------------------------
Query sequence(s):                                               1  (1886 residues searched)
Query sequences re-searched for truncated hits:                  1  (834.7 residues re-searched, avg per model)
Target model(s):                                                 3  (382 consensus positions)
Windows   passing  local HMM SSV           filter:              20  (0.3992); expected (0.35)
Windows   passing  local HMM Viterbi       filter:                  (off)
Windows   passing  local HMM Viterbi  bias filter:                  (off)
Windows   passing  local HMM Forward       filter:               4  (0.07976); expected (0.02)
Windows   passing  local HMM Forward  bias filter:               4  (0.07976); expected (0.02)
Windows   passing glocal HMM Forward       filter:               2  (0.06641); expected (0.02)
Windows   passing glocal HMM Forward  bias filter:               2  (0.06641); expected (0.02)
Envelopes passing glocal HMM envelope defn filter:               3  (0.03357); expected (0.02)
Envelopes passing  local CM  CYK           filter:               3  (0.03357); expected (0.0001)
Total CM hits reported:                                          3  (0.03222); includes 0 truncated hit(s)

# CPU time: 0.02u 0.00s 00:00:00.02 Elapsed: 00:00:00.09
//
\end{widesreoutput}

This report is similar to the one you saw earlier from
\prog{cmsearch}, but not identical. One big difference is that
\prog{cmscan} will report a summary per query sequence, instead of per
query model. In this case, the sequence was 943 residues long, so a
%                                           ^^^
total of 1886 residues were searched, since both strands were
%        ^^^^
examined. The next line reports that the average number of residues
re-searched for truncated hits per model was 834.7. An average is
%                                            ^^^^^
reported here because remember the number of residues re-searched per
model depends on the expected maximum size of a hit, which varies per
model, because only sequence terminii are examined for truncated hits
(see ``Truncated RNA detection'' above and
section~\ref{section:pipeline}).  The remainder of the output is the
same as in \prog{cmsearch} except that the fractional values are
averages per model. For example, for all three models a total of 3
envelopes survived the CYK filter stage, and those surviving envelopes
contained 3.357\% of the target sequence \emph{on average per model}.
%         ^^^^^

You may have noticed that the expected survival fractions are
different than they were in the \prog{cmsearch} example. This is
because the P-value filter thresholds are set differently depending on
the search space. For this example, the search space ($Z$) is roughly
only 6 Kb (1886 residues in the query sequence multiplied by 3 target
models), so the thresholds are set differently than for the
\prog{cmsearch} example which had a search space size of roughly 6
Mb. The exact thresholds used for various $Z$ values can be found in
Table~\ref{tbl:thresholds} in section~\ref{section:pipeline}.

\subsubsection{Truncated hit and local end alignment example}
Next, take a look at the \prog{cmscan} output for the second query
sequence \otext{AAFY01022046.1}. The bottom strand of this sequence
contains one hit to the Cobalamin riboswitch model, which is truncated
at the 5' end. Here is the alignment:

\begin{widesreoutput}
Hit alignments:
>> Cobalamin  Cobalamin riboswitch
 rank     E-value  score  bias mdl mdl from   mdl to       seq from      seq to       acc trunc   gc
 ----   --------- ------ ----- --- -------- --------    ----------- -----------      ---- ----- ----
  (1) !   6.6e-09   29.9   0.0  cm       31      190 ~]         934         832 - ~. 0.92    5' 0.48

                                 ???              v           v      v    v                                     ???? NC
                     ~~~~~~______>>>,,,,,(((,,,<.<<<<_______>>>>>,,<<<____>>>,<<<---<<<<~~~~~~>>>>---->>>,,,,)))]]]] CS
       Cobalamin   1 <[30]*agucggggguuAAaaGGGAAc.ccGGUGcaAauCCgggGCuGcCCCCgCaACUGUAAgcGg*[61]*cCgcgAGCCAGGAGACCuGCCg 173
                           + + G     + AA: GGAA: : GGUG AAAUCC ::+C:G CCC  C:ACUGUAA:C:        :G:+AG+CAG A AC :  C 
  AAFY01022046.1 934 <[ 0]*GUAGGCAAAAGGAAGAGGAAGgAUGGUGGAAAUCCUUCACGGGCCCGGCCACUGUAACCAG*[ 4]*UUGGAAGUCAG-AUACUCUUCU 849
                     ......44455677788899******9899***********************************96...7..69*********.9********* PP

                     ??                NC
                     ]]::::::::::::::: CS
       Cobalamin 174 ucaaaaauaguaucacc 190
                       +AA +  G+A+C+ C
  AAFY01022046.1 848 AUUAAGGCGGAAACUAC 832
                     ***************** PP
\end{widesreoutput}

This alignment has some important differences with the ones we've seen
so far because it is of a truncated hit. First, notice that the
\otext{trunc} column reads \otext{5'} indicating that Infernal
predicts the 5' end (beginning) of this Cobalamin riboswitch is
missing. (Note that this hit is on the bottom (reverse) strand so the
Cobalamin hit is actually predicted to extend past the end of the
input sequence (past residue 934), but on the opposite strand.) The 5'
%                            ^^^
end of the alignment indicates this with special annotation: the
\otext{<[30]*} in the model line indicates that the missing sequence
%        ^^
is expected to align to the 30 5'-most positions of the Cobalamin
%                           ^^
model (i.e. about 30 residues are missing) and the first \otext{a}
%                 ^^
residue in this line corresponds to model position 31; the
%                                                  ^^
\otext{<[0]*} annotation in the sequence line indicates that there are
no observed residues which align to those 31 positions and the first
%                                         ^^
\otext{G} residue is at position 934 of the sequence, which is the
%                                ^^^
first position (on the opposite strand) of the query sequence. If,
alternatively, this sequence was 3' truncated, or both 5' \emph{and}
3' truncated, there would be analogous annotation at the 3' end of the
alignment. Also notice the \verb+~]+ and \verb+~.+ symbols following
the model and sequence coordinates, respectively. The \verb+~+
leftmost symbol in both of these pairs indicate that this hit is
truncated at the 5' end. To make sense of this alignment display, it
may help to look at the \prog{cmalign} alignment of the same sequence
to the same model on page~\pageref{cmalign-cobalamin}, which shows all
model and sequence positions.

Truncated hit alignments also contain different annotation in the
\otext{NC} lines. Instead of only containing blank spaces or \otext{v}
characters indicating negative scoring noncanonical basepairs,
\otext{?} characters are used to denote basepairs for which the other
half is missing due to the truncation. For example, the second
\otext{g} in the first model line corresponds to the right half of a
basepair for which the left half is included in the stretch of 30 5'
%                                                              ^^
truncated model positions, so it is annotated with a \otext{?}.  A
\otext{?} is used because it is impossible to tell if such basepairs
are negative scoring non-canonicals or not since we don't know the
identity of the other half.

This hit alignment also demonstrates another type of annotation not
yet seen in the previous examples, for \emph{local end}
alignments. Notice the stretch of six \verb+~+ characters towards the
end of the first \otext{CS} line, at the same position as the string
\otext{*[61]*} in the model line immediately below. This indicates a a
%        ^^
special type of alignment called a local end. Local ends occur when a
large insertion or deletion is used in the optimal alignment at
reduced penalty \citep{KleinEddy03} and allow Infernal to be tolerant
of the insertion and/or deletion of RNA substructures not modeled by
the CM. An example of how local ends enable remote homology detection
in RNase P can be see in Figure~\ref{fig:bsu-alignment} on
page~\pageref{fig:bsu-alignment}. In this case, 61 model positions are
%                                               ^^                 
deleted and 4 residues are inserted in the sequence (indicated by
%           ^                 
\otext{*[ 4]*} at the corresponding positions in the sequence line).
%         ^                 
It is possible for zero sequence residues to be inserted by a local
end, and for the number of residues inserted in a local end to exceed
the number of model positions. Note that a single \otext{7} annotates
%                                                        ^                 
the posterior probability of the four sequence residues in the
\otext{PP} line. This means that the average posterior probability for
these four residues is between 65 and 75\%. If no sequence residues
%                              ^^     ^^
were in this EL, the \otext{PP} annotation would be a gap (\otext{.})
character.

\subsection{Searching the Rfam CM database with a query sequence}

The Rfam database \url{http://rfam.xfam.org/} is a collection
of RNA families, each represented by a CM and multiple sequence
alignment used to build that CM. As of September 2023, the current release
%                                      ^^^^^^^^^^^^^^
is 14.9, which includes 4108 families. The Rfam website allows
web-based searches using \prog{cmscan} of the Rfam CM database against
query sequences that the user can upload. Alternatively, you can perform
the same searches by running \prog{cmscan} locally, as shown in this
example. By searching all of Rfam with your sequence dataset, you will
be annotating your dataset for most known types of structural RNAs
with a single command.

To complete this step of the tutorial you'll need to download the Rfam
14.9 CM file from here: \\
\url{https://ftp.ebi.ac.uk/pub/databases/Rfam/14.9/Rfam.cm.gz}
and gunzipping it, like this: 

\user{curl -k -L -o Rfam.cm.gz https://ftp.ebi.ac.uk/pub/databases/Rfam/14.9/Rfam.cm.gz}

\user{gunzip Rfam.cm.gz}

Then, as in the previous example, you'll need to run \prog{cmpress} on
this CM database:

\user{cmpress Rfam.cm}

The next step is to run \prog{cmscan}. In order to reproduce how Rfam
searches are performed \citep{Nawrocki15} several command line
options are required. Each of these options is explained below. The
full command is (split up into two lines so it fits on the page): 

\indent\indent\small\verb+> cmscan --rfam --cut_ga --nohmmonly --tblout mrum-genome.tblout --fmt 2 \+\\
\indent\indent\small\verb+> --clanin tutorial/Rfam.14.9.clanin Rfam.cm tutorial/mrum-genome.fa > mrum-genome.cmscan+\\

This command will take at least several minutes and possibly up to
about 30 minutes depending on the number of cores and speed of your
computer. 

The command line options used in the above command are as follows:

\begin{sreitems}{\emprog{--nohmmonly}}
\item[\otext{--rfam}] Specifies that the filter pipeline run in fast
  mode, with the same strict filters that are used for Rfam searches
  and for other sequence databases larger than 20 Gb (see
  section~\ref{section:pipeline}). 
%
\item[\otext{--cut\_ga}] Specifies that the special Rfam \emph{GA}
  (gathering) thresholds be used to determine which hits are
  reported. These thresholds are stored in the \prog{Rfam.cm} file.
  Each model has its own GA bit score threshold, which was determined
  by Rfam curators as the bit score at and above which all hits are
  believed to be true homologs to the model. These determinations were made based on
  observed hit results against the large Rfamseq database used by
  Rfam~\citep{Nawrocki15}.
%
\item[\otext{--nohmmonly}] All models, even those with zero basepairs,
  are run in CM mode (not HMM mode). This ensures all GA cutoffs, which were
  determined in CM mode for each model, are valid.
%
\item[\otext{--tblout}] Specifies that a tabular output
  file should be created, see section~\ref{section:tabular}.
%
\item[\otext{--fmt 2}] The tabular output file will be in format 2,
  which includes annotation of overlapping hits. See
  page~\pageref{tabular-format2} for a complete description of this
  format. 
%
\item[\otext{--clanin}] Clan information should be read
  from the file \prog{tutorial/Rfam.14.9.claninfo}. This file lists
  which models belong to the same clan. Clans are groups of models
  that are homologous and therefore it is expected that some hits to
  these models will overlap. For example, the LSU\_rRNA\_archaea and
  LSU\_rRNA\_bacteria models are both in the same clan.
%
\end{sreitems}

When the \prog{cmscan} command finishes running, the file
\prog{mrum-genome.cmscan} will contain the standard output of the
program. This file will be similar to what we saw in the earlier example of
\prog{cmscan}. The file \prog{mrum-genome.tblout} has also been
created, which is a tabular representation of all hits, one line per
hit.  Take a look at this file. The first two lines are comment lines
(prefixed with \prog{\#} characters) with the labels of each of the 27
columns of data in the file. Each subsequent line has 27 space
delimited tokens. The specific meaning of these tokens is described in
detail in section~\ref{section:tabular}.  Below I'm including the
first 24 lines of the file, with columns 3-5, 7-9 and 13-16 removed
(replaced with \prog{...}) so that the text will fit on this page:

\begin{tinysreoutput}
#idx target name            ... clan name ... seq from   seq to strand ...  score   E-value inc olp anyidx afrct1 afrct2 winidx wfrct1 wfrct2 description of target
#--- ---------------------- ... --------- ... -------- -------- ------ ... ------ --------- --- --- ------ ------ ------ ------ ------ ------ ---------------------
1    LSU_rRNA_archaea       ... CL00112   ...   762872   765862      + ... 2763.5         0  !   ^       -      -      -      -      -      - Archaeal large sub...
2    LSU_rRNA_archaea       ... CL00112   ...  2041329  2038338      - ... 2755.0         0  !   ^       -      -      -      -      -      - Archaeal large sub...
3    LSU_rRNA_bacteria      ... CL00112   ...   762874   765862      + ... 1872.9         0  !   =       1  1.000  0.999      "      "      " Bacterial large su...
4    LSU_rRNA_bacteria      ... CL00112   ...  2041327  2038338      - ... 1865.5         0  !   =       2  1.000  0.999      "      "      " Bacterial large su...
5    LSU_rRNA_eukarya       ... CL00112   ...   763018   765851      + ... 1581.3         0  !   =       1  1.000  0.948      "      "      " Eukaryotic large s...
6    LSU_rRNA_eukarya       ... CL00112   ...  2041183  2038349      - ... 1572.1         0  !   =       2  1.000  0.948      "      "      " Eukaryotic large s...
7    SSU_rRNA_archaea       ... CL00111   ...  2043361  2041888      - ... 1552.0         0  !   ^       -      -      -      -      -      - Archaeal small sub...
8    SSU_rRNA_archaea       ... CL00111   ...   760878   762351      + ... 1546.5         0  !   ^       -      -      -      -      -      - Archaeal small sub...
9    SSU_rRNA_bacteria      ... CL00111   ...  2043366  2041886      - ... 1161.9         0  !   =       7  0.995  1.000      "      "      " Bacterial small su...
10   SSU_rRNA_bacteria      ... CL00111   ...   760873   762353      + ... 1156.4         0  !   =       8  0.995  1.000      "      "      " Bacterial small su...
11   SSU_rRNA_eukarya       ... CL00111   ...  2043361  2041891      - ...  970.4    3e-289  !   =       7  1.000  0.998      "      "      " Eukaryotic small s...
12   SSU_rRNA_eukarya       ... CL00111   ...   760878   762348      + ...  963.8    3e-287  !   =       8  1.000  0.998      "      "      " Eukaryotic small s...
13   SSU_rRNA_microsporidia ... CL00111   ...  2043361  2041891      - ...  919.9  2.3e-277  !   =       7  1.000  0.998      "      "      " Microsporidia small..
14   SSU_rRNA_microsporidia ... CL00111   ...   760878   762348      + ...  917.2  1.6e-276  !   =       8  1.000  0.998      "      "      " Microsporidia small..
15   RNaseP_arch            ... CL00002   ...  2614544  2614262      - ...  184.9   3.4e-50  !   *       -      -      -      -      -      - Archaeal RNase P
16   Archaea_SRP            ... CL00003   ...  1064321  1064634      + ...  197.6   2.1e-45  !   *       -      -      -      -      -      - Archaeal signal re...
17   FMN                    ... -         ...   193975   193837      - ...  115.2   2.1e-24  !   *       -      -      -      -      -      - FMN riboswitch (RF...
18   tRNA                   ... CL00001   ...   735136   735208      + ...   72.1   1.5e-12  !   *       -      -      -      -      -      - tRNA
19   tRNA                   ... CL00001   ...  2350593  2350520      - ...   71.0     3e-12  !   *       -      -      -      -      -      - tRNA
20   tRNA                   ... CL00001   ...  2680310  2680384      + ...   70.9   3.2e-12  !   *       -      -      -      -      -      - tRNA
21   tRNA                   ... CL00001   ...  2351254  2351181      - ...   69.7   6.7e-12  !   *       -      -      -      -      -      - tRNA
22   tRNA                   ... CL00001   ...   361676   361604      - ...   69.5   7.6e-12  !   *       -      -      -     
\end{tinysreoutput}

This tabular format includes the target model name, sequence name (in
column 3, which is omitted above to save space), clan name, sequence
coordinates, bit score, E-value and more. Because the \prog{--fmt 2}
option was used, this file includes information on which hits overlap
with other hits, starting at the column labelled ``olp'' and ending
with ``wfrct2''. Hits with the ``*'' character in the ``olp'' column
do not overlap with any other hits. Those with ``\verb+^+'' do overlap with
at least one other hit, but none of those overlapping hits have a
better score (that occurs higher in the list). Those with ``='' also
overlap with at least one other hit that does have a better score, the
index of which is given in the ``anyidx'' column. For more detailed
explanation of these columns, see page~\pageref{tabular-format2}.

The top two hits are both to the \prog{LSU\_rRNA\_archaea}
model. These are the two copies of LSU rRNA in the
\emph{Methanobrevibacter ruminantium} genome. Hits number 3 and 4
%                                                         ^     ^
are to the \prog{LSU\_rRNA\_bacteria} model and overlap with hits 1
%                                                                 ^
and 2 nearly completely (hit 1 is from sequence positions 762872 to 765862 and
%   ^                        ^                            ^^^^^^    ^^^^^^
hit 3 is from sequence positions 762874 to 765862). This overlap is not
%   ^                            ^^^^^^    ^^^^^^
surprising because the bacterial and archaeal LSU rRNA models are very
similar, and so are assigning high scores to the same
subsequences. Further, hit 5 is to \prog{LSU\_rRNA\_eukarya} and also
%                          ^
overlaps hits 1 and 3. Because these three LSU models are all expected
%             ^     ^
to produce overlapping hits due to their homology, Rfam has grouped
them into the same \emph{clan}, note the ``CL00112'' value in the ``clan
%                                          ^^^^^^^
name'' column for all three hits. This clan information was provided
in the \prog{rfam.14.9.claninfo} input file we provided to \prog{cmscan}
by using the \prog{--clanin} option. 

The ``olp'' column indicates that hit 1 is the highest scoring of the
three overlapping hits because it contains the ``\verb+^+''
character. Hits 3 and 5 both have ``='' in the ``olp'' column
%               ^     ^
indicating that there is another hit to another model which overlaps
these hits and has a better score.

If you were using these results to produce annotations for the 
\emph{Methanobrevibacter ruminantium} genome, you may want to ignore 
any hits that have higher scoring overlaps. To do this you can just
remove any hits with ``='' in the ``olp'' column. Alternatively, you
can have these hits not printed to the tabular output file by
additionally providing the \prog{--oskip} option to \prog{cmscan}.
You can also modify the overlap annotation behavior with
\prog{--oclan} option which restricts the annotation of overlaps to
hits for models within the same clan. Overlapping hits from models
that are not in the same clan will not be marked as
overlaps, instead they will marked as ``*'' in the ``olp'' field.

For more information on using Infernal for genome annotation see a
similar example in the Rfam documentation
(https://rfam.readthedocs.io/en/latest/genome-annotation.html)

\subsection{Creating multiple alignments with cmalign}
The file \otext{tutorial/mrum-tRNAs10.fa} is a FASTA file containing
the 10 tRNA hits above the inclusion threshold (with an E-value
less than $0.01$) found by \prog{cmsearch} in our 
search of \emph{M. ruminantium} genome\footnote{The
  \otext{-A <f>} option to \prog{cmsearch} can be used to save a
  Stockholm formatted multiple alignment of all hits above the
  inclusion threshold to file \otext{<f>}.}
To align these 10 sequences to our example tRNA model and make a multiple sequence
alignment:

\user{cmalign tRNA5.cm tutorial/mrum-tRNAs10.fa}

The output of this is a Stockholm format multiple alignment file:

\begin{tinysreoutput}
# STOCKHOLM 1.0
#=GF AU Infernal 1.1.5

mrum-tRNA.1          GGAGCUAUAGCUCAAU..GGC..AGAGCGUUUGGCUGACAU........................................CCAAAAGGUUAUGGGUUCGAUUCCCUUUAGCCCCA
#=GR mrum-tRNA.1  PP ****************..***..******************........................................***********************************
mrum-tRNA.2          GGGCCCGUAGCUCAGU.uGGG..AGAGCGCUGCCCUUGCAA........................................GGCAGAGGCCCCGGGUUCAAAUCCCGGUGGGUCCA
#=GR mrum-tRNA.2  PP ****************.****..******************........................................***********************************
mrum-tRNA.3          GGGCCCAUAGCUUAGCcaGGU..AGAGCGCUCGGCUCAUAA........................................CCGGGAUGUCAUGGGUUCGAAUCCCAUUGGGCCCA
#=GR mrum-tRNA.3  PP *********************..******************........................................***********************************
mrum-tRNA.4          AGGCUAGUGGCACAGCcuGGU.cAGCGCGCACGGCUGAUAA........................................CCGUGAGGUCCUGGGUUCGAAUCCCAGCUAGCCUA
#=GR mrum-tRNA.4  PP ***************999***.*******************........................................***********************************
mrum-tRNA.5          CCCGACUUAGCUCAAUuuGGC..AGAGCGUUGGACUGUAGA........................................UCCAAAUGUUGCUGGUUCAAGUCCGGCAGUCGGGA
#=GR mrum-tRNA.5  PP *********************..******************........................................***********************************
mrum-tRNA.6          GCUUCUAUGGGGUAAU.cGGC.aAACCCAUCGGACUUUCGA........................................UCCGAUAA-UCCGGGUUCAAAUCCCGGUAGAAGCA
#=GR mrum-tRNA.6  PP ****************.****.*******************........................................********.9*************************
mrum-tRNA.7          GCUCCGAUGGUGUAGUccGGCcaAUCAUUUCGGCCUUUCGA........................................GCCGAAGA-CUCGGGUUCGAAUCCCGGUCGGAGCA
#=GR mrum-tRNA.7  PP ********************9999*****************........................................********.**************************
mrum-tRNA.8          GCGGUGUUAGUCCAGCcuGGU.uAAGACUCUAGCCUGCCAC........................................GUUAGAGA-CCCGGGUUCAAAUCCCGGACGCCGCA
#=GR mrum-tRNA.8  PP ***************999***.*******************........................................********.**************************
mrum-tRNA.9          GCCGGGGUGGCUCAGC.uGGU.uAGAGCGCACGGCUC----auaggguaacuaagcgugcucugacuuuuuuccugggauaCCGUGAGAUCGCGGGUUCGAAUCCCGCCCCCGGCA
#=GR mrum-tRNA.9  PP ****************.****.************995....678************************************************************************
mrum-tRNA.10         GGUUCUAUAGUUUAAC.aGGU..AAAACAACUGGCUGUUAA........................................CCGGCAGA-UAGGAGUUCGAAUCUUCUUAGAACCG
#=GR mrum-tRNA.10 PP ****************.****..******************........................................********.9*************************
#=GC SS_cons         (((((((,,<<<<___..___.._>>>>,<<<<<_______~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~>>>>>,,,,,<<<<<_______>>>>>))))))):
#=GC RF              gCcggcAUAGcgcAgU..GGu..AgcgCgccagccUgucAa~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~gcuggAGgUCCgggGUUCGAUUCcccGUgccgGca
//
\end{tinysreoutput}

The first thing to notice here is that \prog{cmalign} uses both lower
case and upper case residues, and it uses two different characters for
gaps. This is because there are two different kinds of columns:
``match'' columns in which residues are assigned to match states and
gaps are treated as deletions relative to consensus, and ``insert''
columns where residues are assigned to insert states.
In a match column, residues are upper case, and a '\otext{-}'
character means a deletion relative to the consensus. In an insert
column, residues are lower case, and a '\otext{.}' is padding.  A '\otext{-}' deletion
has a cost: transition probabilities were assessed, penalizing the
transition into and out of a deletion. A '\otext{.}' pad has no cost per se;
instead, the sequence(s) with insertions are paying transition
probabilities into and out of their inserted residue.

Actually, there's two types of insert columns: standard insert columns
where residues are assigned to insert states and less common ``local end''
insertion columns where residues are assigned to a special local end
state called the ``EL'' state. There is one EL state per model and it
allows a CM to permit large insertions or deletions in the structure
present in the alignment the model was built from, at a reduced
cost. This can help Infernal detect remote homologs in some cases. We
saw an example of this in our Cobalamin \prog{cmscan} hit
above. Another example is shown in Figure~\ref{fig:bsu-alignment} on
page~\pageref{fig:bsu-alignment}. Columns containing EL insertions are
denoted in the \otext{\#=GC RF} annotation, described next.

Take a look at the final two lines of the alignment. The \otext{\#=GC
  RF} line is Stockholm-format \emph{reference coordinate annotation},
with a residue marking each column that the CM considered to be
consensus, a '\otext{.}' marking insert columns and a '\verb+~+'
marking EL insert columns. For match columns, upper case residues
denote strongly conserved columns, and lower case denotes weakly
conserved ones. The \otext{\#=GC SS\_cons} line gives the consensus
secondary structure of the model. The symbols here have the same
meaning that they did in a pairwise alignment from \prog{cmsearch},
for a detailed description see section~\ref{section:formats}. As in the
RF annotation, EL columns are indicated by '\verb+~+'. In this
alignment, \otext{mrum-tRNA.9} is the only sequence that uses the EL
state.

Important: both standard and local end insertions in a CM are
\emph{unaligned} (the same way that insertions are unaligned in
profile HMM alignments produced by the HMMER package). Suppose one
sequence has an insertion of length 10 and another has an insertion of
length 2 in the same place in the model. The alignment will show ten
insert columns, to accomodate the longest insertion.  The residues of
the shorter insertion are thrown down in an arbitrary order. (If you
must know: by arbitrary Infernal convention, the insertion is divided
in half; half is left-justified, and the other half is
right-justified, leaving '.' characters in the middle.)  Notice that
in the previous paragraph I oh-so-carefully said residues are
``assigned'' to a state, not ``aligned''. For match states, assigned
and aligned are the same thing: a one-to-one correspondence between a
residue and a consensus match state in the model. But there may be one
\emph{or more} residues assigned to the same insert state.

Don't be confused by the unaligned nature of CM insertions. You're
sure to see cases where lower-case inserted residues are ``obviously
misaligned''. This is just because Infernal isn't trying to ``align''
them in the first place: it is assigning them to unaligned
insertions. For example of an obvious misalignment look at sequences
\otext{mrum-tRNA.6} and \otext{mrum-tRNA.8} in the above example
alignment. The first inserted (lowercase) \otext{c} in
\otext{mrum-tRNA.6} would be better aligned with respect to
\otext{mrum-tRNA.8} if it were placed one position to the left.

Enough about the sequences in the alignment. Now notice all those
\otext{PP} annotation lines. That's posterior probability annotation,
as in the single sequence alignments that \prog{cmscan} and
\prog{cmsearch} showed. This essentially represents the confidence
that each residue is aligned where it should be.

Er, I mean, ``assigned'', not ``aligned''. The posterior probability
assigned to an inserted residue is the probability that it is assigned
to the insert state that corresponds to that column, or the EL state
in the case of local end insertions. Because the same insert state
might correspond to more than one column, the probability on an insert
residue is \emph{not} the probability that it belongs in that
particular column; again, where there's a choice of column for
inserted residues, that choice is arbitrary. 

\subsubsection{cmalign assumes sequences may be truncated}
\prog{cmalign} assumes that input sequences may be truncated and uses
specialized dynamic programming algorithms to align them appropriately
based on that assumption. Importantly, output
alignments from \prog{cmalign} do not indicate when sequences have
been truncated, as those from \prog{cmsearch} and \prog{cmscan} do. As
an example, use the \ccode{tutorial/Cobalamin.c.cm} model to align the
truncated Cobalamin hit (in the file \ccode{tutorial/Cobalamin.fa})
from the \prog{cmscan} example above (I've split up the alignment
below so it will fit on the page):

\user{cmalign tutorial/Cobalamin.c.cm tutorial/Cobalamin.fa}

\label{cmalign-cobalamin}
\begin{sreoutput}
# STOCKHOLM 1.0
#=GF AU Infernal 1.1.5

Cobalamin.1         ------------------------------GUAGGCAAAAGGAAGAGGAAGgAUGGUGGAAAUCCUUCACGGGCCCGGCCA
#=GR Cobalamin.1 PP ..............................44455677788899******9899***************************
#=GC SS_cons        :::::::::::::::[[[[[[<<<____________>>>,,,,,(((,,,<.<<<<_______>>>>>,,<<<____>>>,
#=GC RF             uuaaauuauucuggugacGGUcccccuuaaagucagggguuAAaaGGGAAc.ccGGUGcaAauCCgggGCuGcCCCCgCaA


Cobalamin.1         CUGUAACCAG---------------------------------auuu----------------------------UUGGAA
#=GR Cobalamin.1 PP ********96.................................5789............................69****
#=GC SS_cons        <<<---<<<<------<<<<<<-----<<<-<<<<<<______~~~~>>>>>>--->>>>>>>>>---------->>>>--
#=GC RF             CUGUAAgcGggaagcaccccccaaaauGCCACUGgcccguaag~~~~ggcCGGGAAGGCggggggaagcgaugAccCgcgA


Cobalamin.1         GUCAG-AUACUCUUCUAUUAAGGCGGAAACUAC
#=GR Cobalamin.1 PP *****.9**************************
#=GC SS_cons        -->>>,,,,)))]]]]]]:::::::::::::::
#=GC RF             GCCAGGAGACCuGCCgucaaaaauaguaucacc

\end{sreoutput}

If you look back to the \prog{cmscan} hit alignment for this sequence,
you'll notice that this looks a little different. Instead of the
\otext{<[30]*} annotation in the model line indicating 30 model
%        ^^
positions were missing due to a presumed 5' truncation, \prog{cmalign}
includes those 30 positions, and their secondary structure
%              ^^
annotation. The local end alignment annotation in the middle of the
second line is different too. In the \prog{cmscan} hit alignment the
model line included \otext{*[61]*} indicating that 61 model positions
%                            ^^                    ^^
were skipped to a local end insertion. In the \prog{cmalign} output
those 61 positions are included, aligned to gaps in the sequence. 
%     ^^
\prog{cmalign} also shows the identity of the four inserted residues
in the EL state, and annotates each with a posterior probability,
whereas \prog{cmscan} only indicated the number of inserted EL
residues and their average posterior probability. 

If you examine the \prog{cmscan} and \prog{cmalign} alignments closely
you'll notice that they are identical (the same sequence residues are
aligned to the same model positions in each) and include identical
posterior probability annotation. While this will be the case for the
large majority of sequences, it is not guaranteed by Infernal, and
sometimes a \prog{cmsearch} or \prog{cmscan} alignment of a hit will
differ from a \prog{cmalign} alignment of the same hit. There are
several technical reasons for this, including the fact that the HMM
band constraints used by \prog{cmalign} can differ from those used by
\prog{cmsearch} and \prog{cmscan} which in rare cases can lead to a
different alignment or different posterior probabilities. Another
reason is that while \prog{cmalign} always assumes a sequence may be
5' or 3' truncated, \prog{cmsearch} and \prog{cmscan} only allow
certain types of truncation (5', 3' or both) at sequence
terminii. The types of truncated alignment allowed modify the
dynamic programming alignment algorithm, so this difference can also
result in different alignments. However, most of the time the
alignments will be identical, and when they are different they will
usually only be slightly different.

\subsection{Searching a sequence database for RNAs with unknown or no
 secondary structure}

Some functional RNAs, including many types of small nucleolar RNAs
(snoRNAs) do not conserve a secondary structure. It is more
appropriate to model such RNA families with a profile HMM than a CM,
because CMs only differ from profile HMMs in their modeling of
basepairs and because profile HMM scoring algorithms are more
efficient than CM ones.

Likewise, profile HMMs are more appropriate than CMs for modeling RNAs
of unknown structure. For such RNAs, a profile HMM search may reveal
additional homologs that can help elucidate the conserved structure
for the family, if there is one. Given the multiple homologs, you can
use other programs like RNAalifold~\citep{Bernhart08,Hofacker02} or
Locarna~\citep{Will12} to try to predict the conserved secondary
structure of a collection of homologous RNAs. Currently, Infernal
itself does not have the capability of predicting structure, but it's
predecessor COVE did with the \prog{covet} program, still available at
\url{eddylab.org/software/cove/cove.tar.Z}.

Infernal automatically detects when a model has zero basepairs and
uses efficient profile HMM algorithms in \prog{cmsearch} and
\prog{cmscan}. As an example, let's revisit the construction of our
tRNA model but pretend that we do not know the conserved cloverlear
structure of tRNA. The file 
\otext{tutorial/tRNA5-noss.sto} is identical to the file 
\otext{tutorial/tRNA5.sto} except that it does not contain consensus
secondary structure annotation. 
To build a structureless tRNA model from this alignment, execute:

\user{cmbuild --noss tRNA5-noss.cm tutorial/tRNA5-noss.sto}

The \otext{--noss} option is required when the alignment file lacks
structure annotation. By forcing the use of \otext{--noss} we're
forcing the user to realize that having no secondary structure is a
special case for Infernal.

The output will be similar to be the earlier example, but not
identical:

\begin{sreoutput}
# cmbuild :: covariance model construction from multiple sequence alignments
# INFERNAL 1.1.5 (Sep 2023)
# Copyright (C) 2023 Howard Hughes Medical Institute.
# Freely distributed under the BSD open source license.
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# CM file:                                            tRNA5-noss.cm
# alignment file:                                     tutorial/tRNA5-noss.sto
# ignore secondary structure, if any:                 yes
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#                                                                      rel entropy
#                                                                      -----------
# idx    name                     nseq eff_nseq   alen  clen  bps bifs    CM   HMM description
# ------ -------------------- -------- -------- ------ ----- ---- ---- ----- ----- -----------
       1 tRNA5-noss                  5     5.00     74    72    0    0 0.552 0.552 
#
# CPU time: 0.05u 0.00s 00:00:00.05 Elapsed: 00:00:00.06
\end{sreoutput}

The output reports that this model has 0 basepairs (``bps'') (the
earlier example had 21) and that the CM relative entropy is different
from before because basepairs have been ignored. Now the CM and HMM 
relative entropy are identical, because the CM can be mirrored nearly
identically by a profile HMM.

Remember that after we built our tRNA CM with structure above, we
needed to use \prog{cmcalibrate} to calibrate the model before we
could perform searches. Importantly, zero basepair models do not need
to be calibrated prior to running \prog{cmsearch}\footnote{Zero
  basepair models do need to be calibrated before they can be used
  with \prog{cmscan} however.}, so we can skip that step here. 

Now, let's repeat the search of the \emph{M. ruminantium} genome but
using our structureless tRNA model: 

\user{cmsearch tRNA5-noss.cm tutorial/mrum-genome.fa}\\

This search is faster than the first one because only profile HMM
algorithms are used this time around. Take a look at the list of hits:

\begin{sreoutput}
 rank     E-value  score  bias  sequence      start     end   mdl trunc   gc  description
 ----   --------- ------ -----  ----------- ------- -------   --- ----- ----  -----------
  (1) !   3.8e-09   38.2   0.0  NC_013790.1  362022  361960 - hmm     - 0.46  Methanobrevibacter rumi ...
  (2) !   3.1e-07   32.1   0.0  NC_013790.1  762496  762559 + hmm     - 0.64  Methanobrevibacter rumi ...
  (3) !   3.1e-07   32.1   0.0  NC_013790.1 2041698 2041635 - hmm     - 0.64  Methanobrevibacter rumi ...
  (4) !   4.9e-06   28.4   0.0  NC_013790.1 2585264 2585203 - hmm     - 0.56  Methanobrevibacter rumi ...
  (5) !   6.2e-06   28.0   0.0  NC_013790.1  735143  735198 + hmm     - 0.59  Methanobrevibacter rumi ...
  (6) !     7e-06   27.9   0.0  NC_013790.1 2351247 2351188 - hmm     - 0.57  Methanobrevibacter rumi ...
  (7) !   1.2e-05   27.2   0.0  NC_013790.1  662193  662251 + hmm     - 0.64  Methanobrevibacter rumi ...
  (8) !   1.2e-05   27.1   0.0  NC_013790.1 2186009 2185951 - hmm     - 0.51  Methanobrevibacter rumi ...
  (9) !   4.4e-05   25.4   0.0  NC_013790.1 2585183 2585118 - hmm     - 0.56  Methanobrevibacter rumi ...
 (10) !   6.6e-05   24.8   0.0  NC_013790.1 1873882 1873820 - hmm     - 0.63  Methanobrevibacter rumi ...
 (11) !   0.00015   23.6   0.0  NC_013790.1  360882  360824 - hmm     - 0.51  Methanobrevibacter rumi ...
 (12) !   0.00061   21.7   0.0  NC_013790.1  361910  361851 - hmm     - 0.38  Methanobrevibacter rumi ...
 (13) !   0.00095   21.1   0.0  NC_013790.1 2350586 2350528 - hmm     - 0.58  Methanobrevibacter rumi ...
 (14) !    0.0018   20.3   0.0  NC_013790.1  995341  995267 - hmm     - 0.51  Methanobrevibacter rumi ...
 (15) !    0.0027   19.7   0.0  NC_013790.1   97728   97788 + hmm     - 0.49  Methanobrevibacter rumi ...
 (16) !    0.0029   19.6   0.0  NC_013790.1 2186083 2186024 - hmm     - 0.50  Methanobrevibacter rumi ...
 (17) !    0.0032   19.5   0.0  NC_013790.1 2130421 2130351 - hmm     - 0.59  Methanobrevibacter rumi ...
 (18) !    0.0047   19.0   0.0  NC_013790.1  360727  360670 - hmm     - 0.43  Methanobrevibacter rumi ...
 (19) !    0.0058   18.7   0.0  NC_013790.1 1160527 1160608 + hmm     - 0.59  Methanobrevibacter rumi ...
 (20) !    0.0077   18.3   0.0  NC_013790.1  361056  360994 - hmm     - 0.40  Methanobrevibacter rumi ...
 ------ inclusion threshold ------
 (21) ?     0.012   17.7   0.0  NC_013790.1 2151679 2151737 + hmm     - 0.56  Methanobrevibacter rumi ...
 (22) ?     0.019   17.0   0.0  NC_013790.1 2327123 2327043 - hmm     - 0.62  Methanobrevibacter rumi ...
 (23) ?     0.025   16.7   0.0  NC_013790.1  360973  360920 - hmm     - 0.54  Methanobrevibacter rumi ...
 (24) ?     0.038   16.1   0.0  NC_013790.1 2350982 2350919 - hmm     - 0.50  Methanobrevibacter rumi ...
 (25) ?     0.039   16.0   0.0  NC_013790.1  361671  361606 - hmm     - 0.50  Methanobrevibacter rumi ...
 (26) ?     0.039   16.0   0.0  NC_013790.1 2680176 2680227 + hmm     - 0.62  Methanobrevibacter rumi ...
 (27) ?     0.062   15.4   0.0  NC_013790.1 1064778 1064857 + hmm     - 0.62  Methanobrevibacter rumi ...
 (28) ?     0.062   15.4   0.0  NC_013790.1  362793  362738 - hmm     - 0.46  Methanobrevibacter rumi ...
 (29) ?     0.064   15.4   0.0  NC_013790.1 2585067 2585000 - hmm     - 0.59  Methanobrevibacter rumi ...
 (30) ?     0.073   15.2   0.0  NC_013790.1 2749938 2749884 - hmm     - 0.47  Methanobrevibacter rumi ...
 (31) ?     0.074   15.2   0.0  NC_013790.1 2749832 2749778 - hmm     - 0.47  Methanobrevibacter rumi ...
 (32) ?      0.13   14.4   0.0  NC_013790.1  361729  361689 - hmm     - 0.56  Methanobrevibacter rumi ...
 (33) ?       1.3   11.2   0.0  NC_013790.1  361439  361393 - hmm     - 0.49  Methanobrevibacter rumi ...
 (34) ?       4.4    9.6   0.0  NC_013790.1 2583867 2583803 - hmm     - 0.49  Methanobrevibacter rumi ...
 (35) ?         6    9.1   0.0  NC_013790.1  546054  546021 - hmm     - 0.68  Methanobrevibacter rumi ...
\end{sreoutput}

Note that this time we only find 35 hits (20 with E-values less than
the inclusion threshold of 0.01) compared to 56 (with 54 less than
0.01) in the original search with the structure-based CM. The
increased sensitivity in the initial CM search exemplifies the
additional power that comes with knowledge of the conserved secondary
structure. Not all RNAs will show such a dramatic difference
though. In fact, tRNA is a particularly strong example of the power of
CMs versus sequence-only based methods like HMMs because about as much
statistical signal is present in their structure as in their primary
sequence. Many structural RNAs contain significantly less information
in their structure than in sequence conservation, but many include 10
bits or more of information in their structure. An additional 10 bits
roughly translates into lowering the expected statistical significance
of homologs detected in database searches by about 3 orders or
magnitude\footnote{For a graph showing the relative amount of
information in structure for sequence for many different types of RNAs
see Figure 1.9 of \citep{Nawrocki09}}.

HMM hit alignments differ slightly from CM hit alignments. Take a look
at the first alignment:

\begin{sreoutput}
>> NC_013790.1  Methanobrevibacter ruminantium M1 chromosome, complete genome
 rank     E-value  score  bias mdl mdl from   mdl to       seq from      seq to       acc trunc   gc
 ----   --------- ------ ----- --- -------- --------    ----------- -----------      ---- ----- ----
  (1) !   3.8e-09   38.2   0.0 hmm        5       67 ..      362022      361960 - .. 0.93     - 0.46

                     ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: CS
   tRNA5-noss      5 auAUaGcgcAGUGGuAGcGCGgcagccUgucaauguggAGGUCCuggGUUCGAUUCccaGUgu 67    
                      uAUaGc+cA+UGG AG+GCG   g cUg ca  +   AGGU  uggGUUCGAUUCcc  U+ 
  NC_013790.1 362022 CUAUAGCUCAAUGGCAGAGCGUUUGGCUGACAUCCAAAAGGUUAUGGGUUCGAUUCCCUUUAG 361960
                     689******************************************************988876 PP
\end{sreoutput}

First, note that the \otext{mdl} field reads \otext{hmm} instead of
\otext{cm}. Also, because HMMs do not model basepairs, there is no
\otext{NC} annotation pointing out negative scoring noncanonical
basepairs.

If you were to look at more hits you may notice that HMM hits are
more likely to not be full length than CM hits are. This hit, for
example, begins at model position 5 and ends at model position 67,
%                                 ^                            ^^
whereas our earlier CM search included a full length alignment from
model position 1 to 72 for the same region of the genome. The profile
%              ^    ^^
HMMs used in Infernal allow local alignments that begin and end at any
position in the model, with an equal score given for all possible
start and end positions (this is the same local alignment strategy
used in HMMER3~\citep{Eddy08}). In contrast, the CM local alignment
strategy used by Infernal encourages global alignments by enforcing a
score penalty for local ones. Partly because CM and HMM alignments
differ in this way, ``truncated'' HMM hits can start and end anywhere
in the target sequence (instead of being identified only with
specialized algorithms only at sequence ends like a CM), and so the
\otext{trunc} field is invalid for HMM hits and will always read
``\otext{-}''.

In \prog{cmscan}, HMM algorithms will be used to compute alignments
for all models in the CM database that contain zero basepairs, and CM
algorithms will be used for all others. Hits found with the HMM will
be annotated like this example above, while CM hits will be annotated
like the earlier CM examples.

By default, Infernal uses HMM algorithms for models with zero
basepairs, mainly because they are more efficient. If you'd like, you
can force the use of CM algorithms for such models with the
\otext{--nohmmonly} option with \prog{cmsearch} or \prog{cmscan}. Using
\otext{--nohmmonly} will encourage more full length hits, but will
cause the program to run a few-fold slower, and also requires that the
CM be calibrated with \prog{cmcalibrate} first.

\subsection{Forcing global CM alignment with the -g option}

By default, \prog{cmsearch}, \prog{cmscan} and \prog{cmalign} allow
local \emph{or} global alignments between the sequence and the
model. This allows an alignment to use what's called a ``local begin''
and start at any match position in the model and to use what's called
a ``local end'' to prematurely exit a section of the model to handle
large insertions or deletions of the model structure and potentially
insert nonhomologous residues into the alignment (these were the
\verb+~+ annotated residues in the \prog{cmscan} and \prog{cmalign}
examples above). Local alignments are permitted by default, because
our internal benchmarks suggest this strategy is more sensitive for
remote homology detection. (An example of a local RNase P alignment is
demonstrated in Figure~\ref{fig:bsu-alignment} on
page~\pageref{fig:bsu-alignment} of this guide.)
However, some users may wish to turn off local alignment. This can be
done using the \otext{-g} option to \prog{cmsearch}, \prog{cmscan},
and \prog{cmalign}. The use of local mode by default in \prog{cmalign}
is new to version 1.1 - all previous versions used global alignment by
default. Importantly, using \otext{-g} does not turn off truncated hit
detection and alignment; to do that use \otext{--notrunc}.

\subsection{Specifying and annotating match positions with cmbuild --hand}

When \prog{cmbuild} constructs a model from an input alignment it must
decide which columns of the input alignment to define as match (also
called ``consensus'') and insert columns. To do this, it first
computes weights for each sequences using an \emph{ad hoc} algorithm
to downweight closely related sequences and upweight distantly related
ones. Then, for each position, the sum of the weights of the sequences
that include a residue at the position is computed. If this sum is at
least half the total number of sequences then the position is defined
as a match position, else it is an insert position. If you're curious
which positions get annotated as match versus insert in your
alignment, use the \otext{-O <f>} option with \prog{cmbuild}, which
will save a annotated version of your input alignment to file
\otext{<f>}, including the weights of each sequence. In this output
alignment, positions with residues in the \otext{\#=GC RF} lines of
the alignment are match positions, and positions with gap characters
are inserts.

Sometimes you may want to override this automated behavior by
specifying which columns be defined as match/insert, especially if
you've painstakingly created your alignment by hand. You can do this
with the \otext{--hand} option\footnote{This is the same option as
  \otext{--rf} from previous versions of Infernal.}. Using this option
you can also define a character to annotate each match position, which
will be propagated to alignment output. In this section we'll go
through an example of using \otext{--hand} with the 5 sequence tRNA
alignment from earlier in the tutorial.

Take a look at the file \otext{tutorial/tRNA5-hand.sto}:

\vspace{1em}
\begin{minipage}{4.0in}
\begin{sreoutput}[xleftmargin=0em]
# STOCKHOLM 1.0

tRNA1             GCGGAUUUAGCUCAGUUGGG.AGAGCGCCAGACUGAAGAUCUGGAGGUCC
tRNA2             UCCGAUAUAGUGUAAC.GGCUAUCACAUCACGCUUUCACCGUGGAGA.CC
tRNA3             UCCGUGAUAGUUUAAU.GGUCAGAAUGGGCGCUUGUCGCGUGCCAGA.UC
tRNA4             GCUCGUAUGGCGCAGU.GGU.AGCGCAGCAGAUUGCAAAUCUGUUGGUCC
tRNA5             GGGCACAUGGCGCAGUUGGU.AGCGCGCUUCCCUUGCAAGGAAGAGGUCA
#=GC SS_cons      <<<<<<<..<<<<.........>>>>.<<<<<.......>>>>>.....<
#=GC RF           [accep]======[=Dloop=]============acd=======[vlp]=

tRNA1             UGUGUUCGAUCCACAGAAUUCGCA
tRNA2             GGGGUUCGACUCCCCGUAUCGGAG
tRNA3             GGGGUUCAAUUCCCCGUCGCGGAG
tRNA4             UUAGUUCGAUCCUGAGUGCGAGCU
tRNA5             UCGGUUCGAUUCCGGUUGCGUCCA
#=GC SS_cons      <<<<.......>>>>>>>>>>>>.
#=GC RF           ====[Tloop]=====[accep]=
//
\end{sreoutput}
\end{minipage}
\begin{minipage}{1.5in}
\includegraphics[scale=0.4]{Figures/trna1-DF6280-hand}
\end{minipage}
\vspace{1em}

This file is the same as \otext{tutorial/tRNA5.sto} except for the two
additional lines beginning with \otext{\#=GC RF}. This RF (reference)
annotation is required for using \otext{--hand}. When \otext{--hand}
is used, any non-gap character in the reference annotation will be
assigned as a match (consensus) position. Importantly, four different
characters are considered gaps: dashes (\otext{-}), underscores
(\otext{\_}), dots (\otext{.}) and tildes (\verb+~+). In this example
alignment, all columns are non-gap characters, so all columns will be
considered match positions.

Different regions of the secondary structure have been marked up using
abbreviations for the names of the regions in the reference
annotation. For example, \otext{acd} annotates the three positions of
the anticodon, and \otext{[vlp]} annotates the so-called variable
loop. I've used \otext{[} and \otext{]} to indicate region boundaries
in some cases. Crucially, I've avoided the use of any gap characters
for positions between named regions which I still want to be
considered match positions, and opted to use \otext{=} (which is not
considered a gap by \prog{cmbuild}) for these positions.

To build the hand-specified model from this alignment, do:

\user{cmbuild --hand tRNA5-hand.cm tutorial/tRNA5-hand.sto}

\begin{sreoutput}
# cmbuild :: covariance model construction from multiple sequence alignments
# INFERNAL 1.1.5 (Sep 2023)
# Copyright (C) 2023 Howard Hughes Medical Institute.
# Freely distributed under the BSD open source license.
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# CM file:                                            tRNA5-hand.cm
# alignment file:                                     tutorial/tRNA5-hand.sto
# use #=GC RF annotation to define consensus columns: yes
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#                                                                      rel entropy
#                                                                      -----------
# idx    name                     nseq eff_nseq   alen  clen  bps bifs    CM   HMM description
# ------ -------------------- -------- -------- ------ ----- ---- ---- ----- ----- -----------
       1 tRNA5-hand                  5     3.59     74    74   21    2 0.763 0.476 
#
# CPU time: 0.09u 0.00s 00:00:00.09 Elapsed: 00:00:00.10
\end{sreoutput}

The output reports that the model now has 74 match (consensus)
positions in the \otext{clen} column. If we had built this model
without specifying \otext{--hand} (as we did earlier in this tutorial)
the resulting model would have had only 72 consensus positions.
(I've annotated the two extra match positions with three gaps in
\otext{tRNA5.hand.sto} as match solely to demonstrate how
\otext{--hand} works, not because I think it's better to model these
positions as matches than inserts.)

Now, let's use this model to search the \emph{M. ruminantium} genome
again. First, the model must be calibrated. To save time, a calibrated
version of the file is in \otext{tutorial/tRNA5-hand.c.cm}. To do the
search:

\user{cmsearch tutorial/tRNA5-hand.c.cm tutorial/mrum-genome.fa}

The results are very similar to the earlier search with the
tRNA model built with default \prog{cmbuild} parameters (though not
identical since the model now has two additional match positions). The
important difference involves the hit alignments. Take a look at the
alignment for hit number 46 as an illustrative example:

\newpage

\begin{widesreoutput}
>> NC_013790.1  Methanobrevibacter ruminantium M1 chromosome, complete genome
 rank     E-value  score  bias mdl mdl from   mdl to       seq from      seq to       acc trunc   gc
 ----   --------- ------ ----- --- -------- --------    ----------- -----------      ---- ----- ----
 (46) !   1.6e-09   43.2   0.0  cm        1       74 []      995344      995263 - .. 0.90    no 0.49

                                 v          v                                                            NC
                     (((((((,,<<<<________._>>>>,<<<<<_______>>>>>,,,........,,<<<<<_______>>>>>))))))): CS
   tRNA5-hand      1 gCcggcaUAGcgcAgUUGGuu.AgcgCgccagccUgucAagcuggAGg........UCCgggGUUCGAUUCcccGugccgGca 74    
                     :::G:CAUAGCG AG  GGU+ A CGCG:CAG:CU +++A:CUG: G+        UC:GGGGUUCGA UCCCC:UG:C:::A
  NC_013790.1 995344 AGAGACAUAGCGAAGC-GGUCaAACGCGGCAGACUCAAGAUCUGUUGAuuaguucuUCAGGGGUUCGAAUCCCCUUGUCUCUA 995263
                     ************9***.8886258***********************9444444445************************** PP
                     [accep]======[=Dloop=.]============acd=======[vl........p]=====[Tloop]=====[accep]= RF
\end{widesreoutput}

The reference annotation from the training alignment to \prog{cmbuild}
has been propagated to the hit as an extra \otext{RF} line at the
bottom of the alignment. All inserts in the alignment are annotated as
\otext{.} columns in the RF annotation. Note that the variable loop
(annotated as \otext{[vlp]} in the training alignment) contains 8
inserted residues. The RF annotation will also be transferred to
multiple alignments created with \prog{cmalign}.

\subsection{Building models from fragmentary sequences}

Sometimes sequences in the \prog{cmbuild} input alignment are sequence
\emph{fragments} with missing data at the 5' or 3' ends. Gaps at the
ends of fragments are treated as missing data by \prog{cmbuild}
instead of as deletions. The definition of fragments can be important
because if sequences have missing data at the end but are not
defined as fragments then the parameters of the resulting model will
incorrectly reflect that deletions are common at the beginning and/or end of the
sequence. Such a model will be less specific for full length matches
to the model. To build the best model from a given alignment, all sequences
with missing data should be defined as fragments.

By default, \prog{cmbuild} will try to guess which sequences are
fragments by defining any sequence which spans less than 0.5 fraction of the
length of the input alignment as a fragment, but this
fraction can be changed to \prog{<x>} with the \prog{--fragthresh <x>}
option. There are two other ways to define fragments. If the
alignment contains RF annotation like the example above, the 
\prog{--hand} and \prog{--fragnrfpos <n>} options can be used in
combination to define a sequence as a fragment if it has more than
\prog{<n>} terminal gaps in nongap RF positions at the 5' \emph{or} 3'
ends. Finally, fragments can be defined manually in the input alignment
using \verb+~+ characters at the ends of fragments as explained more below.

The following modified version of the tRNA alignment will help
illustrate the different fragment definition options. It can be found in
the file \prog{tutorial/tRNA5-hand.frag.sto}:

\begin{sreoutput}
# STOCKHOLM 1.0

tRNA1             ................UGGG.AGAGCGCCAGACUGAAGAUCUGGAGGUCCUGUGUUCGAUCCACAG........
tRNA2             UCCGAUAUAGUGUAAC.GGCUAUCACAUCACGCUUUCACCGUGGAGA.CCGGGGUUCGACUCCCCGUAUCGGAG
tRNA3             UCCGUGAUAGUUUAAU.GGUCAGAAUGGGCGCUUGUCGCGUGCCAGA.UCGGGGUUCAAUUCCCCGUCGCGGAG
tRNA4             ....................................AAAUCUGUUGGUCCUUAGUUCGAU..............
tRNA5             ..GCACAUGGCGCAGUUGGU.AGCGCGCUUCCCUUGCAAGGAAGAGGUCAUCGGUUCGAUUCCGGUUGCGUCCA
#=GC SS_cons      <<<<<<<..<<<<.........>>>>.<<<<<.......>>>>>.....<<<<<.......>>>>>>>>>>>>.
#=GC RF           [accep]======[=Dloop=]============acd=======[vlp]=====[Tloop]=====[accep]=
//
\end{sreoutput}

If this alignment is used as input to \prog{cmbuild} with the command:

\user{cmbuild -O tRNA-hand.O.stk tRNA-hand.frag.cm tutorial/tRNA5-hand.frag.sto}\\

then only the \prog{tRNA4} sequence will be defined as a fragment,
because it spans positions $37$ to $60$ out of a total of $74$
positions ($(60-37+1)/74 = 0.324$) which is less than the $0.5$
threshold. The \prog{tRNA1} sequence spans positions $17$ to $66$ and
so is above the threshold $((66-17+1)/74 = 0.676)$. Due to the \prog{-O}
option above, \prog{cmbuild} saved the annotated alignment file to 
\prog{tRNA-hand.O.stk}. That file, shown below, indicates that \prog{tRNA4} is
defined as a fragment with the \verb+~+ characters used for terminal gaps (also note that
the RF annotation has been changed to the consensus sequence because
\prog{--hand} was not used):

\newpage

\begin{sreoutput}
# STOCKHOLM 1.0
#=GF ID tRNA5-hand.frag
#=GF AU Infernal 1.1.5

#=GS tRNA1 WT 0.89
#=GS tRNA2 WT 1.00
#=GS tRNA3 WT 1.08
#=GS tRNA4 WT 0.86
#=GS tRNA5 WT 1.17

tRNA1        ----------------uGGG-AGAGCGCCAGACUGAAGAUCUGGAGGUCCUGUGUUCGAUCCACAG--------
tRNA2        UCCGAUAUAGUGUAAC.GGCUAUCACAUCACGCUUUCACCGUGGAGA-CCGGGGUUCGACUCCCCGUAUCGGAG
tRNA3        UCCGUGAUAGUUUAAU.GGUCAGAAUGGGCGCUUGUCGCGUGCCAGA-UCGGGGUUCAAUUCCCCGUCGCGGAG
tRNA4        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~AAAUCUGUUGGUCCUUAGUUCGAU~~~~~~~~~~~~~~
tRNA5        --GCACAUGGCGCAGUuGGU-AGCGCGCUUCCCUUGCAAGGAAGAGGUCAUCGGUUCGAUUCCGGUUGCGUCCA
#=GC SS_cons (((((((,,<<<<___._____>>>>,<<<<<_______>>>>>,,,,,<<<<<_______>>>>>))))))):
#=GC RF      UCcgacAUAGugcAAU.GGuuAgcaCGcccgcCUuucAagcgggAGgUCCgggGUUCGAUUCcccGUgucgGAg
//
\end{sreoutput}

Fragments can be defined differently if \prog{--hand} and
\prog{--fragnrfpos 5} options are used, with the following command:

\begin{widesreoutput}
> cmbuild --hand --fragnrfpos 5 -O tRNA-hand.O.2.stk tRNA-hand.frag2.cm tutorial/tRNA5-hand.frag.sto
\end{widesreoutput}
  
In this case, both \prog{tRNA1} and \prog{tRNA4} sequences are defined as
fragments because they include more than $5$ gaps in nongap RF positions
at the 5' or 3' ends (in this case both ends). Again, the fragment
definitions are indicated in the output alignment file
\prog{tRNA-hand.O.2.stk}:

\begin{sreoutput}
# STOCKHOLM 1.0
#=GF ID tRNA5-hand.frag
#=GF AU Infernal 1.1.5

#=GS tRNA1 WT 0.89
#=GS tRNA2 WT 1.02
#=GS tRNA3 WT 1.09
#=GS tRNA4 WT 0.84
#=GS tRNA5 WT 1.16

tRNA1        ~~~~~~~~~~~~~~~~UGGG-AGAGCGCCAGACUGAAGAUCUGGAGGUCCUGUGUUCGAUCCACAG~~~~~~~~
tRNA2        UCCGAUAUAGUGUAAC-GGCUAUCACAUCACGCUUUCACCGUGGAGA-CCGGGGUUCGACUCCCCGUAUCGGAG
tRNA3        UCCGUGAUAGUUUAAU-GGUCAGAAUGGGCGCUUGUCGCGUGCCAGA-UCGGGGUUCAAUUCCCCGUCGCGGAG
tRNA4        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~AAAUCUGUUGGUCCUUAGUUCGAU~~~~~~~~~~~~~~
tRNA5        --GCACAUGGCGCAGUUGGU-AGCGCGCUUCCCUUGCAAGGAAGAGGUCAUCGGUUCGAUUCCGGUUGCGUCCA
#=GC SS_cons (((((((,,<<<<_________>>>>,<<<<<_______>>>>>,,,,,<<<<<_______>>>>>))))))):
#=GC RF      [accep]======[=Dloop=]============acd=======[vlp]=====[Tloop]=====[accep]=
//
\end{sreoutput}

Alternatively, if fragments are annotated in the input alignment, the 
\prog{--fraggiven} option can be used. The
\prog{tutorial/tRNA5-hand.fraggiven.sto} shows an example with the \prog{tRNA1},
\prog{tRNA2} and \prog{tRNA5} sequences all defined as fragments:

\begin{sreoutput}
# STOCKHOLM 1.0

tRNA1             ~~~~~~~~~~~~~~~~UGGG.AGAGCGCCAGACUGAAGAUCUGGAGGUCCUGUGUUCGAUCCACAG~~~~~~~~
tRNA2             UCCGAUAUAGUGUAAC.GGCUAUCACAUCACGCUUUCACCGUGGAGA.CCGGGGUUCGACUCCCCGUAUCGGAG
tRNA3             UCCGUGAUAGUUUAAU.GGUCAGAAUGGGCGCUUGUCGCGUGCCAGA.UCGGGGUUCAAUUCCCCGUCGCGGAG
tRNA4             ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~AAAUCUGUUGGUCCUUAGUUCGAU~~~~~~~~~~~~~~
tRNA5             ~~GCACAUGGCGCAGUUGGU.AGCGCGCUUCCCUUGCAAGGAAGAGGUCAUCGGUUCGAUUCCGGUUGCGUCCA
#=GC SS_cons      <<<<<<<..<<<<.........>>>>.<<<<<.......>>>>>.....<<<<<.......>>>>>>>>>>>>.
#=GC RF           [accep]======[=Dloop=]============acd=======[vlp]=====[Tloop]=====[accep]=
//
\end{sreoutput}

The command to build a model using the fragment definitions in
that alignment is:

\user{cmbuild --fraggiven tRNA-hand.fraggiven.cm tutorial/tRNA5-hand.fraggiven.sto}

To obtain alignments annotated with fragments with \verb+~+
characters, you can either manually edit an existing alignment file,
or generate alignments with the \prog{cmsearch} program using the
\prog{-A} option, or the \prog{cmalign} program with the \prog{--miss}
option which will have sequences inferred to have missing data
annotated as fragments.

\subsection{Parallelizing calibration of large models by splitting
  into partitions}

For large models, like small and large subunit ribosomal RNA (SSU and
LSU rRNA) models, the \prog{cmcalibrate} step can be extremely slow
and require a large amount of memory. Parallelization with
multithreading can help, but memory can still be limiting. For
example, the eukaryotic LSU rRNA model from Rfam (RF02543) requires
about 14Gb of RAM per thread and roughly 30 hours on 12 threads (so
168Gb RAM total) for calibration\footnote{The \prog{--memreq} and
  \prog{--forecast} options to \prog{cmcalibrate} can be used to
  estimate the required memory and time, respectively.}.  MPI
parallelization is also implemented for \prog{cmcalibrate} but not all
systems are set up to use it. If you have access to a compute farm or
cluster, an alternative way to parallelize calibration is
using the \prog{--split} and \prog{--merge} options introduced in
version 1.1.5 which allow parallelization across different CPUs
without using MPI by partitioning the calibration calculations. To use
this method to parallelize calibration of a Cobalamin riboswitch CM
in the file built from the \prog{tutorial/Cobalamin.sto} alignment, we
first build the model with the command:

\user{cmbuild Cobalamin.cm tutorial/Cobalamin.sto}

Then the first step for the parallelized calibration is to create the
script that will execute the calibration calculations for each partition:

\user{cmcalibrate --split --ptot 16 --cfile Cobalamin.cfile Cobalamin.cm}

This will result in the following output:

\begin{sreoutput}
# cmcalibrate :: fit exponential tails for CM E-values
# INFERNAL 1.1.5 (Sep 2023)
# Copyright (C) 2023 Howard Hughes Medical Institute.
# Freely distributed under the BSD open source license.
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# CM file:                                     Cobalamin.cm
# split mode:                                  on
# saving partition commands to file:           Cobalamin.cfile
# total number of partitions:                  16
# number of worker threads:                    0
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# List of commands for each partition saved to file Cobalamin.cfile.
#
# Command to merge per-partition scores after all commands are finished:
# src/cmcalibrate --merge --ptot 16 Cobalamin.cm
#
# CPU time: 0.00u 0.00s 00:00:00.00 Elapsed: 00:00:00.00
[ok]
\end{sreoutput}

And will have created the file \prog{Cobalamin.cfile} which will look
like:

\begin{sreoutput}
cmcalibrate --ptot 16 --part 1 --pfile Cobalamin.cm.calib.1 Cobalamin.cm
cmcalibrate --ptot 16 --part 2 --pfile Cobalamin.cm.calib.2 Cobalamin.cm
cmcalibrate --ptot 16 --part 3 --pfile Cobalamin.cm.calib.3 Cobalamin.cm
cmcalibrate --ptot 16 --part 4 --pfile Cobalamin.cm.calib.4 Cobalamin.cm
cmcalibrate --ptot 16 --part 5 --pfile Cobalamin.cm.calib.5 Cobalamin.cm
cmcalibrate --ptot 16 --part 6 --pfile Cobalamin.cm.calib.6 Cobalamin.cm
cmcalibrate --ptot 16 --part 7 --pfile Cobalamin.cm.calib.7 Cobalamin.cm
cmcalibrate --ptot 16 --part 8 --pfile Cobalamin.cm.calib.8 Cobalamin.cm
cmcalibrate --ptot 16 --part 9 --pfile Cobalamin.cm.calib.9 Cobalamin.cm
cmcalibrate --ptot 16 --part 10 --pfile Cobalamin.cm.calib.10 Cobalamin.cm
cmcalibrate --ptot 16 --part 11 --pfile Cobalamin.cm.calib.11 Cobalamin.cm
cmcalibrate --ptot 16 --part 12 --pfile Cobalamin.cm.calib.12 Cobalamin.cm
cmcalibrate --ptot 16 --part 13 --pfile Cobalamin.cm.calib.13 Cobalamin.cm
cmcalibrate --ptot 16 --part 14 --pfile Cobalamin.cm.calib.14 Cobalamin.cm
cmcalibrate --ptot 16 --part 15 --pfile Cobalamin.cm.calib.15 Cobalamin.cm
cmcalibrate --ptot 16 --part 16 --pfile Cobalamin.cm.calib.16 Cobalamin.cm
\end{sreoutput}

Alternatively, if the \prog{--cbash} option is added to the above
command the \prog{Cobalamin.cfile} would be a more compact bash
shell script with a for loop:

\begin{sreoutput}
# !/bin/bash
for ((i = 1; i <= 16; i++)); do
  cmcalibrate --ptot 16 --part $i --pfile Cobalamin.cm.calib.$i Cobalamin.cm
done
\end{sreoutput}

For this command, we specified that there by 16 partitions with the
option \prog{--ptot 16}, but this number can be as high as
160. This is because by default there are 160 10Kb sequences (1.6Mb
total) searched during the calibration, so the minimal workunit is 
a search of a single 10Kb sequence. If the \prog{-L} option is used to
modify the 1.6Mb total length, then the maximum allowed \prog{<n>}
with \prog{--ptot <n>} will change accordingly (see the
\prog{cmcalibrate} manual page in section 7 of this guide for more
%                                         ^hard-coded, careful
details on \prog{cmcalibrate} options).

If you'd like to use additional options for the calibration, for
example if you want to specify that each partition use only one thread
with the \prog{--cpu 1} option, or run in non-threaded serial mode
with \prog{--cpu 0}, perhaps to limit the required amount of RAM, you
should specify those additional option in the original command above
that includes \prog{--split} and they will be included in the
partition commands in the \prog{Cobalamin.cfile} script.

The next step is to run these 16 commands. If you execute either of
these scripts without modification, the 16 partition commands will run
in serial, one after the other. This would not save any time or reduce
memory requirement versus a standard \prog{cmcalibrate} command like
those shown earlier in this tutorial. So for this approach to make any
sense, you will need to modify the \prog{Cobalamin.cfile} file so that
it will run on a compute cluster. For example, if you use Sun Grid
Engine (SGE) and the program \prog{qsub} to submit jobs to a cluster,
then you might modify the \prog{Cobalamin.cfile} file to something
like:

\begin{sreoutput}
# !/bin/bash
for ((i = 1; i <= 16; i++)); do
  qsub -n Cobalamin.$i -V -o Cobalamin.$i.out -e Cobalamin.$i.err -m n \
  ``cmcalibrate --ptot 16 --part $i --pfile Cobalamin.cm.calib.$i Cobalamin.cm''
done
\end{sreoutput}

The modified script can then be used to submit all 16 jobs. When they
are all finished the final step is to merge the results and update the
CM file with the calibration statistics:

\user{cmcalibrate --merge --ptot 16 Cobalamin.cm}

After this, you will be able to use the \prog{cmsearch} and
\prog{cmscan} with \prog{Cobalamin.cm}. The E-value-related statistics
that are added to the \prog{Cobalamin.cm} file using this split/merge
approach will be identical to those added with standard
\prog{cmcalibrate} (unless the \prog{--seed <n>} option is used to
change the seeding of the random number generator).

\subsubsection{Reading from a stdin pipe using - (dash) as a filename argument}

Generally, Infernal programs read their sequence and/or profile input
from files. Unix power users often find it convenient to string an
incantation of commands together with pipes (indeed, such wizardly
incantations are a point of pride). For example, you might extract a
subset of query sequences from a larger file using a one-liner
combination of scripting commands (perl, awk, whatever). To facilitate
the use of Infernal programs in such incantations, you can almost
always use an argument of '-' (dash) in place of a filename, and the
program will take its input from a standard input pipe instead of
opening a file.

For example, the following three commands are entirely equivalent, and
give essentially identical output:

\user{cmalign tRNA5.cm mrum-tRNAs10.fa}

\user{cat tRNA5.cm | ../src/cmalign - mrum-tRNAs10.fa}

\user{cat mrum-tRNAs10.fa | ../src/cmalign tRNA5.cm -}

Most Easel ``miniapp'' programs share the same ability of pipe-reading.

Because the programs for CM fetching (\prog{cmfetch}) and
sequence fetching (\prog{esl-sfetch}) can fetch any number of profiles
or (sub)sequences by names/accessions given in a list, \emph{and} these
programs can also read these lists from a stdin pipe, you can craft
incantations that generate subsets of queries or targets on the
fly. For example, you can extract and align all hits found by
\prog{cmsearch} with an E-value below the inclusion threshold as 
follows (using the \textbackslash character twice below to split up the final
command onto three lines):

%note: can't use \user{} here because too many special characters
%(believe me I tried). Only difference between \user{} and the way
%I've done it below is that we're not bold, oh well. 
\indent\indent\small\verb+> cmsearch --tblout tRNA5.mrum-genome.tbl tRNA5.c.cm mrum-genome.fa+ \\
\indent\indent\small\verb+> esl-sfetch --index mrum-genome.fa+ \\
\indent\indent\small\verb+> cat tRNA5.mrum-genome.tbl | grep -v ^# | grep ! \ + \\
\indent\indent\small\verb+> | awk '{ printf(``%s/%s-%s %s %s %s\n'', $1, $8, $9, $8, $9, $1); }' \ + \\
\indent\indent\small\verb+> | esl-sfetch -Cf mrum-genome.fa - | cmalign tRNA5.cm - + \\

The first command performed the search using the CM file
\ccode{tRNA5.c.cm} and sequence file \ccode{mrum-genome.fa} (these
were used in the tutorial), and saved tabular output to
\prog{tRNA5.mrum-genome.tbl}.  The second command indexed the genome
file to prepare it for fast (sub)sequence retrieval. In the third
command we've extracted only those lines from
\prog{tRNA5.mrum-genome.tbl} that do not begin with a \prog{\#} (these
are comment lines) and also include a \prog{!} (these are hits that
have E-values below the inclusion threshold) using the first two
\prog{grep} commands. This output was then sent through \prog{awk} to
reformat the tabular output to the ``GDF'' format that
\prog{esl-sfetch} expects: \otext{<newname> <from> <to> <source
seqname>}.  These lines are then piped into \prog{esl-sfetch} (using
the '-' argument) to retrieve each hit (only the subsequence that
comprises each hit -- not the full target sequence). \prog{esl-sfetch}
will output a FASTA file, which is finally being piped into
\prog{cmalign}, again using the '-' argument. The output that is
actually printed to the screen will be a multiple alignment of all the
included tRNA hits.

You can do similar commands piping subsets of CMs. Supposing you have a copy of Rfam in Rfam.cm:

\user{cmfetch --index Rfam.cm} \\ 
\user{cat myqueries.list | cmfetch -f Rfam.cm - | cmsearch - mrum-genome.fa}

This takes a list of query CM names/accessions in
\prog{myqueries.list}, fetches them one by one from Rfam, and does an
cmsearch with each of them against the sequence file
\prog{mrum-genome.fa}. As above, the \prog{cat myqueries.list} part
can be replaced by any suitable incantation that generates a list of
profile names/accessions.

There are three kinds of cases where using '-' is restricted or
doesn't work. A fairly obvious restriction is that you can only use
one '-' per command; you can't do a \prog{cmalign - -} that tries to
read both a CM and sequences through the same stdin
pipe. Second, another case is when an input file must be obligately
associated with additional, separately generated auxiliary files, so
reading data from a single stream using '-' doesn't work because the
auxiliary files aren't present (in this case, using '-' will be
prohibited by the program). An example is \prog{cmscan}, which needs
its \prog{<cmfile>} argument to be associated with four auxiliary
files named \prog{<cmfile>.i1\{mifp\}} that \prog{cmpress} creates,
so \prog{cmscan} does not permit a '-' for its \prog{<cmfile>}
argument. Finally, when a command would require multiple passes over
an input file the command will generally abort after the first pass
if you are trying to read that file through a standard input pipe
(pipes are nonrewindable in general; a few Easel programs
will buffer input streams to make multiple passes possible, but this
is not usually the case). An important example is trying to search a
database that is streamed into \prog{cmsearch}. This is not allowed
because \prog{cmsearch} first reads the entire sequence file to
determine its size (which dictates the filter thresholds that will be
used for the search), then needs to rewind the file before beginning
the search.

In general, Infernal, HMMER and Easel programs document in their man
page whether (and which) command line arguments can be replaced by
'-'.  You can find the HMMER and Easel program man pages in PDF format
in the HMMER user's guide\cite{hmmer3guide}.  You can always check by
trial and error, too. The worst that can happen is a ``Failed to open
file -'' error message, if the program can't read from pipes.