soapdenovo(1)

soapdenovo πŸ’…βœ¨

okay bestie, gather round, because this is the manual for soapdenovo β€” a short-read assembly method that can build a de novo draft assembly. Yes, that's the whole NAME line, and yes, it's a lot. We're gonna get through it together. 🫢


SYNOPSIS 🌈

soapdenovo_31mer soapdenovo_63mer soapdenovo_127mer

Three flavors, one legend. Pick your kmer girlie and go. πŸ‘‘


Introduction πŸ¦‹

SOAPdenovo is a novel short-read assembly method that can build a de novo draft assembly for the human-sized genomes. It's specially designed to assemble Illumina GA short reads. The way it creates new opportunities for building reference sequences and carrying out accurate analyses of unexplored genomes in a cost effective way?? Iconic. πŸ’–

  1. Support large kmer up to 127 to utilize long reads. Three version are provided.
    • I. The 31mer version support kmer only <=31.
    • II. The 63mer version support kmer only <=63 and doubles the memory consumption than 31mer version, even being used with kmer <=31.
    • III. The 127mer version support kmer only <=127 and double the memory consumption than 63mer version, even being used with kmer <=63.
Please notice that, with longer kmer, the quantity of nodes would decrease significantly, thus the memory consumption is usually smaller than double with shifted version. 🧠✨
  1. New parameter added in "pregraph" module. This parameter initiates the memory assumption to avoid further reallocation. Unit of the parameter is GB. Without further reallocation, SOAPdenovo runs faster and provide the potential to eat up all the memory of the machine. For example, if the workstation provides 50g free memory, use -a 50 in pregraph step, then a static amount of 50g memory would be allocated before processing reads. This can also avoid being interrupted by other users sharing the same machine. (so considerate honestly 😌)
  2. Gap filled bases now represented by lowercase characters in 'scafSeq' file. lowercase = filled, periodt. πŸ”‘
  3. Introduced SIMD instructions to boost the performance. speed. πŸŽοΈπŸ’¨

Configuration file πŸŽ€

For big genome projects with deep sequencing, the data is usually organized as multiple read sequence files generated from multiple libraries. The configuration file tells the assembler where to find these files and the relevant information. β€œexample.config” is an example of such a file. πŸ“„

The configuration file has a section for global information, and then multiple library sections. Right now only β€œmax_rd_len” is included in the global information section. Any read longer than max_rd_len will be cut to this length. βœ‚οΈ

The library information and the information of sequencing data generated from the library should be organized in the corresponding library section. Each library section starts with tag [LIB] and includes the following items:

  • avg_ins β€” This value indicates the average insert size of this library or the peak value position in the insert size distribution figure.
  • reverse_seq β€” This option takes value 0 or 1. It tells the assembler if the read sequences need to be complementarily reversed. Illumima GA produces two types of paired-end libraries: a) forward-reverse, generated from fragmented DNA ends with typical insert size less than 500 bp; b) forward-forward, generated from circularizing libraries with typical insert size greater than 2 Kb. The parameter β€œreverse_seq” should be set to indicate this: 0, forward-reverse; 1, forward-forward.
  • asm_flags=3 β€” This indicator decides in which part(s) the reads are used. It takes value 1 (only contig assembly), 2 (only scaffold assembly), 3 (both contig and scaffold assembly), or 4 (only gap closure).
  • rd_len_cutoff β€” The assembler will cut the reads from the current library to this length.
  • rank β€” It takes integer values and decides in which order the reads are used for scaffold assembly. Libraries with the same β€œrank” are used at the same time during scaffold assembly.
  • pair_num_cutoff β€” This parameter is the cutoff value of pair number for a reliable connection between two contigs or pre-scaffolds.
  • map_len β€” This takes effect in the β€œmap” step and is the minimum alignment length between a read and a contig required for a reliable read location.

The assembler accepts read file in two formats: FASTA or FASTQ. Mate-pair relationship could be indicated in two ways: two sequence files with reads in the same order belonging to a pair, or two adjacent reads in a single file (FASTA only) belonging to a pair. πŸ‘―

In the configuration file single end files are indicated by β€œf=/path/filename” or β€œq=/pah/filename” for fasta or fastq formats separately. Paired reads in two fasta sequence files are indicated by β€œf1=” and β€œf2=”. While paired reads in two fastq sequences files are indicated by β€œq1=” and β€œq2=”. Paired reads in a single fasta sequence file is indicated by β€œp=” item. 🧬

All the above items in each library section are optional. The assembler assigns default values for most of them. If you are not sure how to set a parameter, you can remove it from your configuration file. (no stress, bestie, it's got you πŸ’œ)


Get it started πŸš€

Once the configuration file is available, a typical way to run the assembler is:

${bin} all –s config_file –K 63 –R –o graph_prefix

User can also choose to run the assembly process step by step as:

${bin} pregraph –s config_file –K 63 [–R -d –p -a] –o graph_prefix
${bin} contig –g graph_prefix [–R –M 1 -D]
${bin} map –s config_file –g graph_prefix [-p]
${bin} scaff –g graph_prefix [–F -u -G -p]

✨ Options ✨

  • -a INT β€” Initiate the memory assumption (GB) to avoid further reallocation
  • -s STR β€” configuration file
  • -o STR β€” output graph file prefix
  • -g STR β€” input graph file prefix
  • -K INT β€” K-mer size [default 23, min 13, max 127]
  • -p INT β€” multithreads, n threads [default 8]
  • -R β€” use reads to solve tiny repeats [default no]
  • -d INT β€” remove low-frequency K-mers with frequency no larger than [default 0]
  • -D INT β€” remove edges with coverage no larger that [default 1]
  • -M INT β€” strength of merging similar sequences during contiging [default 1, min 0, max 3]
  • -F β€” intra-scaffold gap closure [default no]
  • -u β€” un-mask high coverage contigs before scaffolding [default mask]
  • -G INT β€” allowed length difference between estimated and filled gap
  • -L β€” minimum contigs length used for scaffolding

Output files πŸ“¦

These files are output as assembly results:

a. *.contig

contig sequences without using mate pair information

b. *.scafSeq

scaffold sequences (final contig sequences can be extracted by breaking down scaffold sequences at gap regions)

There are some other files that provide useful information for advanced users, which are listed in Appendix B. πŸ”


FAQ πŸ™‹β™€οΈ

How to set K-mer size?

The program accepts odd numbers between 13 and 31. Larger K-mers would have higher rate of uniqueness in the genome and would make the graph simpler, but it requires deep sequencing depth and longer read length to guarantee the overlap at any genomic location. (big kmer = big commitment, choose wisely πŸ’…)

How to set library rank?

SOAPdenovo will use the pair-end libraries with insert size from smaller to larger to construct scaffolds. Libraries with the same rank would be used at the same time. For example, in a dataset of a human genome, we set five ranks for five libraries with insert size 200-bp, 500-bp, 2-Kb, 5-Kb and 10-Kb, separately. It is desired that the pairs in each rank provide adequate physical coverage of the genome. πŸ“


APPENDIX A: an example.config πŸ“

#maximal read length
max_rd_len=50
[LIB]
#average insert size
avg_ins=200
#if sequence needs to be reversed
reverse_seq=0
#in which part(s) the reads are used
asm_flags=3
#use only first 50 bps of each read
rd_len_cutoff=50
#in which order the reads are used while scaffolding
rank=1
# cutoff of pair number for a reliable connection (default 3)
pair_num_cutoff=3
#minimum aligned length to contigs for a reliable read location (default
32)
map_len=32
#fastq file for read 1
q1=/path/**LIBNAMEA**/fastq_read_1.fq
#fastq file for read 2 always follows fastq file for read 1
q2=/path/**LIBNAMEA**/fastq_read_2.fq
#fasta file for read 1
f1=/path/**LIBNAMEA**/fasta_read_1.fa
#fastq file for read 2 always follows fastq file for read 1
f2=/path/**LIBNAMEA**/fasta_read_2.fa
#fastq file for single reads
q=/path/**LIBNAMEA**/fastq_read_single.fq
#fasta file for single reads
f=/path/**LIBNAMEA**/fasta_read_single.fa
#a single fasta file for paired reads
p=/path/**LIBNAMEA**/pairs_in_one_file.fa
[LIB]
avg_ins=2000
reverse_seq=1
asm_flags=2
rank=2
# cutoff of pair number for a reliable connection
#(default 5 for large insert size)
pair_num_cutoff=5
#minimum aligned length to contigs for a reliable read location
#(default 35 for large insert size)
map_len=35
q1=/path/**LIBNAMEB**/fastq_read_1.fq
q2=/path/**LIBNAMEB**/fastq_read_2.fq
q=/path/**LIBNAMEB**/fastq_read_single.fq
f=/path/**LIBNAMEB**/fasta_read_single.fa

πŸ¦‹ APPENDIX B: output files πŸ¦‹

ok bestie, this is the part where we open the lil gift bags 🎁 that SOAPdenovo leaves behind after every stage. Every command dumps its own set of files, and knowing what's in them is literally how you debug your assembly instead of just crying at it. Let's go file by file, no cap. πŸ’…


✨ 1. Output files from the command pregraph ✨

  • *.kmerFreq β€” every row is the count of Kmers whose frequency is exactly equal to the row number. So row 5 = "how many kmers showed up 5 times", periodt. πŸ”’
  • *.edge β€” one record per edge in the pre-graph, and it's giving: length, the Kmers sitting on both ends, average kmer coverage, whether the thing is reverse-complementarily identical, and the sequence itself. 🧬
  • *.markOnEdge & *.path β€” these two cuties work as a duo 🫢 and they're for using reads to solve small repeats.
  • *.preArc β€” the connections between edges that get established by the read paths. πŸ”—
  • *.vertex β€” the Kmers chillin' at the ends of edges. πŸ“
  • *.preGraphBasic β€” the vibe check of the whole pre-graph: number of vertex, K value, number of edges, maximum read length, etc. ✨
lowkey note: the letters skip from c to e up there β€” that's just how the original list rolls, we don't question it, we respect it πŸ™ˆ

✨ 2. Output files from the command contig ✨

  • *.contig β€” contig info time: corresponding edge index, length, kmer coverage, whether it's a tip, and the sequence. Heads up tho β€” either a contig or its reverse-complementary counterpart is included, not both. Every reverse complementary contig index gets flagged over in the *.ContigIndex file. 🧡
  • *.Arc β€” arcs coming out of each edge, plus their corresponding coverage by reads. 🏹
  • *.updated.edge β€” the glow-up info for each edge in the graph: length, Kmers at both ends, and the index difference between the reverse-complementary edge and this one. ✨
  • *.ContigIndex β€” each record spills the tea on a contig in *.contig: its edge index, length, and the index difference between its reverse-complementary counterpart and itself. πŸͺž

✨ 3. Output files from the command map ✨

  • *.peGrads β€” info for each clone library: insert-size, read index upper bound, rank, and pair number cutoff for a reliable link. And yes bestie, this file can be revised manually for scaffolding tuning β€” you're allowed to touch it. πŸŽ›οΈ
  • *.readOnContig β€” where reads land on contigs. Contigs get referred to by their edge index here. BUT β€” about half of them are not listed in the *.contig file, because their reverse-complementary counterparts are already in there. Not a bug, just the rules. πŸ’€
  • *.readInGap β€” reads that could be located in the gaps between contigs. This intel gets used to close gaps in scaffolds later. πŸ•³οΈβž‘οΈπŸ©Ή

✨ 4. Output files from the command scaff ✨

  • *.newContigIndex β€” contigs get sorted by their length before scaffolding, and their shiny new indexes live here. Super useful if you wanna match contigs in *.contig with the ones in *.links. πŸ”’
  • *.links β€” the links between contigs established by read pairs. New index are used here. πŸ”—
  • *.scaf_gap β€” contigs sitting in gaps found by the contig graph that the contiging procedure outputted. New index are used here too. 🧩
  • *.scaf β€” for each scaffold: the contigs it's made of β€” contig index (concordant to index in *.contig), approximate start position on the scaffold, orientation, contig length, and its links to others. πŸ‘‘
  • *.gapSeq β€” the gap sequences between contigs. 🫧
  • *.scafSeq β€” the actual sequence of each scaffold. The final main character moment. 🌟

πŸ‘‘ AUTHOR πŸ‘‘

Olivier Sallou (olivier.sallou (at) irisa.fr) β€” Man page and packaging. A legend, fr fr. 🫢


and that's a wrap on the whole page, bestie!! πŸ’– you just read version 1.1.0, dated July 30, 2012 β€” old but iconic, the blueprint that never misses. now go assemble something gorgeous, you've got this. ✨🧬🌸