In this manual it is assumed that "Sibelia" is properly installed and the directory with "Sibelia.py" is in your "PATH" environment variable or input files are in the same folder with "Sibelia" executable.
Directory "examples/Sibelia" contains two sets of bacterial genomes. The easiest way to run "Sibelia" is to type:
Sibelia -s loose <input FASTA file(s)>
For example:
Sibelia -s loose Helicobacter_pylori.fasta
Important note -- "Sibelia" requires some free space on HDD to run. If you experience any problems and see error messages that mention temporary files, try to change directory used for temporary files (see section "Directory for temporary files").
Above commands will run "Sibelia" on the file "Helicobacter_pylori.fasta" with the "loose" simplification parameters. There are two another simplification parameters sets, called "fine" and "far". To run "Sibelia" on "Helicobacter_pylori.fasta" with "fine" parameters set, type:
Sibelia -s fine Helicobacter_pylori.fasta
The difference between "loose" and "fine" set is that "loose" usually produces fewer blocks, but longer. And it may lose some small synteny blocks (shorter than 15 000 BP), while "fine" option produces more blocks (but shorter) and their coverage is worse. Usually "loose" is the best choice, but if you do not want to lose information about small-scale rearrangements, use "fine". See also section "Output description" for detailed depiction of the output format. The "far" parameters set can be used for analysis of distantly-related genomes. However, usage of this set requires more computational time and space and may be not appropriate in case of many genomes.
If you are not satisfied by the results (poor coverage, for example), try to set simplification parameters manually (see section "Fine tuning").
By default, Sibelia filters out synteny blocks shorter than 5 000 BP. You can change this behaviour, see section "Minimum block size".
You also may be interested in blocks that occur exactly once in each input sequence, such blocks are used as input for MGRA algorithm, for example. To get such output use option "-a":
Sibelia -s loose -a Helicobacter_pylori.fasta
Synteny blocks are visualized with an interactive diagram (see section "d3" visualization"). The blocks are also can be visualized with "Circos" (see section "Circos" visualization"). While "Circos" is better for publications, "d3" diagram is more suitable for analysis. Please note that sequences that do not contain any synteny blocks instances are not shown on these diagrams.
Genomes from the "examples/Sibelia" dir were taken from [5, 6]. Note that you can specify multiple FASTA files, just separate them with spaces.
Default directory = "." You can change this by setting cmd parameter:
-o <dir name> or --outdir <dir name>
By default, "Sibelia" places output files in current working directory. Setting this parameter will change output directory.
Default directory is output directory. You can change this by setting cmd parameter:
-t <dir name> or --tempdir <dir name>
"Sibelia" creates some temporary files while running. By default these files are placed in the output directory. If you want to place temporary files in another folder due to some reasons, use this parameter. Although the files exist only for a very short period of time, they can be quite big -- ~20*N bytes, where N is the total size of all input genomes. You can also use switch
-r or --inram
that will force "Sibelia" to not create any temporary files and store all it's data in RAM.
By default, "Sibelia" produces following files:
- Blocks coordinates
- Genomes represented as permutations of the synteny blocks
- Coverage report
- Files for generating a "Circos" picture
- Interactive html-diagram of synteny blocks
There are also optional output files:
- Sequences file
- Dot file with resulting de Bruijn graph
All these files are described below in details.
File name = "blocks_coords.txt". First part of this file lists input sequences, their IDs, sizes and descriptions. IDs are just index numbers of sequences (in the same order as they apper in the input files).
Second part of the file describes synteny blocks in sections separated by dashes. Each section starts with the synteny block ID. Then all instances of this block are listed in tabular format. Each row of a table depicts one instance of this synteny block. Columns of the table designate following:
- Seq_id -- ID of the sequence, that the instance belongs to
- Strand -- strand of the synteny block instance, either '+' or '-'. Input sequences are treated as positive strands
- Start -- one-based index of the first base pair of the instance
- End -- one-based index of the last base pair of the instance
- Length -- length of the instance of the synteny block
Note that all coordinates are given relatively to POSITIVE strand of the sequence. If an instance of a synteny block is located on the positive strand, then start < end, otherwise start > end.
If you wish, you can obtain coordinates in GFF format. If you use the flag:
-gff
Then coordinates of synteny blocks are listed in file "blocks_coords.gff". For description of this format, see:
https://cgwb.nci.nih.gov/FAQ/FAQformat.html#format3
Each record represents different copies of a synteny block. Copies having the same number in the "tag" field (last column) are instances of the same synteny block.
File name = "genomes_permutations.txt".
This file contains input sequences represented as permutations of the synteny block. It has two lines for each input sequence:
- Header line -- FASTA header of the sequence (starting with '>')
- Genome line -- sequence of synteny blocks instances as they appear on the positive strand of the sequence. Each instance is represented as a signed integer, where '+' sign depicts direct version of the block, and '-' depicts reversed block
File name = "coverage_report.txt".
The file describes portion of the genomes, that found synteny block cover. First part of the file describes input sequencess (see "Blocks coordinates" section). Second part of the file is a table with the following columns:
- Degree -- multiplicity of a synteny block. For example, if a synteny block has degree = 3, then the are three instances of this block in the input sequence
- Count -- number of synteny blocks with a given degree
- Total -- portion of all the input sequences that are covered by the blocks with a given degree
- Seq -- portion of the sequence with id that is covered by blocks with a given degree
Here is an example of such table from a report file.
| Degree | Count | Total | Seq 1 | Seq 2 | Seq 3 | Seq 4 |
|---|---|---|---|---|---|---|
| 2 | 11 | 3.82% | 6.59% | 2.41% | 2.96% | 3.30% |
| 3 | 4 | 1.68% | 2.24% | 2.19% | 1.44% | 0.85% |
| 4 | 21 | 91.93% | 91.34% | 94.71% | 87.54% | 94.53% |
| All | 36 | 95.66% | 97.44% | 97.98% | 90.67% | 96.89% |
This table contains one row for each degree (2, 3, 4) and one ("All") row for the overall coverage. It means that there are 11 blocks with degree = 2, i.e. 11 * 2 instances, and they cover 3.82% of all four genomes, 6.59% of Seq 1 and 2.41% of Seq 2. And also there are 21 synteny blocks with degree = 4, i.e. 4 * 21 instances and they cover 91.93% of all genomes. All the blocks cover 95.66% of all the input sequences.
File name = "blocks_sequences.fasta". By default this file is not written. To output this file, set cmd parameter:
-q or --sequencesfile
This FASTA file contains sequences of instances of the synteny block. Each sequence has header in following format:
>Seq="<header>",Strand=<sign>,Block_id=<block id>,Start=<x>,End=<y>
Where "
" is a header of the FASTA sequence where the block instance is located. Other fields are described in section "Coordinates file". Sequences of synteny blocks are also written in SAM format in the file "blocks_sequences.sam".File name = "d3_blocks_diagram.html".
"Sibelia" generates an interactive html diagram that shows found synteny blocks. Coordinates follow the same convention as described in section "Coordinates file".
You can visualize synteny blocks with a colorful circular diagram by using the "Circos" software [3]. Files for generating such diagram are written in directory "circos" inside the output directory. To generate "Circos" diagram do following:
- Download and install "Circos" software
- Run "Sibelia"
- Run "Circos" in "circos" directory
For example of such diagrams (generated from "Helicobacter_pylori.fasta), see "examples/Sibelia/Helicobacter_pylori/circos/circos.png". Also note that such diagrams can become very piled with larger number of genomes. To overcome this, plot only big blocks, see section "Minimum block size". Blocks located on the positive strand are colored green, while blocks from negative strand are red.
By default, "Sibelia" plots only blocks obtained after the last stage. You can also view blocks at the intermediate stages by using switch:
-v or --visualize
On the resulting diagram the outermost circle shows blocks obtained at the first stage, then the second stage and so on. Please note that this option slows down the computation.
File name = "de_bruijn_graph.dot". By default this file is not written. To output this file, set cmd parameter:
-g or --graphfile
If you are a curious person, you can also view condensed de Bruijn graph that is used for generating synteny blocks. To understand the graph, see [1]. Condensed means that only bifurcations in the graph are plotted and non-branching paths are collapsed into a single edge. Blue edges are generated from positive strand and red edges are from negative strand respectively. Note that this graph is generated for K = min(Kn, MinimumBlockSize) or for value of "--lastk" cmd parameter if it is set, where Kn is the value of K used for the last stage (see section "Fine tuning").
If one is interested in graph output only, he or she can use the "--noblocks" option:
--noblocks
In this case Sibelia doesn't compute the synteny blocks and doesn't ouput them, but can output the graph. For example, to get only non-modified compressed de Bruijn grahp for k=25, one can use the following command line:
Sibelia -k run.stage --noblocks -g -m 25 <input_file>
Where "run.stage" contains single number "0".
Here we will describe parameters that can affect computation results.
Default value = 5000. To change this value, set cmd parameter:
-m <integer> or --minblocksize <integer>
If you are interested only in big synteny blocks, like > 100 000 BP, set this parameter to an appropriate value.
Default = not set. Add flag to cmd parameters to set:
-a or --sharedonly
Output only blocks that occur exactly once in each input sequence. This option assumes that all input genomes contain single chromosome.
By default, synteny blocks are postprocessed after computation by gluing "stripes" consisting of the same synteny blocks. For example, if each occurence of synteny block 1 is followed by synteny block 2 and vice a versa, their directions are consistent, then they are "glued" together to form a single synteny block. Postprocessing could be turned off by specifying flag:
--nopostprocess
"Sibelia" performs computations in multiple stages. Every stage produces it's own synteny blocks. You can get blocks from all stages by specifying flag:
--allstages
Files "blocks_coordsN.(txt|gff)" will contain their coordinates, where N is the number of stage. Zero corresponds to blocks obtained without any simplification.
Algorithm of "Sibelia" depends on presence of solid k-mers within syntenic regions. If such regions contain variations close to their borders, they will be truncated. In case of two genomes, some synteny blocks (of multiplicity two) could be corrected using local alinment algorithms. Use flag:
--correctboundaries
This flag is used by "C-Sibelia".
Default value = not set. To select the parameters set, use cmd parameter:
-s <loose|fine|far> or --parameters <loose|fine>
This option is incompatible with "-k", you must specify one of these, not both. Approach used in "Sibelia" is parameter dependent. To understand the details, please see the next section and [1]. The "loose" option produces longer blocks and better coverage, while "fine" can capture small-scale rearrangements, for example, inversions of size < 15 000 BP. The "far" set is for distant genomes.
Default value = not set. To specify the file that contains custom parameters set, use cmd parameter:
-k <file name> or --stagefile <file name>
The algorithm consists of several stages of computations. Each stage has two parameters, K and D. Let's call K-mer a substring of length K. At each stage "Sibelia" constructs so called de Bruijn graph, graph of K-mers that occur in the genome, and simplifies it by removing special type of undirected cycles called "bulges", see [1] for more details.
Graph is a good model for describing the algorithm, but to understand "physical meaning" of the parameters it is useful to consider operations that are actually performed with the genome behind the graph model. Suppose that somewhere in the genome exist two pairs of K-mers K1 and K2:
1st pair: ... K1 ABCD K2 ...
2nd pair: ... K1 FGHE K2 ...
If the distance between K1 and K2 within each pair is less than D, then "Sibelia" replaces FGHE with ABCD to obtain longer "synteny block":
1st pair: ... K1 ABCD K2 ...
2nd pair: ... K1 ABCD K2 ...
More concrete example. Suppose that K = 3, D = 5 and somewhere in the genome we find:
1st pair: ... act gaga ggc ...
2nd pair: ... act gatg ggc ...
As we see, distance between "act" and "ggc" is less than 5 nucleotides so we replace "gatg" by "gaga":
1st pair: ... act gaga ggc ...
2nd pair: ... act gaga ggc ...
"Sibelia" keeps track of all changes so it is able to locate original locations of the synteny blocks obtained by the simplification. This process continues step by step, we start with small values of K to obtain longer K-mers shared between synteny regions and then increase K and D. The "loose" parameters set has 4 stages:
| K | D |
|---|---|
| 30 | 150 |
| 100 | 1000 |
| 1000 | 5000 |
| 5000 | 15000 |
The "fine" set consists of 3 stages and it's final values are less:
| K | D |
|---|---|
| 30 | 150 |
| 100 | 500 |
| 500 | 1500 |
The "far" set is for distant genomes:
| K | D |
|---|---|
| 15 | 120 |
| 100 | 500 |
| 500 | 1500 |
As you can see, "loose" set is more aggressive -- at it's final stage it glues together 5000-mers that are separated from each other by at most 15000 symbols. Although this description is very simplified and lacks many important technical details, it is enough to infer your own parameter set. Stage file that you may use to specify your own parameters has following simple format:
M
K1 D1
K2 D2
...
KM KM
Where M is the number of stages. So, running with the stage file:
4
30 150
100 1000
1000 5000
5000 15000
Is equivalent to running with the -s "loose" cmd option. As you may notice, the algorithm relies on exact K-mers shared between the genomes. If input genomes doesn't have such shared substrings, then "Sibelia" won't be able to locate the synteny blocks.
If you cannot find synteny blocks with the default parameters, try to start with smaller values of K (~20), increase D values or vary number of stages.
In "Sibelia" de Bruijn graph is constructed (N + 1) times, where N is the number of simplification stages. De Bruijn graph which is constructed last is used for inferring synteny blocks. Value of K for this grahp is determined by min(Kn, MinimumBlockSize) where Kn is the value of K used for the last stage. You can override this value by setting the "--lastk" parameter:
--lastk 1>
Default value = 4. Tho change this value, set cmd parameter:
-i <integer> or --maxiterations <integer>
Maximum number of iterations during a stage of simplification. Increasing this parameter may slightly improve coverage.
- Ilya Minkin, Nikolay Vyahhi, Son Pham. "SyntenyFinder: A Synteny Blocks Generation and Genome Comparison Tool" (poster), WABI 2012 http://bioinf.spbau.ru/sites/default/files/SyntenyFinder.pdf
- Max A. Alekseyev and Pavel A. Pevzner. "Breakpoint graphs and ancestral genome reconstructions", Genome Res. 2009. 19: 943-957.
- Circos. http://circos.ca
- D3. http://d3js.org/
- Helicobacter pylori. http://www.ncbi.nlm.nih.gov/genome/169
- Staphylococcus aureus. http://www.ncbi.nlm.nih.gov/genome/154