bio-clinical-databases-tumor-mutational-burden

Version Compatibility

Reference examples tested with: Ensembl VEP 111+, SnpEff 5.2+, pandas 2.2+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed

package and adapt the example to match the actual API rather than retrying.

Tumor Mutational Burden

"Calculate TMB from my tumor sequencing data" → Compute tumor mutational burden as nonsynonymous coding mutations per megabase with proper panel normalization for immunotherapy eligibility assessment.

• Python: cyvcf2 for VCF parsing + variant counting per panel region

TMB Calculation from VCF (Ensembl VEP 111+)

Goal: Calculate tumor mutational burden as nonsynonymous coding mutations per megabase from a somatic VCF.

Approach: Iterate through VCF variants, filter for coding nonsynonymous consequences via VEP/SnpEff annotations, and divide count by panel size.

from cyvcf2 import VCF

def calculate_tmb(vcf_path, panel_size_mb):
    '''Calculate TMB (mutations per megabase)

    Args:
        vcf_path: Path to somatic VCF
        panel_size_mb: Capture region size in megabases

    Returns:
        TMB value (mutations/Mb)
    '''
    vcf = VCF(vcf_path)
    mutation_count = 0

    for variant in vcf:
        # Count nonsynonymous coding mutations
        # Adjust filters based on VCF annotation format
        if is_coding_nonsynonymous(variant):
            mutation_count += 1

    tmb = mutation_count / panel_size_mb
    return tmb

def is_coding_nonsynonymous(variant):
    '''Check if variant is coding nonsynonymous

    Adjust logic based on your VCF annotation tool:
    - VEP: CSQ field
    - SnpEff: ANN field
    - Funcotator: FUNCOTATION field
    '''
    # Example for VEP annotation
    csq = variant.INFO.get('CSQ', '')
    if not csq:
        return False

    # Check consequence types
    nonsynonymous = ['missense_variant', 'nonsense', 'frameshift',
                     'inframe_insertion', 'inframe_deletion', 'stop_gained',
                     'stop_lost', 'start_lost']

    for transcript in csq.split(','):
        fields = transcript.split('|')
        consequence = fields[1] if len(fields) > 1 else ''
        if any(ns in consequence for ns in nonsynonymous):
            return True
    return False

Panel-Specific TMB (Ensembl VEP 111+)

Goal: Calculate TMB normalized to known gene panel capture region sizes.

Approach: Look up the panel's megabase coverage from a reference table and pass to the TMB calculator.

# Common panel sizes (in megabases)
# Check your specific panel's capture region size
PANEL_SIZES_MB = {
    'FoundationOne CDx': 0.8,
    'MSK-IMPACT': 1.14,
    'TruSight Oncology 500': 1.94,
    'Oncomine Comprehensive': 1.5,
    'WES (exome)': 30.0,  # Approximate coding region
    'WGS': 3000.0,        # Approximate
}

def calculate_tmb_panel(vcf_path, panel_name):
    '''Calculate TMB for known panel'''
    if panel_name not in PANEL_SIZES_MB:
        raise ValueError(f'Unknown panel: {panel_name}')
    return calculate_tmb(vcf_path, PANEL_SIZES_MB[panel_name])

TMB with Variant Filtering (Ensembl VEP 111+)

Goal: Calculate TMB with quality and germline filters to reduce false positives.

Approach: Apply VAF, depth, and gnomAD population frequency filters before counting coding nonsynonymous variants.

def calculate_tmb_filtered(vcf_path, panel_size_mb, min_vaf=0.05, min_depth=100):
    '''Calculate TMB with quality filters

    Args:
        vcf_path: Path to somatic VCF
        panel_size_mb: Panel size in Mb
        min_vaf: Minimum variant allele frequency (default 5%)
        min_depth: Minimum read depth (default 100)

    Filters:
    - VAF >= 5%: Reduce false positives from sequencing errors
    - Depth >= 100: Ensure reliable variant calls
    - Exclude known polymorphisms (gnomAD AF > 1%)
    - Include only coding nonsynonymous
    '''
    vcf = VCF(vcf_path)
    mutation_count = 0

    for variant in vcf:
        # Quality filters
        depth = variant.INFO.get('DP', 0)
        vaf = get_vaf(variant)

        if depth < min_depth:
            continue
        if vaf < min_vaf:
            continue

        # Exclude germline polymorphisms
        gnomad_af = variant.INFO.get('gnomAD_AF', 0)
        if gnomad_af > 0.01:
            continue

        # Count coding nonsynonymous
        if is_coding_nonsynonymous(variant):
            mutation_count += 1

    return mutation_count / panel_size_mb

def get_vaf(variant):
    '''Extract variant allele frequency from variant'''
    # Format depends on caller (e.g., Mutect2, Strelka)
    # Mutect2 format: AD field in genotype
    try:
        ad = variant.format('AD')[0]  # First sample
        if sum(ad) > 0:
            return ad[1] / sum(ad)
    except:
        pass
    return 0

Clinical TMB Thresholds (Ensembl VEP 111+)

Goal: Classify a TMB value as TMB-High or TMB-Low based on clinical cutoffs.

Approach: Compare the TMB value against FDA-approved or study-specific thresholds (10, 16, or 20 mut/Mb).

def classify_tmb(tmb_value, threshold='FDA'):
    '''Classify TMB as high or low

    Clinical thresholds:
    - FDA (pembrolizumab): 10 mut/Mb
    - ESMO: 10 mut/Mb
    - Some studies use 16, 20 mut/Mb for specific cancers

    Note: Panel-specific thresholds may differ
    '''
    thresholds = {
        'FDA': 10,
        'conservative': 16,
        'strict': 20
    }

    cutoff = thresholds.get(threshold, 10)

    if tmb_value >= cutoff:
        return 'TMB-High'
    else:
        return 'TMB-Low'

# Example
tmb = 12.5
status = classify_tmb(tmb)
print(f'TMB: {tmb} mut/Mb -> {status}')

TMB by Variant Type (Ensembl VEP 111+)

Goal: Break down TMB by mutation type (missense, nonsense, frameshift, etc.) for detailed characterization.

Approach: Classify each variant by consequence type, count per category, and compute TMB from nonsynonymous subtotal.

def detailed_tmb_analysis(vcf_path, panel_size_mb):
    '''Calculate TMB broken down by variant type'''
    vcf = VCF(vcf_path)

    counts = {
        'missense': 0,
        'nonsense': 0,
        'frameshift': 0,
        'inframe_indel': 0,
        'splice': 0,
        'synonymous': 0,
        'other': 0
    }

    for variant in vcf:
        vtype = classify_variant_type(variant)
        counts[vtype] = counts.get(vtype, 0) + 1

    # TMB typically excludes synonymous
    nonsynonymous_count = sum(v for k, v in counts.items()
                               if k != 'synonymous' and k != 'other')

    results = {
        'counts': counts,
        'total_nonsynonymous': nonsynonymous_count,
        'tmb': nonsynonymous_count / panel_size_mb,
        'panel_size_mb': panel_size_mb
    }
    return results

TMB vs MSI Comparison (Ensembl VEP 111+)

Goal: Assess concordance between TMB status and microsatellite instability for immunotherapy biomarker evaluation.

Approach: Cross-tabulate TMB-High/Low with MSI-H/MSS to identify concordant and discordant cases.

def tmb_msi_concordance(tmb_value, msi_status):
    '''Compare TMB with MSI status

    MSI-H tumors typically have high TMB (>10-20 mut/Mb)
    TMB-H and MSI-H are correlated but not identical:
    - ~80% MSI-H are TMB-H
    - Many TMB-H are MSS (especially smoking-related)

    Both predict immunotherapy response
    '''
    tmb_high = tmb_value >= 10

    if msi_status == 'MSI-H' and tmb_high:
        return 'Concordant TMB-H/MSI-H'
    elif msi_status == 'MSI-H' and not tmb_high:
        return 'Discordant MSI-H/TMB-L (uncommon)'
    elif msi_status == 'MSS' and tmb_high:
        return 'TMB-H/MSS (e.g., smoking-related)'
    else:
        return 'TMB-L/MSS'

Batch TMB Calculation

Goal: Calculate TMB for an entire cohort of samples and export results with clinical classification.

Approach: Iterate over VCF files in a directory, compute filtered TMB for each, and collect into a summary DataFrame.

import pandas as pd
from pathlib import Path

def batch_tmb(vcf_dir, panel_size_mb, output_file):
    '''Calculate TMB for multiple samples'''
    results = []

    for vcf_path in Path(vcf_dir).glob('*.vcf.gz'):
        sample_id = vcf_path.stem.replace('.vcf', '')
        tmb = calculate_tmb_filtered(str(vcf_path), panel_size_mb)
        status = classify_tmb(tmb)

        results.append({
            'sample': sample_id,
            'tmb': round(tmb, 2),
            'status': status
        })

    df = pd.DataFrame(results)
    df.to_csv(output_file, index=False)
    return df

Related Skills

• variant-calling/somatic-variant-calling - Input variants
• variant-calling/clinical-interpretation - ACMG/AMP classification
• variant-calling/variant-annotation - VEP/SnpEff annotation