Part 2 — MOFA tools over the Model Context Protocol (MCP)#
In Part 1, the agent loop called the MOFA functions directly from Python. In this part, those same functions are provided by a separate program, an MCP server, and the agent accesses them by sending requests through it.
We start with our own MCP server, built around the same MOFA functions, and check that it gives the same results as Part 1. We then connect the agent to an MCP server developed by someone else, so that it can use those additional tools alongside our own within the same conversation. Finally, we look at how an MCP server like ours is actually built.
Throughout, the analysis itself stays fixed: the same MOFA functions, the same fitted model, and the same biological questions. What changes is only how the agent reaches the tools.
Learning objectives#
By the end of this notebook you should be able to:
Explain what an MCP server is, and why moving a tool onto one changes where it runs and who can reach it, rather than what it does.
Connect an agent to an MCP server and use its tools exactly as you used local tools in Part 1.
Connect the same agent to a server written and run by someone else, and use both sets of tools in one conversation.
Describe what information leaves your machine when an external server is involved.
What is MCP?#
An MCP server is a program that makes a set of tools available to other programs. It runs on its own, holds whatever data those tools need, and waits for requests. An agent that wants to use one of those tools does not import any code: it connects to the server, asks what tools are available, and asks for one to be run.
For that exchange to work, the agent and the server have to agree on how to ask those questions. The Model Context Protocol (MCP) is that agreement. It defines how a client asks a server what it offers, how it requests that a tool be run, and how results come back. Because the protocol is the same everywhere, any MCP-compatible agent can use any MCP server.
You may ask: why add this extra layer at all? It is useful because it separates the tool provider from the application using the tools. The same set of tools can then be reused by different notebooks or agents without each one having to integrate the underlying code itself, and access to the data and functions can stay behind a single controlled interface. Just as importantly, the same standard lets your agent connect to tools built by other people without needing a new custom integration for each one.
The server we use here is server/mofa_mcp_server.py. We wrote it for this practical: it holds the same eight analysis functions from Part 1, backed by the same cached MOFA model. This notebook is about using it. The last section shows how it was built.
0. Setup#
Load the API key and locate the MCP server script.
from pathlib import Path
import os, sys, json
from dotenv import load_dotenv
SESSION_DIR = Path(os.environ.get(
"ECCB_SESSION3_DIR",
Path.home() / "ECCB2026_TEST/sessions/session-3-agentic-llm-workflows"))
assert SESSION_DIR.exists(), f"Session directory not found: {SESSION_DIR}"
DOTENV_PATH = "/.env"
DATA_DIR = Path("/data") if Path("/data").exists() else SESSION_DIR / "data" # shared mount on workshop servers
PROJECT_ROOT = SESSION_DIR # anchors imports (src/...) and paths below
sys.path.insert(0, str(PROJECT_ROOT))
load_dotenv(DOTENV_PATH)
if not os.environ.get("ANTHROPIC_API_KEY"):
raise RuntimeError("ANTHROPIC_API_KEY not found. Add it to a .env in the project root.")
SERVER_PATH = PROJECT_ROOT / "server" / "mofa_mcp_server.py"
assert SERVER_PATH.exists(), f"MCP server not found at {SERVER_PATH}"
print("API key loaded:", bool(os.environ.get("ANTHROPIC_API_KEY")))
print("MCP server :", SERVER_PATH.name)
API key loaded: True
MCP server : mofa_mcp_server.py
1. Connecting to the server#
This section sets up the connection between the notebook and the MCP server, then rebuilds the same agent loop used in Part 1.
MultiServerMCPClient is an MCP client: it launches the server, performs the exchange described above, and hands back the server’s tools in the form LangChain expects.
from langchain_mcp_adapters.client import MultiServerMCPClient
client = MultiServerMCPClient({
"mofa": {
"command": sys.executable,
"args": [str(SERVER_PATH)],
"transport": "stdio", # the server runs as a subprocess of this notebook
}
})
# Ask the server what it offers. Each tool arrives with its name, its arguments
# and its description, the same three things Claude saw in Part 1.
tools = await client.get_tools()
tools_by_name = {t.name: t for t in tools}
print(f"{len(tools)} tools from the MOFA server:")
for t in tools:
print(f" {t.name}")
8 tools from the MOFA server:
data_summary
split_summary
active_factors
factor_view_r2
factor_subtype_association
top_features_for_factor
classify_subtype_from_factors
train_vs_test_subtype_association
Binding and the agent loop are unchanged from Part 1. The only difference is that the tools now live in another program, so calling one means sending a request and waiting for the reply — which is what await and ainvoke mark on the tool call.
from langchain_anthropic import ChatAnthropic
from langchain_core.messages import AIMessage, HumanMessage, SystemMessage, ToolMessage
MODEL = "claude-haiku-4-5"
llm = ChatAnthropic(model=MODEL, temperature=0)
llm_with_tools = llm.bind_tools(tools)
# Identical to the system prompt in Part 1.
SYSTEM = ("You are a computational-biology assistant analysing a fitted MOFA model "
"of TCGA breast-cancer multi-omics data. Use the tools to gather evidence; "
"ground every quantitative claim in tool results and name which tool you used. "
"Factors are named 'Factor1'..'Factor10'; subtypes are PAM50 (LumA, LumB, "
"Basal, Her2, Normal).")
async def run_agent(question: str, max_steps: int = 6, verbose: bool = True) -> AIMessage:
"""
Identical to Part 1's run_agent, except that it waits for the server to reply.
"""
messages = [SystemMessage(content=SYSTEM), HumanMessage(content=question)]
for step in range(max_steps):
# Ask the model what to do next.
model_response = await llm_with_tools.ainvoke(messages)
messages.append(model_response)
# If no tool was requested, the model has finished its answer. Return it.
if not model_response.tool_calls:
return model_response
# Otherwise, run each requested tool
for call in model_response.tool_calls:
if verbose:
print(f"[step {step}] -> {call['name']}({call['args']})")
# ...except that the tool does not run here: the request goes to the
# server, which runs the function and sends the result back.
result = await tools_by_name[call["name"]].ainvoke(call["args"])
messages.append(ToolMessage(content=json.dumps(result, default=str),
tool_call_id=call["id"]))
raise RuntimeError(f"Agent did not finish within {max_steps} steps.")
print("Bound", len(tools), "MCP tools to", MODEL)
Bound 8 MCP tools to claude-haiku-4-5
Note — why does run_agent() need async, await and ainvoke now?
In Part 1, tools were called like any other Python function:
result = tools_by_name[call["name"]].invoke(call["args"])In Part 2, the same line reads:
result = await tools_by_name[call["name"]].ainvoke(call["args"])This is because the agent is no longer calling a function in this notebook. An MCP tool may be running in another program, so the notebook sends it a request and waits for the result to come back.
Python has a way of handling that kind of waiting without stopping everything else: asynchronous programming (async). An asynchronous function can pause while it waits for a result and let Python get on with other work in the meantime — useful whenever a program is talking to several things at once.
That is what the two new keywords are. The a in ainvoke means asynchronous: it starts the operation rather than completing it. await then means, roughly, pause here until the result is ready, then continue.
Python only allows await inside an asynchronous function, so run_agent() has to change too — from def run_agent(...) to async def run_agent(...).
We gain nothing from this here: the notebook asks one question at a time and has no other work to get on with while it waits. The library that connects us to MCP servers is simply written this way, because it is built for programs handling several servers and many requests at once.
2. The same questions, answered through the server#
Below are two of the questions from Part 1: the first one, which needs a single tool, and the flagship question, which needs several. The tool names in the trace and the answers themselves should match what you saw in Part 1.
# Q1 from Part 1: one tool
query = "How many patients and features are in each omics view, and how many patients per PAM50 subtype?"
answer = await run_agent(query, max_steps=12)
print(answer.content)
[step 0] -> data_summary({})
Based on the **data_summary** tool results, here's the breakdown:
## Patients and Features per Omics View
| View | Patients | Features |
|------|----------|----------|
| **Transcriptomics** | 500 | 29,995 |
| **Proteomics** | 500 | 464 |
| **Methylation** | 500 | 200,000 |
## Patients per PAM50 Subtype
| Subtype | Count |
|---------|-------|
| **LumA** | 237 |
| **LumB** | 100 |
| **Basal** | 97 |
| **Her2** | 41 |
| **Normal** | 25 |
| **Total** | **500** |
The dataset is dominated by Luminal A (LumA) subtype (47.4% of patients), followed by Luminal B and Basal subtypes (~20% each), with Her2 and Normal subtypes being less common. Methylation is by far the largest omics view with 200,000 features, while proteomics is the smallest with 464 features.
# Q5 from Part 1, the flagship: several tools, the second chosen after seeing the first result
query = "Which MOFA factor is most associated with breast-cancer subtype, and which transcriptomic features most strongly drive it?"
answer = await run_agent(query, max_steps=12)
print(answer.content)
[step 0] -> factor_subtype_association({})
[step 1] -> top_features_for_factor({'factor': 'Factor2', 'view': 'transcriptomics', 'n': 5})
## Summary
**Factor2** is the MOFA factor most strongly associated with breast-cancer PAM50 subtype, with an eta-squared of **0.824** (based on `factor_subtype_association`). This is substantially higher than all other factors, indicating it captures subtype-defining variation.
### Top Transcriptomic Drivers of Factor2:
**Positive weights (strongest drivers):**
1. ENSG00000160182.3 (weight: 0.589)
2. ENSG00000173467.9 (weight: 0.568)
3. ENSG00000160180.15 (weight: 0.532)
4. ENSG00000106541.12 (weight: 0.493)
5. ENSG00000235687.9 (weight: 0.492)
**Negative weights:**
1. ENSG00000166535.20 (weight: -0.365)
2. ENSG00000185686.18 (weight: -0.36)
3. ENSG00000164434.12 (weight: -0.352)
4. ENSG00000186832.9 (weight: -0.349)
5. ENSG00000102243.13 (weight: -0.344)
These genes represent the primary transcriptomic signatures that distinguish breast-cancer subtypes in this MOFA model.
Same tools, same order, same answer. The functions ran in another process, and neither Claude nor the agent loop had to be told anything about that.
Now write your own question. Anything the eight tools can answer between them — a good one needs more than a single tool, so that you can watch Claude decide what to ask for next.
#################################################
# Write a question that needs more than one tool.
# Additional clue: look at Part 1 notebook for ideas
query = "Is the factor that explains the most variance also the one that best separates subtypes? " \
"Compare the R2 ranking with the subtype association ranking and explain any mismatch."
answer = await run_agent(query, max_steps=12)
print(answer.content)
#################################################
[step 0] -> active_factors({})
[step 0] -> factor_subtype_association({})
Excellent! Now I can provide a detailed comparison. Here's what the data shows:
## R² Ranking (Variance Explained)
1. **Factor1**: 34.64%
2. **Factor2**: 29.12%
3. **Factor3**: 17.22%
4. **Factor4**: 10.65%
5. **Factor5**: 8.87%
## Subtype Association Ranking (η² with PAM50)
1. **Factor2**: 0.824
2. **Factor6**: 0.373
3. **Factor1**: 0.281
4. **Factor4**: 0.181
5. **Factor8**: 0.154
## Key Finding: **Significant Mismatch**
**No, the factor explaining the most variance is NOT the one that best separates subtypes.** In fact, there's a striking reversal:
- **Factor1** ranks **#1 in variance explained** (34.64% R²) but only **#3 in subtype association** (η² = 0.281)
- **Factor2** ranks **#2 in variance explained** (29.12% R²) but is **#1 in subtype association** (η² = 0.824)
## Explanation of the Mismatch
This mismatch reveals an important principle in multi-omics analysis:
1. **Factor1 captures general biological variation** that is not specifically related to PAM50 subtypes. It explains the largest amount of overall variance across the omics views, but this variance is driven by features that vary across all samples regardless of subtype.
2. **Factor2 is subtype-specific** despite explaining less total variance. Its 0.824 η² indicates that Factor2 is highly discriminative for PAM50 subtypes—the variation it captures is strongly structured by subtype membership, even though it accounts for less total variance.
3. **Factor6** (η² = 0.373) is the second-best subtype separator but ranks only 6th in total variance (6.50% R²), further illustrating that subtype-discriminative power is independent of total variance explained.
This demonstrates that **variance and biological relevance are orthogonal properties** in MOFA: a factor can be biologically meaningful for a specific phenotype (subtypes) without being the largest source of overall variation in the data.
3. Connecting to somebody else’s server#
Everything so far used a server we wrote, fronting our own model. The more useful case is a server somebody else runs, exposing data we do not have.
BioMCP is one: an open-source MCP server covering around fifteen public biomedical sources — PubMed, ClinicalTrials.gov, ClinVar, MyGene.info and others — behind one set of tools. Because it speaks the same protocol, connecting to it takes the same few lines as connecting to our own server, and the agent can use both at once.
It is already installed in this environment. Elsewhere, you would install it with:
pip install biomcp-python
The next cell finds the BioMCP command and adds it as a second entry alongside our own server.
import shutil, subprocess
# biomcp is installed in the same env as this kernel, but the kernel's PATH may
# not include that env's bin/ (e.g. Jupyter server running from another env),
# so fall back to the directory of the kernel's own interpreter.
BIOMCP_BIN = shutil.which("biomcp") or str(Path(sys.executable).with_name("biomcp"))
assert Path(BIOMCP_BIN).exists(), "biomcp not found -- pip install biomcp-python into this kernel's env"
subcommand = "run" # stdio-transport server on current biomcp versions
client = MultiServerMCPClient({
"mofa": {"command": sys.executable, "args": [str(SERVER_PATH)], "transport": "stdio"},
"biomcp": {"command": BIOMCP_BIN, "args": [subcommand], "transport": "stdio"},
})
mofa_tool_names = set(tools_by_name)
all_tools = await client.get_tools()
all_tools_by_name = {t.name: t for t in all_tools}
print(f"{len(all_tools)} tools in total")
print(" ours :", sorted(mofa_tool_names))
print(" theirs :", sorted(set(all_tools_by_name) - mofa_tool_names))
44 tools in total
ours : ['active_factors', 'classify_subtype_from_factors', 'data_summary', 'factor_subtype_association', 'factor_view_r2', 'split_summary', 'top_features_for_factor', 'train_vs_test_subtype_association']
theirs : ['alphagenome_predictor', 'article_getter', 'article_searcher', 'disease_getter', 'drug_getter', 'enrichr_analyzer', 'fetch', 'gene_getter', 'nci_biomarker_searcher', 'nci_disease_searcher', 'nci_intervention_getter', 'nci_intervention_searcher', 'nci_organization_getter', 'nci_organization_searcher', 'openfda_adverse_getter', 'openfda_adverse_searcher', 'openfda_approval_getter', 'openfda_approval_searcher', 'openfda_device_getter', 'openfda_device_searcher', 'openfda_label_getter', 'openfda_label_searcher', 'openfda_recall_getter', 'openfda_recall_searcher', 'openfda_shortage_getter', 'openfda_shortage_searcher', 'search', 'think', 'trial_getter', 'trial_locations_getter', 'trial_outcomes_getter', 'trial_protocol_getter', 'trial_references_getter', 'trial_searcher', 'variant_getter', 'variant_searcher']
With two sets of tools available, the system prompt has to say which is for what. It also has to mention one practical detail: our MOFA tools return Ensembl IDs with a version suffix, and BioMCP does not recognise those.
llm_with_all_tools = llm.bind_tools(all_tools)
SYSTEM_MULTI = (
"You are a computational-biology assistant. You have two kinds of tools: "
"(1) MOFA tools analysing a fitted multi-omics model of TCGA breast-cancer "
"data (factors 'Factor1'..'Factor10', PAM50 subtypes), and (2) BioMCP tools "
"for general biomedical knowledge (genes, drugs, diseases, literature). "
"Use MOFA tools for questions about our fitted model/factors/patients; use "
"BioMCP tools for general biomedical facts about specific genes, drugs, or "
"diseases. Ground every claim in tool output and name which tool you used. "
"NOTE: MOFA tools return Ensembl gene IDs with a version suffix (e.g. "
"'ENSG00000160180.15'). BioMCP tools do not recognise the version suffix -- "
"strip it (to 'ENSG00000160180') before passing an ID to any BioMCP tool. "
"Resolve one gene at a time rather than issuing many lookups in parallel. "
"If a question asks for something no available tool result actually supports, "
"say so explicitly rather than inferring an answer from general biomedical knowledge."
)
async def run_agent_multi(question: str, max_steps: int = 12, verbose: bool = True) -> AIMessage:
"""Same loop as before, over the combined tool set."""
messages = [SystemMessage(content=SYSTEM_MULTI), HumanMessage(content=question)]
for step in range(max_steps):
model_response = await llm_with_all_tools.ainvoke(messages)
messages.append(model_response)
if not model_response.tool_calls:
return model_response
for call in model_response.tool_calls:
if verbose:
print(f"[step {step}] -> {call['name']}({call['args']})")
result = await all_tools_by_name[call["name"]].ainvoke(call["args"])
if verbose:
# Worth watching: the arguments above are what leaves your machine.
preview = json.dumps(result, default=str)
print(f" <- {preview[:200]}"
+ (f"... [{len(preview)} chars]" if len(preview) > 200 else ""))
messages.append(ToolMessage(content=json.dumps(result, default=str),
tool_call_id=call["id"]))
raise RuntimeError(f"Agent did not finish within {max_steps} steps.")
print("Bound", len(all_tools), "tools to", MODEL)
Bound 44 tools to claude-haiku-4-5
The question below cannot be answered by either server alone. Finding the genes that drive a factor needs our model; saying what those genes do needs the literature.
query = "What is the top negative-weight transcriptomic driver of Factor2? " \
"Search the literature for what is known about its role in breast cancer."
answer = await run_agent_multi(query, max_steps=12)
print(answer.content)
[step 0] -> top_features_for_factor({'factor': 'Factor2', 'n': 50, 'view': 'transcriptomics'})
<- [{"type": "text", "text": "{\n \"factor\": \"Factor2\",\n \"view\": \"transcriptomics\",\n \"top_negative\": {\n \"ENSG00000166535.20\": -0.365,\n \"ENSG00000185686.18\": -0.36,\n \"ENSG00... [3807 chars]
[step 1] -> enrichr_analyzer({'genes': ['ENSG00000166535', 'ENSG00000185686', 'ENSG00000164434', 'ENSG00000186832', 'ENSG00000102243', 'ENSG00000135069', 'ENSG00000198729', 'ENSG00000156219', 'ENSG00000094755', 'ENSG00000280916', 'ENSG00000143546', 'ENSG00000163064', 'ENSG00000136928', 'ENSG00000163220', 'ENSG00000135374', 'ENSG00000205420', 'ENSG00000258910', 'ENSG00000102854', 'ENSG00000261175', 'ENSG00000019169', 'ENSG00000143452', 'ENSG00000229544', 'ENSG00000185479', 'ENSG00000176887', 'ENSG00000114805', 'ENSG00000108176', 'ENSG00000158258', 'ENSG00000286379', 'ENSG00000225362', 'ENSG00000224559', 'ENSG00000124664', 'ENSG00000124935', 'ENSG00000170743', 'ENSG00000259459', 'ENSG00000115648', 'ENSG00000163491', 'ENSG00000255774', 'ENSG00000171428', 'ENSG00000134240', 'ENSG00000074410', 'ENSG00000163879', 'ENSG00000138696', 'ENSG00000136457', 'ENSG00000286208', 'ENSG00000162989', 'ENSG00000147041', 'ENSG00000139865', 'ENSG00000144218', 'ENSG00000183888', 'ENSG00000251141', 'ENSG00000006071'], 'database': 'go_process'})
<- [{"type": "text", "text": "{\"error\": \"Failed to retrieve enrichment results from Enrichr API\", \"genes\": [\"ENSG00000166535\", \"ENSG00000185686\", \"ENSG00000164434\", \"ENSG00000186832\", \"ENS... [1264 chars]
[step 2] -> gene_getter({'gene_id_or_symbol': 'ENSG00000160182'})
<- [{"type": "text", "text": "# Record 1\nError: Gene 'ENSG00000160182' not found\nSuggestion: Please check the gene symbol or ID\n", "id": "lc_a36424a5-9bc2-4115-9e85-0fcf22da59da"}]
[step 2] -> gene_getter({'gene_id_or_symbol': 'ENSG00000173467'})
<- [{"type": "text", "text": "# Record 1\nError: Gene 'ENSG00000173467' not found\nSuggestion: Please check the gene symbol or ID\n", "id": "lc_1216eb50-9d2b-4d19-b810-f212f0a8d2be"}]
[step 2] -> gene_getter({'gene_id_or_symbol': 'ENSG00000160180'})
<- [{"type": "text", "text": "# Record 1\nError: Gene 'ENSG00000160180' not found\nSuggestion: Please check the gene symbol or ID\n", "id": "lc_939b6a55-5a7a-415b-a839-2b9c0248f159"}]
[step 2] -> gene_getter({'gene_id_or_symbol': 'ENSG00000106541'})
<- [{"type": "text", "text": "# Record 1\nError: Gene 'ENSG00000106541' not found\nSuggestion: Please check the gene symbol or ID\n", "id": "lc_44407d89-eb7e-4a16-8e20-6beffbaf133a"}]
[step 2] -> gene_getter({'gene_id_or_symbol': 'ENSG00000235687'})
<- [{"type": "text", "text": "# Record 1\nError: Gene 'ENSG00000235687' not found\nSuggestion: Please check the gene symbol or ID\n", "id": "lc_3ada252e-2001-4901-84e0-256b057d9bcf"}]
[step 3] -> search({'query': 'gene:ENSG00000160182 OR gene:ENSG00000173467 OR gene:ENSG00000160180', 'domain': 'gene'})
<- [{"type": "text", "text": "{\n \"results\": [\n {\n \"id\": 37786776,\n \"title\": \"Establishment of two oxaliplatin-resistant gallbladder cancer cell lines and comprehensive analysis o... [5553 chars]
[step 4] -> think({'thought': "The user wants to know what biological processes are enriched among Factor2's top 50 genes. I've retrieved the top genes from the MOFA model, but I'm having difficulty converting Ensembl IDs to gene symbols for enrichment analysis. The enrichr_analyzer tool requires gene symbols, not Ensembl IDs. Let me try a different strategy - I can search PubMed for articles that mention these Ensembl IDs to identify what genes they are, or I can try to use the fetch tool with the Ensembl IDs directly to get gene information.", 'thoughtNumber': 1, 'totalThoughts': 3, 'nextThoughtNeeded': True})
<- [{"type": "text", "text": "{\n \"domain\": \"thinking\",\n \"result\": \"Added thought 1 to main sequence. Progress: 1/3 thoughts. Next thought needed.\",\n \"thoughtNumber\": 1,\n \"nextThoughtNe... [267 chars]
[step 5] -> fetch({'id': 'ENSG00000160182', 'domain': 'gene'})
<- [{"type": "text", "text": "{\n \"error\": \"Gene ENSG00000160182 not found\"\n}", "id": "lc_a47f185c-3756-4c62-abe1-5d50dce3048a"}]
[step 5] -> fetch({'id': 'ENSG00000173467', 'domain': 'gene'})
<- [{"type": "text", "text": "{\n \"error\": \"Gene ENSG00000173467 not found\"\n}", "id": "lc_e4b6f47b-fb35-49b7-8083-2fcfeda814fe"}]
[step 5] -> fetch({'id': 'ENSG00000160180', 'domain': 'gene'})
<- [{"type": "text", "text": "{\n \"error\": \"Gene ENSG00000160180 not found\"\n}", "id": "lc_31b2d2cf-67de-483d-a97e-bf10a8395c54"}]
[step 5] -> fetch({'id': 'ENSG00000106541', 'domain': 'gene'})
<- [{"type": "text", "text": "{\n \"error\": \"Gene ENSG00000106541 not found\"\n}", "id": "lc_6367abfd-8534-4571-ac90-e3befefd749a"}]
[step 5] -> fetch({'id': 'ENSG00000235687', 'domain': 'gene'})
<- [{"type": "text", "text": "{\n \"error\": \"Gene ENSG00000235687 not found\"\n}", "id": "lc_8fdff01f-5dab-4d96-a809-557c51e1c90f"}]
[step 6] -> think({'thought': "The BioMCP tools are not recognizing these Ensembl IDs. This could be because: 1) The IDs are outdated or from a different genome build, 2) The tools don't support Ensembl ID lookups directly, or 3) These are non-coding RNAs or other features that aren't in the standard gene databases. Let me try a different approach - I'll search PubMed for articles that mention these specific Ensembl IDs to see if I can identify what genes they correspond to.", 'thoughtNumber': 2, 'totalThoughts': 3, 'nextThoughtNeeded': True})
<- [{"type": "text", "text": "{\n \"domain\": \"thinking\",\n \"result\": \"Added thought 2 to main sequence. Progress: 2/3 thoughts. Next thought needed.\",\n \"thoughtNumber\": 2,\n \"nextThoughtNe... [267 chars]
[step 7] -> article_searcher({'keywords': 'ENSG00000160182 ENSG00000173467 ENSG00000160180', 'page_size': 5})
<- [{"type": "text", "text": "[]", "id": "lc_0135a1d0-5a45-4dcd-bfa8-0cc2c39a333e"}]
[step 8] -> think({'thought': "The BioMCP tools are unable to resolve these Ensembl IDs. This is a limitation - the tools don't have a direct mapping from Ensembl IDs to gene symbols. However, I can still provide useful information based on what I know from the MOFA output. The top genes show strong positive and negative loadings on Factor2. Rather than trying to convert all 50 genes, let me acknowledge this limitation to the user and explain what I can determine from the MOFA output itself, or suggest an alternative approach.", 'thoughtNumber': 3, 'totalThoughts': 3, 'nextThoughtNeeded': False})
<- [{"type": "text", "text": "{\n \"domain\": \"thinking\",\n \"result\": \"Added thought 3 to main sequence. Progress: 3/3 thoughts. Thinking sequence complete.\",\n \"thoughtNumber\": 3,\n \"nextTh... [275 chars]
I apologize, but I've encountered a limitation in performing the enrichment analysis you requested. Here's what happened:
**What I retrieved:**
I successfully obtained Factor2's top 50 genes from the MOFA model (25 with the highest positive loadings and 25 with the highest negative loadings). The top positive genes include ENSG00000160182, ENSG00000173467, and ENSG00000160180, while the top negative genes include ENSG00000166535, ENSG00000185686, and ENSG00000164434.
**The problem:**
The BioMCP tools (enrichr_analyzer, gene_getter, and fetch) do not recognize these Ensembl IDs. This appears to be because:
1. The Ensembl IDs from the MOFA model may be from a specific genome build that the BioMCP tools don't have current mappings for
2. The tools require gene symbols (like "TP53", "BRCA1") rather than Ensembl IDs for enrichment analysis
3. Some of these IDs may correspond to non-coding RNAs or other genomic features not in standard gene databases
**What I recommend:**
To complete this analysis, you would need to:
1. Convert the Ensembl IDs to gene symbols using a tool like Ensembl's BioMart or a local annotation file
2. Then provide those gene symbols to me for enrichment analysis
Alternatively, if you have access to the gene symbols for Factor2's top genes, I can perform the biological process enrichment analysis directly.
Watch the trace: it starts in our tools, then crosses into theirs, using the gene IDs the first set returned.
Now write your own. Something that starts in the MOFA model and ends in the public data.
#################################################
# Write a question that has to use both sets of tools.
# Additional clue: get a gene or a set of genes out of the MOFA model first, then ask something
# about them that only public data can answer -- clinical trials, variants, or published literature.
query = "Run BioMCP's gene-set enrichment on Factor2's top 5 genes and report only the " \
"pathways/processes it returns."
answer = await run_agent_multi(query, max_steps=12)
print(answer.content)
#################################################
[step 0] -> top_features_for_factor({'factor': 'Factor2', 'n': 5, 'view': 'transcriptomics'})
<- [{"type": "text", "text": "{\n \"factor\": \"Factor2\",\n \"view\": \"transcriptomics\",\n \"top_negative\": {\n \"ENSG00000166535.20\": -0.365,\n \"ENSG00000185686.18\": -0.36,\n \"ENSG00... [560 chars]
[step 1] -> enrichr_analyzer({'genes': ['ENSG00000160182', 'ENSG00000173467', 'ENSG00000160180', 'ENSG00000106541', 'ENSG00000235687'], 'database': 'pathway'})
<- [{"type": "text", "text": "{\"error\": \"Failed to retrieve enrichment results from Enrichr API\", \"genes\": [\"ENSG00000160182\", \"ENSG00000173467\", \"ENSG00000160180\", \"ENSG00000106541\", \"ENS... [295 chars]
[step 2] -> gene_getter({'gene_id_or_symbol': 'ENSG00000160182'})
<- [{"type": "text", "text": "# Record 1\nError: Gene 'ENSG00000160182' not found\nSuggestion: Please check the gene symbol or ID\n", "id": "lc_b91f3721-22c2-46ba-b6d3-d13357da6ae5"}]
[step 2] -> gene_getter({'gene_id_or_symbol': 'ENSG00000173467'})
<- [{"type": "text", "text": "# Record 1\nError: Gene 'ENSG00000173467' not found\nSuggestion: Please check the gene symbol or ID\n", "id": "lc_3e0a48ec-c122-4bbf-9d0b-602a1af6be5f"}]
[step 2] -> gene_getter({'gene_id_or_symbol': 'ENSG00000160180'})
<- [{"type": "text", "text": "# Record 1\nError: Gene 'ENSG00000160180' not found\nSuggestion: Please check the gene symbol or ID\n", "id": "lc_5dd459cc-8d48-452f-89f4-536b14d87dac"}]
[step 2] -> gene_getter({'gene_id_or_symbol': 'ENSG00000106541'})
<- [{"type": "text", "text": "# Record 1\nError: Gene 'ENSG00000106541' not found\nSuggestion: Please check the gene symbol or ID\n", "id": "lc_edfb74eb-816d-433c-b7d7-8de35a9ac8df"}]
[step 2] -> gene_getter({'gene_id_or_symbol': 'ENSG00000235687'})
<- [{"type": "text", "text": "# Record 1\nError: Gene 'ENSG00000235687' not found\nSuggestion: Please check the gene symbol or ID\n", "id": "lc_47929459-6bef-4028-b667-53d8417a6198"}]
Unfortunately, I'm unable to complete this request. The BioMCP tools have limitations:
1. **MOFA tool output**: The top 5 genes for Factor2 are identified by Ensembl IDs (ENSG00000160182, ENSG00000173467, ENSG00000160180, ENSG00000106541, ENSG00000235687)
2. **Gene symbol conversion failed**: The `gene_getter` tool cannot resolve these Ensembl IDs to gene symbols, which are required by the `enrichr_analyzer` tool for pathway enrichment
3. **Enrichr limitation**: The enrichment analysis tool requires standard gene symbols (like "TP53", "BRCA1") rather than Ensembl IDs
To proceed, you would need either:
- A mapping of these Ensembl IDs to gene symbols from an external source
- Or access to the original MOFA data with gene symbols already annotated
Would you be able to provide the gene symbols for these Ensembl IDs, or do you have access to the MOFA model's gene annotation file?
4. What this bought us#
The first half of the notebook changed nothing about what the agent could do. That was the point: the same tools, reached differently, gave the same answers. What changed is that the tools and the data behind them now live in one place, loaded once, and any number of clients can use them without a copy of the code or the omics tables.
The second half is where we see the usefulness of MCPs. Connecting to BioMCP took the same handful of lines as connecting to our own server, and the agent gained tools over fifteen public databases that we did not write and do not maintain.
What leaves your machine. The omics tables and the fitted model stay inside the MOFA server’s process; BioMCP never sees them. What reaches BioMCP is whatever Claude puts in a tool call’s arguments — in the traces above, gene identifiers — which BioMCP then passes on to the public APIs behind it. Separately, and true of Part 1 as well, the whole conversation goes to Anthropic as part of the agent loop.
For this dataset that is low-stakes: TCGA is public and de-identified. On data that is not, note that the system prompt is a guideline and not a boundary — nothing prevents Claude from putting more context into a tool call than you intended. If you need a guarantee, validate the arguments in code before they reach the server.
Reflection#
Compare the traces in Section 2 with the ones from Part 1. Did Claude call the same tools, in the same order?
The agent loop never mentions MCP. Why did it not need changing when the tools moved onto a server?
When you asked a question needing both servers, how did Claude get from a MOFA result to a BioMCP query?
Our server exposes only read-only tools. What would you want to change if one of them wrote to a file?
5. How the server works#
For the interested reader. Nothing here is needed to run anything above.
Building the server#
MCP has official software development kits in several languages. Ours is written with the Python one, using its high-level interface, FastMCP, which turns ordinary functions into MCP tools. The whole idea fits in a few lines:
from mcp.server.fastmcp import FastMCP
mcp = FastMCP("my-server")
@mcp.tool()
def add(a: int, b: int) -> int:
"""Add two numbers."""
return a + b
if __name__ == "__main__":
mcp.run()
That is all MCP asks of you. Everything else in server/mofa_mcp_server.py is the MOFA analysis from Part 1, unchanged.
The decorator is what does the work. @mcp.tool() takes the function’s name, its type-hinted arguments, and its docstring, and makes those three things visible to any client that connects. They are the same three things Claude saw in Part 1, and the same job LangChain’s @tool did there — one layer further out. As in Part 1, the docstring is the interface: it is what the model reads when deciding whether the tool is worth calling, and the body is what it never sees.
Everything above the decorators in our file — loading the omics tables, the train/test split, opening the fitted model, projecting the test patients — runs once, when the server starts. That is what makes the reuse described earlier practical: one copy of the data, loaded once, shared by every client that connects, rather than a copy loaded inside every notebook.
Two other decorators publish things a plain function list cannot. @mcp.resource("mofa://summary") publishes data addressed by a URI — here a compact description of the fitted model — which the application reads into context rather than the model choosing to call. @mcp.prompt() publishes a reusable template that the user picks; ours is interpret_factor(factor), four lines that ask for a full interpretation of one factor in a fixed form. A prompt supplies instructions rather than capability, which is an idea Part 3 takes much further.
Each tool also carries ToolAnnotations(readOnlyHint=True, ...), saying it only reads and never changes anything. A client can use that to run a tool without pausing to ask a human first — worth having precisely because the caller here is a model rather than a fixed script.
Finally, mcp.run() at the bottom is what lets MultiServerMCPClient start the server with python server/mofa_mcp_server.py. Run that command yourself in a terminal and it will sit there apparently doing nothing, waiting for input — which is what a server of this kind looks like when nothing is talking to it.
The exchange underneath#
Now that you have seen how a server publishes things, here is what the client and the server actually say to each other.
Below we connect with MCP’s own client library instead of the adapter. The pattern is always the same: open a session, initialize it — the server replies with what it can do — and then ask for things. No model is involved, so this costs nothing to run.
from mcp import ClientSession, StdioServerParameters
from mcp.client.stdio import stdio_client
server_params = StdioServerParameters(command=sys.executable, args=[str(SERVER_PATH)])
async with stdio_client(server_params) as (read, write):
async with ClientSession(read, write) as session:
init = await session.initialize()
print("connected to:", init.serverInfo.name, "\n")
# TOOLS -- the model may choose to call these
for t in (await session.list_tools()).tools:
args = ", ".join(t.inputSchema.get("properties", {}))
print(f" tool {t.name}({args})")
# RESOURCES -- data to read into context, not called by the model
for r in (await session.list_resources()).resources:
print(f" resource {r.uri}")
# PROMPTS -- reusable templates a user can pick
for p in (await session.list_prompts()).prompts:
args = ", ".join(a.name for a in (p.arguments or []))
print(f" prompt {p.name}({args})")
print("\nmofa://summary ->")
summary = await session.read_resource("mofa://summary")
print(" ", summary.contents[0].text)
print("\ncalling factor_subtype_association ->")
called = await session.call_tool("factor_subtype_association", {})
print(" ", called.content[0].text[:150], "...")
connected to: eccb2026-mofa
tool data_summary()
tool split_summary()
tool active_factors()
tool factor_view_r2(factor)
tool factor_subtype_association()
tool top_features_for_factor(factor, view, n)
tool classify_subtype_from_factors()
tool train_vs_test_subtype_association()
resource mofa://summary
prompt interpret_factor(factor)
mofa://summary ->
{
"n_patients": 500,
"views": [
"transcriptomics",
"proteomics",
"methylation"
],
"n_factors": 10,
"active_factors": [
"Factor1",
"Factor2",
"Factor3",
"Factor4",
"Factor5",
"Factor6",
"Factor7",
"Factor8",
"Factor9",
"Factor10"
],
"subtypes": {
"LumA": 237,
"LumB": 100,
"Basal": 97,
"Her2": 41,
"Normal": 25
}
}
calling factor_subtype_association ->
{
"factor": "Factor2",
"eta_squared": 0.824
} ...
Two things that makes visible.
A server offers more than tools. Tools are chosen by the model. Resources are data the application reads into context, like mofa://summary. Prompts are templates the user picks, like interpret_factor. In Section 1 we called get_tools(), which pulled in only the first of the three, because that is what a tool-calling loop needs; get_resources() and get_prompt() fetch the others.
The messages are ordinary. They travel as JSON-RPC — a plain request/response format — over a transport. Ours is stdio: the server runs as a subprocess and the two talk through its standard input and output, the same way you would pipe one command-line program into another. Over a network it would be HTTP instead, and nothing above would change.