PDB, PubMed & Enrichment

Fetches an entry from the RCSB Protein Data Bank. Returns a record with structural metadata including the title, experimental method, resolution, source organism, and release date.

# requires: internet connection (PDB API)
let hb = pdb_entry("4HHB")
println(hb.title)         # "THE CRYSTAL STRUCTURE OF HUMAN DEOXYHAEMOGLOBIN..."
println(hb.method)        # "X-RAY DIFFRACTION"
println(hb.resolution)    # 1.74
println(hb.organism)      # "Homo sapiens"
println(hb.release_date)  # "1984-07-17"

Searches the PDB by keyword. Returns a list of PDB IDs matching the query. Use the results with pdb_entry() to fetch full metadata.

# requires: internet connection (PDB API)
let ids = pdb_search("insulin receptor")
println(len(ids))   # number of matching structures
println(ids[0..5])  # first 5 PDB IDs

# Fetch details for the top hits
ids[0..3]
  |> map(|id| pdb_entry(id))
  |> each(|e| println(e.title + " — " + str(e.resolution) + " A"))

Lists the polymer chains in a structure. Each entry is a record with chain ID, entity ID, and description.

# requires: internet connection (PDB API)
let chains = pdb_chains("4HHB")
chains |> each(|c| println(str(c.entity_id) + ": " + c.description))
# 1: Hemoglobin subunit alpha
# 2: Hemoglobin subunit beta

Fetches detailed information about a specific entity within a structure. Returns a record with the entity description, type, molecular weight, chain assignments, and organism source.

# requires: internet connection (PDB API)
let alpha = pdb_entity("4HHB", 1)
println(alpha.description)   # "Hemoglobin subunit alpha"
println(alpha.entity_type)   # "polypeptide(L)"
println(len(alpha.sequence)) # 141

Retrieves the amino acid sequence for a given entity. Returns a plain protein sequence string suitable for alignment or further analysis.

# requires: internet connection (PDB API)
let seq = pdb_sequence("4HHB", 1)
println(len(seq))       # 141 residues
println(seq[0..20])     # "VLSPADKTNVKAAWGKVGAH"

Searches PubMed via the NCBI E-utilities. Returns a record with count (total matches in the database) and ids (list of PMIDs for the returned page). The optional second argument limits how many IDs are returned (default 20).

# requires: internet connection (NCBI PubMed API)
let results = pubmed_search("CRISPR Cas9 therapy", 5)
# pubmed_search returns a list of PMID strings directly
println(len(results))    # number of IDs returned (up to limit)
println(results)         # ["39847123", "39841056", "39838901", "39835444", "39831200"]

Fetches the abstract and metadata for a single PubMed article. Returns a record with title, abstract, authors, journal, year, and doi.

# requires: internet connection (NCBI PubMed API)
# pubmed_abstract returns the MEDLINE record as text, not a parsed record.
let paper = pubmed_abstract("32553272")
println(paper |> lines() |> first())
# 1. Cell. 2020 Jul 9;182(1):145-161.e23. doi: 10.1016/j.cell.2020.05.021. Epub

# Pull out a field with a plain text search.
let title_line = paper |> lines() |> find(|l| starts_with(l, "Search-and-replace"))
println(title_line)

Reads a Gene Matrix Transposed (GMT) file. Each line defines one gene set as: set_name\tdescription\tgene1\tgene2\t.... Returns a record mapping set names to their gene lists.

# GMT format (tab-separated):
#   GO:0006915  Apoptotic process  BAX  BCL2  CASP3  CASP9  TP53
#   GO:0007049  Cell cycle         CDK1 CDK2  CCNB1  CCND1  RB1

let gene_sets = read_gmt("c5.go.bp.v2024.1.symbols.gmt")
println(len(keys(gene_sets)))  # number of gene sets loaded
println(gene_sets["GO:0006915"])  # ["BAX", "BCL2", "CASP3", "CASP9", "TP53"]

Performs over-representation analysis (ORA) using the hypergeometric test. Tests whether each gene set is enriched in your gene list relative to the background. Returns a table with columns: term, overlap, set_size, p_value, fdr (Benjamini-Hochberg adjusted), and genes (the overlapping gene symbols). Rows are sorted by p-value ascending.

let de_genes = ["TP53", "BRCA1", "CDK2", "CASP3", "BAX", "BCL2", "MYC"]
let gene_sets = read_gmt("c5.go.bp.v2024.1.symbols.gmt")
let background = read_lines("all_expressed_genes.txt")

let ora = enrich(de_genes, gene_sets, background)
ora
  |> filter(|r| r.fdr < 0.05)
  |> println()

#  term           | overlap | set_size | p_value  | fdr     | genes
# ----------------+---------+----------+----------+---------+-----------------------
#  GO:0006915     | 4       | 312      | 1.2e-06  | 3.4e-04 | TP53,BAX,BCL2,CASP3
#  GO:0008219     | 3       | 198      | 8.7e-05  | 0.012   | TP53,BAX,CASP3
#  ...

Performs Gene Set Enrichment Analysis on a pre-ranked gene list. The ranked_genes table must have gene and score columns (e.g., log2 fold change or signed -log10(p)). Uses a permutation test (1000 permutations) to compute enrichment scores and significance. Returns a table with: term, es (enrichment score), nes (normalized ES), p_value, fdr, and leading_edge (the core enrichment genes).

# Prepare ranked gene list (e.g., from differential expression)
let de_results = tsv("de_results.tsv")
let ranked = de_results
  |> mutate("score", |r| r.log2fc * -log10(r.pvalue))
  |> select("gene", "score")
  |> sort_by(|r| -r.score)

let gene_sets = read_gmt("h.all.v2024.1.symbols.gmt")

let gsea_results = gsea(ranked, gene_sets)
gsea_results
  |> filter(|r| r.fdr < 0.25)
  |> println()

#  term                      | es    | nes   | p_value | fdr   | leading_edge
# ---------------------------+-------+-------+---------+-------+-------------------
#  HALLMARK_APOPTOSIS        | 0.62  | 2.14  | 0.0     | 0.002 | BAX,CASP3,TP53...
#  HALLMARK_P53_PATHWAY      | 0.58  | 2.01  | 0.001   | 0.008 | CDKN1A,MDM2,BAX...
#  HALLMARK_MYC_TARGETS_V1   | -0.45 | -1.87 | 0.002   | 0.015 | NOP56,RRP1,NPM1...