For genome assembly, the version 5 series of MIRA have been reported to work on projects with something like a million Sanger reads (~80 to 100 megabases at 10x coverage), five to ten million 454 Titanium reads (~100 megabases at 20x coverage) and 20 to 40 million Illumina reads (enough for de-novo of a bacterium or a small eukaryote with 76mers to 300mers).

Provided you have the memory, MIRA is expected to work in de-novo mode with

and "normal" coverages, whereas "normal" would be at no more than 50x to 70x for genome projects. Higher coverages will also work, but may create somewhat larger temporary files without heavy parametrisation. Lower coverages (<4x for Sanger, <10x for 454, < 10x for IonTorrent) also need special attention in the parameter settings.

As the complexity of mapping is a lot lower than de-novo, one can basically double (perhaps even triple) the number of reads compared to 'de-novo'. The limiting factor will be the amount of RAM though, and MIRA will also need lots of it if you go into eukaryotes.

The main limiting factor regarding time will be the number of reference sequences (backbones) you are using. MIRA being pedantic during the mapping process, it might be a rather long wait if you have more than 40 megabase of reference sequences.

The default values for MIRA should allow it to work with many EST and RNASeq data sets, sometimes even from non-normalised libraries. For extreme coverage cases however (like, something with a lot of cases at and above 10k coverage), one would perhaps want to resort to data reduction routines before feeding the sequences to MIRA.

On the other hand, recent developments of MIRA were targeted at making de-novo RNASeq assembly of non-normalised libraries liveable, and indeed I now regularly use MIRA for data sets with up to 50 million Illumina 250bp paired-end reads.

Mapping 50 to 60 million sequences against 20k to 30k gene reference sequences is doable, but slow and memory intensive.

	Note
	For mapping against an eukaryotic genome: MIRA does not splice reads in mapping mode. While it is able to more or less handle insertion / deletions of up to 15% of a read length, splices larger than that will not be made: MIRA will still map the reads, but call a lot of SNPs in the parts which should have been spliced.

MIRA is an iterative assembler (it works in several passes) and acts a bit like a child when exploring the world: it explores the assembly space and is specifically parameterised to allow a couple of assembly errors during the first passes. But after each pass some routines (the "parents", if you like) check the result, searching for assembly errors and deduce knowledge about specific assemblies MIRA should not have ventured into. MIRA will then prevent these errors to re-occur in subsequent passes.

As an example, consider the following multiple alignment:

Figure 1.1. How MIRA learns from misassemblies (1). Multiple alignment after 1st pass with an obvious assembly error, notice the clustered columns discrepancies. Two slightly different repeats were assembled together.

These kind of errors can be easily spotted by a human, but are hard to prevent by normal alignment algorithms as sometimes there's only one single base column difference between repeats (and not several as in this example).

MIRA spots these things (even if it's only a single column), tags the base positions in the reads with additional information and then will use that information in subsequent passes. The net effect is shown in the next two figures:

Figure 1.2.  Multiple alignment after last pass where assembly errors from previous passes have been resolved (1st repeat site)

Figure 1.3.  Multiple alignment after last pass where assembly errors from previous passes have been resolved (2nd repeat site)

The ability of MIRA to learn and discern non-identical repeats from each other through column discrepancies is nothing new. Here's the link to a paper from a talk I had at the German Conference on Bioinformatics in 1999: http://www.bioinfo.de/isb/gcb99/talks/chevreux/

I'm sure you'll recognise the basic principle in figures 8 and 9. The slides from the corresponding talk also look very similar to the screen shots above:

Figure 1.4.  Slides presenting the repeat locator at the GCB 99

You can get the talk with these slides here: http://chevreux.org/dkfzold/gcb99/bachvortrag_gcb99.ppt

In the first versions in 1999, the EdIt automatic Sanger sequence editor from Thomas Pfisterer was integrated into MIRA.

Figure 1.5.  Slides presenting the Edit automatic Sanger editor at the GCB 99

The routines used a combination of hypothesis generation/testing together with neural networks (trained on ABI and ALF traces) for signal recognition to discern between base calling errors and true multiple alignment differences. They wen back to the trace data to resolve potential conflicts and eventually recall bases using the additional information gained in a multiple alignment of reads.

Figure 1.6.  Sanger assembly without EdIt automatic editing routines. The bases with blue background are base calling errors.

Figure 1.7.  Sanger assembly with EdIt automatic editing routines. Bases with pink background are corrections made by EdIt after assessing the underlying trace files (SCF files in this case). Bases with blue background are base calling errors where the evidence in the trace files did not show enough evidence to allow an editing correction.

With the introduction of 454 and Illumina reads, MIRA also got in 2007 specialised editors to search and correct for typical 454 and Illumina sequencing problems, like the homopolymer run over-/undercalls or the Illumina GGCxG motif problem. These editors are now integrated into MIRA itself and are not part of EdIt anymore.

While not being paramount to the assembly quality, both editors provide additional layers of safety for the MIRA learning algorithm to discern non-perfect repeats even on a single base discrepancy. Furthermore, the multiple alignments generated by these two editors are way more pleasant to look at (or automatically analyse) than the ones containing all kind of gaps, insertions, deletions etc.pp.

Figure 1.8.  454 assembly without 454 automatic editing routines

Figure 1.9.  454 assembly with 454 automatic editing routines

A very useful feature for finishing are kmer (hash) frequency tags which MIRA sets in the assembly. Provided your finishing editor understands those tags (gap4, gap5 and consed are fine but there may be others), they'll give you precious insight where you might want to be cautious when joining to contigs or where you would need to perform some primer walking. MIRA colourises the assembly with the hash frequency (HAF) tags to show repetitiveness.

You will need to read about the HAF tags in the reference manual, but in a nutshell: the HAF5, HAF6 and HAF7 tags tell you potentially have repetitive to very repetitive read areas in the genome, while HAF2 tags will tell you that these areas in the genome have not been covered as well as they should have been.

As an example, the following figure shows the coverage of a contig.

Figure 1.10.  Coverage of a contig.

The question is now: why did MIRA stop building this contig on the left end (left oval) and why on the right end (right oval).

Looking at the HAF tags in the contig, the answer becomes quickly clear: the left contig end has HAF5 tags in the reads (shown in bright red in the following figure). This tells you that MIRA stopped because it probably could not unambiguously continue building this contig. Indeed, if you BLAST the sequence at the NCBI, you will find out that this is an rRNA area of a bacterium, of which bacteria normally have several copies in the genome:

Figure 1.11.  HAF5 tags (reads shown with red background) covering a contig end show repetitiveness as reason for stopping a contig build.

The right end of the same contig however ends in HAF3 tags (normal coverage, bright green in the next figure) and even HAF2 tags (below average coverage, pale green in the next image). This tells you MIRA stopped building the contig at this place simply because there were no more reads to continue. This is a perfect target for primer walking if you want to finish a genome.

Figure 1.12.  HAF2 tags covering a contig end show that no more reads were available for assembly at this position.

Many people combine Sanger & 454 or 454 & Illumina to improve the sequencing quality of their project through two (or more) sequencing technologies. To reduce time spent in finishing, MIRA automatically tags those bases in a consensus of a hybrid assembly where reads from different sequencing technologies severely contradict each other.

The following example shows a hybrid 454 / Illumina assembly where reads from 454 (highlighted read names in following figure) were not sure whether to have one or two "G" at a certain position. The consensus algorithm would have chosen "two Gs" for 454, obviously a wrong decision as all Illumina reads at the same spot (the reads which are not highlighted) show only one "G" for the given position. While MIRA chose to believe Illumina in this case, it tagged the position anyway in case someone chooses to check these kind of things.

Figure 1.13.  A "STMS" tag (Sequencing Technology Mismatch Solved, the black square base in the consensus) showing a potentially difficult decision in a hybrid 454 / Illumina de-novo assembly.

While MIRA is not really suited to assemble PacBio or Oxford Nanopore reads de-novo, it is very well suited to polish the results of such assemblies with Illumina reads in mapping projects.

TODO: example

Quality control is paramount when you do mutation analysis for biologists: I know they'll be on my doorstep the very next minute they found out one of the SNPs in the resequencing data wasn't a SNP, but a sequencing artefact. And I can understand them: why should they invest -- per SNP -- hours in the wet lab if I can invest a couple of minutes to get them data false negative rates (and false discovery rates) way below 1%? So, finishing and quality control for any mapping project is a must.

Both gap4 and consed start to have a couple of problems when projects have millions of reads: you need lots of RAM and scrolling around the assembly gets a test to your patience. Still, these two assembly finishing programs are amongst the better ones out there, although gap5 starts to quickly arrive in a state in which it allows itself to substitute to gap4.

So, MIRA reduces the number of reads in Illumina mapping projects without sacrificing information on coverage. The principle is pretty simple: for 100% matching reads, MIRA tracks coverage of every reference base and creates long synthetic, coverage equivalent reads (CERs) in exchange for the Illumina reads. Reads that do not match 100% are kept as own entities, so that no information gets lost. The following figure illustrates this:

Figure 1.14.  Coverage equivalent reads (CERs) explained.

Left side of the figure: a conventional mapping with eleven reads of size 4 against a consensus (in uppercase). The inversed base in the lowest read depicts a sequencing error.

Right side of the figure: the same situation, but with coverage equivalent reads (CERs). Note that there are less reads, but no information is lost: the coverage of each reference base is equivalent to the left side of the figure and reads with differences to the reference are still present.

This strategy is very effective in reducing the size of a project. As an example, in a mapping project with 9 million Illumina 36mers, MIRA created a project with 1.7m reads: 700k CER reads representing ~8 million 100% matching Illumina reads, and it kept ~950k mapped reads as they had ≥ mismatch (be it sequencing error or true SNP) to the reference. A reduction of 80%, and numbers for mapping projects with Illumina 100bp reads are in a similar range.

Also, mutations of the resequenced strain now really stand out in the assembly viewer as the following figure shows:

Figure 1.15.  Coverage equivalent reads let SNPs become very visible in assembly viewers

Want to assemble two or several very closely related genomes without reference, but finding SNPs or differences between them?

Tired of looking at some text output from mapping programs and guessing whether a SNP is really a SNP or just some random junk?

MIRA tags all SNPs (and other features like missing coverage etc.) it finds so that -- when using a finishing viewer like gap4 or consed -- one can quickly jump from tag to tag and perform quality control. This works both in de-novo assembly and in mapping assembly, all MIRA needs is the information which read comes from which strain.

The following figure shows a mapping assembly of Illumina 36mers against a bacterial reference sequence, where a mutant has an indel position in an gene:

Figure 1.16.  "SROc" tag (Snp inteR Organism on Consensus) showing a SNP position in a Illumina mapping assembly.

Other interesting places like deletions of whole genome parts are also directly tagged by MIRA and noted in diverse result files (and searchable in assembly viewers):

Figure 1.17.  "MCVc" tag (Missing CoVerage in Consensus, dark red stretch in figure) showing a genome deletion in Illumina mapping assembly.

	Note
	For bacteria -- and if you use annotated GenBank files as reference sequence -- MIRA will also output some nice lists directly usable (in Excel) by biologists, telling them which gene was affected by what kind of SNP, whether it changes the protein, the original and the mutated protein sequence etc.pp.

There are two kind of versions for MIRA that can be compiled form source files: production and development.

Production versions are from the stable branch of the source code. These versions are available for download from SourceForge.

Development versions are from the development branch of the source tree. These are also made available to the public and should be compiled by users who want to test out new functionality or to track down bugs or errors that might arise at a given location. Release candidates (rc) also fall into the development versions: they are usually the last versions of a given development branch before being folded back into the production branch.

© 1997-2000 Deutsches Krebsforschungszentrum Heidelberg -- Dept. of Molecular Biophysics and Bastien Chevreux (for MIRA) and Thomas Pfisterer (for EdIt)

© 2001 (and later) Bastien Chevreux.

All rights reserved.

MIRA uses the excellent Expat library to parse XML files. Expat is Copyright © 1998, 1999, 2000 Thai Open Source Software Center Ltd and Clark Cooper as well as Copyright © 2001, 2002 Expat maintainers.

See http://www.libexpat.org/ and http://sourceforge.net/projects/expat/ for more information on Expat.

Please use these citations:

If you find this software useful, please send the author a postcard. If postcards are not available, a treasure chest full of Spanish doubloons, gold and jewellery will do nicely, thank you.

Source packages are usually named

mira-miraversion.tar.bz2

Note

For miraversion, the stable versions of MIRA with the general public as audience usually have a version number in three parts, like 3.0.5, sometimes also followed by some postfix like in 3.2.0rc1 to denote release candidate 1 of the 3.2.0 version of MIRA. On very rare occasions, stable versions of MIRA can have four part like in, e.g., 3.4.0_1: these versions create identical binaries to their parent version (in this case: 3.4.0) and just contains fixes to the source build machinery or documentation. The version string sometimes can have a different format: sometext_somenumber+gsomehexnumber like in, e.g., ft_fastercontig_10+g4a27c91. These versions of MIRA are snapshots from the development tree of MIRA and usually contain new functionality which may not be as well tested as the rest of MIRA, hence contains more checks and more debugging output to catch potential errors

Precompiled binary packages are named in the following way:

mira_miraversion_OS-and-binarytype.tar.bz2

	Note
	The OS-and-binarytype part defines for which operating system and which processor class the package is destined. E.g., linux-gnu_x86_64_static contains static binaries for Linux running a X86 64 bit processor (most Intel/AMD nowadays).

SourceForge: https://sourceforge.net/projects/mira-assembler/files/MIRA/

There you will normally find two directories (stable and development), both with a couple of precompiled binaries -- usually for Linux and Mac OSX -- or the source package for compiling yourself.

Examples for packages at SourceForge:

GitHub main page: https://github.com/bachev/mira

GitHub releases (go here for binary packages): https://github.com/bachev/mira/releases

GitHub contains the live source of MIRA as well as source packages and binary packages for selected versions. Two main branches are generally interesting to the public: "master" and "development".

To get a clone of the repository on on GitHub, do:

git clone https://github.com/bachev/mira.git
cd mira
./bootstrap.sh

Once bootstrap has run through (see below for necessary prerequisites), the system is ready for a ./configure call.

Compiling the 5.x series of MIRA needs a C++14 compatible tool chain, i.e., systems starting from 2013/2014 should be OK. The requisites for compiling MIRA are:

For building the documentation, additional prerequisites are:

When building from a GitHub checkout, additional prerequisites are a couple of additional tools:

If you're lucky, after fetching the source (and eventually calling ./bootstrap.sh if you worked from a GitHub checkout), all you need to do is basically:

arcadia:/path/to/mira-5.0.0$ ./configure
arcadia:/path/to/mira-5.0.0$ make
arcadia:/path/to/mira-5.0.0$ make install

This should install the following programs:

Should the ./configure step fail for some reason or another, you should get a message telling you at which step this happens and and either install missing packages or tell configure where it should search the packages it did not find. See also next section or the system specific walkthroughs in this chapter.

MIRA understands all standard autoconf configure arguments, the most important one probably being --prefix= to tell where you want to have MIRA installed. Please consult ./configure --help to get a full list of currently supported arguments.

You will need to install a couple of tools and libraries before compiling MIRA. Here's the recipe:

sudo apt-get install gcc make flex
sudo apt-get install libexpat1-dev libboost-all-dev

Once this is done, you can unpack and compile MIRA. For a dynamically linked version, use:

tar xvjf mira-5.0.0.tar.bz2
cd mira-5.0.0
./configure
make && make install

For a statically linked version, just change the configure line from above into:

./configure --enable-mirastatic

In case you also want to build documentation yourself, you will need to install this in addition:

sudo apt-get install xsltproc docbook-xsl dblatex

And building the documentation in PDF and HTML format is done like this:

make docs

	Note
	People working on git checkouts of the MIRA source code will obviously need some more tools. Get them with this: sudo apt-get install automake libtool xutils-dev

	Note
	This guide for OpenSUSE is outdated, but may give you clues on how to perform this on current systems.

You will need to install a couple of tools and libraries before compiling MIRA. Here's the recipe:

sudo zypper install gcc-c++ boost-devel
sudo zypper install flex libexpat-devel zlib-devel

Once this is done, you can unpack and compile MIRA. For a dynamically linked version, use:

tar xvjf mira-5.0.0.tar.bz2
cd mira-5.0.0
./configure
make && make install

In case you also want to build documentation yourself, you will need this in addition:

sudo zypper install docbook-xsl-stylesheets dblatex

	Note
	People working on git checkouts of the MIRA source code will obviously need some more tools. Get them with this: sudo zypper install automake libtool xutils-dev

	Note
	This guide for Fedora is outdated, but may give you clues on how to perform this on current systems.

You will need to install a couple of tools and libraries before compiling MIRA. Here's the recipe:

sudo yum -y install gcc-c++ boost-devel
sudo yum install flex expat-devel vim-common zlib-devel

Once this is done, you can unpack and compile MIRA. For a dynamically linked version, use:

tar xvjf mira-5.0.0.tar.bz2
cd mira-5.0.0
./configure
make && make install

In case you also want to build documentation yourself, you will need this in addition:

sudo yum -y install docbook-xsl dblatex

	Note
	People working on git checkouts of the MIRA source code will obviously need some more tools. Get them with this: sudo yum -y install automake libtool xorg-x1-util-devel

These instructions are for OSX 10.14 (Mojave) and use Homebrew. There are other ways to do this (e.g., see the "compile everything from scratch"), but they are definetly more painful.

First, you will need to install the Apple XCode command line tools. This may be as easy as simply open a terminal and typing:

xcode-select --install

However, on the machine I have, I got weird compilation errors and found https://github.com/frida/frida/issues/338#issuecomment-426777849 which told me to do this for things to work (don't ask me why, but this resolved my problems):

cd /Library/Developer/CommandLineTools/Packages/
open macOS_SDK_headers_for_macOS_10.14.pkg

Then, if you do not already have it, install Homebrew. See https://brew.sh/.

Then go on and install with homebrew the LLVM clang compiler suite and a couple of other things needed:

brew install llvm boost flex expat

To compile from a git checkout, you will also need these:

brew install autoconf automake libtool git

To build the HTML documentation, this:

brew install docbook-xsl

And to build the PDF documentation, we need MacTex (huge download, around 1 GiB, 6 GiB on disk) and the dblatex package:

brew cask install mactex
pip install dblatex

Finally, if you want to build small, rudimentary man pages:

brew install help2man

Now unpack MIRA, configure it and compile:

tar xvf mira-5.0.0.tar.bz2
cd mira-5.0.0
./configure --with-brew --enable-debug
make -j 2

That's it for the dynamic version.

For building an almost static version, we need some trickery: configure with the --enable-mirastatic argument), then ask the build system to create a distribution. This will launch a toolchain at the end of which you will get a bundle which contains everything needed to run, both on OSX and other Unix systems.

./configure --with-brew --enable-mirastatic --enable-debug
make
make distrib

This lets you build a self-contained static MIRA binary. The only prerequisite here is that you have a working gcc with the minimum version described above. Please download all necessary files (expat, flex, etc.pp) and then simply follow the script below. The only things that you will want to change are the path used and, maybe, the name of some packages in case they were bumped up a version or revision.

Contributed by Sven Klages.

## whatever path is appropriate
cd /home/gls/SvenTemp/install

## expat
tar zxvf expat-2.0.1.tar.gz
cd expat-2.0.1
./configure --prefix=/home/gls/SvenTemp/expat
make && make install

## flex
cd /home/gls/SvenTemp/install
tar zxvf flex-2.5.35.tar.gz
cd flex-2.5.35
./configure --prefix=/home/gls/SvenTemp/flex
make && make install
cd /home/gls/SvenTemp/flex/bin
ln -s flex flex++
export PATH=/home/gls/SvenTemp/flex/bin:$PATH

## boost
cd /home/gls/SvenTemp/install
tar zxvf boost_1_48_0.tar.gz
cd boost_1_48_0
./bootstrap.sh --prefix=/home/gls/SvenTemp/boost
./b2 install

## mira itself
export CXXFLAGS="-I/home/gls/SvenTemp/flex/include"

cd /home/gls/SvenTemp/install
tar zxvf mira-3.4.0.1.tar.gz
cd mira-3.4.0.1
./configure --prefix=/home/gls/SvenTemp/mira \
--with-boost=/home/gls/SvenTemp/boost \
--with-expat=/home/gls/SvenTemp/expat \
--enable-mirastatic
make && make install

In case you do not want a static binary of MIRA, but a dynamically linked version, the following script by Robert Bruccoleri will give you an idea on how to do this.

Note that he, having root rights, puts all additional software in /usr/local, and in particular, he keeps updated versions of Boost and Flex there.

#!/bin/sh -x

make distclean
oze=`find . -name "*.o" -print`
if [[ -n "$oze" ]]
then
   echo "Not clean."
   exit 1

fi

export prefix=${BUILD_PREFIX:-/usr/local}
export LDFLAGS="-Wl,-rpath,$prefix/lib"

./configure --prefix=$prefix \
           --enable-debug=yes \
           --enable-mirastatic=no \
           --with-boost-libdir=$prefix/lib \
           --enable-optimisations \
           --enable-boundtracking=yes \
           --enable-bugtracking=yes \
           --enable-extendedbugtracking=no
make
make install

Contributed by Thomas Vaughan

The system flex (/usr/bin/flex) is too old, but the devel/flex package from a recent pkgsrc works fine. BSD make doesn't like one of the lines in src/progs/Makefile, so use GNU make instead (available from pkgsrc as devel/gmake). Other relevant pkgsrc packages: devel/boost-libs, devel/boost-headers and textproc/expat. The configure script has to be told about these pkgsrc prerequisites (they are usually rooted at /usr/pkg but other locations are possible):

FLEX=/usr/pkg/bin/flex ./configure --with-expat=/usr/pkg --with-boost=/usr/pkg

If attempting to build a pkgsrc package of MIRA, note that the LDFLAGS passed by the pkgsrc mk files don't remove the need for the --with-boost option. The configure script complains about flex being too old, but this is harmless because it honours the $FLEX variable when writing out makefiles.

Depending on options/paramaters, the MIRA/mirabait binary may need to load some additional data during the run. By default this data will always be searched at this location: LOCATION_OF_BINARY/../share/mira/...

That is: If the binary is, e.g., /opt/mira5/bin/mira with a softlink pointing from /usr/local/bin/mira -> /opt/mira5/bin/mira (because, e.g., /usr/local/bin may be by default in your PATH variable), then the additional data will be searched in /opt/mira5/share/mira/... and NOT in /usr/local/share/mira/....

	Note
	In short: since MIRA 4.9.6, moving the binary is not enough anymore. Take care to have the share directory in the right place, i.e., adjacent to the directory the MIRA binary lives in.

You need to create (or get from your sequencing provider) the sequencing data in any supported file format. Amongst these, FASTQ and FASTA + FASTA-quality will be the most common, although the latter is well on the way out nowadays. The following walkthrough uses what most people nowadays get: FASTQ.

Create a new project directory (e.g. myProject) and a subdirectory of this which will hold the sequencing data (e.g. data).

arcadia:/path/to mkdir myProject
arcadia:/path/to cd myProject
arcadia:/path/to/myProject$ mkdir data

Put the FASTQ data into that data directory so that it now looks perhaps like this:

arcadia:/path/to/myProject$ ls -l data
-rw-r--r-- 1 bach users 263985896 2008-03-28 21:49 bchocwt_lane6.solexa.fastq

	Note
	I completely made up the file names above. You can name them anyway you want. And you can have them live anywhere on the hard-disk, you do not need to put them in this data directory. It's just the way I do it ... and it's where the example manifest files a bit further down in this chapter will look for the data files.

We're almost finished with the setup. As I like to have things neatly separated, I always create a directory called assemblies which will hold my assemblies (or different trials) together. Let's quickly do that:

arcadia:/path/to/myProject$ mkdir assemblies
arcadia:/path/to/myProject$ mkdir assemblies/1sttrial
arcadia:/path/to/myProject$ cd assemblies/1sttrial

A manifest file is a configuration file for MIRA which tells it what type of assembly it should do and which data it should load. In this case we'll make a simple assembly of a genome with unpaired Illumina data

# Example for a manifest describing a genome de-novo assembly with
# unpaired Illumina data

# First part: defining some basic things
# In this example, we just give a name to the assembly
#  and tell MIRA it should map a genome in accurate mode

project = MyFirstAssembly
job = genome,denovo,accurate

# The second part defines the sequencing data MIRA should load and assemble
# The data is logically divided into "readgroups"

# here comes the unpaired Illumina data

readgroup = SomeUnpairedIlluminaReadsIGotFromTheLab
data = ../../data/bchocwt_lane6.solexa.fastq
technology = solexa

	Note
	Please look up the parameters of the manifest file in the main manual or the example manifest files in the following section. The ones above basically say: make an accurate denovo assembly of unpaired Illumina reads.

Starting the assembly is now just a matter of a simple command line:

arcadia:/path/to/myProject/assemblies/1sttrial$ mira manifest.conf >&log_assembly.txt

For this example - if you followed the walk-through on how to prepare the data - everything you might want to adapt in the first time are the following thing in the manifest file: options:

Of course, you are free to change any option via the extended parameters, but this is the topic of another part of this manual.

Well, we've seen that already in the section above, but here it is again ... but this time with 454 data.

# Example for a manifest describing a denovo assembly with
# unpaired 454 data

# First part: defining some basic things
# In this example, we just give a name to the assembly
#  and tell MIRA it should map a genome in accurate mode

project = MyFirstAssembly
job = genome,denovo,accurate

# The second part defines the sequencing data MIRA should load and assemble
# The data is logically divided into "readgroups"

# here's the 454 data

readgroup = SomeUnpaired454ReadsIGotFromTheLab
data = ../../data/some454data.fastq
technology = 454

Hybrid de-novo assemblies follow the general manifest scheme: tell what you want in the first part, then simply add as separate readgroup the information MIRA needs to know to find the data and off you go. Just for laughs, here's a manifest for 454 shotgun with Illumina shotgun

# Example for a manifest describing a denovo assembly with
# shotgun 454 and shotgun Illumina data

# First part: defining some basic things
# In this example, we just give a name to the assembly
#  and tell MIRA it should map a genome in accurate mode

project = MyFirstAssembly
job = genome,denovo,accurate

# The second part defines the sequencing data MIRA should load and assemble
# The data is logically divided into "readgroups"

# now the shotgun 454 data
readgroup = DataForShotgun454
data = ../../data/project454data.fastq
technology = 454

# now the shotgun Illumina data

readgroup = DataForShotgunIllumina
data = ../../data/someillumina.fastq
technology = solexa

When using paired-end data, you should know

In case you do not know one (or any) of the above, don't panic! MIRA is able to estimate the needed values during the assembly if you tell it to.

The following manifest shows you the most laziest way to define a paired data set by simply adding autopairing as keyword to a readgroup (using Illumina just as example):

# Example for a lazy manifest describing a denovo assembly with
# one library of paired reads

# First part: defining some basic things
# In this example, we just give a name to the assembly
#  and tell MIRA it should map a genome in accurate mode

project = MyFirstAssembly
job = genome,denovo,accurate

# The second part defines the sequencing data MIRA should load and assemble
# The data is logically divided into "readgroups"

# now the Illumina paired-end data

readgroup = DataIlluminaPairedLib
autopairing
data = ../../data/project_1.fastq ../../data/project_2.fastq
technology = solexa

If you know the orientation of the reads and/or the library size, you can tell this MIRA the following way (just showing the readgroup definition here):

readgroup = DataIlluminaPairedEnd500Lib
data = ../../data/project_1.fastq ../../data/project_2.fastq
technology = solexa
template_size = 250 750
segment_placement = ---> <---

In cases you are not 100% sure about, e.g., the size of the DNA template, you can also give a (generous) expected range and then tell MIRA to automatically refine this range during the assembly based on real, observed distances of read pairs. Do this with autorefine modifier like this:

template_size = 50 1000 autorefine

The following manifest file is an example for assembling with several different libraries from different technologies. Do not forget you can use autopairing or autorefine :-)

# Example for a manifest describing a denovo assembly with
# several kinds of sequencing libraries

# First part: defining some basic things
# In this example, we just give a name to the assembly
#  and tell MIRA it should map a genome in accurate mode

project = MyFirstAssembly
job = genome,denovo,accurate

# The second part defines the sequencing data MIRA should load and assemble
# The data is logically divided into "readgroups"

# now the Illumina paired-end data

readgroup = DataIlluminaForPairedEnd500bpLib
data = ../../data/project500bp-1.fastq ../../data/project500bp-2.fastq
technology = solexa
strain = bchoc_se1
template_size = 250 750
segment_placement = ---> <---

# now the Illumina mate-pair data

readgroup = DataIlluminaForMatePair3kbLib
data = ../../data/project3kb-1.fastq ../../data/project3kb-2.fastq
technology = solexa
strain = bchoc_se1
template_size = 2500 3500
segment_placement = <--- --->

# some Sanger data (6kb library)

readgroup = DataForSanger6kbLib
data = ../../data/sangerdata.fastq
technology = sanger
template_size = 5500 6500
segment_placement = ---> <---

# some 454 data

readgroup = DataFo454Pairs
data = ../../data/454data.fastq
technology = 454
template_size = 8000 1200
segment_placement = 2---> 1--->

# some Ion Torrent data

readgroup = DataFoIonPairs
data = ../../data/iondata.fastq
technology = iontor
template_size = 1000 300
segment_placement = 2---> 1--->

MIRA will make use of ancillary information present in the manifest file. One of these is the information to which strain (or organism or cell line etc.pp) the generated data belongs.

You just need to tell in the manifest file which data comes from which strain. Let's assume that in the example from above, the "lane6" data were from a first mutant named bchoc_se1 and the "lane7" data were from a second mutant named bchoc_se2. Here's the manifest file you would write then:

# Example for a manifest describing a de-novo assembly with
# unpaired Illumina data, but from multiple strains

# First part: defining some basic things
# In this example, we just give a name to the assembly
#  and tell MIRA it should map a genome in accurate mode

project = MyFirstAssembly
job = genome,denovo,accurate

# The second part defines the sequencing data MIRA should load and assemble
# The data is logically divided into "readgroups"

# now the Illumina data

readgroup = DataForSE1
data = ../../data/bchocse_lane6.solexa.fastq
technology = solexa
strain = bchoc_se1

readgroup = DataForSE2
data = ../../data/bchocse_lane7.solexa.fastq
technology = solexa
strain = bchoc_se2

	Note
	While assembling de-novo with multiple strains is possible, the interpretation of results may become a bit daunting in some cases. For many scenarios it might therefore be preferable to successively use the data sets in own assemblies or mappings.

This strain information for each readgroup is really the only change you need to perform to tell MIRA everything it needs for handling strains.

You need to create (or get from your sequencing provider) the sequencing data in either FASTQ or FASTA + FASTA quality format. The following walkthrough uses what most people nowadays get: FASTQ.

Create a new project directory (e.g. myProject) and a subdirectory of this which will hold the sequencing data (e.g. data).

arcadia:/path/to mkdir myProject
arcadia:/path/to cd myProject
arcadia:/path/to/myProject$ mkdir data

Put the FASTQ data into that data directory so that it now looks perhaps like this:

arcadia:/path/to/myProject$ ls -l data
-rw-r--r-- 1 bach users 263985896 2008-03-28 21:49 bchocse_lane6.solexa.fastq
-rw-r--r-- 1 bach users 264823645 2008-03-28 21:51 bchocse_lane7.solexa.fastq

The reference sequence (the backbone) can be in a number of different formats: GFF3, GenBank (.gb, .gbk, .gbf, .gbff), MAF, CAF, FASTA. The first three have the advantage of being able to carry additional information like, e.g., annotation. In this example, we will use a GFF3 file like the ones one can download from the NCBI.

	Note
	TODO: Write why GFF3 is better and where to get them at the NCBI.

So, let's assume that our wild type strain is in the following file: NC_someNCBInumber.gff3.

You do not need to copy the reference sequence to your directory, but I normally copy also the reference file into the directory with my data as I want to have, at the end of my work, a nice little self-sufficient directory which I can archive away and still be sure that in 10 years time I have all data I need together.

arcadia:/path/to/myProject$ cp /somewhere/NC_someNCBInumber.gff3 data
arcadia:/path/to/myProject$ ls -l data
-rw-r--r-- 1 bach users   6543511 2008-04-08 23:53 NC_someNCBInumber.gff3
-rw-r--r-- 1 bach users 263985896 2008-03-28 21:49 bchocse_lane6.solexa.fastq
-rw-r--r-- 1 bach users 264823645 2008-03-28 21:51 bchocse_lane7.solexa.fastq

We're almost finished with the setup. As I like to have things neatly separated, I always create a directory called assemblies which will hold my assemblies (or different trials) together. Let's quickly do that:

arcadia:/path/to/myProject$ mkdir assemblies
arcadia:/path/to/myProject$ mkdir assemblies/1sttrial
arcadia:/path/to/myProject$ cd assemblies/1sttrial

A manifest file is a configuration file for MIRA which tells it what type of assembly it should do and which data it should load. In this case we have unpaired sequencing data which we want to map to a reference sequence, the manifest file for that is pretty simple:

# Example for a manifest describing a mapping assembly with
# unpaired Illumina data

# First part: defining some basic things
# In this example, we just give a name to the assembly
#  and tell MIRA it should map a genome in accurate mode

project = MyFirstMapping
job = genome,mapping,accurate

# The second part defines the sequencing data MIRA should load and assemble
# The data is logically divided into "readgroups"

# first, the reference sequence
readgroup
is_reference
data = ../../data/NC_someNCBInumber.gff3
strain = bchoc_wt

# now the Illumina data

readgroup = SomeUnpairedIlluminaReadsIGotFromTheLab
data = ../../data/*fastq
technology = solexa
strain = bchoc_se

Note

Please look up the parameters of the manifest file in the main manual or the example manifest files in the following section. The ones above basically say: make an accurate mapping of Solexa reads against a genome; in one pass; the name of the backbone strain is 'bchoc_wt'; the data with the backbone sequence (and maybe annotations) is in a specified GFF3 file; for Solexa data: assign default strain names for reads which have not loaded ancillary data with strain info and that default strain name should be 'bchoc_se'.

Starting the assembly is now just a matter of a simple command line:

arcadia:/path/to/myProject/assemblies/1sttrial$ mira manifest.conf >&log_assembly.txt

For this example - if you followed the walk-through on how to prepare the data - everything you might want to adapt in the first time are the following thing in the manifest file: options:

Of course, you are free to change any option via the extended parameters, but this is the topic of another part of this manual.

Well, we've seen that already in the section above, but here it is again ... this time with Ion Torrent data though.

# Example for a manifest describing a mapping assembly with
# unpaired Ion data

# First part: defining some basic things
# In this example, we just give a name to the assembly
#  and tell MIRA it should map a genome in accurate mode

project = MyFirstMapping
job = genome,mapping,accurate

# The second part defines the sequencing data MIRA should load and assemble
# The data is logically divided into "readgroups"

# first, the reference sequence
readgroup
is_reference
data = ../../data/NC_someNCBInumber.gff3
strain = bchoc_wt

# now the Ion Torrent data

readgroup = SomeUnpairedIonReadsIGotFromTheLab
data = ../../data/someiondata.fastq
technology = iontor
strain = bchoc_se

When using paired-end data in mapping, you must decide whether you want

The Illumina pipeline generally normally gives you two files for paired-end data: a project-1.fastq and project-2.fastq. The first file containing the first read of a read-pair, the second file the second read. Depending on the preprocessing pipeline of your sequencing provider, the names of the reads are either the very same in both files or already have a /1 or /2 appended. Also, your sequencing provider may give you one big file where the reads from both ends are present.

	Note
	MIRA can read all FASTQ variants produced by various Illumina pipelines, be they with or without the /1 and /2 already appended to the names. You generally do not need to do any name mangling before feeding the data to MIRA. However, MIRA will shell out a warning if read names are longer than 40 characters.

When using paired-end data, you should know

In case you do not know one (or any) of the above, don't panic! MIRA is able to estimate the needed values during the assembly if you tell it to.

The following manifest shows you the most laziest way to define a paired data set by simply adding autopairing as keyword to a readgroup (using Illumina just as example):

# Example for a lazy manifest describing a mapping assembly with
# one library of paired reads

# First part: defining some basic things
# In this example, we just give a name to the assembly
#  and tell MIRA it should map a genome in accurate mode

project = MyFirstMapping
job = genome,mapping,accurate

# The second part defines the sequencing data MIRA should load and assemble
# The data is logically divided into "readgroups"

# first the reference sequence
readgroup
is_reference
data = ../../data/NC_someNCBInumber.gff3
strain = bchoc_wt

# now the Illumina paired-end data

readgroup = DataIlluminaPairedLib
autopairing
data = ../../data/project_1.fastq ../../data/project_2.fastq
technology = solexa
strain = bchoc_se1

See? Wasn't hard and it did not hurt, did it? One just needs to tell MIRA it should expect paired reads via the autopairing keyword and that is everything you need.

If you are lazy, or want to efficiently include multiple FASTQ files, you can use regular expressions in the line which defines the data to load. E.g., if there were only two FASTQ files in the directory named _1 and _2, you can write this:

data = ../../data/project_?.fastq

Notice the question mark in the above line? Read more about how to use this and other shortcuts in Section 9.3.3.2: “ Defining data files to load ”

If you know the orientation of the reads and/or the library size, you can tell this MIRA the following way (just showing the readgroup definition here):

readgroup = DataIlluminaPairedEnd500Lib
data = ../../data/project_1.fastq ../../data/project_2.fastq
technology = solexa
template_size = 250 750
segment_placement = ---> <---

In cases you are not 100% sure about, e.g., the size of the DNA template, you can also give a (generous) expected range and then tell MIRA to automatically refine this range during the assembly based on real, observed distances of read pairs. Do this with autorefine modifier like this:

template_size = 50 1000 autorefine

The following manifest file is an example for mapping a 500 bp paired-end and a 3kb mate-pair library of a strain called bchoc_se1 against a GenBank reference file containing a strain called bchoc_wt:

# Example for a manifest describing a mapping assembly with
# paired Illumina data, not merging reads and therefore keeping
# all pair information

# First part: defining some basic things
# In this example, we just give a name to the assembly
#  and tell MIRA it should map a genome in accurate mode
# As special parameter, we want to switch off merging of Solexa reads

project = MyFirstMapping
job = genome,mapping,accurate
parameters = SOLEXA_SETTINGS -CO:msr=no

# The second part defines the sequencing data MIRA should load and assemble
# The data is logically divided into "readgroups"

# first, the reference sequence
readgroup
is_reference
data = ../../data/NC_someNCBInumber.gff3
strain = bchoc_wt

# now the Illumina data

readgroup = DataForPairedEnd500bpLib
autopairing
data = ../../data/project500bp-1.fastq ../../data/project500bp-2.fastq
technology = solexa
strain = bchoc_se1

readgroup = DataForMatePair3kbLib
data = ../../data/project3kb-1.fastq ../../data/project3kb-2.fastq
technology = solexa
strain = bchoc_se1
template_size = 2000 4000 autorefine
segment_placement = <--- --->

Please look up the parameters used in the main manual. The ones above basically say: make an accurate mapping of Solexa reads against a genome. Additionally do not merge short short Solexa reads to the contig.

For the paired-end library, be lazy and let MIRA find out everything it needs. However, that information should be treated as "information only" by MIRA, i.e., it is not used for deciding whether a pair is well mapped.

For the mate-pair library, assume a DNA template template size of 2000 to 4000 bp (but let MIRA automatically refine this using observed distances) and the segment orientation of the read pairs follows the reverse / forward scheme. That information should be treated as "information only" by MIRA, i.e., it is not used for deciding whether a pair is well mapped.

Comparing this manifest with a manifest for unpaired-data, two parameters were added in the section for Solexa data:

	Note
	Note that in mapping assemblies, these template_distance and segment_placement parameters are normally treated as information only, i.e., MIRA will map the reads regardless whether the distance and orientation criterions are met or not. This enables post-mapping analysis programs to hunt for genome rearrangements or larger insertions/deletion.

Warning

If template size and segment placement checking were on, the following would happen at, e.g. sites of re-arrangement: MIRA would map the first read of a read-pair without problem. However, it would very probably reject the second read because it would not map at the specified distance or orientation from its partner. Therefore, in mapping assemblies with paired-end data, checking of the template size must be switched off to give post-processing programs a chance to spot re-arrangements.

I'm sure you'll have picked up the general scheme of manifest files by now. Hybrid mapping assemblies follow the general scheme: simply add as separate readgroup the information MIRA needs to know to find the data and off you go. Just for laughs, here's a manifest for 454 shotgun with Illumina paired-end

# Example for a manifest describing a mapping assembly with
# shotgun 454 and paired-end Illumina data, not merging reads and therefore keeping
# all pair information

# First part: defining some basic things
# In this example, we just give a name to the assembly
#  and tell MIRA it should map a genome in accurate mode
# As special parameter, we want to switch off merging of Solexa reads

project = MyFirstMapping
job = genome,mapping,accurate
parameters = SOLEXA_SETTINGS -CO:msr=no

# The second part defines the sequencing data MIRA should load and assemble
# The data is logically divided into "readgroups"

# first, the reference sequence
readgroup
is_reference
data = ../../data/NC_someNCBInumber.gff3
strain = bchoc_wt

# now the shotgun 454 data
readgroup = DataForShotgun454
data = ../../data/project454data.fastq
technology = 454
strain = bchoc_se1

# now the paired-end Illumina data

readgroup = DataForPairedEnd500bpLib
data = ../../data/project500bp-1.fastq ../../data/project500bp-2.fastq
technology = solexa
strain = bchoc_se1
template_size = 250 750
segment_placement = ---> <---

MIRA will make use of ancillary information present in the manifest file. One of these is the information to which strain (or organism or cell line etc.pp) the generated data belongs.

You just need to tell in the manifest file which data comes from which strain. Let's assume that in the example from above, the "lane6" data were from a first mutant named bchoc_se1 and the "lane7" data were from a second mutant named bchoc_se2. Here's the manifest file you would write then:

# Example for a manifest describing a mapping assembly with
# unpaired Illumina data

# First part: defining some basic things
# In this example, we just give a name to the assembly
#  and tell MIRA it should map a genome in accurate mode

project = MyFirstMapping
job = genome,mapping,accurate

# The second part defines the sequencing data MIRA should load and assemble
# The data is logically divided into "readgroups"

# first, the reference sequence
readgroup
is_reference
data = ../../data/NC_someNCBInumber.gff3
strain = bchoc_wt

# now the Illumina data

readgroup = DataForSE1
data = ../../data/bchocse_lane6.solexa.fastq
technology = solexa
strain = bchoc_se1

readgroup = DataForSE2
data = ../../data/bchocse_lane7.solexa.fastq
technology = solexa
strain = bchoc_se2

	Note
	While mapping (or even assembling de-novo) with multiple strains is possible, the interpretation of results may become a bit daunting in some cases. For many scenarios it might therefore be preferable to successively use the data sets in own mappings or assemblies.

This strain information for each readgroup is really the only change you need to perform to tell MIRA everything it needs for handling strains.

For bacteria you can map EST / RNASeq data against any reference sequence(s), be it a genome reference or transcriptome reference.

For eukaryotes, I recommend mapping against transcriptome sequences only as MIRA does not handle splice sites very well. You can get away mapping RNASeq data against 'simple' eukaryotes which have very few and small introns, like S. cerevisiae or simple fungi; but just note that the expression level for exons near larger introns will be slightly lower than the rest of the gene.

Note

The standard MIRA parameters are geared towards handling difficult mapping cases where the reference sequence can be quite different to what is actually mapped against it (10% average difference, up to 15% locally). If the data you are mapping is actually from a strain which is very similar to your reference (give or take a couple of SNPs, genome reorganisations and / or insertions and deletions), you can save a lot of time in the mapping process by relaxing a few parameters like [-SK:kms=(high-value)], [-SB:bnb=no], [-FM:mte=(low-value)].

Here's a very simple, albeit powerful manifest which will allow you to retrieve both SNPs and feature coverage data from the mapping results:

# Example for a lazy manifest describing a mapping assembly with
# one library of paired reads

# First part: defining some basic things
# In this example, we just give a name to the assembly
#  and tell MIRA it should map EST / RNASeq data in accurate mode
# Accurate mode may be "too good" for mapping, see next example

project = MyFirstMapping
job = est,mapping,accurate

# The second part defines the sequencing data MIRA should load and assemble
# The data is logically divided into "readgroups"

# first the reference sequence
readgroup
is_reference
data = ../../data/NC_someNCBInumber.gff3
strain = bchoc_wt

# now the Illumina data

readgroup
data = ../../data/project_1.fastq ../../data/project_2.fastq
technology = solexa
strain = bchoc_se1

If you compare the above manifest to a manifest for a simple genome mapping, the only line which really changed is the job=... line, which changed from genome,mapping,accurate to est,mapping,accurate

As already alluded to, the standard MIRA mapping parameters are made to find almost anything that maps, and that takes time. In case you are mapping data which is very closely related to the reference sequence you are using, you can save time by changing the manifest slightly. Start it like this:

# Example for a lazy manifest describing a mapping assembly with
# one library of paired reads which is close to your reference sequence

# First part: defining some basic things

project = MyFirstMapping
job = est,mapping,accurate

# Accurate mode may be "too good" for mapping, therefore these parameters to
# speed up things
parameters = -SK:kms=20 SOLEXA_SETTINGS -SB:bnb=no -FM:mte=6

# remaining manifest as above
# ...

TODO: Sorry, needs to be re-written for the relatively new SRR format distributed at the NCBI ... and changes in MIRA 3.9.x. Please come back later.

Poly-A tails are part of the mRNA and therefore also part of sequenced data. They can occur as poly-A or poly-T, depending from which direction and which part of the mRNA was sequenced. Having poly-A/T tails in the data is a something of a double edged sword. More specifically., if the 3' poly-A tail is kept unmasked in the data, transcripts having this tail will very probably not align with similar transcripts from different splice variants (which is basically good). On the other hand, homopolymers (multiple consecutive bases of the same type) like poly-As are features that are pretty difficult to get correct with today's sequencing technologies, be it Sanger, Solexa or, with even more problems problems, 454. So slight errors in the poly-A tail could lead to wrongly assigned splice sites ... and wrongly split contigs.

This is the reason why many people cut off the poly-A tails. Which in turn may lead to transcripts from different splice variants being assembled together.

Either way, it's not pretty.

Single transcripts (or very lowly expressed transcripts) containing SNPs, splice variants or similar differences to other, more highly expressed transcripts are a problem: it's basically impossible for an assembler to distinguish them from reads containing junky data (e.g. read with a high error rate or chimeras). The standard setting of many EST assemblers and clusterers is therefore to remove these reads from the assembly set. MIRA handles things a bit differently: depending on the settings, single transcripts with sufficiently large differences are either treated as debris or can be saved as singlet.

Another interesting problem for de-novo assemblers are non-normalised libraries. In each cell, the number of mRNA copies per gene may differ by several orders of magnitude, from a single transcripts to several tens of thousands. Pre-sequencing normalisation is a wet-lab procedure to approximately equalise those copy numbers. This can however, introduce other artifacts.

If an assembler is fed with non-normalised EST data, it may very well be that an overwhelming number of the reads comes only from a few genes (house-keeping genes). In Sanger sequencing projects this could mean a couple of thousand reads per gene. In 454 sequencing projects, this can mean several tens of thousands of reads per genes. With Solexa data, this number can grow to something close to a million.

Several effects then hit a de-novo assembler, the three most annoying being (in ascending order of annoyance): a) non-random sequencing errors then look like valid SNPs, b) sequencing and library construction artefacts start to look like valid sequences if the data set was not cleaned "enough" and more importantly, c) an explosion in time and memory requirements when attempting to deliver a "good" assembly. While MIRA has methods to deal with this kind of data (e.g. via digital normalisation), a sure sign of the latter are messages from MIRA about megahubs in the data set.

	Note
	The guide on how to tackle hard projects with MIRA gives an overview on how to hunt down sequences which can lead to the assembler getting confused, be it sequencing artefacts or highly expressed genes.

Chimeras are sequences containing adjacent base stretches which are not occurring in an organism as sequenced, neither as DNA nor as (m)RNA. Chimeras can be created through recombination effects during library construction or sequencing. Chimeras can, and often do, lead to misassemblies of sequence stretches into one contig although they do not belong together. Have a look at the following example where two stretches (denoted by x and o are joined by a chimeric read r4 containing both stretches:

r1 xxxxxxxxxxxxxxxx
r2 xxxxxxxxxxxxxxxxx
r3 xxxxxxxxxxxxxxxxx
r4 xxxxxxxxxxxxxxxxxxx|oooooooooooooo
r5                        ooooooooooo
r6                        ooooooooooo
r7                          ooooooooo

The site of the recombination event is denoted by x|o in read r4.

MIRA does have a chimera detection -- which works very well in genome assemblies due to high enough coverage -- by searching for sequence stretches which are not covered by overlaps. In the above example, the chimera detection routine will almost certainly flag read r4 as chimera and only use a part of it: either the x or o part, depending on which part is longer. There is always a chance that r4 is a valid read though, but that's a risk to take.

Now, that strategy would also work totally fine in EST projects if one would not have to account for lowly expressed genes. Imagine the following situation:

s1 xxxxxxxxxxxxxxxxx
s2         xxxxxxxxxxxxxxxxxxxxxxxxx
s3                          xxxxxxxxxxxxxxx
    

Look at read s2; from an overlap coverage perspective, s2 could also very well be a chimera, leading to a break of an otherwise perfectly valid contig if s2 were cut back accordingly. This is why chimera detection is switched off by default in MIRA.

Warning

When starting an EST assembly via the --job=est,... switch, chimera detection is switched off by default. It is absolutely possible to switch on the SKIM chimera detection afterwards via [-CL:ascdc]. However, this will have exactly the effects described above: chimeras in higher coverage contigs will be detected, but perfectly valid low coverage contigs will be torn apart. It is up to you to decide what you want or need.

Imagine this simple case: a gene has two slightly different alleles and you've sequenced this:

A1-1  ...........T...........
A1-2  ...........T...........
A1-3  ...........T...........
A1-4  ...........T...........
A1-5  ...........T...........
B2-1  ...........G...........
B2-2  ...........G...........
B2-3  ...........G...........
B2-4  ...........G...........
      

Depending on base qualities and settings used during the assembly like, e.g., [-CO:mr:mrpg:mnq:mgqrt:emea:amgb] MIRA will recognise that there's enough evidence for a T and also enough evidence for a G at that position and create two contigs, one containing the "T" allele, one the "G". The consensus will be >99% identical, but not 100%.

Things become complicated if one has to account for errors in sequencing. Imagine you sequenced the following case:

A1-1  ...........T...........
A1-2  ...........T...........
A1-3  ...........T...........
A1-4  ...........T...........
A1-5  ...........T...........
B2-1  ...........G...........
      

It shows very much the same like the one from above, except that there's only one read with a "G" instead of 4 reads. MIRA will, when using standard settings, treat this as erroneous base and leave all these reads in a contig. It will likewise also not mark it as SNP in the results. However, this could also very well be a lowly expressed transcript with a single base mutation. It's virtually impossible to tell which of the possibilities is right.

	Note
	You can of course force MIRA to mark situations like the one depicted above by, e.g., changing the parameters for [-CO:mrpg:mnq:mgqrt]. But this may have the side-effect that sequencing errors get an increased chance of getting flagged as SNP.

Further complications arise when SNPs and potential sequencing errors meet at the same place. consider the following case:

A1-1  ...........T...........
A1-2  ...........T...........
A1-3  ...........T...........
A1-4  ...........T...........
B1-5  ...........T...........
B2-1  ...........G...........
B2-2  ...........G...........
B2-3  ...........G...........
B2-4  ...........G...........
E1-1  ...........A...........
      

This example is exactly like the first one, except an additional read E1-1 has made it's appearance and has an "A" instead of a "G" or "T". Again it is impossible to tell whether this is a sequencing error or a real SNP. MIRA handles these cases in the following way: it will recognise two valid read groups (one having a "T", the other a "G") and, in assembly mode, split these two groups into different contigs. It will also play safe and define that the single read E1-1 will not be attributed to either one of the contigs but, if it cannot be assembled to other reads, form an own contig ... if need to be even only as single read (a singlet).

	Note
	Depending on some settings, singlets may either appear in the regular results or end up in the debris file.

Gaps in alignments of transcripts are handled very cautiously by MIRA. The standard settings will lead to the creation of different contigs if three or more consecutive gaps are introduced in an alignment. Consider the following example:

A1-1  ..........CGA..........
A1-2  ..........*GA..........
A1-3  ..........**A..........
B2-1  ..........***..........
B2-2  ..........***..........
      

Under normal circumstances, MIRA will use the reads A1-1, A1-2 and A1-3 to form one contig and put B2-1 and B2-2 into a separate contig. MIRA would do this also if there were only one of the B2 reads.

The reason behind this is that the probability for having gaps of three or more bases only due to sequencing errors is pretty low. MIRA will therefore treat reads with such attributes as coming from different transcripts and not assemble them together, though this can be changed using the [-AL:egp:egpl] parameters of MIRA if wanted.

Problems with homopolymers, especially in 454, Ion Torrent and high coverage Illumina

As 454 and Ion Torrent sequencing has a general problem with homopolymers, this rule of MIRA will sometimes lead formation of more contigs than expected due to sequencing errors at "long" homopolymer sites ... where long starts at ~6-7 bases. Though MIRA does know about the problem in 454 homopolymers and has some routines which try to mitigate the problem. this is not always successful. The same applies for Illumina data with long homopolymers (~ 8-9 bp) and high coverage (≥ 100x).

The following files in projectname_d_results contain results of the assembly in different formats. Depending on the output options you defined for MIRA, some files may or may not be there. As long as the CAF or MAF format are present, you can translate your assembly later on to about any supported format with the miraconvert program supplied with the MIRA distribution:

The following files in projectname_info contain statistics and other information files of the assembly:

miraconvert is tool in the MIRA package which reads and writes a number of formats, ranging from full assembly formats like CAF and MAF to simple output view formats like HTML or plain text.

Figure 7.1. miraconvert supports a wide range of format conversions to simplify export / import of results to and from other programs

The question "How Do I convert to / from other tools?" is complicated by the plethora of file formats and tools available. This section gives an overview on what is needed to reach the most important ones.

Figure 7.2.  Conversion steps, formats and programs needed to reach some tools like assembly viewers, editors or scaffolders.

Please also read the chapter on MIRA utilities in this manual to learn more on miraconvert and have a look at miraconvert -h which lists all possible formats and other command line options.

$ miraconvert maf -t sam INPUT.maf OUTPUT

MIRA sets a number of different tags in resulting assemblies. They can be set in reads (in which case they mostly end with a r) or in the consensus.(then ending with a c).

Note

If you use the Staden gap4, gap5 or consed assembly editor to tidy up the assembly, you can directly jump to places of interest that MIRA marked for further analysis by using the search functionality of these programs. However, you need to tell these programs that these tags exist. For that you must change some configuration files. More information on how to do this can be found in the support/README file of the MIRA distribution.

You should search for the following "consensus" tags for finding places of importance (in this order).

of lesser importance are the "read" versions of the tags above:

In normal assemblies (only one sequencing technology, just one strain), search for the IUPc, UNSc, SRMc and WRMc tags.

In hybrid assemblies, searching for the IUPc, UNSc, SRMc, WRMc, and STMU tags and correcting only those places will allow you to have a qualitatively good assembly in no time at all.

Columns with SRMr tags (SRM in Reads) in an assembly without a SRMc tag at the same consensus position show where mira was able to resolve a repeat during the different passes of the assembly ... you don't need to look at these. SRMc and WRMc tags however mean that there may be unresolved trouble ahead, you should take a look at these.

Especially in mapping assemblies, columns with the MCVc, SROx, SIOx and SAOx tags are extremely helpful in finding places of interest. As they are only set if you gave strain information to MIRA, you should always do that.

For more information on tags set/used by MIRA and what they exactly mean, please look up the according section in the reference chapter.

The read coverage histogram as well as the template display of gap4 will help you to spot other places of potential interest. Please consult the gap4 documentation.

I recommend to invest a couple of minutes (in the best case) to a few hours in joining contigs, especially if the uniform read distribution option of MIRA was used (but first filter for large contigs). This way, you will reduce the number of "false repeats" in improve the overall quality of your assembly.

7.5.3.1.  Joining contigs at true repetitive sites

Joining contigs at repetitive sites of a genome is always a difficult decision. There are, however, two rules which can help:

1. If the sequencing was done without a paired-end library, don't join.
2. If the sequencing was done with a paired-end library, but no pair (or template) span the join site, don't join.

The following screen shot shows a case where one should not join as the finishing program (in this case gap4) warns that no template (read-pair) span the join site:

Figure 7.3.  Join at a repetitive site which should not be performed due to missing spanning templates.

The next screen shot shows a case where one should join as the finishing program (in this case gap4) finds templates spanning the join site and all of them are good:

Figure 7.4.  Join at a repetitive site which should be performed due to spanning templates being good.

Remember that MIRA takes a very cautious approach in contig building, and sometimes creates two contigs when it could have created one. Three main reasons can be the cause for this:

When working with resequencing data and a mapping assembly, I always load finished projects into an assembly editor and perform a quick cleanup of the results. SNP or small indels normally do not need cleanups, but every mapper will get larger indels mostly wrong, and MIRA is no exception to this.

For close relatives of the reference strain this doesn't take long as MIRA will have set tags (see section earlier in this document) at all sites you should have a look at. For example, very close mutant bacteria with just SNPs or simple deletions and no genome reorganisation, I usually clean up in 10 to 15 minutes. That gives the last boost to data quality and your users (biologists etc.) will thank you for that as it reduces their work in analysing the data (be it looking at data or performing wet-lab experiments).

The general workflow I use is to convert the CAF file to a gap4 or gap5 database. Then, in gap4 or gap5, I

After this, I convert the gap4 or gap5 database back to CAF format. But beware: gap4 does not have the same consensus calling routines as MIRA and will have saved it's own consensus in the new CAF. In fact, gap4 performs rather badly in projects with multiple sequencing technologies. So I use miraconvert from the MIRA package to recall a good consensus (and save it in MAF as it's more compact and a lot faster in handling than CAF):

And from this MAF file I can then convert with miraconvert to any other format I or my users need: CAF, FASTA, ACE, WIG (for coverage analysis), SNP and coverage analysis (see below), HTML etc.pp.

Biologists are not really interested in SNPs coordinates, and why should they? They're more interested where SNPs are, how good they are, which genes or other elements they hit, whether they have an effect on a protein sequence, whether they may be important etc. For organisms without intron/exon structure or splice variants, MIRA can generate pretty comprehensive tables and files if an annotated GenBank file was used as reference and strain information was given to MIRA during the assembly.

Well, MIRA does all that automatically for you if the reference sequence you gave was annotated.

For this, miraconvert should be used with the asnp format as target and a MAF (or CAF) file as input:

$ miraconvert -t asnp input.maf output

Note that it is strongly suggested to perform a quick manual cleanup of the assembly prior to this: for rare cases (mainly at site of small indels of one or two bases), MIRA will not tag SNPs with a SNP tag (SROc, SAOc or SIOc) but will be fooled into a tag denoting unsure positions (UNSc). This can be quickly corrected manually. See further down in this manual in the section on post-processing.

After conversion, you will have four files in the directory which you can all drag-and-drop into spreadsheet applications like OpenOffice Calc or Excel.

The files should be pretty self-explanatory, here's just a short overview:

I've come to realise that people who don't handle data from NextGen sequencing technologies on a regular basis (e.g., many biologists) don't want to be bothered with learning to handle specialised programs to have a look at their resequenced strains. Be it because they don't have time to learn how to use a new program or because their desktop is not strong enough (CPU, memory) to handle the data sets.

Something even biologist know to operate are browsers. Therefore, miraconvert has the option to load a MAF (or CAF) file of a mapping assembly at output to HTML those areas which are interesting to biologists. It uses the tags SROc, SAOc, SIOc and MCVc and outputs the surrounding alignment of these areas together with a nice overview and links to jump from one position to the previous or next.

This is done with the '-t hsnp' option of miraconvert:

$ miraconvert -t hsnp input.maf output

	Note
	I recommend doing this only if the resequenced strain is a very close relative to the reference genome, else the HTML gets pretty big. But for a couple of hundred SNPs it works great.

miraconvert can also dump a coverage file in WIG format (using '-t wig') or a WIG file for GC content (using '-t gcwig'). This comes pretty handy for searching genome deletions or duplications in programs like the Affymetrix Integrated Genome Browser (IGB, see http://igb.bioviz.org/) or when looking for foreign sequence in a genome.

When having data mapped against a reference with annotations (either from GenBank formats or GFF3 formats), miraconvert can generate tables depicting either expression values (in RNASeq/EST data mappings) or probable genome multiplication and deletion factors (in genome mappings). For this to work, you must use a MAF or CAF file as input, specify fcov as output format and the reference sequence must have had annotations during the mapping with MIRA.

TODO: add example

miraconvert -t fcov mira_out.maf myfeaturetable

miraconvert is a tool to convert, extract and sometimes recalculate all kinds of data related to sequence assembly files.

More specifically, miraconvert can

…

…

In the following examples, the CAF and MAF files used are expected to contain full assembly data like the files created by MIRA during an assembly or by the gap2caf program. CAF and MAF could be used interchangeably in these examples, depending on which format currently is available. In general though, MAF is faster to process and smaller on disk.

mirabait selects reads from a read collection which are partly similar or equal to sequences defined as target baits. Similarity is defined by finding a user-adjustable number of common k-mers (sequences of k consecutive bases) which are the same in the bait sequences and the screened sequences to be selected, either in forward or reverse complement direction.

When used on paired files (-p or -P), selects read pairs where at least one read matches.

One can use mirabait to do targeted assembly by fishing out reads belonging to a gene and just assemble these; or to clean out rRNA sequences from data sets; or to fish out and iteratively reconstruct mitochondria from metagenomic data; or, or, or ... whenever one has to take in or take out subsets of reads based on kmer equality, this tool should come in quite handy.

	Note
	The search performed is exact, that is, sequences selected are guaranteed to have the required number of matching k-mers to the bait sequences while sequences not selected are guaranteed not have these.

	Note
	The examples below, together with the manual above, should be enough to get you going. If there's a typical use case you are missing, feel free to suggest it on the MIRA talk mailing list.

Baiting unpaired sequences, bait sequences in FASTA, sequences in FASTQ:

mirabait -b b.fasta file.fastq

Same as above, but baits in two files (FASTA and GenBank):

mirabait -b b1.fasta -b b2.gbk file.fastq

Baiting paired sequences, read pairs are in two files:

mirabait -b b.fasta -p file_1.fastq file_2.fastq

Baiting paired sequences, pairs are interleaved in one file:

mirabait -b b.fasta -P file.fastq

Like above, but selecting sequences which do not match the baits:

mirabait -i -b b.fasta -P file.fastq

Baiting paired sequences (file_1.fastq, file_2.fastq and file3.fastq) and unpaired sequences (file4.fastq), all at once and different file types:

mirabait -b b.fasta -p file_1.fastq file_2.fastq -P file3.fasta file4.caf

Like above, but writing sequences matching baits and sequences not matching baits to different files:

mirabait -I -b b.fasta -p file_1.fastq file_2.fastq -P file3.fasta file4.caf

Change bait criterion to need 10 kmers of size 27:

mirabait -k 27 -n 10 -b b.fasta file.fastq

Change bait criterion to baiting only reads which have all kmers present in the bait:

mirabait -n 0 -b b.fasta file.fastq

Change bait criterion to baiting all reads having almost all kmers present in the bait, but allowing for up to 40 kmers not in the bait:

mirabait -n -40 -b b.fasta file.fastq

Force bait sequences to load as FASTA, force sequences to be baited to be loaded as FASTQ:

mirabait -b fasta::b.dat fastq::file.dat

Write result files to directory /dev/shm/:

mirabait -o /dev/shm/ -b b.fasta -p file_1.fastq file_2.fastq

Merge all result files containing sequences hitting baits to file /dev/shm/match:

mirabait -o /dev/shm/match -b b.fasta -p file_1.fastq file_2.fastq

Like above, but also merge all result files containing sequences not hitting baits to file /dev/shm/nomatch:

mirabait -o /dev/shm/match -O /dev/shm/nomatch -b b.fasta -p file_1.fastq file_2.fastq

Fetch all reads having rRNA motifs in a paired FASTQ files:

mirabait -j rrna -p file1.fastq file2.fastq

Fetch all reads not having rRNA motifs in a paired FASTQ files:

mirabait -j rrna -i -p file1.fastq file2.fastq

Split a paired FASTQ file into two sets of files (4 files total), one containing rRNA reads and one containing non-rRNA reads:

mirabait -j rrna -I -p file1.fastq file2.fastq

Assuming the file human_genome.fasta contains the human genome: bait all read pairs matching the human genome. Also, save the compute kmer statistics for later re-use to file HG_kmerstats.mhs.gz:

mirabait -b human_genome.fasta -K HG_kmerstats.mhs.gz -p file1.fastq file2.fastq

The same as above, but just precompute and save the kmer statistics, no actual baiting done.

mirabait -b human_genome.fasta -K HG_kmerstats.mhs.gz

Using the precomputed kmer statistics from the command above: bait files with read pairs for human reads:

mirabait -B HG_kmerstats.mhs.gz -p file_1.fastq file_2.fastq

The standard database for rRNA baiting supplied with the MIRA source code and binary packages is called rfam_rrna-21-12.sls.gz which will get installed as MHS (MiraHashStatistics) file into $BINDIR/../share/mira/mhs/rfam_rrna-21-12.mhs.gz (where $BINDIR is the directory where the mira/mirabait binary resides) and a soft link pointing from filter_default_rrna.mhs.gz to rfam_rrna-21-12.mhs.gz like so:

arcadia:~$ which mira
/usr/local/bin/mira
arcadia:~$ ls -l /usr/local/share/mira/mhs
lrwxrwxrwx 1 bach bach        22 Mar 24 23:58 filter_default_rrna.mhs.gz -> rfam_rrna-21-12.mhs.gz
-rw-rw-r-- 1 bach bach 148985059 Mar 24 23:58 rfam_rrna-21-12.mhs.gz

The file naming scheme for the database is as following: dbidentifier-kmerlength-kmerfreqcutoff. The standard database is therefore: rfam_rrna as identifier for the RFAM rRNA sequences (currently RFAM 12), then 21 defining a kmer length of 21 and finally a kmer cutoff frequency of 12, meaning that kmers must have been seen at least 12 times in the RFAM database to be stored in the subset.

	Note
	The value of 12 as frequency cutoff for the standard mirabait rRNA database was chosen as a compromise between sensitivity and database size.

Although rRNA are pretty well conserved overall, the cutoff frequency also implies that kmers from rare rRNA variants will not be present in the database, eventually losing some sensitivity for rRNA from rarely sequenced organisms. It follows that more sensitive versions of the rRNA database can be installed by downloading a file from the MIRA repository at SourceForge and calling a script provided by MIRA. To install a version with a kmer size of 21 and a cutoff frequency of, e.g., 3, download rfam_rrna-21-3.sls.gz and install it like this:

arcadia:~/tmp$ ls
arcadia:~/tmp$ wget https://sourceforge.net/projects/mira-assembler/files/MIRA/slsfiles/rfam_rrna-21-3.sls.gz
...

TODO: continue docs here.

A manifest file can be seen as a two part configuration file for an assembly: the first part contains some general information and tell MIRA how exactly you want the data to be assembled while the second part contains information about the sequencing data to be loaded. Examples being always easier to follow than long texts, here's an example for a de-novo assembly with single-end (also called shotgun) 454 data:

# Example for a manifest describing a simple 454 de-novo assembly

# A manifest file can contain comment lines, these start with the #-character

# First part: defining some basic things
# In this example, we just give a name to the assembly
#  and tell MIRA it should assemble a genome de-novo in accurate mode
# As special parameter, we want to use 4 threads in parallel (where possible)


project = MyFirstAssembly
job = genome,denovo,accurate
parameters = -GE:not=4

# The second part defines the sequencing data MIRA should load and assemble
# The data is logically divided into "readgroups": this reflects the
#  ... that read sequences ...

readgroup = SomeUnpaired454ReadsIGotFromTheLab
data = TCMFS456ZH345.fastq TQF92GT7H34.fastq
technology = 454

To make things a bit more interesting, here's an example using a couple more technologies and showing some more options of the manifest file like wild cards in file names, different paired-end/mate-pair libraries and how to let MIRA refine pairing information (or even find out everything by itself):

# Example for a manifest describing a de-novo assembly with
# unpaired 454, paired-end Illumina, a mate-pair Illumina
# and a paired Ion Torrent

# First part: defining some basic things
# In this example, we just give a name to the assembly
#  and tell MIRA it should assemble a genome de-novo in accurate mode
# As special parameter, we want to use 4 passes with kmer sizes of
# 17, 31, 63 and 127 nucleotides. Obviously, read lengths of the
# libraries should be greater than 127 bp.
# Note: usually MIRA will choose sensible options for number of
#  passes and kmer sizes to be used by itself.

project = MyFirstAssembly
job = genome,denovo,accurate
parameters = -AS:kms=17,31,63,127

# The second part defines the sequencing data MIRA should load and assemble
# The data is logically divided into "readgroups": this reflects the
#  ... that read sequences ...

# defining the shotgun (i.e. unpaired) 454 reads
readgroup = SomeUnpaired454ReadsIGotFromTheLab
data = TCMFS456ZH345.fastq TQF92GT7H34.fastq
technology = 454

# defining the paired-end Illumina reads, fixing all needed pair information
readgroup = SomePairedEndIlluminaReadsIGotFromTheLab
data = datape*.fastq
technology = solexa
template_size = 100 300
segment_placement = ---> <---
segment_naming = solexa

# defining the mate-pair Illumina reads, fixing most needed pair information
#  but letting MIRA refine the template_size via "autorefine"
readgroup = SomeMatePairIlluminaReadsIGotFromTheLab
data = datamp*.fastq
technology = solexa
template_size = 2000 4000 autorefine
segment_placement = <--- --->
segment_naming = solexa

# defining paired Ion Torrent reads
# example to show how lazy one can be and simply let MIRA estimate by itself
#  all needed pairing information via "autopairing"
#  Hint: it usually does a better job at it than we do ;-)
readgroup = SomePairedIonReadsIGotFromTheLab
autopairing
data = dataion*.fastq
technology = iontor

The first part of an assembly manifest contains the very basic information the assembler needs to have to know what you want it to do. This part consists of exactly three entries:

The second part of an assembly manifest tells MIRA which files it needs to load, which sequencing technology generated the data, whether there are DNA template constraints it can use during the assembly process and a couple of other things.

9.3.3.1.  Starting a new readgroup

readgroup [= group name] is the keyword which tells MIRA that you are going to define a new read group. You can optionally name that group.

Understanding readgroups and DNA templates

When you send away your DNA for sequencing, it is going to be prepared for sequencing according to your wishes. Sequencing providers call this "constructing a library" and regardless whether you sequence with Sanger, 454, Illumina, Ion Torrent, Pacific Biosciences or other technologies, the "library prep" is always there.

With most library preps, your DNA is first amplified and then cut into small pieces. These pieces are called templates and their length can be anywhere between a few dozen bases, a few hundred bases or even a couple of dozen or even hundred kilobases. The important thing is that these templates can be much bigger in size than the actual read length. While this is a wet lab step, protocols and providers have gotten pretty good at constructing libraries where the DNA templates are all in a given range of bases like, e.g., having a library with template size 500bp (+/- 100bp) and another library with template size around 7kb (+/- 500bp).

Depending on the technology and sequencing strategy used, the DNA templates are used to create either one single read or - and that's important - two or more reads.

Libraries with "single reads" are often called "single read libraries" or "shotgun libraries". They can be found for every sequencing technology and are most of the time easy to construct (therefore cheap) and are often used to provide a decent amount of bases as basic coverage for your project.

Libraries with two reads per DNA template are often called "mate-pair" or "paired-end" libraries. They are harder to construct and sometime have less yield, therefore they are often more expensive. But the sequencing approach using several reads per DNA template allows assembly and scaffolding algorithms to resolve repetitive regions of a genome which are longer than the average read length. Note that Pacific Biosciences has a sequencing mode called "strobed sequencing" which is different from "paired-end/mate-pair" but also creates multiple reads per DNA template.

Long story short: an assembler must know afterwards what kind of reads it has to expect: the sequencing technology, library preparation strategy etc. For this, the notion of read groups has emerged: reads coming from the same technology and same library preparation are pooled together in a read group to tell the assembler: in the assembly, if you see two reads coming from a same DNA template, you should expect them to be at a certain distance from each other and they should be oriented in a certain way.

	Note
	The above was a very simplified view on the whole area of DNA templates, readgroups, shotgun and paired end sequencing. Enough to hopefully understand the concepts, but you might want to read more about it.

9.3.3.2.  Defining data files to load

data = filepath [filepath ...] defines the file paths from which sequences should be loaded. A file path can contain just the name of one (or several) files or it can contain the path, i.e., the directory (absolute or relative) including the file name.

MIRA automatically recognises what type the sequence data is by looking at the postfix of files. For postfixes not adhering widely used naming schemes for file types, there's additionally a way of explicitly defining the type (see further down at the end of this item on how this is done). Currently allowed file types are:

• .fasta for sequences formatted in FASTA format where there exists an additional .fasta.qual file which contains quality data. If the file with quality data is missing, this is interpreted as error and MIRA will abort.

• .fna and .fa also for sequences formatted in FASTA format. The difference to .fasta lies in the way MIRA treats a missing quality file (called .fna.qual or .fa.qual): it does not see that as critical error and continues.

• .fastq or .fq for files in FASTQ format

• .gff3 or .gff for files in GFF3 format. Note that MIRA will load all sequences and annotations contained in this file.

• .gbk, .gbf, .gbff or .gb for files formatted in GenBank format. Note that the MIRA GenBank loader does not understand intron/exon or other multiple-locus structures in this format, use GFF3 instead!

• .caf for files in the CAF format (from Sanger Centre)

• .maf for files in the MIRA MAF format

• .exp for files in the Staden EXP format.

• .fofnexp for a file of EXP filenames which all point to files in the Staden EXP format.

• .xml, .ssaha2 and .smalt for ancillary data in NCBI TRACEINFO, SSAHA2 or SMALT format respectively.

Multiple 'data' lines and multiple entries per line (even different formats) are allowed, as in, e.g.,

data = file1.fastq file2.fastq file3.fasta file4.gbk
data = myscreenings.smalt

You can also use wildcards and/or directory names. E.g., loading all file types MIRA understand from a given directory mydir:

data = mydir

or loading all files starting with mydata and ending with fastq:

data = mydata*fastq

or loading all files in directory mydir starting with mydata and ending with fastq:

data = mydir/mydata*fastq

or loading all FASTQ files in all directories starting with mydir:

data = mydir*/*fastq

or ... well, you get the gist.

	Note
	Giving a directory like in mydir is equivalent to mydir/* (saying: give me all files in the directory mydir), however the first version should be preferred when the directory contains thousands of files.

Note

GenBank and GFF3 files may or may not contain embedded sequences. If annotations are present in these files for which no sequence is present in the same file, MIRA will look for reads of the same name which it already loaded in this or previously defined read groups and add the annotations there. As security measure, annotations in GenBank and GFF3 files for which absolutely no sequence or read has been defined are treated as error.

Explicit definition of file types. It is possible to explicitly tell MIRA the type of a file even if said file does not have a 'standard' naming scheme. For this, the EMBOSS double-colon notation has been adapted to work also for MIRA, i.e., you prepend the type of a file and separate it from the file name by a double colon. E.g., the .dat postfix is not anything MIRA will recognise, but you can define it should be loaded as FASTQ file like this:

data = fastq::myfile.dat

Another frequent usage is forcing MIRA to load FASTA files named .fasta without complaining in case quality files (which MIRA wants you to provide) are not present:

data = fna::myfile.fasta

This does (of course) work also with directories or wildcard characters. In the following example, the first line will load all files from mydirectory as FASTQ while the second line loads just .dat files in a given path as FASTA:

data = fastq::mydirectory
data = fasta::/path/to/somewhere/*.dat

It is entirely possible (although not really sensible), to give contradicting information to MIRA by using a different explicit file type than one would guess from the standard postfix. In this case, the explicit type takes precedence over the automatic type. E.g.: to force MIRA to load a file as FASTA although it is named .fastq, one could use this:

data = fasta::file.fastq

Note that the above does not make any kind of file conversion, file.fastq needs to be already in FASTA format or else MIRA will fail loading that data.

9.3.3.7.  Autopairing: letting MIRA find out pair info by itself

autopairing This keyword is used to tell MIRA it should estimate values for template_size and segment_placement (see below).

This is basically the lazy way to tell MIRA that the data in the corresponding readgroup consists of paired reads and that you trust it will find out the correct values.

	Note
	autopairing usually works quite well for small and mid-sized libraries (up to, say, 10 kb). For larger libraries it might be a good thing to tell MIRA some rough boundaries via template_size / segment_placement and let MIRA refine the values for the template size via autorefine (see below).

	Note
	autopairing is a feature new to MIRA 4.0rc5, it may contain bugs for some corner cases. Feedback appreciated.

9.3.3.8.  Setting size of read templates

template_size = min_size max_size [infoonly|exclusion_criterion] [autorefine]. Defines the minimum and maximum size of "good" DNA templates in the library prep for this read group. This defines at which distance the two reads of a pair are to be expected in a contig, a very useful information for an assembler to resolve repeats in a genome or different splice variants in transcriptome data.

If the term infoonly is present, then MIRA will pass the information on template sizes in result files, but will not use it for any decision making during de-novo or mapping assembly. The term exclusion_criterion makes MIRA use the information for decision making.

If infoonly or exclusion_criterion are missing, then MIRA assumes exclusion_criterion for de-novo assemblies and infoonly for mapping assemblies.

If the term autorefine is present, MIRA will start the assembly with the given size information but switch to refined value computed from observed distances in an assembly. However, please note that the size values can never be expanded, only shrunk. It is therefore advisable to use generous bounds when using the autorefine feature.

	Note
	The template_size line in the manifest file replaces the parameters -GE:uti:tismin:tismax of earlier versions of MIRA (3.4.x and below).

	Note
	The minimum or the maximum size (or both) can be set to a negative value for "don't care and don't check". This allows constructs like template_size= 500 -1 exclusion_criterion which would check only the minimum distance but not the maximum distance.

Note

For mapping assemblies with MIRA, you usually will want to use infoonly as else - in case of genome re-arrangements, larger deletions or insertions - MIRA would probably reject one read of every read pair in the corresponding areas as it would not be at the expected distance and/or orientation ... and you would not be able to simply find the re-arrangement in downstream analysis. For de-novo assemblies however you should not use infoonly except in very rare cases where you know what you do.

Understanding the size of DNA templates

When using a paired-end or mate-pair sequencing strategy, two sequences are generated for the ends of each DNA template (see sidebar above: "understanding readgroups and DNA templates"). That is, if one has a library with 6kb fragments, one knows that the outer ends of the two reads will be approximately 6kb apart, like so:

DNA template    ##############################################################
read 1          .......
read 2                                                                  ......
<------------------------- ~6 kb ---------------------------->

Sequencing labs will try their best to get these two sequences from DNA templates which comply to a given length specification. But as this is chemistry and wet lab, things must be seen with a certain uncertainty and therefore the DNA templates generated are not exactly of the specified size (e.g. 6kb), but the size distribution will vary in a given range, e.g., 5.5kb to 6.5 kb.

9.3.3.9.  Read segment placement

	Note
	You do not need to use this when using 'autopairing' (see above).

segment_placement = placementcode [infoonly|exclusion_criterion]. Allowed placement codes are:

• ? or unknown which are place-holders for "well, in the end: don't care." Segments of a template can be reads in any direction and in any relationship to each other.

  This is typically used for unpaired libraries (sometimes called shotgun libraries), but may be also useful for, e.g., primer walking with Sanger.

• ---> <--- or FR or INNIES. The forward / reverse scheme as used in traditional Sanger sequencing as well as Illumina paired-end sequencing,

  This is the usual placement code for Sanger paired-end protocols as well as Illumina paired-end. Less frequently used in IonTorrent paired-end sequencing.

• <--- ---> or RF or OUTIES. The reverse / forward scheme as used in Illumina mate-pair sequencing.

  This is the usual placement code for Illumina mate-pair protocols.

• 1---> 2---> or samedir forward or SF or LEFTIES. The forward / forward scheme. Segments of a template are all placed in the same direction, the segment order in the contig follows segment ordering of the reads.

• 2---> 1---> samedir backward or SB or RIGHTIES. Segments of a template are all placed in the same direction, the segment order in the contig is reversed compared to segment ordering of the reads.

  This is the usual placement code for 454 "paired-end" and IonTorrent long-mate protocols.

• samedir Segments of a template are all placed in the same direction, the spatial relationship however is not cared of.

• >>> (reserved for sequencing of several equidistant fragments per template like in PacBio strobe sequencing, not implemented yet)

If the term infoonly is present, then MIRA will pass the information on segment placement in result files, but will not use it for any decision making during de-novo assembly or mapping assembly. The term exclusion_criterion makes MIRA use the information for decision making.

If infoonly or exclusion_criterion are missing, then MIRA assumes exclusion_criterion for de-novo assemblies and infoonly for mapping assemblies.

Note

For mapping assemblies with MIRA, you usually will want to use infoonly as else - in case of genome re-arrangements, larger deletions or insertions - MIRA would probably reject one read of every read pair (as it would not be at the expected distance and/or orientation) and you would not be able to simply find the re-arrangement in downstream analysis. For de-novo assemblies however you should not use infoonly except in very rare cases where you know what you do.

Note

As soon as you tell MIRA that a readgroup contains paired reads (via one of the other typical readgroup parameters like template_size, segment_naming etc.), the segment_placement line becomes mandatory in the manifest. This is because different sequencing technologies and/or library preparations result in different read orientations. E.g., Illumina libraries come in paired-end flavour which have FR (forward/reverse) placements, but there are also mate-pair libraries which have reverse/forward (RF) placements.

Understanding read segment placement on DNA templates

bla

9.3.3.10.  Read segment naming

segment_naming = naming_scheme [rollcomment]. Defines the naming scheme reads are following to indicate the DNA template they belong to. Allowed naming schemes are: sanger, stlouis, tigr, FR, solexa, sra.

If not defined, the defaults are sanger for Sanger sequencing data, while solexa for Solexa, 454 and Ion Torrent.

For FASTQ files, the modifier rollcomment can be used to let MIRA take the first token in the comment as name of a read instead of the orginal name. E.g.: for a read

@DRR014327.1.1 HWUSI-EAS547_0013:1:1:1106:4597.1 length=91
TTAGAAGGAGATCTGGAGAACATTTTAAACCGGATTGAACAACGCGGCCGTGAGATGGAGCTTCAGACAAGCCGGTCTTATTGGGACGAAC
+
bbb`bbbbabbR`\_bb_bba`b`bb_bb_`\^\^Y^`\Zb^b``]]\S^a`]]a``bbbb_bbbb]bbb\`^^^]\aaY\`\\^aa__aB

the rollcomment modifier will lead to the read being named HWUSI-EAS547_0013:1:1:1106:4597.1 (which is almost the original instrument read name) instead of DRR014327.1.1 (which is the SRA read name).

	Warning
	For data from the short read archive (SRA), one will usually need to explicitly specify the 'sra' naming scheme or use the 'rollcomment' modifier in FASTQ files.

	Note
	This has changed with MIRA 3.9.1 and sff_extract 0.3.0. Before that, 454 and Ion Torrent were given fr as naming scheme.

Understanding read naming schemes

Read naming is a long story with lots of historical gotchas: it needs to be clear and simple, but still people sometimes wanted to convey additional meta-information with it. Unsurprisingly, several "standards" emerged over time. In short: it's a mess. See also XKCD entry on proliferating standards.

How to choose: please read the documentation available at the different centres or ask your sequence provider. In a nutshell (and probably over-simplified):

Sanger scheme

"somename.[pqsfrw][12][bckdeflmnpt][a|b|c|..." (e.g. U13a08f10.p1ca), but the length of the postfix must be at least 4 characters, i.e., ".p" alone will not be recognised.

Usually, ".p" + 3 characters or "f" + 3 characters are used for forwards reads, while reverse complement reads take either ".q" or ".r" (+ 3 characters in both cases).

TIGR scheme

"somenameTF*|TR*|TA*" (e.g. GCPBN02TF or GCPDL68TABRPT103A58B),

Forward reads take "TF*", reverse reads "TR*".

St. Louis scheme

"somename.[sfrxzyingtpedca]*"

Forward/Reverse scheme

"somename.[fr]*" (e.g. E0K6C4E01DIGEW.f or E0K6C4E01BNDXN.r2nd),

".f*" for forward, ".r*" for reverse.

Solexa scheme

Even simpler than the forward/reverse scheme, it allows only for one two reads per template: "somename/[12]"

SRA scheme

The Short Read Archive (SRA) finally settled on a naming scheme and renames each and every read within its database. When you download sequences from the archive, all reads will be named XXX000000.Y[.Z] (where X's are characters A-Z, 0 are digits from 0 to 9, Y is a counter and Z is a number denoting the segment (usually 1,2 or 3)). This naming scheme is applied to reads from all technologies, therefore the MIRA technology dependent defaults will not apply and one must specify the 'sra' naming scheme in the command line.

Any wildcard in the forward/reverse suffix must be consistent for a read pair, and is treated as part of the template name. This is to allow multiple sequencing of a fragment, particularly common with Sanger capillary data (e.g. given somename.f and somename.r, resequenced as somename.f2 and somename.r2, this would be treated as two pairs, with template names somename and somename_2 respectively).

9.3.3.13.  Renaming read name prefixes

rename_prefix= prefix replacement. Allows to rename reads on the fly while loading data by searching each read name for a given prefix string and, if found, replace it with a given replacement string.

This is most useful for systems like Illumina or PacBio which generate quite long read names which, in the end, are either utterly useless for an end user or are even breaking older programs which have a length restriction on read names. E.g.:

rename_prefix = DQT9AAQ4:436:H371HABMM: Sample1_

will rename reads like DQT9AAQ4:436:H371HABMM:5:1101:9154:3062 into Sample1_5:1101:9154:3062

	Note
	rename_prefix entries are valid per readgroup. I.e., an entry for a readgroup will not rename reads of another readgroup.

Note

Multiple rename_prefix entries are allowed per readgroup. E.g.: rename_prefix = DQT9AAQ4:436:H371HABMM: S1sxa_ rename_prefix = m140328_002546_42149_c100624422550000001823118308061414_s1_ S1pb_ will rename a read called DQT9AAQ4:436:H371HABMM:1:1101:3099:2186 into S1sxa_1:1101:3099:2186 while renaming another read called m140328_002546_42149_c100624422550000001823118308061414_s1_p0/100084/10792_20790/0_9573 into S1pb_p0/100084/10792_20790/0_9573

The parameters= line in the manifest file opens up the full panoply of possibilities the MIRA assembler offers. This ranges from fine-tuning assemblies to setting parameters in a way so that MIRA is suited also for very special assembly cases.

parameters =  -GENERAL:number_of_threads=4 \
              -ALIGN:min_relative_score=70 -ASSEMBLY:minimum_read_length=150 \
              -OUTPUT:output_result_caf=no

9.3.4.2.  Technology sections

With the introduction of new sequencing technologies, MIRA also had to be able to set values that allow technology specific behaviour of algorithms. One simple example for this could be the minimum length a read must have to be used in the assembly. For Sanger sequences, having this value to be 150 (meaning a read should have at least 150 unclipped bases) would be a very valid, albeit conservative choice. For 454 reads and especially Solexa and ABI SOLiD reads however, this value would be ridiculously high.

To allow very fine grained behaviour, especially in hybrid assemblies, and to prevent the explosion of parameter names, MIRA knows two categories of parameters:

1. technology independent parameters which control general behaviour of MIRA like, e.g., the number of assembly passes or file names etc.

2. technology dependent parameters which control behaviour of algorithms where the sequencing technology plays a role. Example for this would be the minimum length of a read (like 200 for Sanger reads and 120 for 454 FLX reads).

More on this a bit further down in this documentation.

As example, a manifest using technology dependent and independent parameters could look like this:

parameters = COMMON_SETTINGS -GENERAL:number_of_threads=4 \
              SANGER_SETTINGS -ALIGN:min_relative_score=70 -ASSEMBLY:minimum_read_length=150 \
              454_SETTINGS -ALIGN:min_relative_score=75 -ASSEMBLY:minimum_read_length=100 \
              SANGER_SETTINGS -ALIGN:min_relative_score=90 -ASSEMBLY:minimum_read_length=75

Now, assume the following read group descriptions in a manifest:

...

readgroup
technology=454
...

readgroup
technology=solexa
...

For MIRA, this means a number of parameters should apply to the assembly as whole, while others apply to the sequencing data itself ... and some parameters might need to be different depending on the technology they apply to. MIRA dumps the parameters it is running with at the beginning of an assembly and it makes it clear there which parameters are "global" and which parameters apply to single technologies.

Here is as example a part of the output of used parameters that MIRA will show when started with 454 and Illumina (Solexa) data:

...

Assembly options (-AS):
    Number of passes (nop)                      : 1
    Skim each pass (sep)                        : yes
    Maximum number of RMB break loops (rbl)     : 1
    Spoiler detection (sd)                      : no
    Last pass only (sdlpo)                      : yes

    Minimum read length (mrl)                   :  [454]  40
                                                   [sxa]  20
    Enforce presence of qualities (epoq)        :  [454]  no
                                                   [sxa]  yes

...

You can see the two different kind of settings that MIRA uses: common settings (like [-AS:nop]) which allows only one value and technology dependent settings (like [-AS:mrl]), where for each sequencing technology used in the project, the setting can be different.

How would one set a minimum read length of 40 and not enforce presence of base qualities for Sanger reads, but for 454 reads a minimum read length of 30 and enforce base qualities? The answer:

job=denovo,genome,draft
parameters= SANGER_SETTINGS -AS:mrl=40:epoq=mo 454_SETTINGS -AS:mrl=40:epoq=yes

Notice the ..._SETTINGS section in the command line (or parameter file): these tell MIRA that all the following parameters until the advent of another switch are to be set specifically for the said technology.

Note

For improved readability, you can distribute parameters across several lines either by pre-fixing every line with parameter=, like so: job=denovo,genome,draft parameters= SANGER_SETTINGS -AS:mrl=80:epoq=no parameters= 454_SETTINGS -AS:mrl=30:epoq=yes Alternatively you can use a backslash at the end of a parameter line to indicate that the next line is a continuing line, like so: job=denovo,genome,draft parameters= SANGER_SETTINGS -AS:mrl=80:epoq=no \ 454_SETTINGS -AS:mrl=30:epoq=yes Note that the very last line of the parameters settings MUST NOT end with a backslash.

Beside COMMON_SETTINGS there are currently 6 technology settings available:

1. SANGER_SETTINGS

2. 454_SETTINGS

3. IONTOR_SETTINGS

4. PCBIOLQ_SETTINGS (currently not supported)

5. PCBIOHQ_SETTINGS

6. SOLEXA_SETTINGS

7. TEXT_SETTINGS

Some settings of MIRA are influencing global behaviour and are not related to a specific sequencing technology, these must be set in the COMMON_SETTINGS environment. For example, it would not make sense to try and set different number of assembly passes for each technology like in

parameters= 454_SETTINGS -AS:nop=4 SOLEXA_SETTINGS -AS:nop=3

Beside being contradictory, this makes not really sense. MIRA will complain about cases like these. Simply set those common settings in an area prefixed with the COMMON_SETTINGS switch like in

parameters= COMMON_SETTINGS -AS:nop=4 454_SETTINGS ... SOLEXA_SETTINGS ...

Since MIRA 3rc3, the parameter parser will help you by checking whether parameters are correctly defined as COMMON_SETTINGS or technology dependent setting.

parameters = COMMON_SETTINGS -GENERAL:number_of_threads=4 \
              SANGER_SETTINGS -ALIGN:min_relative_score=70 -ASSEMBLY:minimum_read_length=150 \
              454_SETTINGS -ALIGN:min_relative_score=75 -ASSEMBLY:minimum_read_length=100 \
              SOLEXA_SETTINGS -ALIGN:min_relative_score=90 -ASSEMBLY:minimum_read_length=75

parameters = COMMON_SETTINGS -GE:not=4 \
              SANGER_SETTINGS -AL:mrs=70 -AS:mrl=150 \
              454_SETTINGS -AL:mrs=75 -AS:mrl=100 \
              SOLEXA_SETTINGS -AL:mrs=90 -AS:mrl=75

parameters = COMMON_SETTINGS -GE:not=4 \
              SANGER_SETTINGS \
                -AL:mrs=70 \
		-AL:mrl=150 \
              454_SETTINGS -AL:mrs=75:mrl=100 \
              SOLEXA_SETTINGS \
	        -AL:mrs=90 \
                -AL:mrl=75