====================================================================================
This software is suitable only in linux-like systems (Unfortunately not Windows or MAC)
- Purpose
- Installation
- Before-launching-the-workflow
- How to use
- Input data
- Details of the worfklow and outputs
- Working examples
- SLURM integration
- CPU and memory requirements
I - Perform TE and gene prediction
II - Identify synteny blocks and rearragements
III - Plot dS along the genome
IV - Perform changepoint analysis to identify evolutionary strata
We recommand mamba over conda as the later can be very slow. Alternatively edit the script in INSTALL/dependencies.sh to replace calls to mamba by conda.
Clone the workflow, then please work from within it to preserve the architecture.
We recommend that you clone the pipeline for each of your new project and work within it, and to keep all projects separated otherwise it will be difficult to recover your results.
All options and full paths to input files must be set in the config file provided in : config/config .
PLEASE, carefully read the user guide below before any attempt at running the workflow.
Note that this workflow uses several softwares and packages, notably BREAKER, TSEBRA, GeneSpace, PAML, R mcp. We recommend that you read their corresponding manuals before launching the workflow.
The first step is to provide the path to your input files and choose settings in the config file. An example config file is provided here
To launch the workflow, simply run:
./master.sh -o X #with X an option from 1 to 8.
There are several options which allow you to choose which steps of the workflow you wish to run. This allows the workflow run from any step in the process. In case of bug you may restart it from whenever it crashes (after fixing the bug) and it should work smoothly.
./master.sh --help #to see all options
./master.sh -o 1 2>&1 |tee log
All steps: performs all steps of the workflow, i.e. gene prediction, synteny analysis with GeneSpace including single copy orthologs inference between sex/mating type chromosomes, synonymous divergence (dS) computation, evolutionary strata inference and production of various plots
The following options allow you to run only certain parts of the workflow.
./master.sh -o 2 2>&1 |tee log
Steps I and II: performs only gene prediction and synteny analysis with GeneSpace (no dS computation or evolutionary strata inference)
./master.sh -o 3 2>&1 |tee log
Steps II to IV: performs synteny analysis with GeneSpace and subsequent analyses : useful if you already have annotated your genome (either from running this pipeline or any other annotation tools)
./master.sh -o 4 2>&1 |tee log
Steps III to IV: performs dS computation and subsequent analysis : useful if you already ran the synteny analysis with GeneSpace, and for customizing the plots produced at step III
./master.sh -o 5 2>&1 |tee log
Step II: performs only the synteny analysis with GeneSpace
./master.sh -o 6 2>&1 |tee log
Step I: performs only gene prediction
./master.sh -o 7 2>&1 |tee log
Step IV(G and H): performs only evolutionary strata inference and the production of various plots: useful if you already ran the synteny analysis with GeneSpace and the dS computation with PAML. This option is useful and recommanded to explore various parameter settings in the MCP analysis, for instance adding priors or tweaking the order of the scaffolds.
./master.sh -o 8 2>&1 |tee log
Step IV(H): performs only the plots subsequent to dS computation: useful if you already ran the rest of the workflow and want to customize your plots
| Option: | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
|---|---|---|---|---|---|---|---|---|
| I.Gene prediction | X | X | X | |||||
| II.Orthology and synteny (GeneSpace & Minimap2) | X | X | X | X | ||||
| II.Synteny plots | X | X | X | |||||
| III.dS computation + plots | X | X | X | |||||
| IV.Evolutionary strata inference | X | X | X | X | ||||
| IV.Evolutionary strata plots | X | X | X | X | X |
The different steps of the worflow are detailed below
/!\ The input required will vary strongly based on which steps of the workflow you want to perform. Several files are compulsory
Again, all input data, including full path to input files, should be provided in the config file
all options
-
Input genome(s) - compulsory: This may be one genome assembly containing both sex/mating type chromosomes, or ideally two separate haplotype assemblies containing each one of the sex/mating-type chromosomes.
-
list of scaffolds - compulsory: names of the contigs/scaffolds/chromosomes composing the sex/mating-type chromosomes.
a tab delimited three column file containing: -
column1: name of the haplotype/species
-
column2: name of the sex/mating type chromosome(s)
-
column3: whether the chromosome(s)/contig(s) is(are) kept in Normal order (N) or reversed (R)
see example here
⚠️ All entries must belong to haplotype 1 only. Do not include chromosomes from haplotype 2. The plotting scripts assign all entries to haplotype 1 internally, so mixing both haplotypes will cause an error inplot_riparian.
- ancestral genome - optional but highly recommended: The genome assembly of a species used as a proxy for the ancestral state. This will allow to plot dS along 'ancestral' gene order, and to infer more accurately single copy orthologs.
- ancestral gene prediction - compulsory with ancestral genome: gene prediction associated with the ancestral genome. format: gtf/gff(.gz)
names of fasta and contigs/scaffolds/chromosomes:
We recommend short names for genome assemblies and NO SPECIAL CHARACTERS apart from underscore.
example: species-1.fasta will not be valid in GeneSpace. => Use species1.fasta instead.
❗ chromosome/scaffold IDs: ❗
you MUST use standardized IDs including the species/individual name NO SPECIAL CHARACTERS apart from underscore.
example: species1_chr1 or species1_contigX or species1_scaffoldZ otherwise the code will failed during renaming steps
gene_id MUST follow this structure:
[individualID][chromosomeID][geneID]" :
-
avoid any special character in the ID
-
use only two underscore as above.
-
[individualID] : any ID for you species/strain/individual of interest
-
[chromosomeID] : ID of the chromosome should be like "chrX", "contigZ", "chrW" etc
-
[geneID] : anyID avoid complex characters
options 1,2,6 if your input genomes are not already softmasked
- TE database - compulsory: the name of the TE database (some are available online depending on your taxon)
- NCBI taxon - compulsory: a taxon name for NCBI (used with repeatmasker)
- TE bed files - optional: a pair of bed files containing TE for your region of interest if already available (will be displayed on the circos plots)
options 1,2,6 if your input genomes are not already annotated
- BUSCO lineage name - compulsory: name of the BUSCO lineage corresponding to your species (the list of busco lineages is available with busco --list-lineage)
- RNAseq - optional: RNAseq data for each genome, will improve BRAKER annotation
- Protein database - optional: a database of proteins from related species. Alternatively, orthoDB12 can be used (downloaded automatically)
- orthoDB12 lineage name - optional: one of "Metazoa" "Vertebrata" "Viridiplantae" "Arthropoda" "Eukaryota" "Fungi" "Alveolata"
Full details on the options in config file are listed in config folder For an example of input files, we provide an example data folder
| Operation | Tools | data type |
|---|---|---|
| I.A.read trimming | Trimmomatic | RNAseq |
| I.A.read mapping | gmap/gsnap | RNAseq |
| I.A.sorting read | samtools | RNAseq |
| I.A.mapping quality assement | samtools + R | RNAseq |
| I.B.TE detection and softmasking | RepeatModeler + RepeadMasker | genome assembly |
| I.C.genome annotation | BRAKER + tsebra | genome assembly + protein + database |
| I.C.quality assessment | BUSCO + Blast + Inter Pro | genome prediction |
| II.D1.whole genome alignement | minimap2 | genome assemblies |
| II.D1.gene microsynteny | R | single copy orthologs |
| II.D2.Synteny and orthogroups | GeneSpace (including OrthoFinder/MCScan) | gene prediction and proteins |
| III.E.cds alignement | muscle + translatorX | gene prediction (single copy orthologs) |
| III.F.dS computation | paml | CDS alignment |
| III.F. dS plot/CIRCOS plot | R | Ds and genome information |
| IV.G.changepoint analysis | R | dS values and gene order |
📝 Corresponding script: 00_scripts/launch_rnaseq.sh
- Reads trimming using trimmomatic, automatic detecting whether the data is Single-End or Paired-End
- Creation of database for gsnap using gmap
- Alignment using gsnap
📝 Corresponding scripts:
00_scripts/03_gsnap_SE.shfor Single-End or00_scripts/03_gsnap_SE.shfor Paired-End - Mapping quality assessment and plot (sequencing depth and MAPQ)
plots and output:
haplo1/04_mapped/Depth/andhaplo2/04_mapped/Depth
Figure 1: example RNAseq depth
📝 Corresponding script: 00_scripts/launch_step05_to_08.sh
- De novo repeat identification using repeatmodeler on the input genome(s)
- Genome masking using repeatmasker with ncbi TE Library, and custom dataset if provided
📝 Corresponding script: 00_scripts/06_braker.sh
- Five successive rounds of gene prediction based on OrthoDB (and custom protein database if provided) with BRAKER
- If RNA-seq data was provided: one round of gene prediction based on RNA-seq data with BRAKER
- Quality assessment and reports production for each round of gene prediction
This uses the raw BRAKER hintsfile, and BUSCO (:pencil: corresponding script:
00_scripts/07_busco_after_braker.sh)
The report includes number of genes, number of introns per gene, gene support, number of complete genes and various histograms useful for error checking. - If RNA-seq data was provided: combination of the best protein-based gene model and the RNA-seq-based gene model using TSEBRA
Warning: TSEBRA parameters intron_support and stasto_support are set to 0 in this workflow (default in TSEBRA: 1 and 2 respectively). This means that only overlapping genes between the two gene models will be filtered. You can change this parameter and others to adjust to your desired level of stringency in the TSEBRA config file:
config/default.cfgPlease read the TSEBRA manual before running the script. - If no RNA-seq data was provided: genome annotation is set as the best protein-based gene model, as evaluated with BUSCO
- Reshaping of genoma annotation: scaffold name is inserted in gene names (facilitating dowonstream analyses) and the longest transcript of each gene is kept (necessary for single copy ortholog identification)
:pencil: Corresponding script:
00_scripts/08_braker_reshaping.sh - Final genome annotation quality assessment using Blast against Uniprot database, BUSCO, and optionally InterProScan Warning running InterProScan is time-consuming and can take up to several days. By default the option is turned off in this workflow. It can be turned on in the config file.
📝 Corresponding script: 00_scripts/11_run_genesSpace_paml_ideogram.sh
- Alignment between the two haplotypes using minimap2 If the scaffolds of interest were provided in the scaffold table, only those will be aligned. Otherwise, the whole haplotypes will be aligned.
- **If ancestral genome was provided: alignment between the two haplotypes and ancestral genome using minimap2
- Construction of synteny plot on the focal scaffolds using pafR
- If no scaffold table was provided: construction of whole genome dotplot using pafR
Figure 2: pafr synteny plot based on minimap between the focal scaffolds
Figure 3: dotplot between the two haplotypes
- Launching of GeneSpace
- Identification of single copy orthologs with OrthoFinder
- If no scaffold table was provided: Production of dotplot and a riparian plot for the whole genome - Production of a riparian plot on focal scaffolds
For more information, consult the GeneSpace readme. - Circos plots production using the R package circlize, between the two haplotypes, and the ancestral genome (if provided).
- Micro-synteny plot production using the R package Rideogram, examples below
📝 Corresponding script:
00_scripts/Rscripts/04.ideogram.R
Figure 4: Synteny plot by GeneSpace showing gene synteny between the ancestral genome (ancestral_sp) and the two mating types of Microbotryum lychnidis-dioiciae
- Align all coding sequences from the focal scaffolds using muscle from within TranslatorX
- Calculation of dS (and dN) using PAML Note: PAML fails if special characters occur in the input fasta file, or if the length of a gene name is above 32 characters. To prevent this, we implemented an automatic renaming procedure to shorten character names and remove special characters. After PAML, genes are converted back to their original names, so this should be transparent to users.
- Plotting of dS values on the focal scaffolds (and on the whole genome if no scaffold table was provided), using a custom R script.
:pencil: Corresponding script:
00_scripts/Rscripts/03_plot_paml.RIf an ancestral genome was provided, dS is plotted along that genome, otherwise it will be plotted along one of the two sex/mating type chromosomes.
Note: It is possible to modify this R script to adapt the plotting options to your needs (for instance position and direction of scaffolds).
Figure 5: A) dS values along the position in the ancestral chromosomes. B) dS values along the ancestral gene order after inverting the chromosomes and removing the large autosomal part on contig 8.
C) and D) gene rank as a function of ancestral gene order, in genome1 and genome2 respectively.
- Plotting of circos plots of the focal scaffolds, with links between single copy orthologs genes, using the R package circlize
dS Corresponding script:
00_scripts/Rscripts/05_plot_circos.RTE and gene density bed files can be provided as arguments. By default fused autosome will be plotted but these can be removed from the contig list. Note: It is possible to modify this R script to adapt the plotting options to your needs (for instance position and direction of scaffolds).
put figure4 panel B here Figure 6: Circos plot between the ancestral genome and one of the focal haplotype (left part) and circos plot between the two focal haplotypes (mating types). The most external track displays the position of gene of interest (red, green and light blue) and the centromeres (in purple). The middle track in lightblue displays gene density. The most interior track in green displays TE density. Red and darkblue interior links display single copy orthologs links. Note that the same figure is produced, with inner links colored according to discrete quantile values of dS example here
📝 Corresponding script: '00_scripts/Rscripts/06.MCP_model_comp.R`
** Warning: This step is automatically performed with default gene order. However, it is unlikely that the default gene order is correct, so we strongly suggest that you consult the results of the workflow, especially the dS plot. Once you have deciphered clear hypotheses regarding scaffold order and orientation, you may change them manually in the "scaffold" table, then relaunch the workflow at step IV:
bash ./master.sh -o 7
- Launch of multiple change point (MCP) analysis based on the dS values, with the R package mcp This allows the automatic inference of evolutionary strata through the identification of regions with distinctively different levels of dS. Note that regions with low or zero dS on the leftmost and rigthmost extremities of the chromosomes of interest usuall correspond to the pseudo-autosomal regions and are therefore not true evolutionary strata. - Running MCP analyses from 1 to 8 changepoints - Plotting changepoints posterior distribution and models convergence
Figure 7: Results from the changepoint analysis for 3 (panel A) to 8 changepoints (panel F). Each changepoint panel displays the distribution of dS values (black dots) along with 25 draws from the joint posterior distribution (grey lines) and 95% highest density interval (red lines). Posterior distributions of the changepoints are shown in blue with one line for each chain.
Figure 8: posterior fit of the model parameter and mixing of the chains. Here only 7-changepoint model posterior fit is shown.
- Export of summary tables The MCP analysis produces many useful information which are extracted and automatically exported in .txt tables:
- modelX.chpt.txt (with X being the number of changepoints tested, from 1 to 8): Contains the output of the summary function from MCP.
The columns correspond to:
1. name: name of the parameter (changepont and interval)
2. mean: mean value of dS and interval (based on gene order)
3. lower/upper: lower and upper boundaries
4. Rhat: Gelman-Rubin convergence diagnostic which is often taken to be acceptable if <1.1
5. n.eff: effective sample size computed using effectiveSize. Low effective sample sizes can also be revealed by poor mixing in trace plots.
- modelchoice.txt : Contains information from the loo model choice operation.
The columns correspond to:
1. elpd_diff
2. se_diff
3. elpd_loo
4. se_elpd_loo
5. p_loo
6. se_p_loo looic
7. se_looic
- weights.txt : Contains the weights of each tested models. Higher weights indicates higher supports.
- HypothesisXstrata.txt (with X being the number of changepoint tested, from 1 to 8): Contains results from hypothesis testing (BayesFactor and posterior probabilities aiming at testing differences among strata).
Differences in dS values among adjacent strata are tested, from left to right in the gene order. Both directionalities are tested, i.e. both the hypothesis that stratum 1 has higher dS than stratum 2, and the hypothesis that stratum 2 has higher dS than stratum 1 are tested.
- classif.sX.$haplo1.$haplo2 (with X being the number of changepoint tested, from 1 to 8): Contains the assignment of single copy orthologs to strata
1. column1: gene name in the first haplotype
2. column2: gene name in the second haplotype
3. column3: stratum of appartenance
- df.txt : Contains a summary of all informations
These figures are produced for each number of changepoints tested
- Violin plot of the dS for each evolutionary strata inferred by the model with results of statistical tests
- dS plot along the ancestral gene order, colored according to strata as inferred by the model
- dS plot along the ancestral gene position, colored according to strata as inferred by the model
- Ideograms between the ancestral genome and the haplotypes, colored according to strata as inferred by the model
- Circos plots adding to the exterior track the genes colored according to the strata as inferred by the model
Figure 9: example violinboxplot for the two "most likely" models inferred by the loo analysis. By default plots are constructed for all models (from 1 to 8 changepoints). Default statiscal test from the ggstats plot package are used assuming parametric tests.
Figure 10: dS values plotted along the ancestral gene order for all possible models from three to eight changepoints
each point is a gene dS value colored according to the strata of assignation.
Figure 11: dS values plotted along the ancestral genome for all possible models from three to eight changepoints
each point is a gene dS value colored according to the strata of assignation
Figure 12: example ideograms infered for the most likely models here. Links are colored according to their strata of appartenance.
Figure 13: circos plot with external links based on inferred strata
Use case: when a run of orthofinder was performed and included many species or sample.
this can reduce the number of false positive single copy orthologs
(to do):
beta testing
the workflow is being deployed on SLURM for testing.
each major step can be run separately see the scripts here
- arrays are implemented for read mapping
- arrays are used for parallel BRAKER run
ideally, once the config/config file is set properly and once you have chosen which option to use you can simply launch:
./00_scripts/slurm_code/runall_slurm_workflow.sh <OPTION>
with OPTION being a number between 1 to 8 depending on the chosen configuration for your data.
see the summary of options here
this will submit all job and their dependencies, making appropriate check depending on whether you have RNAseq data or not, if you have already a softmasked genome, etc.
NOTE: for now only option 1 (full workflow), option 6 (only TE and gene annotation) and option 3 have been tested!
this workflows requires very different CPU and memory settings depending along the way.
Limits: unfortunately EASYstrata was not developped on a cluster implementing a scheduler like slurm so memory and CPU usage is not always optimized. For instance job array could be used for efficient and faster RNAseq mapping in a multisample cases.
Similarly job array could be used to parallelise gene alignment along with dS computation.
with this limit in mind some hidden variable enable the control of memory and cpu number in a non-slurm cluster:
in the config/cpu_mem file you'll find the following parameters:
NCPUS_TRIMMO=8
NCPUS_GSNAP=24
NCPUS_BRAKER=24
NCPUS_REPEATEMODELER=20
MEM_TRIMMO=XX
MEM_GSNAP=XX
MEM_BRAKER=XX
MEM_REPEATMODELER=XX