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
v1.2.2 (released on 2023-09-03) - initial release
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.
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.
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:
-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).
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.
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).
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:
For all the queries, COBRA assigns them to different categories based on their joining status (detailed in the COBRA_joining_status.txt file), i.e.,
“self_circular” - the query contig itself is a circular genome.
“extended_circular” - the query contig was joined and extended into a circular genome.
“extended_partial” - the query contig was joined and extended, but not circular.
“extended_failed” - the query contig was not able to be extended due to COBRA rules.
“orphan_end” - neither end of a given contig shares “expected overlap length” with others.
“complex_end” - the query contig has multiple possible joining paths, so COBRA could not resolve them.
For the joined and extended queries in each category, only the unique ones (*.fa) will be saved for users’ following analyses, and the sequence information (e.g., length, coverage, GC, num of Ns) is summarized in the *fa.summary.tsv files. For categories of “extended_circular”, and “extended_partial”, the joining details of each query are included in the corresponding folder and summarized in the COBRA_joining_summary.tsv file, an example shown below:
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.
Versions
v1.2.2 (released on 2023-09-03) - initial release
v1.2.3 (released on 2024-02-26)
v1.3.0 (released on 2025-06-26)
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-2Introduction
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).
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:
The only third-party software that COBRA will use is BLASTn.
Installation
COBRA could now be installed in different ways.
git clone https://github.com/linxingchen/cobra.gitcd cobrapython cobra.py -hpip install cobra-metaTo confirm the installment,
cobra-meta -hWhich shows something like this
conda create -n cobra python=3.8conda activate cobraconda install bioconda::cobra-metaorconda install linxingchen1987::cobra-metaTo confirm the installment,
cobra-meta -hWhich shows something like this
Update
pip install --upgrade cobra-metaconda activate cobra(if cobra is the conda environment name)conda update cobra-metaInput 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/--fastafile, example below:-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/--fastafile; otherwise, COBRA may have problems extending them.-m/--mapping: the paired-end reads mapping file of all contigs in the-f/--fastafile, 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.fastabowtie2 -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.bamsamtools sort -o sorted_output.bam output.bamHow 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 *samjgi_summarize_bam_contig_depths --outputDepth original.coverage.txt *bamThe 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:
(2) The users could also include the optional parameters like output name (-o), mismatch of mapped reads for linkage identification (-lm)
Output files
Below is a general list of output files in the output folder:
For all the queries, COBRA assigns them to different categories based on their joining status (detailed in the
COBRA_joining_status.txtfile), i.e.,For the joined and extended queries in each category, only the unique ones (
*.fa) will be saved for users’ following analyses, and the sequence information (e.g., length, coverage, GC, num of Ns) is summarized in the*fa.summary.tsvfiles. For categories of “extended_circular”, and “extended_partial”, the joining details of each query are included in the corresponding folder and summarized in theCOBRA_joining_summary.tsvfile, an example shown below:logfile includes the content of each processing step, an example shown below: