An ultra fast, heuristic approach to detect conserved signals in extremely large pairwise genome comparisons.
GCC compiler (any version that is not completely outdated should do), the R programming language (base installation) and python3 (tested on 3.5 and 3.8). Please make sure that on your linux CC resolves to GCC, otherwise it might not work.
Simply download the .zip and unzip it, or clone the repository. Then issue the following command:
cd chromeister && make all -C src/ && python3 -m venv chromeisterenv && source chromeisterenv/bin/activate && pip install -r src/requirements.txt
This will compile CHROMEISTER and create a virtualenv where the python libraries will be installed (see src/requirements.txt
)
If the installation finished without errors, you are ready to go! If you encounter the following problem ImportError: No module named 'skbuild'
you might need to do source chromeisterenv/bin/activate && pip install --upgrade pip
and then run pip install -r src/requirements.txt
to finish the installation.
NOTE: python and its libraries are only used for the detection of events. If the binaries compile (i.e. the make all
) then you can still run CHROMEISTER and plot the results, even if the python installation did not work.
There are several ways in which CHROMEISTER can be used. The simplest one is to run a 1-vs-1 comparison and then compute the score and the plot. To do so, use the binaries at the bin folder:
You can run CHROMEISTER directly by issuing:
CHROMEISTER -query seqX -db seqY -out dotplot.mat && Rscript compute_score.R dotplot.mat 1000
If you do not want a grid on the output dotplot (which is recommended when running comparisons with a lot of scaffolds for instance) then run the same command but replace compute_score by compute_score-nogrid, see below:
CHROMEISTER -query seqX -db seqY -out dotplot.mat && Rscript compute_score-nogrid.R dotplot.mat 1000
The 1000 value is the default size of dotplot.mat, i.e. the resolution of the matrix -- if you want to change this (for example to generate a larger image (if you use 2000 it will generate a plot of 2000x2000, so be careful) include also the parameter -dimension in CHROMEISTER. Example command with larger resolution:
CHROMEISTER -query seqX -db seqY -out dotplot.mat -dimension 2000 && Rscript compute_score.R dotplot.mat 2000
And if you want to run the events detection, use (make sure that your virtualenv chromeisterenv
is in the chromeister root folder:
source ../chromeisterenv/bin/activate && python3 bin/detect_events.py dotplot.mat.raw.txt
This will generate a dotplot.mat.events.txt
file containing the detected events and classified. If you want to get a plot of the signal with the overlapped detected events, issue the same command but add at the end the parameter png
(separated by a space).
You can also use the script that is in the bin folder (which will do all of the above for you):
run_and_plot_chromeister.sh (input sequence A) (input sequence B) (KMER size) (DIMENSION of plot) (inexactitude level) [optional: grid]
(see parameters at the end) (the grid keyword at the end can be included/omitted depending if you want grid in the output dotplot)
This will generate the following items:
You can run massive all versus all comparisons in two diferent ways:
Comparing all the sequences in one folder. This accounts for 1/2 n (n+1) comparisons, hence it will not compare sequence B to sequence A if the comparison for sequence A to sequence B already existed.
allVsAll.sh <sequences folder> <extension (e.g. fasta)> <matrix size (1000 for chromosomes, 2000 for full genomes)> <kmer size 32 (32 is best)> <inexactitude level (4 is recommended)>
Comparing two folders containing sequences. This accounts for n * m comparisons, therefore it will compare ALL to ALL. Use this for instance to compare all chromosomes of one genome to all chromosomes of another genome.
allVsAll_incremental.sh <sequences folder 1> <sequences folder 2> <extension (e.g. fasta)> <matrix size (1000 for chromosomes, 2000 for full genomes)> <kmer size 32 (32 is best)> <inexactitude level (4 is recommended)>
At the end of both comparisons, an index will be created summarizing the scores per each comparison. This index has the following format (see header and example below): header: <SpX, SpY, IDX, IDY, IMG, CHNumberX, CHNumberY, Score, LengthX, LengthY> example: BRAOL.Chr.C1,BRAOL.Chr.C2,>C1 dna:chromosome chromosome:v2.1:C1:1:43764888:1 REF,>C2 dna:chromosome chromosome:v2.1:C2:1:52886895:1 REF,BRAOL.C hr.C1.fasta-BRAOL.Chr.C2.fasta.mat.filt.png,C1,C2, 0.996,43764888,52886895
Notice that you can easily run this in parallel by just re-issuing the command (i.e. execute same command as many times as you want, each time another core will help in the processing).
First of all, consider whether it is interesting or not to use CHROMEISTER for "fine-grained" results. CHROMEISTER is recommended for VERY coarse-grained and full-genome comparisons in order to quickly assess similarity between genomes. Thus it does NOT produce alignments. However, if you find yourself in a situation where you want to convert the signal of CHROMEISTER into alignments (e.g. two large genomes), this can be done. The following tutorial shows how to do it, with human chromosome X and mouse chromosome X as example:
First, run CHROMEISTER like this:
./CHROMEISTER -query HOMSA.Chr.X.fasta -db MUSMU.Chr.X.fasta -out dotplot.mat -dimension 1000 && Rscript compute_score.R dotplot.mat 1000
Check the "dotplot.mat.filt.png" corresponding to the dotplot between both chromosomes to see if there is any similarity. If so, proceed to next step.
Clone the following repository: https://github.com/estebanpw/gecko
git clone https://github.com/estebanpw/gecko
Switch branch to the one named "inmemory_guided_chrom" and compile it. To do so, issue the following command:
cd gecko && git checkout inmemory_guided_chrom && make all -C src
Now run the script "guidefastas" in the bin folder. See below:
bin/guidefastas.sh HOMSA.Chr.X.fasta MUSMU.Chr.X.fasta hits-XY-dotplot.mat.hits 1000 100 60 32
You can add the following arguments to the execution:
--alignments
: This one will extract the alignments and write them to a file with extension .alignments
in the all-results
folder.--names
: This one will output the names of the sequences to which each fragment belongs instead of their sequence number (e.g. if comparing chromosomes 1 and 2 of Homo sapiens vs Mus musculus then 0,1
will now be Homo s. chr1, Mus m. chr2
)--sort
: This will sort output frags in the csv
file according first to the comparison they belong to and secondly by coordinates. --local
: This one will convert the coordinates in the csv
file (which are global in respect to the file) to local in respect to each sequence in the fasta file. Useful for chromosome comparison when we do not want all coordinates to accumulate.Note (1): remember to include the full path to the sequences. Note (2): the "hits-XY-dotplot.mat.hits" file is produced by CHROMEISTER in step 1. Copy it to the folder or include full path. Note (3): the parameters following in the command "1000 200 75 32" are namely (1) size of dotplot, (2) minimum length that an alignment must have to be reported, (3) minimum similarity from 0-100, (4) k-mer seed size (use 32 for chromomsome-like sequences).
This step can take several minutes, e.g. using 1 CPU this execution took around 9-10 minutes.
A CSV file containing the alignments coordinates can be found in the folder all-results
. You can download it here if you wish to do so: http://mango.ac.uma.es/compartir/HOMSA_X-MUSMU_X.csv
If you also wish to visually contrast annotations to the alignments, you can use our genomic browser at https://pistacho.ac.uma.es/. To do so just follow the user guide available at https://pistacho.ac.uma.es/static/data/GeckoMGV-UserGuide.pdf
An example alignment file is shown below:
AAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAAAGAAAGAAAGAAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAGAAAGAAAA
||||||||||||||||||||||||||||||||||| | ||| || || ||| ||| |||||||||||||||||||||||||||||||||| | | | ||
AAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAAAAGAAAGAAGGAAGGAAGGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAA
@ FORWARD STRAND x1: 10501385 y1: 169863509 x2: 10501485 y2: 169863609 Identity: 88/101 (87.1287%)
TTTTCCCATTGATTAATATTTTTCCTGTTGAGCAGATGAGAGAAAGCCAAAAAAAGCACAGCTGGGCCATTTCCCCTCACTGGGAACGTCATTTCCAGGCACTTTGTGCTTACTTGAT
|||||||||||||| ||||||||| | |||| | |||||||||||||||||||||||||||||||||||||||| ||||||| |||||||||||||| |||| |||||| |||||||
TTTTCCCATTGATTGATATTTTTCTTATTGAACTGATGAGAGAAAGCCAAAAAAAGCACAGCTGGGCCATTTCCTCTCACTGTAAACGTCATTTCCAGTCACTCTGTGCTCACTTGAT
@ REVERSE STRAND x1: 10451224 y1: 169981322 x2: 10451341 y2: 169981205 Identity: 107/118 (90.678%)
Note that each fragment begins with an alignment, and has a corresponding @
line where information is included such as strand of the alignments, coordinates, identity, etc.
USAGE:
You can test CHROMEISTER with the two mycoplasma sequences provided in the 'test-data' folder. You can do so by running the following commands (from within the test-data folder):
../bin/CHROMEISTER -query mycoplasma-232.fasta -db mycoplasma-7422.fasta -out mycoplasma-232-7422.mat -dimension 500
Rscript ../bin/compute_score.R mycoplasma-232-7422.mat 500
Note: in this example we used size 500 since the two sequences are quite small.
Comparing two chromosomes (Homo sapiens Chr X vs Mus musculus Chr X) in a minute:
run_and_plot_chromeister.sh HOMSA.Chr.X.fasta MUSMU.Chr.X.fasta 32 1000 4
Comparing the full genome of the Gallus gallus against the Meleagris gallopavo with all their chromosomes including auto-generated grid:
run_and_plot_chromeister.sh GALGA.Chr.complete.fasta MELGA.Chr.complete.fasta 32 1000 4 grid
From the Homo sapiens Chr X vs Mus musculus Chr X
comparison you can inspect the file HOMSA.Chr.X.fasta-MUSMU.Chr.X.fasta.mat.events.txt
which contains the classified events. You can also plot the overlap between detected events and original plot by doing:
python3 chromeister/bin/detect_events.py HOMSA.Chr.X.fasta-MUSMU.Chr.X.fasta.mat.raw.txt png
which will generate the following plot:
Notice that the orange color represents detected events whereas green represents undetected blocks. The detection is based on the Houghs transform and therefore there are some parameters that can be used to change the minimum length to detect an event (e.g. if we want single points as well) or the minimum gap to join two blocks. These parameters can be changed in the detect_events.py
file. The current configuration is mostly tailored to detect larger events and no single points.
Comparing some mycoplasma hyopneumoniae genomes (not all has to be mammalian or plant chromosomes!).
This is afterwards ran with the inmemory_guided_chrom
branch of GECKO (see here) using the following command line (supply the CHROMEISTER hits file):
gecko/bin/guidefastas.sh NC_014.fasta NC_017.fasta hits-XY-NC_014.fasta-NC_017.fasta.mat.hits 1000 50 60 32 names
and we get a csv NC_014-NC_017.csv
which includes:
Frag,1218913,3918445,1221077,3916281,r,0,2165,8644,2163,99.82,1.00,_gi|304372805|ref|NC_014448.1|,_gi|385858114|ref|NC_017519.1|
[...]
Frag,6618091,5316746,6620403,5319058,f,0,2313,8524,2222,92.13,0.96,_gi|321309518|ref|NC_014970.1|,_gi|385858893|ref|NC_017520.1|
[...]
Frag,3062371,6264699,3064642,6266970,f,0,2272,9080,2271,99.91,1.00,_gi|313664890|ref|NC_014751.1|,_gi|392388518|ref|NC_017521.1|
Where the last two columns indicate the origin of the fragment. For instance, in this case, the first fragment belongs to gi|304372805|ref|NC_014448.1|
and gi|385858114|ref|NC_017519.1|
which are sequences 2
and 5
in the x
and y
axis in the previous plot, respectively.
Remember that your csv
file can be uploaded and interactively inspected here!
Hanging output and program does not finish If you experience this kind of output:
[INFO] Generating a 1000x1000 matrix
[INFO] Loading database
100%...[INFO] Database loaded and of length 70039485.
[INFO] Ratios: Q [6.658880e+04] D [7.003949e+04]. Lenghts: Q [66588797] D [70039485]
[INFO] Pixel size: Q [1.501754e-05] D [1.427766e-05].
[INFO] Computing absolute hit numbers.
100%...Scanning hits table.
100%...
[INFO] Query length 66588797.
[INFO] Writing matrix.
[INFO] Found 238693 unique hits for z = 4.
But the program doesnt finish (it "hangs"), then your input sequences probably contain a lot of sequences (i.e. a multifasta with hundreds of contigs). To fix this, simply run the same command but instead of using the "compute_score.R" use the "compute_score-nogrid.R" script. This will remove the drawing of the grid which can get overflown when using too many sequences.
If you use or have used CHROMEISTER in your research, please cite the following article:
Esteban Pérez-Wohlfeil, Sergio Diaz-del-Pino and Oswaldo Trelles. "Ultra-fast genome comparison for large-scale genomic experiments." Scientific reports 9, no. 1 (2019): 1-10.