brankovics/grabb repository overview

⁠GRAbB

GRAbB (Genome Region Assembly by Baiting) is program designed to assemble selected regions of the genome or transcriptome using reference sequences and NGS data.

⁠Table of contents

1. Usage⁠
   1. Installation⁠
   2. Documentation⁠
   3. Examples⁠
2. Prerequisites⁠
   1. mirabait⁠
   2. Edena⁠
   3. Velvet⁠
   4. Seqtk⁠
   5. exonerate⁠
   6. PRINSEQ lite⁠
3. Helper programs⁠
   1. fastq2fasta⁠
   2. get_overlaps⁠
   3. interleaved2pairs⁠
   4. merge_contigs⁠
   5. rename_fastq⁠
   6. single2pairs⁠
   7. uniform_length⁠
   8. fasta_shift⁠
   9. pairwise_alignment_score⁠
   10. reverse_complement⁠
4. Algorithm overview⁠
   1. Main loop⁠
   2. Creating the bait⁠
   3. Baiting⁠
   4. Collecting reads⁠
   5. Assembly⁠
   6. Testing completion⁠
   7. Modes⁠
5. Arguments⁠
   1. ref⁠
   2. bait⁠
   3. reads⁠
   4. folder⁠
   5. prefix⁠
   6. single⁠
   7. min_length⁠
   8. type⁠
   9. arg1⁠
   10. arg2⁠
   11. assembler⁠
   12. clean⁠
6. Using custom assembler program⁠
   1. Adding to the source code of GRAbB⁠
   2. Using external_skeleton⁠
   3. Using SPAdes as assembler⁠
7. Citation⁠
8. Contact⁠

⁠Usage

⁠Installation

Asciinema casts:

• Ubuntu:
• Centos:
• Fedora:

Steps:

1. Install prerequisites (If this step is skipped, then configure_GRAbB.pl tries to use prerequisites included in the package)

   • Minimal set:

     • Baiting program: mirabait (recommended) OR kmer_bait.pl (no installation needed)
     • Read collecting program: seqtk (recommended) OR create_readpool.pl (no installation needed)
     • Assembler: Edena OR Vevlet OR external_scaffold.pl (needs to be modified and requires a working installation of an assembler)

   • Assemblers:

     • Edena: default assembler for GRAbB.pl
     • Vevlet
     • Other assembler: external_scaffold.pl has to be also edited

   • Exonerate:

     • Required for running GRAbB.pl in exonerate mode

2. Configure GRAbB.pl

    ./configure_GRAbB.pl
   

   Copy

   Recommended to add prerequisites to the Path Or add the absolute path of the executables to the source code of GRAbB.pl before running configure_GRAbB.pl

   Configured GRAbB.pl can be found in bin directory

   Bug: On some systems the exonerate binary included in the package runs extremely slow. Configure can get stuck at 'Testing exonerate' block, then issue Ctrl + c. On these systems exonerate has to be installed or build from source code (See Prerequisites::exonerate⁠, before rerunning the configuration script.

3. To test installation run the following (assembler has to be adjusted unless GRAbB.pl is configured without Edena)

    bin/GRAbB.pl --ref for_testing/assembly.fas --reads for_testing/read* --folder test --prefix test
   

   Copy

OR

Use Docker (See Docker.md⁠ for more detailed instructions)

• Either download the docker repository via docker interface docker pull brankovics/grabb

• Or create a local docker image:

    git clone https://github.com/b-brankovics/grabb
    cd grabb/docker
    docker build -t localhost:5000/$USER/grabb .
  

  Copy

⁠Documentation

Run GRAbB.pl without any arguments and it prints the Usage information

The documentation is this file and the files mentioned at the examples⁠.

⁠Examples

See the wiki or the files Docker.md⁠, Examples.md⁠ and Tutorial.md⁠.

⁠Prerequisites

⁠mirabait

1. Download MIRA (4.0) assembler (http://sourceforge.net/projects/mira-assembler/⁠)

2. Extract it. The executable files are in the bin folder.

3. Copy/move or symlink 'mirabait' into somewhere in the path or add the folder to the path (This program uses only mirabait)

   Warning: the name of the executable file has to be mirabait!

   Bug: mira 4.0.2 does not work properly, but mira 4.0 does

OR you may also use kmer_bait.pl⁠, which is less efficient, but uses only perl and standard Unix commands

⁠Edena

1. Download EDENA⁠ and extract it or use the copy in the 3rd_party_programs
2. Change to the directory
3. Type make on the command line (g++ needs to be installed, on ubuntu type sudo apt-get install g++)
4. Copy/move or symlink 'edena' into somewhere in the path or add the folder to the path (The files in the bin folder)

⁠Velvet

1. Download Velvet⁠ and extract it or use the copy in the 3rd_party_programs
2. Change to the directory. First zlib needs to be installed
3. Change to 'third-party/zlib-1.2.3/'
4. Type make on the command line
5. Type sudo make install on the command line
6. Go back to the parent directory (cd ../..)
7. Type make on the command line
8. Copy/move or symlink 'velveth' and 'velvetg' into somewhere in the path or add the folder to the path

⁠Seqtk

1. Download Seqtk⁠ from github and uncompress it or git clone https://github.com/lh3/seqtk.git
2. Change to the directory
3. Type make on the command line (zlib needs to be installed, see 1.3) for instructions)
4. Copy/move or symlink 'seqtk' into somewhere in the path or add the folder to the path (The files in the bin folder)

OR you may also use create_readpool.pl⁠, which is less efficient, but uses only perl and standard Unix commands

⁠exonerate

For Ubuntu run sudo apt-get install exonerate

Else:

1. Download exonerate⁠ from the EBI website and uncompress it or use the version included in the GRAbB package (3rd_party_programs)

2. Change to the directory

3. Type the following commands (The following packages have to be installed on the system before running ./configure: gcc, make and glib2)

    ./configure
    make
    make check
    make install
   

   Copy

4. The executable is found at src/program/exonerate. Copy/move or symlink 'exonerate' into somewhere in the path or add the folder to the path (The files in the src/program folder)

⁠PRINSEQ lite

1. Download PRINSEQ lite⁠ and extract it or use the copy in the 3rd_party_programs
2. Copy/move or symlink 'prinseq-lite.pl' into somewhere in the path or add the folder to the path

⁠Helper programs

⁠fastq2fasta

This program creates a FASTA format read file for each FASTQ read file specified.

Usage:

./fastq2fasta <reads_1.fastq> <reads_2.fastq>

⁠get_overlaps

This program reads the contigs from a fasta file and checks if they are overlapping with each other by using a minimal overlap size that is specified at invocation. Finally, prints all the overlaps found.

Usage:

./get_overlaps <contigs.fasta> <overlap>

⁠interleaved2pairs

This program creates a forward and reverse read file from an interleaved file

Usage:

./interleaved2pairs <reads.fastq>

⁠merge_contigs

This program reads the contigs from a fasta file and checks if they are overlapping with each other by using a minimal overlap size that is specified at invocation. Afterwards it loops through all the contigs and merges contig pairs that only overlap with each other at the given side. In the end it saves the contigs that were created into a file.

Usage:

./merge_contigs <contigs.fasta> <overlap> <output.fas>

⁠rename_fastq

This program creates a FASTA format read file for each FASTQ read file

Usage:

./rename_fastq <reads_1.fastq> <reads_2.fastq>

⁠single2pairs

This program creates paired-end read files from single-end files

1. Selects the reads that are at least as long as the specified length (<int>)
2. Gets the first <int> characters to be used as forward read
3. Gets the last <int> characters to be used as reverse read (reverse complement)

All the produced reads are $len long.

The output file will be created in the current working directory

Usage:

./single2pairs <int> <reads.fastq>

⁠uniform_length

This program creates a read file with reads with uniform length

1. Selects the reads that are at least as long as the specified length (<int>)
2. Trims the reads to the specified length

All the produced reads are $len long

The output file will be created in the current working directory

Usage:

./uniform_length <int> <reads.fastq>

⁠fasta_shift

This program takes a fasta file and shifts the sequence in it according to a position value or a reference file the output is printed to STDOUT

1. Using position value n is the position value

   Usage:

    ./fasta_shift -i <input.fas> -p <int> ><output.fas>
   

   Copy

2. Using a reference file: the output will start with the first occurrence of the sequence in the reference file

   Usage:

    ./fasta_shift -i <input.fas> -r <ref.fas> ><output.fas>
   

   Copy

⁠pairwise_alignment_score

This program compares two sequences and prints a few metrics

Output: <# of identical bases> <# of identical bases, not counting '-'> <percentage of identity w/o '-'>

Usage:

 ./pairwise_alignment_score <input1.fas> <input2.fas>

⁠reverse_complement

This program prints the reverse complement of the input to STDOUT

The input may be files or STDIN

Usage:

 ./reverse_complement <input.fas> ><output.fas>

⁠Algorithm overview

GRAbB is written in Perl and it uses only modules that are part of the core distribution of Perl. In addition to basic UNIX commands, the following third-party programs are used by GRAbB: mirabait (from the MIRA⁠ package), Seqtk⁠, EDENA⁠, Velvet⁠ and Exonerate⁠.

The program is designed to be versatile and flexible with the following functionalities:

• Use pairing information
• Use additional bait sequences
• Assemble multiple regions separately in a single run
• Use any of a range of completion criteria (may use different ones for each region)

These functionalities are detailed below at the appropriate steps of the algorithm.

⁠Main loop

The main loop of GRAbB can be summarized as follows:

1. creating the bait file
2. finding reads by baiting
3. collecting reads
4. de novo assembly of selected reads
5. testing completion

If multi-mode⁠ is selected, then the general baiting step is followed by specific baiting, de novo assembling and completion testing for each of the threads. The threads are generated by splitting the reference file into single-entry FASTA files. These newly created reference files are used as bait for the initial specific baiting steps for the given thread. At the end of each cycle of the main loop, the program checks whether there is any thread that is not completed yet, then it continues or stops accordingly.

⁠Creating the bait

At the invocation of GRAbB it is possible to specify a length filter⁠ that excludes all contigs from the assembly that are shorter than the specified length from being used for generating the bait file for the next iteration.

• Initial bait file

  At the start of the run an initial (general) bait file is created by concatenating the reference⁠ and bait⁠ files. This bait files is used for the first general baiting step.

• General bait file

  At the start of each iteration a general bait file is created by concatenating the latest assembly file(s)

• Specific bait file

  The specific bait file is only created in multi-mode⁠. At the start of the run the reference file is split into single-entry FASTA files and these are used as initial specific bait files. In latter iterations the latest assembly for the give thread is used as the new bait file.

⁠Baiting

Reads belonging to the specified sequence are identified by using exact k-mer (31 bp) matching that is implemented by mirabait (from the MIRA⁠ package). The names of the reads thus identified, are collected and added to the list of read names from previous iterations (There is a separate list for the general baiting and for each thread.) If there are no new reads identified, then the program stops the iteration, either for the given thread (specific) or for all the threads (general).

By using a general baiting step before the specific baiting it is possible to reduce the required run time, because the large input read file(s) is/are only screened once per iteration and the specific baiting is confined to screening the reads that are already identified to be specific during the general baiting.

⁠Collecting reads

The identified reads are collected from the read files using Seqtk⁠ into internal read files. The program identifies reads based on the read names, thus the first word of the identifier line should be unique for each read if they are single-end or should be the same for both reads of the pair.

⁠Assembly

The program can use two assemblers, EDENA⁠ and Velvet⁠, by default, but there is a skeleton code to add a new assembler to the source code of the program. Also it is also possible to write an external perl script⁠ that is used by GRAbB for the assembly. By default EDENA⁠ is used for assembly, but using command line options the other assemblers can be selected, as well.

The fact that single⁠- or paired-mode⁠ is selected is passed on to the assembler program that assembles the specific reads de novo. Also, it is possible to pass additional arguments, such as overlap size (EDENA⁠), to the assembler program at the invocation of the main program.

⁠Testing completion

There are multiple completion criteria that can be specified for the program. These can be specified for each thread separately (multi-diff-mode⁠ or for all of the threads at once. Also, it is possible to specify multiple criteria for the same thread or run. In this case the program stops when any one of these criteria is met.

1. Exhaustive run

   The first completion criterion is implicit, the program stops if there is no new information found. This means that either there are no new reads found or that the new assembly is identical to the bait used for the current iteration. If no other criterion is specified then the run is an exhaustive run, since it iterates until it cannot find any new information.

2. Length criteria (total length, longest contig or N50)

   This is the simplest explicit completion criterion. There are three options that can be used for this setting: total assembly size, the length of the longest contig or the N50 value of the assembly. This criterion is tested independently for each of the threads in multi-mode⁠ or, otherwise, for the single thread. As mentioned before, multiple criteria can be used in a single run, this also applies for the different size criteria. These settings are useful when exploring the vicinity of a specified sequence region.

3. Matching homologous sequence

   In this case the specific reference sequence is used to identify the homologous region within the assembly. To identify the matching region, GRAbB uses Exonerate with settings that ensure that the whole reference sequence is aligned to the assembly contigs. This makes it possible to match sequences that are somewhat dissimilar to the reference and may also contain indels, causing gaps in the alignment. The possible results of the matching can be divided into two groups: the whole sequence is matched—also if there are internal regions that correspond to gaps in the other sequence—or it is not. In the first case the completion criterion is met, and the matched region is extracted from the assembly in the same orientation as the reference sequence and saved to an output file. In the latter case there are two possibilities. If the matched region is larger than the one for the previous iteration, then the thread or run will continue, but if the size of the match has not improved than the thread or run will stop.

⁠Modes

1. Single-⁠ or paired-mode⁠

   This information is passed on to the assembler program.

   • Paired-mode

     It is the default if two read files are specified.

     In paired-mode the read files are tested whether GRAbB can identify pairs as expected⁠.

   • Single-mode

     It is the default mode if there is only one or more than two read files specified. If two read files are specified than single-mode can be chosen by using the --single⁠ argument.

2. Multi-mode

   Multi-mode is selected by using --type multi⁠ option at invocation. In multi-mode the reference file is split into single-entry FASTA files, these are referred to as specific references files. Also for each entry a separate thread is created. The individual threads are independent from each other, thus multiple regions can be assembled in a single run without interference from each other. The specific reference file is used as initial specific bait file. Also in exonerate-mode⁠ the specific reference file is used for the homology matching.

3. Multi-diff-mode

   Multi-diff-mode is selected by using --type multi-diff⁠ option at invocation. It also belongs to the multi-mode with all its properties. The difference is that when the specific reference files are created if there is completion criterion specified in the identification line of the FASTA entry, than this is added to the completion criteria to be used for the given thread.

4. Exonerate-mode

   Exonerate-mode is selected by using --type exonerate⁠ option at invocation. In this case the specific reference sequence is used to identify the homologous region within the assembly. To identify the matching region, GRAbB uses Exonerate with settings that ensure that the whole reference sequence is aligned to the assembly contigs. This makes it possible to match sequences that are somewhat dissimilar to the reference and may also contain indels, causing gaps in the alignment. The possible results of the matching can be divided into two groups: the whole sequence is matched—also if there are internal regions that correspond to gaps in the other sequence—or it is not. In the first case the completion criterion is met, and the matched region is extracted from the assembly in the same orientation as the reference sequence and saved to an output file. In the latter case there are two possibilities. If the matched region is larger than the one for the previous iteration, then the thread or run will continue, but if the size of the match has not improved than the thread or run will stop.

5. Clean-mode

   • Clean-mode

     It is selected by --clean⁠. GRAbB will remove some internal files to save disk space. But there is no information lost, because all the deleted files can be reconstructed using the remaining files.

   • Double clean-mode

     It is selected by --clean --clean⁠. GRAbB will remove some internal files to save disk space. At the end of the run all the output files and folders are deleted except for the result files.

⁠Arguments

GRAbB.pl --ref <reference file> --reads <read file 1> [<read file 2>] --folder <directory> --prefix <prefix> [options]

There are four mandatory arguments: reference file⁠, read file(s)⁠, output folder⁠ and prefix for the log and results⁠

The order of the arguments is not important.

⁠ref

--ref <ref.fas>

The reference file is a FASTA formatted file that contains one or more sequences. The sequence IDs have to be unique for each sequence (as required by mirabait). If the file contains multiple sequences and the program is run in multi-mode⁠ then the reference file is split into separate reference files that contain only a single sequence, the handling of these files is discussed in the segment on the main loop⁠. Furthermore, the description lines may contain specification for the completion criterion to be used for the given sequence that is used if multi-diff mode⁠ is selected. Because the read selection is based on exact k-mer (31 bp) matching, the reference sequence does not have to be highly similar to the target sequence.

⁠bait

--bait <bait.fas>

A separate bait file can be specified besides the reference file, this file together with the reference file will be used as first bait⁠. Useful when using special criterion for the assembly, such as homology⁠.

⁠reads

--reads	<r1.fastq> [<r2.fastq> ...]

Multiple read files can be specified as input. If two read files are given, then it is assumed that reads are paired⁠, but in single-mode⁠, reads are considered as single reads. The program identifies read pairs based on the read names, thus the first word of the identifier line should be the same for both sequences. The read files may be in FASTA or FASTQ format and may be compressed (using gnuzip). If EDENA⁠ is selected as assembler program, then all the reads should be of the same length.

⁠folder

--folder <folder_name>

The directory where all the output will be saved. If the directory is non-empty then the files it contains can be used like internal files. In this manner previous runs can be continued, make sure to remove or replace files that would suggest completion:

Folder structure:

<folder>/

• reference.fas

• bait.fas

• extra_bait.fas

• <prefix>_assembly_thread_<int>.fas

• <prefix>_result_thread_<int>.fas (if exonerate-mode is selected and the sequence was matched)

• <prefix>.log

• old_collection.list

• reads*<int>.fastq or reads<int>*.fasta

• Round*<int>*/

  • hashstat.bin

  • mirabait.log

  • new_collection.list

  • positive_<int>.txt

  • readpool*<int>*.fastq

  • reads*<int>*.fastq

  • thread_<int>/

    • assembly.fas
    • new_collection.list
    • readpool*<int>*.fastq
    • (files or folders generated by the assembler)
    • exonerate.log (if exonerate-mode is selected)
    • result.fas (if exonerate-mode is selected and the sequence was matched)

• thread_<int>/

  • assembly.fas
  • assembly_<int>.fas
  • bait.fas
  • final_assembly.fas
  • old_collection.list
  • reference.fas
  • reference.fas.exonerate (if exonerate-mode is selected)
  • result.fas (if exonerate-mode is selected and the sequence was matched)

⁠prefix

--prefix <prefix_of_output>

The prefix for the output files:

• log file
• assembly file
• result file

⁠single

--single

Treat reads as unpaired reads even if two read files are specified

⁠min_length

--min_length=<int>

Minimum size required for a contig to be included for completion testing and baiting