Previous module:
Next module:

BamToCov is a coverage analysis toolkit developed at the Quadram Institute Bioscience.

logo

BamToCov main tool is inspired by the UNIX Phylosophy where the tools are designed for efficient computation of a very specific task. Integration of multiple samples and specific tasks can be achieved with scripts.

  • bamtocov will produce a coverage BED from a single BAM file, or a count matrix from a set of alignments and a target (in BED, GTF or GFF format). Used without a target, it is a drop-in replacement for covtobed, but discarding invalid alignments by default. When providing the target, it can produce coverage statistics for each region in the target, also with multiple BAM files.

There are two further utilities designed to simplify specific tasks:

  • bamtocounts will count the number of reads covering each target region, rather than the nucleotidic coverage. This can count reads in exons, in RNA-Seq experiments, or insertion events in Tn-Seq or TraDIS Experiments

  • bamcountrefs is a shortcut to count the number of reads per chromosome, with filters on the read flags, length and quality. Yes, this is the one for us today

:question: Why Bamtocov? It supports a new metrics (DisCov) and in our test was 60% faster than CoverM on this dataset

Workflow

BamCountRefs is part of BamToCov.

It builds count tables (and other per-contig metrics) from multiple BAM files aligned against the same reference, which is exactly our situation: several samples mapped against a single co-assembly.

Installing BamCountRefs

We installed BamToCov inside the “mapping” conda environment, otherwise it’s easily installable from BioConda.

# Activate the environment
conda activate mapping

# Check installation
bamcountrefs --help

Running BamCountRefs

We have three BAM files in bams/, each corresponding to the mapping of one sample against the co-assembly. BamCountRefs requires sorted and indexed BAM files: if the .bai files are missing, create them first:

# Index the BAM files (only if not already done)
for BAM in bams/*.bam; do
  samtools index "$BAM"
done

Then we can calculate the metrics for all contigs:

# Calculate coverage metrics for all contigs
bamcountrefs \
  --output abundance \
  --mean --covered-ratio --tpm \
  --workers 3 \
  bams/*.bam

Some explanations about the parameters used:

  • --output: specifies the basename of the output files. Unlike other tools, BamCountRefs writes one table per metric (see below). Without this option, a single table is printed to the screen
  • --workers: number of BAM files to process in parallel (here 3, one per sample)
  • bams/*.bam: the input BAM files, given as positional arguments (wildcards are fine)

We can specify one or more metrices, in our example:

  • --mean: mean coverage depth across the contig (approximate method, no extra memory)
  • --covered-ratio: proportion of the contig covered by reads (0-1)
  • --tpm: Transcripts Per Million (normalised abundance)

The raw read counts are always calculated, so there is no need to request them.

Supplementary alignments do not need to be excluded explicitly: by default BamCountRefs ignores unmapped, secondary, supplementary, duplicate and QC-failed alignments (-F 3844).

Other useful filters are -Q to set a minimum mapping quality and -P to only count properly paired reads.

:bulb: You can use --all-metrics to calculate every available metric in one go, but note that some of them (covered bases and ratio, variance, trimmed mean, DisCov) require memory proportional to the total length of the contigs.

Output Files

Our command produces four TSV files, one per metric:

File Description
abundance_counts.tsv Number of reads mapped to each contig
abundance_mean.tsv Mean coverage depth across the contig
abundance_covered_fraction.tsv Proportion of contig covered by reads (0-1)
abundance_tpm.tsv Transcripts Per Million (normalised abundance)

Each table has the same structure: the first column (Contig) has the contig identifier, followed by one column per sample, named after the BAM file. In this case, with 3 samples, each table has 4 columns (1 for contig names + 3 samples).

# Inspect the first lines of the count table
head abundance_counts.tsv | column -t

Having one metric per file means each table is ready to be loaded as a matrix in R or Python without having to select columns first.

Note: Raw counts are the most useful value to perform statistical analyses, while TPM values are useful for a first “eyeballing” comparison of abundances between samples, as they account for both sequencing depth and contig length.

Previous submodule:
Next submodule: