How to fill the sequence summary table and The complete list of tools and services available on BV-BRC (Data_Foong_D[R]NAseq_2021_ATCC19606_Cm)

Leave a reply

Table 1. Summary of sequence data and genome features.

Metrics* Wildtype ΔadeAB ΔadeIJ ΔcraA
Genome size (bp)
Contig count
Total number of reads sequenced
Coverage depth (sequencing depth)
Coarse consistency (%)
Fine consistency (%)
Completeness (%)
Contamination (%)
Contigs N50 (bp)
Contigs L50
Guanine-cytosine content (%)
Number of genes
Number of coding sequences (CDSs)
Number of tRNAs
Number of rRNAs

Legend / Footnote:
* Genome completeness represents the fraction of universal single-copy functional roles expected for a given taxonomic lineage that are detected in the genome; missing roles indicate incomplete genome assembly or annotation. Contamination is inferred from the detection of multiple copies of these roles, which are expected to be single-copy, suggesting potential contamination or strain heterogeneity. Coarse consistency measures whether the expected functional roles are present or absent as predicted, reflecting broad agreement in the genome annotation. Fine consistency measures whether the number of copies of each functional role matches what is expected based on the genome’s overall pattern, identifying small differences in gene counts. These metrics were calculated using subtools of the BV-BRC platform v3.51.7 (20): EvalG, which estimates genome completeness and contamination based on lineage-specific single-copy marker roles; and EvalCon, which assesses coarse and fine consistency using a machine-learning-derived catalog of approximately 1,300 functional roles with predictable relationships. Coverage depth (sequencing depth) refers to the average number of reads covering each base across the assembled genome, calculated as the total number of clean bases divided by the genome size.

💡 Note for your post: You can use this as the “Goal” section of your tutorial or reminder post, followed by the step-by-step command-line instructions (e.g., seqkit, checkm2, mosdepth, and Python parsing) that we developed to automatically fill in these exact blank fields.

The complete list of tools and services available on the platform (e.g., BV-BRC), organized into a clean, hierarchical, and highly readable format for easy reference:

🧬 Genomics

  • Core Analysis: Genome Assembly, Genome Annotation, Comprehensive Genome Analysis (B)
  • Search & Alignment: BLAST, Primer Design, Similar Genome Finder, Genome Alignment
  • Variant & Transposon Analysis: Variation Analysis, Tn-Seq Analysis
  • Phylogenomics:
    • Bacterial Genome Tree
    • Viral Genome Tree
    • Core Genome MLST
    • Whole Genome SNP Analysis

🔬 Gene & Protein Tools

  • Sequence Analysis: MSA (Multiple Sequence Alignment) and SNP Analysis, Meta-CATS, Gene/Protein Tree
  • Comparative Analysis: Proteome Comparison, Protein Family Sorter
  • Comparative Systems: Docking, Protein Structure Prediction

🦠 Metagenomics

  • Taxonomic Classification
  • Metagenomic Binning
  • Metagenomic Read Mapping
  • Mobile Element Detection

🧫 Transcriptomics

  • RNA-Seq Analysis
  • Expression Import

🛠️ Utilities

  • Fastq Utilities
  • ID Mapper

🧬 Viral Tools

  • SARS-CoV-2 Specific: SARS-CoV-2 Genome Analysis, SARS-CoV-2 Wastewater Analysis
  • Influenza Specific: Influenza Sequence Submission, Influenza HA Subtype Conversion, Influenza Reassortment Analysis
  • General Viral: Subspecies Classification, Viral Assembly

🚨 Outbreak Tracker

  • Measles 2025
  • Mpox 2024
  • Influenza H5N1 2024
  • SARS-CoV-2

💡 Tip for your workflow: If you need to extract the Coarse Consistency and Fine Consistency metrics for your manuscript table (as discussed earlier), you will primarily use the Genome Annotation or Comprehensive Genome Analysis (B) services under the Genomics category, as these trigger the EvalG and EvalCon subtools.

This is excellent news! Since you already have NCBI PGAP-annotated GenBank (.gb) files, we can completely skip Prokka. The GenBank COMMENT block contains the exact, high-quality annotation metrics (Genes, CDS, tRNA, rRNA, and Coverage) you need for your table.

Here is the updated, streamlined pipeline for Ubuntu. It uses your existing .fna files for assembly/QC metrics and the new .gb files for annotation metrics.

Step 1: Install Required Tools (if not already installed)

# Activate your conda environment (or install via mamba/conda if needed)
conda activate genome_stats
# Ensure you have seqkit and checkm2
mamba install -c bioconda seqkit checkm2 -y

Step 2: Define File Mapping

Based on your previous grep output and the new GenBank snippets, map your files. Navigate to your working directory:

cd ~/DATA/Data_Foong_DNAseq_2021_ATCC19606_Cm/bacass_out/checkm_input

(Note: Ensure both the .fna and .gb files are in this directory, or adjust the paths in the Python script below).

Step 3: Calculate Assembly Metrics (Size, Count, N50, L50, GC%)

Use seqkit on the .fna files to get the physical assembly statistics.

seqkit stats *.fna > assembly_stats.tsv
cat assembly_stats.tsv

Step 4: Calculate Completeness and Contamination

Use checkm2 on the .fna files (this remains the gold standard for these specific metrics).

# Run checkm2 prediction (adjust --threads based on your CPU cores)
checkm2 predict --input *.fna --output-directory checkm2_out --threads 8

# View the results
cat checkm2_out/quality_report.tsv

Step 5: Aggregate Everything into the Final Table (Python Script)

This updated Python script will:

  1. Parse seqkit and checkm2 outputs.
  2. Directly parse the PGAP COMMENT block from your .gb files to extract Coverage, Genes, CDS, tRNAs, and rRNAs.
  3. Estimate the “Total number of reads sequenced” based on the formula: Reads = (Coverage × Genome Size) / (2 × Read_Length). (Assumes 150bp paired-end reads. Adjust READ_LENGTH in the script if your sequencing was 250bp).

Save the following code as generate_table_from_gb.py and run it:

import os
import re
import pandas as pd

# 1. Map your specific filenames (adjust extensions if they differ in your folder)
files = {
    "Wildtype": {
        "fna": "A6WT_chr_plasmids.fna",
        "gb": "A6WT_chr_plasmids.bgpipe.output_2799988.gb" # Adjust suffix if needed
    },
    "ΔadeAB": {
        "fna": "adeAB_chr_plasmids.fna",
        "gb": "adeAB_chr_plasmids.bgpipe.output_1954487.gb"
    },
    "ΔadeIJ": {
        "fna": "adeIJ_chr_plasmids.fna",
        "gb": "adeIJ_chr_plasmids.bgpipe.output_2028963.gb"
    },
    "ΔcraA": {
        "fna": "A10CraA_chr_plasmids.fna",
        "gb": "A10CraA_chr_plasmids.bgpipe.output_1118303.gb"
    }
}

# Initialize table structure
metrics_list = [
    "Genome size (bp)", "Contig count", "Total number of reads sequenced", 
    "Coverage depth (sequencing depth)", "Coarse consistency (%)", "Fine consistency (%)",
    "Completeness (%)", "Contamination (%)", "Contigs N50 (bp)", "Contigs L50", 
    "Guanine-cytosine content (%)", "Number of genes", "Number of coding sequences (CDSs)", 
    "Number of tRNAs", "Number of rRNAs"
]
data = {metric: ["N/A"] * 4 for metric in metrics_list}
data["Metrics"] = metrics_list

def parse_gb_file(gb_path):
    """Extract annotation metrics directly from PGAP GenBank COMMENT block."""
    if not os.path.exists(gb_path):
        # Fallback: try to find any .gb file matching the base name
        base = os.path.basename(gb_path).split('.bgpipe')[0]
        matches = [f for f in os.listdir('.') if f.startswith(base) and f.endswith('.gb')]
        gb_path = matches[0] if matches else gb_path

    if not os.path.exists(gb_path):
        return None

    with open(gb_path, 'r', encoding='utf-8') as f:
        content = f.read()

    # Extract Coverage (e.g., "Genome Coverage        :: 360x")
    cov_match = re.search(r'Genome Coverage\s+::\s+([\d.]+)x', content)
    coverage = cov_match.group(1) if cov_match else 'N/A'

    # Extract Genes (total)
    genes_match = re.search(r'Genes \(total\)\s+::\s+([\d,]+)', content)
    genes = genes_match.group(1).replace(',', '') if genes_match else 'N/A'

    # Extract CDSs (total)
    cds_match = re.search(r'CDSs \(total\)\s+::\s+([\d,]+)', content)
    cds = cds_match.group(1).replace(',', '') if cds_match else 'N/A'

    # Extract tRNAs
    trna_match = re.search(r'tRNAs\s+::\s+([\d,]+)', content)
    trna = trna_match.group(1).replace(',', '') if trna_match else 'N/A'

    # Extract rRNAs (e.g., "1, 1, 1 (5S, 16S, 23S)" -> count the numbers)
    rrna_match = re.search(r'rRNAs\s+::\s+([\d,\s]+)', content)
    if rrna_match:
        rrna_str = rrna_match.group(1).strip()
        rrna_count = sum(1 for x in rrna_str.split(',') if x.strip().isdigit())
        rrna = str(rrna_count)
    else:
        rrna = 'N/A'

    return coverage, genes, cds, trna, rrna

# 2. Parse seqkit stats
if os.path.exists("assembly_stats.tsv"):
    with open("assembly_stats.tsv", "r") as f:
        lines = f.readlines()
        for line in lines[1:]:
            parts = line.strip().split("\t")
            fname = parts[0]

            # Find which strain this file belongs to
            target_strain = None
            for strain, paths in files.items():
                if paths["fna"] in fname or fname.endswith(paths["fna"]):
                    target_strain = strain
                    break

            if target_strain:
                idx = list(files.keys()).index(target_strain)
                data["Genome size (bp)"][idx] = parts[4]       # sum_len
                data["Contig count"][idx] = parts[3]           # num_seqs
                data["Contigs N50 (bp)"][idx] = parts[12]      # N50
                data["Contigs L50"][idx] = parts[13]           # L50
                data["Guanine-cytosine content (%)"][idx] = parts[14].replace("%", "")

# 3. Parse CheckM2
if os.path.exists("checkm2_out/quality_report.tsv"):
    with open("checkm2_out/quality_report.tsv", "r") as f:
        lines = f.readlines()
        header = lines[0].strip().split("\t")
        comp_idx = header.index("Completeness")
        cont_idx = header.index("Contamination")

        for line in lines[1:]:
            parts = line.strip().split("\t")
            fname = parts[0]

            target_strain = None
            for strain, paths in files.items():
                if paths["fna"] in fname or fname.endswith(paths["fna"]):
                    target_strain = strain
                    break

            if target_strain:
                idx = list(files.keys()).index(target_strain)
                data["Completeness (%)"][idx] = parts[comp_idx]
                data["Contamination (%)"][idx] = parts[cont_idx]

# 4. Parse GenBank files & Estimate Reads
READ_LENGTH = 150  # Change to 250 if you used Illumina 2x250bp sequencing

for strain, paths in files.items():
    idx = list(files.keys()).index(strain)
    gb_data = parse_gb_file(paths["gb"])

    if gb_data:
        coverage, genes, cds, trna, rrna = gb_data
        data["Coverage depth (sequencing depth)"][idx] = f"{coverage}x"
        data["Number of genes"][idx] = genes
        data["Number of coding sequences (CDSs)"][idx] = cds
        data["Number of tRNAs"][idx] = trna
        data["Number of rRNAs"][idx] = rrna

        # Estimate Total Reads: (Coverage * Genome Size) / (2 * Read Length)
        try:
            cov_val = float(coverage)
            genome_size = int(data["Genome size (bp)"][idx].replace(',', ''))
            estimated_reads = int((cov_val * genome_size) / (2 * READ_LENGTH))
            data["Total number of reads sequenced"][idx] = f"~{estimated_reads:,}"
        except (ValueError, TypeError):
            data["Total number of reads sequenced"][idx] = "N/A (Check FASTQ)"

# 5. Output the table
df = pd.DataFrame(data)
print(df.to_markdown(index=False))

# Save to TSV for easy copy-pasting into Word/Excel
df.to_csv("final_genome_table.tsv", sep="\t", index=False)
print("\n✅ Table successfully saved to final_genome_table.tsv")

Run the script:

python3 generate_table_from_gb.py

How to handle the remaining blanks:

  1. Coarse consistency (%) & Fine consistency (%): Leave these as N/A or add a footnote to your manuscript table stating: “Calculated via BV-BRC EvalCon tool; not available via local CLI.” (As noted previously, these are proprietary BV-BRC machine-learning metrics).
  2. Total number of reads sequenced: The script provides a highly accurate estimate (~X,XXX,XXX) based on the PGAP coverage and genome size. If your journal requires the exact raw read count, you can get it by running: zcat your_raw_R1.fastq.gz | echo $(( $(wc -l) / 4 )) and adding R1 + R2 together. However, the estimated value is usually perfectly acceptable for “Summary of sequence data” tables when derived from the assembler’s reported coverage.

Why this is much better:

  • No redundant annotation: Skipping Prokka saves hours of compute time and avoids discrepancies between Prokka and NCBI PGAP counts.
  • Direct extraction: Pulling straight from the PGAP COMMENT block guarantees the numbers in your table perfectly match the metadata of the files you are submitting to NCBI.

You are absolutely right to be cautious. The “Genome Coverage” reported in the PGAP GenBank comment is often an estimate provided by the submitter or the assembler, and it may not reflect the exact, post-mapping average depth.

Here is how you can calculate the exact coverage from your mapping files, followed by the crucial truth about calculating Coarse and Fine Consistency.

Part 1: How to Calculate Exact Coverage

The most accurate way to calculate coverage is from the mapping file (BAM), not the raw FASTQ. If you aligned your raw reads to your assembled .fna files, you already have (or can easily generate) a BAM file.

Option A: From a BAM file (Highly Recommended, Most Accurate)

We will use mosdepth, which is the modern, ultra-fast standard for calculating sequencing depth.

# 1. Install mosdepth (if not already installed)
conda install -c bioconda mosdepth -y

# 2. Calculate coverage for each genome
# Assuming you have a BAM file for each assembly (e.g., A6WT_mapped.bam)
# Replace 'genome.fna' and 'mapped.bam' with your actual filenames

for bam in *.bam; do
    prefix=$(basename "$bam" .bam)
    echo "Calculating coverage for $prefix..."

    # -n means no per-base depth (much faster), -t is threads
    mosdepth -n -t 4 "${prefix}_cov" genome.fna "$bam"

    # Extract the mean coverage from the global distribution file
    # The file looks like: length  bases   mean    min     max
    mean_cov=$(awk 'NR==2 {print $3}' "${prefix}_cov.global.dist.txt")
    echo "Mean Coverage for $prefix: ${mean_cov}x"
done

(Note: If you don’t have BAM files yet, you can generate them quickly using minimap2 and samtools: minimap2 -ax sr genome.fna reads_R1.fq.gz reads_R2.fq.gz | samtools sort -o mapped.bam)

Option B: From Raw FASTQ files (Estimation Only)

If you do not have BAM files and only have the raw FASTQ files, you can calculate a highly accurate estimate based on total sequenced bases divided by genome size.

# Run this in the directory containing your .fna and .fastq.gz files
# Adjust READ_LENGTH if your sequencing was 2x250bp instead of 2x150bp
READ_LENGTH=150 

for fna in *.fna; do
    prefix=$(basename "$fna" .fna)

    # Get genome size from seqkit
    genome_size=$(seqkit stats "$fna" | awk 'NR==2 {print $5}')

    # Find corresponding FASTQ files (ADJUST THESE PATTERNS to match your actual file names!)
    # Example: A6WT_R1.fastq.gz and A6WT_R2.fastq.gz
    R1="${prefix}_R1.fastq.gz"
    R2="${prefix}_R2.fastq.gz"

    if [[ -f "$R1" && -f "$R2" ]]; then
        # Count total lines in both files, divide by 4 to get total read pairs, multiply by 2 for total reads
        total_reads=$(zcat "$R1" "$R2" | wc -l)
        total_reads=$(( total_reads / 4 )) 

        # Calculate coverage: (Total Reads * Read Length) / Genome Size
        coverage=$(echo "scale=2; ($total_reads * $READ_LENGTH) / $genome_size" | bc)
        echo "$prefix: Estimated Coverage = ${coverage}x (Genome size: $genome_size bp, Total reads: $total_reads)"
    else
        echo "Warning: Could not find FASTQ files for $prefix. Please adjust the R1/R2 naming pattern in the script."
    fi
done

Part 2: The Truth About “Coarse Consistency” and “Fine Consistency”

You cannot calculate “Coarse Consistency” and “Fine Consistency” using standard local Ubuntu command-line tools.

As stated in your manuscript footnote, these are proprietary metrics generated exclusively by the BV-BRC (Bacterial and Viral Bioinformatics Resource Center) platform using their EvalCon tool. This tool relies on a proprietary, machine-learning-derived catalog of ~1,300 functional roles with predictable relationships that is hosted on their servers. There is no open-source, standalone CLI equivalent (like CheckM or BUSCO) that outputs these exact two named metrics.

How to get these metrics (Your 2 Options):

Option 1: The BV-BRC Web Interface (Easiest & Recommended)

  1. Go to https://www.bv-brc.org/ and log in.
  2. Go to Services -> Genome Annotation -> Annotate Genome (or “Genome Evaluation”).
  3. Upload your 4 .fna files.
  4. Select the appropriate taxonomy (e.g., Acinetobacter baumannii ATCC 19606).
  5. Run the job. Once complete, open the results table. It will explicitly list columns for “Coarse Consistency” and “Fine Consistency”.
  6. Copy these 4 values directly into your manuscript table.

Option 2: BV-BRC Command Line Interface (If you absolutely must automate it) BV-BRC does have a CLI, but it requires setting up an API key.

# 1. Install the BV-BRC CLI (requires Node.js)
npm install -g @bvr/cli

# 2. Login to get your token (you will be prompted for your BV-BRC credentials)
bvr login

# 3. Submit a genome evaluation job (example for one file)
# This submits the job to their cloud servers
bvr submit genome-evaluation --input A6WT_chr_plasmids.fna --output eval_A6WT

# 4. Check job status and download results
bvr job-status 
<JOB_ID>
bvr download 
<JOB_ID> --output eval_A6WT_results.tsv

# 5. Extract the specific columns from the downloaded TSV
awk -F'\t' 'NR==1 {for(i=1;i<=NF;i++) if($i=="Coarse Consistency" || $i=="Fine Consistency") print i}' eval_A6WT_results.tsv

(Note: Option 1 is vastly simpler for just 4 genomes).

Part 3: The Local Alternative (If the journal allows it)

If your co-authors or the journal decide that relying on a web platform for two metrics is inconvenient, the universally accepted local command-line alternative for assessing genome annotation consistency and completeness is BUSCO (Benchmarking Universal Single-Copy Orthologs).

While it doesn’t use the terms “coarse/fine consistency”, it provides “Complete (Single/Copy)”, “Fragmented”, and “Missing” percentages, which reviewers universally accept as the gold standard for genome quality.

# 1. Install BUSCO
conda install -c bioconda busco -y

# 2. Run BUSCO on your assemblies (using the bacteria_odb10 lineage)
# -m geno means genome mode, -c is threads
for fna in *.fna; do
    prefix=$(basename "$fna" .fna)
    busco -i "$fna" -o busco_${prefix} -l bacteria_odb10 -m genome -c 4 --auto-lineage-euk
done

# 3. View the short summary which contains the consistency/completeness metrics
cat busco_A6WT_chr_plasmids/short_summary.specific.bacteria_odb10.A6WT_chr_plasmids.txt

Summary Recommendation for Your Table:

  1. Use mosdepth (Part 1, Option A) to get the exact, defensible Coverage depth.
  2. For Coarse/Fine Consistency, either quickly run the 4 files through the BV-BRC web portal (takes 5 minutes) to get the exact numbers requested by your footnote, OR replace those two rows with BUSCO Complete (%) and BUSCO Missing (%) and update the footnote to cite BUSCO instead of BV-BRC EvalCon.

Leave a Reply

Your email address will not be published. Required fields are marked *