Custom Function Definitions

Session 1 helper functions.

Most of these are adapted from an earlier gene co-expression network workshop, keeping the same names so they stay recognisable. The originals were written for co-expression networks - one node type, every edge carrying a correlation weight, always undirected. A curated knowledge graph breaks all three of those assumptions, so the ported versions here additionally handle:

• several node types (gene / disease / icd10) rather than just genes

• edges that carry no weight at all (is_a, maps_to)

• directed graphs, where connected components and clustering differ

Deliberately lightweight: matplotlib, pandas, networkx, numpy and seaborn only (numpy arrives with pandas anyway). The original module pulled in torch, dgl and astropy at import time, which is a slow and fragile thing to ask of a workshop room.

s1_helpers.ancestors_of(G: Graph, node: str, with_depth: bool = False) → list

All ancestors reachable by following is_a edges upwards, nearest first.

Breadth-first and alphabetically ordered within each level, so the result is identical on every machine and every run. That matters more than it sounds: a disease can have several parents, and iterating a set gives a different order in each Python process, so an un-sorted version prints a different answer every time the notebook is run.

With with_depth, yields (ancestor, depth) where depth is the number of is_a hops from node - a real distance up the tree, not a position in a traversal.

s1_helpers.clean_graph(G: Graph, degree_threshold: int = 1, keep_largest_component: bool = True) → Graph

Remove self-loops, isolates and low-degree nodes; optionally keep only the largest component. Directed graphs use weakly connected components.

s1_helpers.coexpression_network(expression: DataFrame, threshold: float, absolute: bool = True, correlations: DataFrame | None = None) → Graph

Build a computed network: an edge wherever two genes correlate above threshold.

Every node is a gene - one node type, unlike the knowledge graph - and every edge carries the correlation as its weight. Nodes are added for all genes, including ones that end up with no edges, because “which genes dropped out at this threshold?” is a question worth being able to ask.

absolute thresholds on r, so strong anti-correlation counts as a relationship. Pass correlations to reuse a matrix already computed - it is by far the expensive part, and the threshold sweep calls this repeatedly.

s1_helpers.coexpression_partners(correlations: DataFrame, gene_id: str, G: Graph | None = None, top_n: int = 10) → DataFrame

The genes most strongly co-expressed with one gene, strongest r first.

Pass the knowledge graph as G to get gene symbols alongside the Ensembl ids - the ranking is unreadable without them.

s1_helpers.coexpression_threshold_sweep(expression: DataFrame, thresholds, correlations: DataFrame | None = None) → DataFrame

Rebuild the network at several thresholds and report what survives each one.

The point of the table is that there is no principled place to stop. Every row is a defensible network built from identical data, and the choice of row is the analyst’s, not the data’s.

s1_helpers.correlation_matrix(expression: DataFrame) → DataFrame

Gene x gene Pearson correlation, with the diagonal zeroed.

The diagonal is every gene’s perfect correlation with itself. Left in, it dominates any “strongest partners” ranking, so it goes to zero here once rather than being special-cased at every call site.

s1_helpers.datatype_summary(evidence: DataFrame) → DataFrame

How much of each kind of evidence the graph rests on.

s1_helpers.diseases_for_genes(G: Graph, gene_ids, top_n: int = 10) → DataFrame

Given a gene list, which diseases does it touch?

This is the smallest useful knowledge-graph query: one hop out from a set of genes, counting where we land. It is the same shape as the queries an LLM agent will be asked to plan in Sessions 3 and 4.

s1_helpers.draw_network_with_node_attrs(G, node_attributes=None, communities=None, title='Network Visualization', color_attr=None, shape_attr=None, figsize=(20, 10), layout='spring', cmap_name='tab20', with_labels=False, node_size=400, seed=0)

Draw a graph with nodes coloured and/or shaped by attribute.

Adapted from the co-expression original with two fixes and one signature change:

• The original grouped nodes for drawing by shape, but built that grouping from shape_attr. Called with a color_attr and no shape_attr the grouping came out empty and no nodes were drawn at all - only edges. Colour-only calls now work.

• Attributes are read from the graph itself by default. The original required a separate node_attributes dict of the shape {attr_name: {node: value}}, which had to be kept in sync by hand. Pass one to override to colour by something not on the graph.

• An empty attribute no longer divides by zero when building the colormap.

s1_helpers.edges_of_type(G: Graph, edge_type: str) → list

All edges of a given type, as (u, v) pairs.

s1_helpers.evidence_for_pair(evidence: DataFrame, gene_id: str, disease_id: str, G: Graph | None = None) → DataFrame

Break a single gene-disease edge down into the evidence behind it.

s1_helpers.filter_by_datatype(G: Graph, evidence: DataFrame, datatypes, min_score: float = 0.0) → Graph

Keep only the associated_with edges supported by particular evidence types.

Structural edges (is_a, maps_to) are always kept - they are not the sort of claim evidence types apply to, and dropping them would tear out the hierarchy exactly as an over-eager threshold does.

Pass CAUSAL_DATATYPES to keep the edges that assert the gene has something to do with causing the disease, and watch how many disappear.

s1_helpers.gen_graph_legend(G: Graph, attr: str = 'type') → list

Legend patches for a graph coloured by attr.

The original took a parallel series of colours and zipped it against the attribute values, which relied on the two being in the same order. This derives both from the graph, so they cannot disagree.

s1_helpers.gene_by_symbol(G: Graph, symbol: str) → str

Find a gene node’s Ensembl id from its symbol, e.g. “BRCA1” -> ENSG00000012048.

Participants think in symbols and the graph is keyed by Ensembl id, so this saves a manual lookup every time. Raises rather than returning None: a typo should stop the cell, not silently produce an empty result further down.

s1_helpers.genes_for_disease(G: Graph, disease_id: str) → set

The set of gene nodes linked to a disease by an associated_with edge.

s1_helpers.get_highest_degree_nodes(G: Graph, top_n: int = 10) → DataFrame

The top_n highest-degree nodes.

Returns a DataFrame rather than the original’s list of tuples, so that the readable name sits next to the opaque id - ENSG00000012048 means little, BRCA1 means a lot.

s1_helpers.icd10_for_disease(G: Graph, disease_id: str) → dict

Find an ICD-10 code for a disease, climbing the ontology if necessary.

Returns a dict with the code, the node it was actually found on, and how many is_a steps that took. steps == 0 means the disease carries the code itself; anything higher means we inherited it from an ancestor, which is a weaker claim and should be reported as such.

s1_helpers.load_evidence(data_dir: Path | str = PosixPath('/data/session-1-data')) → DataFrame

The gene-disease edges split by kind of evidence.

One row per (gene, disease, datatype). The weight column is the association score contributed by that evidence type alone, so the same gene-disease pair appears once per datatype supporting it.

s1_helpers.load_expression(data_dir: Path | str = PosixPath('/data/session-1-data')) → DataFrame

Load the committed expression matrix: patients (rows) x genes (columns).

This is the TCGA-BRCA transcriptomics view used in Session 2, cut down to the genes that are also nodes in the knowledge graph (737 of 760) so that both networks describe the same genes. Values are log2-scale and already library-size normalised, so correlations can be taken directly.

Generated ahead of the workshop by build_coexpression_data.py, which lives in the instructors’ repository.

s1_helpers.load_kg(data_dir: Path | str = PosixPath('/data/session-1-data'), directed: bool = False) → Graph

Build the Session 1 knowledge graph from the generated CSVs.

Node attributes: type, name, extra. Edge attributes: type, and weight / evidence on associated_with only.

The underlying edges are directional (gene -> disease, disease -> parent, disease -> icd10). Undirected is the default because degree, hubs and connected components all behave more intuitively that way; pass directed=True when the direction is the point.

s1_helpers.load_omics(path, layers=('transcriptomics',))

Load the Session 2 omics pickle, keeping only the layers asked for.

The full file is ~900 MB, most of it the 200,000-probe methylation matrix, so the default keeps transcriptomics only. meta (the PAM50 subtype per patient) is always returned.

Returns (dict_of_layers, meta_series).

s1_helpers.map_genes_to_kg(G: Graph, gene_ids) → DataFrame

Look a list of gene ids up in the knowledge graph.

Accepts versioned or unversioned Ensembl ids. Returns one row per input id with whether it was found and, if so, its symbol and degree - so an omics gene list becomes an entry point into the graph.

s1_helpers.nodes_of_type(G: Graph, node_type: str) → list

All node ids of a given type.

s1_helpers.plot_degree_distribution(G: Graph, bins: int = 30, by_type: bool = False)

Degree distribution, optionally split by node type.

s1_helpers.print_graph_info(G: Graph) → None

Print basic information about a graph.

Unlike the co-expression original this is safe on directed graphs, and it breaks the counts down by node and edge type - on a typed graph the totals alone hide most of what is interesting.

s1_helpers.remove_by_degree(G: Graph, min_degree: int) → Graph

Drop nodes whose degree is below min_degree.

s1_helpers.shared_gene_projection(G: Graph, disease_ids=None, min_shared: int = 1) → Graph

Project the bipartite gene-disease graph down onto diseases alone.

Two diseases are joined if they share at least min_shared genes. The edge carries both the raw count (n_shared) and the Jaccard index (weight), because the two tell different stories: a well-studied disease shares many genes with everything simply by having many genes.

This is the standard way to turn a two-mode network into a one-mode one, and it is where “which diseases resemble each other?” becomes a question about graph structure rather than about biology directly.

s1_helpers.shared_genes_between(G: Graph, a: str, b: str) → list

Gene symbols shared by two diseases, for eyeballing a single pair.

s1_helpers.strip_ensembl_version(ids) → list

Turn versioned Ensembl ids into bare ones: ENSG00000012048.23 -> ENSG00000012048.

The TCGA matrices carry the version suffix; Open Targets does not. Joining the two without this step matches exactly nothing, silently.

s1_helpers.threshold_sparsification(G: Graph, threshold: float, keep_unweighted: bool = True) → Graph

Drop edges whose weight is below threshold.

keep_unweighted is the knowledge-graph-specific part: is_a and maps_to edges carry no weight, and silently deleting them (which the co-expression version would, by treating a missing weight as 0) would tear the disease hierarchy out of the graph. Keep them by default.

s1_helpers.top_percentage_sparsification(G: Graph, top_percentage: float, keep_unweighted: bool = True) → Graph

Keep only the top top_percentage % of weighted edges.

s1_helpers.visualise_edge_weight_distribution(G: Graph, bins: int = 30)

Distribution of edge weights.

Only edges that actually carry a weight are plotted; the original indexed G[u][v][‘weight’] directly and would raise a KeyError on the first is_a edge it met.

s1_helpers.visualise_graph(G: Graph, title: str = 'Knowledge Graph', colour_by: str | None = 'type', figsize=(10, 10), seed: int = 0)

Quick look at a graph, nodes coloured by an attribute.

The co-expression version drew every node the same blue, which is the right call for a co-expression network and the wrong one here - node type is the first thing we want to see.