Assembling a plant genome from short reads

In development · manuscript in preparation

One paired-end library of 110 bp reads, a repeat-rich legume genome near a gigabase, and no long-read instrument. Below is the sequence of decisions we made for wild soybean accession KBL-107, in the order we made them.

Profile the reads before choosing an assembler

Reference-free profiling costs an afternoon and changes the parameters of the run that follows. Order of operations:

• FastQC per lane, then MultiQC to merge. Read the duplication rate and the GC distribution. Identical GC across forward and reverse lanes indicates no contaminating nucleic acid; a shoulder means a second organism to strip before assembly. Our run gave 35% GC on both lanes, duplication of 28.0% and 24.4%, and coverage above 72× against the ungapped length of a published reference.
• Count k-mers with KMC or Jellyfish, then fit the histogram in GenomeScope 2.0 for estimated genome length, repeat content, heterozygosity and error rate.
• Run smudgeplot on the k-mer dump for ploidy structure. Dumping the full KMC set needs more memory than a workstation holds, so this step runs on a Galaxy server.
• Run kmergenie in haploid mode for the graph-construction k. This answers a different question from the profiling k and there is no reason for the two to agree; ours came out at 25 for profiling and 69 for assembly.

Two results from that step changed the plan. The k-mer size recommended as a default in the tooling produced an implausible genome length and an implausible ploidy call, and the profile only became coherent at a larger k. Treat the documented default as a starting guess and vary k until GenomeScope and smudgeplot agree with each other.

Second, the ploidy call was mixed. Most of the genome behaved as a diploid, but a minority showed a third allele, which points to incomplete diploidisation. Soybean's ancient polyploidy is well documented. Seeing it in our own k-mer spectrum mattered because duplicated regions inflate the assembly graph without adding assembled length, so we read the contiguity metrics cautiously from then on.

Three assembler architectures, one library

We compared three assemblers with different underlying methods: ABySS 2.0 (single-k de Bruijn), MEGAHIT (multi-k de Bruijn) and MaSuRCA (super-reads plus overlap–layout–consensus). We left each at sensible defaults, since the comparison was about method, not tuning.

MaSuRCA gave the shortest total assembly of the three and the best contiguity. That combination is normal for fragmented short-read data: the assembler emits the same repetitive region many times as separate small contigs, so total length grows with repeat content. We ranked the assemblies on contig count and longest contig. Length can be recovered later by gap closing. Split coding sequence cannot.

Spending borrowed long reads

Neither the budget nor the instrument for our own long-read library existed. Public nanopore data for a second wild accession of the same species did. The placement of that data in the pipeline determines what the final assembly represents.

Route A, mix at assembly. Feed short and long reads to MaSuRCA together. Contiguity improves by roughly two orders of magnitude, as the last row of the table shows. Contigs are built from two accessions at once, including across coding sequence, so the product is a chimera of both.

Route B, spend after scaffolding. Assemble from the short reads alone, scaffold against a published reference with RagTag, then close the remaining gaps with TGSGapCloser using the borrowed long reads. Gaps are by construction the positions the short-read data failed to resolve, so coding sequence assembled from your own reads stays unmixed.

We ran both and compared them, which is the part worth reporting. Route B produces a weaker contiguity table and an assembly whose provenance you can state per region.

Testing whether the hybrid is chimeric

Contiguity metrics cannot detect chimerism, so use alignment. Align the two assemblies against each other with D-GENIES for a whole-genome dot plot and identity bands, then draw JupiterPlot circos figures for synteny and rearrangement. Read four things off the result.

• The share of total length with no counterpart in the other assembly. That fraction estimates how much sequence the borrowed library contributed outright.
• The share sitting in intermediate identity bands. A large middle band means the contigs are blended, which is the chimeric signature.
• Where low-identity regions concentrate. Ours clustered around the centromeres. Neither library has the long-range information those regions need, so this is a data problem. Hi-C or ultra-long reads would fix it. More short-read depth would not.
• Organellar sequence. Chloroplast and mitochondrial contigs agreed poorly between the two assemblies, so we treated organellar DNA as accession-specific during validation.

One chromosome showed an inverted segment between the two routes. Inversions at that scale are worth checking against a third assembly or against optical mapping before writing them up as biology, since assembly graphs resolve inverted repeats inconsistently.

Commands

Profiling, with the k-mer size varied until the two ploidy estimates agreed:

kmc -k25 -m64 -ci1 -cs10000 @reads.lst kmc_db tmp/
kmc_tools transform kmc_db histogram kmer_k25.hist -cx10000
genomescope2 -i kmer_k25.hist -k 25 -p 2 -o gscope_k25/
kmergenie reads.lst --diploid -l 21 -k 91 -s 4 -o kmergenie/

The three de novo assemblies. ABySS takes one k, MEGAHIT takes a list, and MaSuRCA is configured through a file:

abyss-pe name=kbl107 k=69 in='R1.fq.gz R2.fq.gz' j=32

megahit -1 R1.fq.gz -2 R2.fq.gz --k-list 33,55,77,99,121 \
        -t 32 -o megahit_out/

# masurca config: PE= pe 350 50 R1.fq.gz R2.fq.gz
#                 NANOPORE= ont.fq.gz        (route A only)
#                 USE_LINKING_MATES=1
masurca config.txt && ./assemble.sh

Scaffolding and gap closing for route B, where the borrowed reads enter after assembly:

ragtag.py scaffold reference.fna masurca_contigs.fasta -o ragtag_out/ -t 32
tgsgapcloser --scaff ragtag_out/ragtag.scaffold.fasta \
             --reads ont.fq.gz --output kbl107_gapclosed \
             --ne --thread 32

Comparison of the two routes:

minimap2 -x asm5 -t 32 routeA.fasta routeB.fasta > A_vs_B.paf
# D-GENIES for the dot plot and identity bands; JupiterPlot for synteny
jupiter name=A_vs_B ref=routeA.fasta fa=routeB.fasta ng=90

Stack

Read QC

FastQC 0.11.8, MultiQC 1.19

Profiling

KMC 3, Jellyfish 2, GenomeScope 2.0, smudgeplot, kmergenie

Assembly

ABySS 2.0, MEGAHIT, MaSuRCA 3.2.8

Scaffolding

RagTag, TGSGapCloser

Annotation

Liftoff

Comparison

D-GENIES, JupiterPlot, Fasta Statistics, sequence-stats

Compute

Intel Xeon Platinum 8368, 512 GB RAM; smudgeplot run on a Galaxy server for the k-mer dump

Current state

Structural and functional annotation runs by lift-over from an annotated reference with Liftoff. Lift-over inherits the reference's gene models, so it cannot support claims about novel or diverged loci. The draft says this explicitly and does not treat the transferred annotation as an independent result.

The comparison between the two routes is the part worth publishing. A non-model plant, a small budget and no long-read instrument is the ordinary situation for most groups. The same choice comes up whenever public long reads exist for a relative of your sample.

Presented in part as Inside the northernmost soybean, February 2025.

References

1. Jackman S. D. et al. ABySS 2.0: resource-efficient assembly of large genomes using a Bloom filter. doi:10.1101/gr.214346.116
2. Li D. et al. MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly. doi:10.1093/bioinformatics/btv033
3. Zimin A. V. et al. The MaSuRCA genome assembler. doi:10.1093/bioinformatics/btt476
4. Ranallo-Benavidez T. R. et al. GenomeScope 2.0 and Smudgeplot for reference-free profiling of polyploid genomes. doi:10.1038/s41467-020-14998-3
5. Chikhi R., Medvedev P. Informed and automated k-mer size selection for genome assembly. doi:10.1093/bioinformatics/btt310
6. Kokot M. et al. KMC 3: counting and manipulating k-mer statistics. doi:10.1093/bioinformatics/btx304
7. Alonge M. et al. Automated assembly scaffolding with RagTag. doi:10.1186/s13059-022-02823-7
8. Xu M. et al. TGS-GapCloser: fast and accurate gap closing with error-prone long reads. doi:10.1093/gigascience/giaa094
9. Shumate A., Salzberg S. L. Liftoff: accurate mapping of gene annotations. doi:10.1093/bioinformatics/btaa1016
10. Cabanettes F., Klopp C. D-GENIES: dot plot large genomes in an interactive, efficient and simple way. doi:10.7717/peerj.4958