Coverage Plots

Plot Nanopore read coverage from BAM files for genes and transcripts.

Compact Layout Monkey-Patch

This monkey-patch improves trackplot layouts by: - Making the x-axis ruler thinner - Reducing vertical spacing between tracks - Pulling annotation tracks closer to data

Helper Functions

Main Coverage Plot Function

API reference

---

plot_gene_coverage

def plot_gene_coverage(
    adata, transcript_data, gene:str | None=None, chrom:str | None=None, start:int | None=None, end:int | None=None,
    strand:str | None=None, groupby:str | Sequence[str]=('cell_type',), bc_col:str='barcode_raw',
    bam_paths:Mapping[str, str] | None=None, bam_column:str='batch', bc_tag:str='BC', padding:int=1000,
    figsize:tuple[float, float]=(6, 3), font_size:int=6, n_y_ticks:int=3, threshold:int=1, same_y:bool=False,
    raster:bool=True, show_site_plot:bool=False, transcripts:List[str] | None=None, top_n:int | None=None,
    choose_primary:bool=False, annotation_scale:float=0.25, intron_scale:float=0.15, exon_width:float=0.3,
    show_cell_labels:bool=True, title:str | None=None, # kept for API compatibility (unused)
    gtf_file:str | None=None, color_key:Optional[str]='cell_type', palette:str | list='tab10',
    group_color_map:Optional[Dict[str, str]]=None,
    group_order:Optional[List[str]]=None, # deprecated: use visible_groups
    visible_groups:Optional[List[str]]=None, show_title:bool=False, # kept for API compatibility (ignored)
    distance_between_label_axis:float=0.2, xaxis_height_scale:float=0.2, annotation_height_scale:float=0.6,
    grid_hspace:float=0.18, verbose:bool=False
): # The figure object (so the caller can save it via fig.savefig()).

Plot Nanopore read coverage from BAM files with transcript annotations.

Example Usage

Load test data and plot coverage for the Myl6 gene.

# Setup barcodes
import pathlib, re, collections

BAMDIR = pathlib.Path("/data/analysis/data_mcandrew/00_allos_0.6_TAB/data/bams")

# Strip any "-1/-2 ..." suffix, keep clean ID
adata.obs["barcode_prefixed"] = adata.obs_names.str.replace(r"-\d+$", "", regex=True)

# Extract 16-mer barcode
if "barcode_raw" not in adata.obs:
    adata.obs["barcode_raw"] = adata.obs["barcode_prefixed"].str.extract(r"([ACGT]{16})$")[0]

# Define BAM file paths
bam_paths = {
    "0": "/data/analysis/data_mcandrew/00_allos_0.6_TAB/data/bams/s190_isobam.bam",
    "1": "/data/analysis/data_mcandrew/00_allos_0.6_TAB/data/bams/s951_isobam.bam",
}

# Plot coverage for Myl6 gene
fig = plot_gene_coverage(
    adata=adata,
    transcript_data=td,
    gene="Myl6",
    groupby=["cell_type"],
    bam_paths=bam_paths,
    bam_column="batch",
    gtf_file=str(gtf_local),
    figsize=(6, 4),
    font_size=7,
    threshold=10,
    top_n = 2,
    group_order = celltype_order,
    group_color_map = celltype_palette,
    show_title = True
)
plt.show()

# Plot coverage for Myl6 gene with top_n parameter
# This automatically selects the 2 most abundant transcripts
p = plot_gene_coverage(
    adata=adata,
    transcript_data=td,
    gene="Myl6",
    groupby=["cell_type"],
    bam_paths=bam_paths,
    bam_column="batch",
    gtf_file=str(gtf_local),
    figsize=(6, 3),
    font_size=7,
    threshold=10,
    top_n=2,  # Automatically select top 2 most abundant transcripts
)
plt.show()