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")
print(hb.title)         # "THE CRYSTAL STRUCTURE OF HUMAN DEOXYHAEMOGLOBIN..."
print(hb.method)        # "X-RAY DIFFRACTION"
print(hb.resolution)    # 1.74
print(hb.organism)      # "Homo sapiens"
print(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")
print(len(ids))   # number of matching structures
print(ids[0..5])  # first 5 PDB IDs

# Fetch details for the top hits
ids[0..3]
  |> map(|id| pdb_entry(id))
  |> each(|e| print(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| print(c.chain_id + ": " + c.description))
# A: HEMOGLOBIN (DEOXY) (ALPHA CHAIN)
# B: HEMOGLOBIN (DEOXY) (BETA CHAIN)
# C: HEMOGLOBIN (DEOXY) (ALPHA CHAIN)
# D: HEMOGLOBIN (DEOXY) (BETA CHAIN)

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)
print(alpha.description)   # "HEMOGLOBIN ALPHA CHAIN"
print(alpha.type)          # "polypeptide(L)"
print(alpha.weight)        # 15126.4
print(alpha.chains)        # ["A", "C"]

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)
print(len(seq))       # 141 residues
print(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
print(len(results))    # number of IDs returned (up to limit)
print(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)
let paper = pubmed_abstract("32553272")
print(paper.title)     # "Search-and-replace genome editing without..."
print(paper.journal)   # "Nature"
print(paper.year)      # 2020
print(paper.authors[0]) # "Anzalone AV"

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")
print(len(keys(gene_sets)))  # number of gene sets loaded
print(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)
  |> print()

#  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)
  |> print()

#  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...