Python for Genomics

intermediate

Python Genomics Toolkit

Essential Python Libraries for Genomics

🧬BiopythonFASTA/GenBank I/O; BLAST; sequence alignment; phylogenetics

📂pysamBAM/SAM/CRAM/VCF read/write; Python wrapper for htslib

🔬cyvcf2Fast VCF parsing; 4x faster than PyVCF; all VCF4.3 features

📊pandasDataFrame for variant tables; groupby, merge, pivot operations

📈matplotlib/seabornVolcano plots, coverage plots, heatmaps, Manhattan plots

🐍snakemakePython-based workflow; rules with input/output; cluster/cloud

Essential Libraries

Biopython - sequence I/O, alignment, BLAST, phylogenetics; Bio.SeqIO, Bio.AlignIO
pysam - BAM/SAM/CRAM/BCF/VCF reading/writing; Python wrapper for htslib
cyvcf2 - fast VCF/BCF parsing; 4× faster than PyVCF; supports VCF4.3
pandas - dataframe operations; essential for variant table manipulation
numpy/scipy - numerical computation; statistics for QC metrics
matplotlib/seaborn - publication-quality plots; heatmaps, volcano plots, coverage plots
plotly - interactive visualisations; dashboards; widely used in genomics tools
pyranges - fast genomic interval operations; BED-like operations in Python

pysam - BAM Operations

code

import pysam

# Open BAM file
bam = pysam.AlignmentFile("sample.bam", "rb")

# Count reads in a region
count = bam.count("chr17", 43044294, 43125364)
print(f"BRCA1 region reads: {count}")

# Iterate over reads
for read in bam.fetch("chr17", 43044294, 43125364):
    if read.is_unmapped or read.is_duplicate:
        continue
    print(read.query_name, read.mapping_quality, read.cigarstring)

# Per-base coverage
for col in bam.pileup("chr17", 43044294, 43044300, min_base_quality=20):
    print(f"pos {col.pos}: depth {col.nsegments}")

bam.close()

cyvcf2 - VCF Parsing

code

from cyvcf2 import VCF, Writer
import pandas as pd

# Parse and filter variants
records = []
for v in VCF("annotated.vcf.gz"):
    # Skip non-PASS
    if v.FILTER and v.FILTER != "PASS":
        continue
    # Get annotation fields
    gene = v.INFO.get("Gene_refGene", ".")
    func = v.INFO.get("Func_refGene", ".")
    gnomad_af = v.INFO.get("gnomAD_exome_AF", 0.0) or 0.0
    cadd = v.INFO.get("CADD_phred", 0.0) or 0.0
    clinvar = v.INFO.get("CLNSIG", ".")

    # Filter: rare + potentially damaging
    if gnomad_af < 0.001 and cadd > 15:
        records.append({
            "CHROM": v.CHROM, "POS": v.POS,
            "REF": v.REF, "ALT": v.ALT[0],
            "Gene": gene, "Function": func,
            "gnomAD_AF": gnomad_af, "CADD": cadd,
            "ClinVar": clinvar,
            "GT": v.genotypes[0][:2]  # first sample
        })

df = pd.DataFrame(records)
df.to_csv("filtered_variants.csv", index=False)
print(f"Retained {len(df)} variants")

Automation with Snakemake

code

# Snakefile - simple WES pipeline
SAMPLES = ["sample1", "sample2", "sample3"]

rule all:
    input:
        expand("vcf/{sample}.vcf.gz", sample=SAMPLES)

rule align:
    input:
        r1 = "fastq/{sample}_R1.fastq.gz",
        r2 = "fastq/{sample}_R2.fastq.gz"
    output: "bam/{sample}.sorted.bam"
    threads: 16
    shell:
        "bwa-mem2 mem -t {threads} hg38.fa {input.r1} {input.r2} "
        "| samtools sort -@ 8 -o {output}"

rule call_variants:
    input: "bam/{sample}.recal.bam"
    output: "vcf/{sample}.vcf.gz"
    shell:
        "gatk HaplotypeCaller -R hg38.fa -I {input} "
        "-O {output} -ERC GVCF"

# Run with 4 parallel jobs
# snakemake --cores 32 --jobs 4 --use-conda

Long-read Sequencing R for Genomics