Cobra

COBRA (Contig Overlap Based Re-Assembly) is a bioinformatics tool to get higher quality viral genomes assembled from metagenomes of short paired-end reads. COBRA was written in Python. COBRA has so far only been tested on assembled contigs and scaffolds from metaSPAdes, IDBA_UD, and MEGAHIT.

# Developed by Dr. LinXing Chen
# University of California, Berkeley
# The Banfield Lab
# Email: linkingchan@gmail.com

Versions

1. v1.2.2 (released on 2023-09-03) - initial release

2. v1.2.3 (released on 2024-02-26)

   • The GC function issue due to the update of Biopython.
   • The abnormal exit in the middle of processing some samples.
   • If none of the queries was extended, the process will break. If your runs do not have the expected output files, see the log file.

3. v1.3.0 (released on 2025-06-26)

   • refactor+feature: format code and add trim_readno
   • fix the handle of "6" shape path

Citation

The paper is out at Nature Microbiology (https://www.nature.com/articles/s41564-023-01598-2). Please cite as follows if you find COBRA is helpful for your analyses. Chen, L., Banfield, J.F. COBRA improves the completeness and contiguity of viral genomes assembled from metagenomes. Nat Microbiol (2024). https://doi.org/10.1038/s41564-023-01598-2

Introduction

1. Why are metagenomic contigs fragmented?

The genomes assembled from short paired-end reads based metagenomes are usually fragmented due to (1) intra-genome repeats, (2) inter-genome shared region, and (3) within-population variations, as the widely utilized assemblers based on de Bruijn graphs, e.g., metaSPAdes, IDBA_UD and MEGAHIT, tend to have a breaking point when multiple paths are available instead of making risky extension (see example in Figure 1).

Figure 1. An example of how assemblers break in assembly when within-population occurs.

2. Contigs may be joined with expected end overlap.

According to the principles of the abovementioned assemblers, the broken contigs have an end overlap with a determined length, that is, the max-kmer (maxK hereafter) used in de nono assembly for metaSPAdes and MEGAHIT, and the maxK-1 for IDBA_UD, which we termed as "expected overlap length" (Figures 1 and 2).

• Note: as COBRA will use the information provided by paired-end reads, only those samples sequenced by paired-end technology should work.

Figure 2. The "expected overlap length" has been documented in manual genome curation, see Chen et al. 2020 for details.

How COBRA works

COBRA determines the "expected overlap length" (both the forward direction and reverse complement direction) for all the contigs from an assembly, then looks for the valid joining path for each query that users provide (should be a fraction of the whole assembly) based on a list of features including contig coverage, contig overlap relationships, and contig continuity (based on paired-end reads mapping) (Figure 3).

Note that scaffolds (for example, metaSPAdes assembly) could be used as input for COBRA extension; however, we suggest not using scaffolds from IDBA_UD as the potential errors in the scaffolding step (see Chen et al. 2020 for details). Thus, for IDBA_UD and MEGAHIT assembly, the contigs should be used.

Given that COBRA has only tested for contigs/scaffolds from IDBA_UD, metaSPAdes, and MEGAHIT, it will be risky to use it on contigs/scaffolds from any other assemblers.

Figure 3. The workflow of COBRA.

Dependencies

• COBRA is a Python script (tested for version 3.7 or higher) that uses a list of frequently used Python packages, including:

Bio
Bio.Seq
collections
argparse
math
pysam
time

• The only third-party software that COBRA will use is BLASTn.

Installation

COBRA could now be installed in different ways.

• (1) git

git clone https://github.com/linxingchen/cobra.git

cd cobra

python cobra.py -h

• (2) pip

pip install cobra-meta

To confirm the installment,

cobra-meta -h

Which shows something like this

usage: cobra-meta [-h] -q QUERY -f FASTA -a {idba,megahit,metaspades} -mink MINK -maxk MAXK -m MAPPING -c COVERAGE [-lm LINKAGE_MISMATCH] [-o OUTPUT] [-t THREADS] [-v]

...

• (3) conda

conda create -n cobra python=3.8

conda activate cobra

conda install bioconda::cobra-meta or conda install linxingchen1987::cobra-meta

To confirm the installment,

cobra-meta -h

Which shows something like this

usage: cobra.py [-h] -q QUERY [-i IGNORE] -f FASTA -a {idba,megahit,metaspades} -mink MINK -maxk MAXK -m MAPPING
                [--mapping-link-cache [MAPPING_LINK_CACHE ...]] -c COVERAGE [-lm LINKAGE_MISMATCH] [-tr {no,trim,auto}]
                [--skip_new_assembly] [-o OUTPUT] [-t THREADS] [-v]

...

Update

• pip

pip install --upgrade cobra-meta

• conda

conda activate cobra (if cobra is the conda environment name)

conda update cobra-meta

Input files

(1) COBRA needs four files as inputs, i.e.,

• -f/--fasta: A fasta format file containing all the contigs from a single assembly. Note that IDBA_UD and MEGAHIT usually save contigs with a minimum length of 200 bp.

• -c/--coverage: a two-column (separated by tab) file of the sequencing coverage of all contigs in the -f/--fasta file, example below:

contig-140_0    25.552
contig-140_1    42.1388
contig-140_2    14.6023
contig-140_3    15.4817
contig-140_4    41.2746
...

• -q/--query: the query contigs that the user wants COBRA to extend, could be provided in a fasta format file, or a one-column text file with the names of the query contigs. Please make sure the names are exactly the same format as in the -f/--fasta file; otherwise, COBRA may have problems extending them.

• -m/--mapping: the paired-end reads mapping file of all contigs in the -f/--fasta file, could be sam or bam file.

(2) and three parameters

• -a/--assembler: the name of the de novo assembler used, currently only 'idba' (for IDBA_UD), 'metaspades' (for metaSPAdes), and 'megahit' (for MEGAHIT).
• -maxk/--maxk: the largest kmer used in de novo assembly.
• -mink/--mink: the smallest kmer used in de novo assembly.

(3) Optional flags

• -lm/--linkage_mismatch: the number of read mapping mismatches allowed when determining if paired reads spanned two contigs.
• -o/--output: the name of the output folder, otherwise it will be "{-q/--query}.COBRA" if not provided.
• -t/--threads: the number of threads used for BLASTn search.

How to obtain the mapping file

The mapping file could be obtained with tools like Bowtie2 and BBMap. Please refer to the manual descriptions for details of the tools. Below is the general way to get the sorted sam/bam file; you thus need to be available to samtools (which could be downloaded here - https://github.com/samtools/samtools).

For example,

• contig file = "contigs.fasta"

• first read file = "R1.fastq.gz"

• second read file = "R2.fastq.gz"

(1) with Bowtie2 (https://github.com/BenLangmead/bowtie2)

• bowtie2-build contigs.fasta contigs.fasta

• bowtie2 -p 16 -x contigs.fasta -1 R1.fastq.gz -2 R2.fastq.gz -S output.sam && samtools view -bS output.sam | samtools sort -o sorted_output.bam -

(2) with BBMap (https://github.com/BioInfoTools/BBMap)

• bbmap.sh ref=contigs.fasta in1=R1.fastq.gz in2=R2.fastq.gz threads=16 out=output.sam (good)
• samtools view -bS output.sam > output.bam
• samtools sort -o sorted_output.bam output.bam

How to obtain the coverage file

(1) with jgi_summarize_bam_contig_depths

Once the sorted sam or bam file is ready, the tool of "jgi_summarize_bam_contig_depths" from MetaBAT (https://bitbucket.org/berkeleylab/metabat/src/master/), or could be used to obtain the coverage file, the resulting profile should be transferred to get a two-column file divided by tab.

• jgi_summarize_bam_contig_depths --outputDepth original.coverage.txt *sam

• jgi_summarize_bam_contig_depths --outputDepth original.coverage.txt *bam

The output file from jgi_summarize_bam_contig_depths could be converted to a two-column file divided by tab using the script provided in this study (coverage.transfer.py).

• python coverage.transfer.py -i original.coverage.txt -o coverage.txt

(2) CoverM

CoverM is a fast DNA read coverage and relative abundance calculator focused on metagenomics applications. Usage could be found here (https://github.com/wwood/CoverM).

(3) pyCoverM

pyCoverM is a simple Python interface to CoverM's fast coverage estimation functions, which could be found here (https://github.com/apcamargo/pycoverm).

How to run

(1) The users can only specify the required parameters:

cobra-meta -f input.fasta -q query.fasta -c coverage.txt -m mapping.sam -a idba -mink 20 -maxk 140

(2) The users could also include the optional parameters like output name (-o), mismatch of mapped reads for linkage identification (-lm)

cobra-meta -f all.contigs.fasta -q query.fasta -o query.fasta.COBRA.out -c coverage.txt -m mapping.sam -a idba -mink 20 -maxk 140 -lm 2

cobra-meta -f all.contigs.fasta -q query.fasta -o query.fasta.COBRA.out -c coverage.txt -m mapping.sam -a metaspades -mink 21 -maxk 127 -lm 2

cobra-meta -f all.contigs.fasta -q query.fasta -o query.fasta.C