Bio Copy Number Cnv Visualization

cnvkit.py scatter sample.cnr -s sample.cns -o scatter.png            # genome-wide
cnvkit.py scatter sample.cnr -s sample.cns -c chr17 -o chr17.png      # one chromosome
cnvkit.py scatter sample.cnr -s sample.cns -v sample.vcf.gz -o baf.png  # with BAF panel
cnvkit.py diagram sample.cnr -s sample.cns -o diagram.pdf             # ideogram
cnvkit.py heatmap cohort/*.cns -d -o cohort_heatmap.pdf              # cohort, desaturated

import pandas as pd
import matplotlib.pyplot as plt

def plot_genome_profile(cnr_file, cns_file, output=None, gain=0.3, loss=-0.3):
    '''Genome-wide log2 scatter with segment overlay.'''
    cnr = pd.read_csv(cnr_file, sep='\t')
    cns = pd.read_csv(cns_file, sep='\t')
    chroms = [f'chr{i}' for i in range(1, 23)] + ['chrX', 'chrY']

    offsets, cum = {}, 0
    for c in chroms:
        sub = cnr[cnr['chromosome'] == c]
        if sub.empty:
            continue
        offsets[c] = cum
        cum += sub['end'].max()

    cnr = cnr[cnr['chromosome'].isin(offsets)].copy()
    cnr['x'] = cnr.apply(lambda r: offsets[r['chromosome']] + r['start'], axis=1)

    fig, ax = plt.subplots(figsize=(16, 4))
    ax.scatter(cnr['x'], cnr['log2'], s=1, c='0.7', alpha=0.5, rasterized=True)
    for _, seg in cns.iterrows():
        if seg['chromosome'] not in offsets:
            continue
        x0 = offsets[seg['chromosome']] + seg['start']
        x1 = offsets[seg['chromosome']] + seg['end']
        color = 'red' if seg['log2'] > gain else 'blue' if seg['log2'] < loss else '0.3'
        ax.hlines(seg['log2'], x0, x1, colors=color, linewidth=2.5)
    for c, x in offsets.items():
        ax.axvline(x, color='0.9', linewidth=0.5)
    ax.axhline(0, color='black', linewidth=0.6)
    ax.set_ylim(-2, 2)
    ax.set_ylabel('log2 copy ratio')
    ax.set_xlabel('genomic position')
    fig.tight_layout()
    if output:
        fig.savefig(output, dpi=200)
    return fig, ax

import numpy as np

def plot_log2_baf(cnr_file, baf_df, output=None):
    '''Two-panel plot: log2 copy ratio above, B-allele frequency below.
    baf_df: columns chromosome, position, baf (germline-het sites only).'''
    cnr = pd.read_csv(cnr_file, sep='\t')
    fig, (ax_cn, ax_baf) = plt.subplots(2, 1, figsize=(16, 6), sharex=True)

    ax_cn.scatter(range(len(cnr)), cnr['log2'], s=1, c='0.6', alpha=0.5)
    ax_cn.axhline(0, color='black', linewidth=0.6)
    ax_cn.set_ylabel('log2 ratio')
    ax_cn.set_ylim(-2, 2)

    # Plot BAF and its mirror; a tight band at 0.5 = balanced, split bands = imbalance/LOH
    ax_baf.scatter(range(len(baf_df)), baf_df['baf'], s=2, c='0.4', alpha=0.5)
    ax_baf.scatter(range(len(baf_df)), 1 - baf_df['baf'], s=2, c='0.4', alpha=0.5)
    ax_baf.axhline(0.5, color='black', linewidth=0.6)
    ax_baf.set_ylabel('B-allele frequency')
    ax_baf.set_ylim(0, 1)
    ax_baf.set_xlabel('het SNP index')
    fig.tight_layout()
    if output:
        fig.savefig(output, dpi=200)
    return fig

import seaborn as sns

def plot_cohort_heatmap(cns_files, bin_size=1_000_000, output=None):
    '''Recurrent-CNV heatmap across a cohort on a uniform bin grid.'''
    chroms = [f'chr{i}' for i in range(1, 23)]
    columns = []
    for c in chroms:
        columns += [(c, b) for b in range(0, 250_000_000, bin_size)]
    matrix = {}
    for f in cns_files:
        name = f.split('/')[-1].replace('.cns', '')
        cns = pd.read_csv(f, sep='\t')
        row = {}
        for c, b in columns:
            hits = cns[(cns['chromosome'] == c) &
                       (cns['start'] < b + bin_size) & (cns['end'] > b)]
            row[(c, b)] = hits['log2'].mean() if not hits.empty else 0.0
        matrix[name] = row
    df = pd.DataFrame(matrix).T
    fig, ax = plt.subplots(figsize=(14, max(4, 0.3 * len(cns_files))))
    sns.heatmap(df, cmap='RdBu_r', center=0, vmin=-1.5, vmax=1.5,
                xticklabels=False, ax=ax)
    ax.set_xlabel('genomic bin')
    ax.set_ylabel('sample')
    if output:
        fig.savefig(output, dpi=200, bbox_inches='tight')
    return fig

# GATK: denoised ratios + modeled segments with allelic info
gatk PlotModeledSegments --denoised-copy-ratios tumor.denoisedCR.tsv \
    --allelic-counts tumor.hets.tsv --segments tumor.modelFinal.seg \
    --sequence-dictionary reference.dict --output-prefix tumor -O plots/