iMAGine is a metagenomic workflow which includes filtering, assembling, and binning.
This workflow includes the following tools which need to be installed in the system.
You can get the latest stable release here.
git clone https://github.com/avishekdutta14/iMAGine.git
chmod a+x imagine.sh
chmod a+x dependency-check.shMost of the tools have their own set of installation protocols and dependencies. They can be obtained from the links mentioned above. This pipeline requires an output of each tool to be piped into the next tool (except for QUAST and the CheckM- since this is the last tool). For checking the dependencies are installed and declared to the $PATH, you can run the following command
./dependency-check.shYou can install all of the dependencies of the dependecies from imagine.yml. You'll still need to install the specific dependencies listed above.
conda env create -f imagine.ymlThen, to enter the environment:
conda activate imaginedocker pull avishekdutta14/imagine:v0.2
docker run --rm -v "$(pwd)":/dir -w /dir imagine /root/iMAGine/./imagine.sh -r sample-R1.fastq.gz -R sample-R2.fastq.gz -s sample_name |& tee -a imagine.txtAll the tools can also be run separately inside the imagine docker container. The standalone version of emapper (for annotation) is integrated inside the image, but due to large sizes of the databses they are not present inside the image. Please download the database and declare it in the PATH before using emapper.py Details.
Important consideration: Create a new folder for each sample and run iMAGine from inside the folder. Modify the imagine.txt filename with a unique modifier to avoid overlaps
./imagine.sh -r sample-R1.fastq.gz -R sample-R2.fastq.gz -s sample_name |& tee -a imagine.txtsample_name is the name used for all the intermediate files and also serves as a prefix for the final bin names
./imagine.sh -hAll the parameters set in the shell script are optimized for specific marine system metagenomes. These parameters should be tested and tried while working with metagenomes from other systems/environments. Some thoughts about selecting certain parameters are shared in the following section How it works? and a detailed description for parameter selection can be obtained from the manuals of each aforesaid tool. Please be cautious about selecting and changing the number of threads selected for each process, as this will vary among systems.
Many of the steps/processes/tools use multi-threads for parallelization, and the user can modify it accordingly. This increases resource usage and also decreases the time needed to run processes. This is highly desirable, but it comes with a sacrifice. Many processes used in this pipeline uses probabilistic models. Probabilistic models and many other machine learning models are impacted by the number of threads used since the probabilities and predictions can vary while computing in different threads/CPUs. But, the good news is that these variations are not significant. The results generated will be different (unfortunately), but the difference is small (fortunately). To understand and identify the steps in which you can expect variations are mentioned in the How it works? section. Morover, the default parameters worked well for specific marine environments and should be modified for other environments. Using default parameters will not gurantee the best result for other environmnets or system (or even other marine systems). To know more on the parameters please consider reading How it works? section.
Bioinformatics tools are evolving fast, so future versions may not be compatible with iMAGine. This tool works with Python version: 3.6.8, fastp v0.21.0, SPAdes version: 3.15.2, QUAST 5.0.2, bwa 0.7.17-r1188, samtools version 1.12, metabat2 (latest version as of May-21), and CheckM v1.1.3. If you feel that any future version of these tools didn't had any major modification with terms of parameters used you can use the future versions.
Quality assessment and filtering of raw reads are done using fastp. The user can play around with a lot of parameters and also modify iMAGine accordingly. iMAGine uses default parameters along with -e 30 for removing all the reads whose average quality score is less than 30. This is done to remove low-quality reads since the assembling using metaSPAdes is very sensitive to low-quality data.
:warning: The -w (number of threads used) is kept to 1 to make the run deterministic (reproducible). The user can change -w to higher values but should be aware of non-reproducible results (the changes in result will be very low and not significant) :zebra:. Since this process does not take a huge time, setting the number of threads to 1 will not increase the run time hugely.
Assembling will be done using metaSPAdes. This is the most time-taking process. The selection of k-mers can be an important thing, but I wouldn't waste too much time on it. The default k-mer values are set to 21,33,55. Selection of k-mers should be governed by sample type (natural or engineered environment), the complexity of microbiome (high or low diverse :confused: - I know these words are relative), size of the sequenced reads, and other parameters (I might have missed some :worried:). My perspective (might vary among researchers and bioinformaticians), for high complex metagenomes from a natural environment having 150 bp paired-end sequencing, the default parameters (21, 33, 55) should work well; and if you think that the system/environment is defined and are low-complex and the depth of sequencing is very high, you can increase the k-mer (I even went up to 127 with 21, 33, 55, 77, 99, 127 for an engineered environment and got good results :sunglasses:). You can play around with other parameters, and also learn more about spades here.
:warning: The resource usage should be set depending on the resource available and the data size. The maximum threads (-t) and RAM (-m) are set to 20 and 300 GB. The good thing is that the actual amount of RAM consumed will be lower for many steps. There are two main steps done by SPAdes assembler viz. 1. Read correction and 2. Assembly. For assembly, similar number of threads should give reproducible results, but for the Read corrections (done using BayesHammer), I found that even for a similar number of threads, the results might vary. I feel that this variation is governed by resource usage by other processes on the computer/server system that is being used. One can stop the read correction step and only use the assembler to get reproducible results, but using read correction results in high-quality assemblies (developer recommendation). I feel that one should use the read correction step since the results are not very different.
The qualities of the assemblies are checked using QUAST. There are different outputs in QUAST that can be used for assembly assessment. N50 values, L50 values, longest contigs, the size distribution of contigs are some of the parameters to look for to understand the assembly quality. One can also use metaQUAST, which has some enhanced features, but for the sake of resource usage and run time, the default is set to QUAST. Details of QUAST can be found here
Binning is done using metabat2, and the mapping/alignment of reads is done using bwa to generate alignment files, which can be further converted to an abundance profile using jgi_summarize_bam_contig_depths (a tool present in the metabat2 repository). bwa mem (details here) is used for alignment. The user can also use bowtie2 for this. If someone wants to use bowtie2 for the alignments, please contact me. I can modify the script accordingly.
Samtools is used convert .sam file to .bam file and to sort the .bam file since the input of jgi_summarize_bam_contig_depths should be a sorted .bam file. There are also many other piped samtools syntaxes available which can be replaced here. For the clarity of the script and breaking things clearly, conversion of sam to bam and sorting bam file is done in two separate steps.
Then comes metabat2 (finally!!!!), which will help to bin based on abundance and tetranucleotide frequency of the contigs. Though the abundance profile is optional, I would recommend using the abundance profile for better binning. I have changed the minimum size of a contig for binning (-m) from 2500 (default) to 1500 since I wanted to consider the smaller contigs (but if you want to stick with the default value, you can remove the -m flag). For more details please refer the metabat2 manual. Here comes the good part :grinning:, thanks to the developers :pray: that metabat2 can be reproducible, though it has some sections to pick things (again probabilities/predictabilities). You can set seed values to make it reproducible. I chose a random number 1234 (it came to my mind when writing the script) as a seed for the script. You can modify it if you want.
Last but not least is checking the bin qualities. Bin qualities are checked using the CheckM tool. Bins can be referred to as high-quality draft (>90% completeness, <5% contamination), medium-quality draft (>50% completeness, <10% contamination) or low-quality draft (<50% completeness, <10% contamination) MAGs as suggested by Bowers et al., 2017.
mag_extract.py helps in segregating the bins based on quality, as mentioned in Bowers et al., 2017. But it only segregates the bins based on completeness and contaminations, and not on basis of the presence of the 23S, 16S, and 5S rRNA genes and at least 18 tRNAs, which is one of the criteria mentioned in Bowers et al., 2017 for high quality drafts. For running it, place the script in the folder where the bin output folder are present. Rename the bin containing folder to bins_dir (If iMAGine is used, the folder will be already named as bins_dir/). Also, place the checkm quality (i.e., bin_stats_ext.tsv) information out to the folder in which mag_extract.py is present. The required checkm input can be obtained from the storage folder (checkm_ouptut/storage/bin_stats_ext.tsv) under the checkm output (i.e., bin_stats_ext.tsv).
This script requires python3, pandas, and numpy (python libraries).
chmod a+x mag_extract.py
./mag_extract.pymag_extract.py can also be added to PATH, but the script should be implemented in a folder where bins_dir/ and bin_stats_ext.tsv is present
chmod a+x mag_extract.py
./mag_extract.py -cs [% of completeness] -cn [% of contamination]
./mag_extract.py -cs 90 -cn 5mag_extract.py can also be added to PATH, but the script should be implemented in a folder where bins_dir/ and bin_stats_ext.tsv is present. Now the script also gives you average coverage information for each bins.
Kaiju can be used for this. For offline use, please visit the Kaiju GitHub link. Do not forget to download the reference database for Kaiju if you are using the offline version. In this step, you will get an idea of the taxonomic distribution of different domains of life (depending on the database you are using you can also get a rough idea of the eukaryotic population distribution).
If you are using any other downstream annotation tools mentioned in this documentation, please cite them in the literature.
This tool and the suggested workflow is strictly intended for educational and academic purpose. The results for different tools might vary depending on the input dataset.