Inside the Omics Enclave
A simulated multi-omics project inside a trusted research environment: what the analyst sees, which tools run on each data type, how the omics files connect to phenopackets and OMOP clinical data, and what a finding looks like on screen.
A project you can walk through end to end
When I onboard a new analyst into an enclave, the first question is never about statistics. It is "what is in the folder, and what can I run on it?" This tool answers that question for a project that holds genomic, transcriptomic, proteomic, and metabolomic data next to clinical phenotypes.
Every participant, sample, file, value, and result on these pages is invented for teaching. The file layouts, the tools, the commands, and the ways the clinical and omics layers connect are the real ones. Where a tool version or a threshold is shown, treat it as a snapshot from September 2026 and check the current release before you rely on it.
The sample project
The enclave holds a rare epilepsy cohort recruited through a registry. Participants consented to research use of their clinical records and biosamples. For each participant the project holds a GA4GH phenopacket (a structured description of the person's phenotypes, disease, and samples, using the Human Phenotype Ontology), an extract of their longitudinal records in the OMOP Common Data Model (the OHDSI standard for observational health data), and up to four omics layers:
- Whole genome sequencing (WGS) from blood, as a parent-child trio where both parents enrolled
- RNA sequencing (RNA-seq) from cultured skin fibroblasts
- Plasma proteomics by data-independent acquisition (DIA) mass spectrometry
- Plasma metabolomics by liquid chromatography mass spectrometry (LC-MS)
The research question that drives the walkthrough is a common one in rare disease work: can the omics layers explain why a subgroup of participants shows developmental regression, and does any layer support a candidate genetic cause in participants who are still undiagnosed?
How to read each stage
Each stage shows a simulated window with tabs. Click the tabs to move between the file browser, the terminal, the notebook, and the results the analyst would see at that point. Below each window, a shaded block explains how the omics step connects back to the phenopackets or the OMOP tables, and an amber block records what the analyst found. Gene names in the findings use the placeholder GENEX1 so that no real gene community is implied.
Use the jump bar at the top to move between stages, or the previous and next buttons at the bottom of each page.
Entering the enclave
The analyst signs in through a browser or a remote desktop client, reaches a launcher, and finds a Linux workspace with no route to the internet.
A trusted research environment (TRE), also called a data enclave, is a computing environment where approved analysts work on sensitive data that cannot leave. Everything the analyst needs has to be inside already: the pipelines, the container images, the reference genomes and databases, and the R and Python packages. Anything that is not there has to be requested through an import process, and anything that needs to leave goes through an export review.
The launcher page the analyst sees after signing in. Each tile opens inside the enclave; nothing opens a connection outward.
analyst@enclave-login:~$ whoami dboyce_r analyst@enclave-login:~$ ls /project/EPI-OMICS-2026 README.md MANIFEST.tsv data_dictionary clinical genomics transcriptomics proteomics metabolomics references work outputs analyst@enclave-login:~$ curl -I https://nf-co.re curl: (6) Could not resolve host: nf-co.re analyst@enclave-login:~$ module avail 2>&1 | head -20 --------------------- /opt/modules --------------------- nextflow/24.10 apptainer/1.3 R/4.4 python/3.12 samtools/1.20 bcftools/1.20 bwa-mem2/2.2 gatk/4.6 deepvariant/1.6 vep/113 exomiser/14 star/2.7 salmon/1.10 diann/1.9 openms/3.2 duckdb/1.1 analyst@enclave-login:~$ ls /project/EPI-OMICS-2026/references/nf-core rnaseq sarek quantms metaboigniter # offline bundles from nf-core download analyst@enclave-login:~$ sinfo PARTITION AVAIL TIMELIMIT NODES STATE NODELIST compute* up 7-00:00:00 8 idle c[01-08] highmem up 7-00:00:00 2 idle hm[01-02] gpu up 2-00:00:00 1 idle g01
What the analyst will not find, and what to do instead
No internet Pipelines, containers, and references are mirrored under
references/. New items go through the airlock as an import.
No Docker Containers run through Apptainer (rootless). Nextflow
profiles in the enclave already point at it.
No pip install A local package index mirrors approved packages. Anything
no install.packages() else is an import request with a justification.
No AI assistant Code is written and read by the analyst. The enclave logs
every session.
No copy and paste out Screenshots are blocked. Results leave only through the
airlock after output review.
No write access to Only work/ and outputs/ are writable. Raw data are
raw data read-only, which also protects the analyst from mistakes.
The tools the analyst will use across the project
Nextflow is the workflow engine. It runs each pipeline as a graph of processes, submits them to the Slurm scheduler, and records every input, output, and container in a log that becomes the provenance record. The nf-core pipelines are community-maintained Nextflow workflows for common assays; the enclave stages them offline with nf-core download so that the analyst runs them with the -offline flag and a pinned release. JupyterLab and RStudio Server are where the analyst reads pipeline output, joins it to the clinical layer, and fits models. DuckDB reads the OMOP Parquet extracts directly with SQL and needs no database server.
Where the clinical layer sits
The phenopackets and the OMOP tables live in clinical/ next to the omics folders, not in a separate system. A single linkage table maps participant identifiers to sample identifiers. That table is the join key for the whole project, and every stage that follows uses it.
The project folder
Before running anything, the analyst reads the README, the manifest, and the data dictionary. The folder layout tells you which assays exist, which participants have which samples, and where the pipelines will write.
/project/EPI-OMICS-2026 ├── README.md ├── MANIFEST.tsv one row per file: participant, sample, assay, path, checksum, batch ├── data_dictionary/ │ ├── manifest_columns.md │ ├── omop_tables_in_scope.md │ └── phenopacket_conventions.md ├── clinical/ │ ├── phenopackets/ one JSON per participant, GA4GH Phenopacket schema v2 │ │ ├── P-0001.json ... P-0060.json │ ├── omop/ Parquet extracts, OMOP CDM v5.4, tables in scope only │ │ ├── person.parquet │ │ ├── condition_occurrence.parquet │ │ ├── drug_exposure.parquet │ │ ├── measurement.parquet │ │ ├── observation.parquet │ │ ├── procedure_occurrence.parquet │ │ └── visit_occurrence.parquet │ └── linkage/ │ ├── participant_sample_map.tsv │ └── pedigrees.ped trio structure for de novo analysis ├── genomics/ │ ├── fastq/ paired-end WGS, ~30x, gzip, ~25 GB per file │ ├── samplesheet_sarek.csv │ └── README_genomics.md ├── transcriptomics/ │ ├── fastq/ fibroblast RNA-seq, paired-end, stranded, ~3 GB per file │ ├── samplesheet_rnaseq.csv │ └── README_transcriptomics.md ├── proteomics/ │ ├── mzml/ plasma DIA runs converted from vendor raw, ~2 GB per file │ ├── sdrf/experiment.sdrf.tsv sample and data relationship format for proteomics │ ├── spectral_library/ │ └── README_proteomics.md ├── metabolomics/ │ ├── mzml/ plasma LC-MS, positive and negative mode, QC pools, blanks │ ├── sample_metadata.tsv injection order, batch, QC flag │ └── README_metabolomics.md ├── references/ mirrored, read-only │ ├── GRCh38/ genome FASTA, index, known sites │ ├── gencode_v46/ gene annotation GTF, transcript FASTA │ ├── vep_cache/ Ensembl VEP cache and plugins │ ├── exomiser_data/ phenotype and variant data release │ ├── uniprot_human.fasta │ ├── hmdb/ metabolite reference, MS2 spectra │ ├── vocab/ OMOP vocabulary tables as Parquet │ └── nf-core/ offline pipeline bundles ├── work/ writable scratch, pipeline work directories, notebooks └── outputs/ the only folder eligible for export review
The manifest is the analyst's map from participant to file. Filtered here to one participant.
participant_id sample_id tissue assay file_path md5 batch
P-0007 S-0007-DNA blood WGS genomics/fastq/S-0007-DNA_L001_R1.fastq.gz 9f3a1c… G2
P-0007 S-0007-DNA blood WGS genomics/fastq/S-0007-DNA_L001_R2.fastq.gz 0b77e4… G2
P-0007F S-0007F-DNA blood WGS genomics/fastq/S-0007F-DNA_L001_R1.fastq.gz c21d90… G2
P-0007F S-0007F-DNA blood WGS genomics/fastq/S-0007F-DNA_L001_R2.fastq.gz 4e0a52… G2
P-0007M S-0007M-DNA blood WGS genomics/fastq/S-0007M-DNA_L001_R1.fastq.gz 77b3f1… G2
P-0007M S-0007M-DNA blood WGS genomics/fastq/S-0007M-DNA_L001_R2.fastq.gz a9c8d3… G2
P-0007 S-0007-RNA fibroblast RNA-seq transcriptomics/fastq/S-0007-RNA_R1.fastq.gz e13b6f… T1
P-0007 S-0007-RNA fibroblast RNA-seq transcriptomics/fastq/S-0007-RNA_R2.fastq.gz 58d2a0… T1
P-0007 S-0007-PL plasma DIA proteomics proteomics/mzml/S-0007-PL.mzML b4f9c7… PR1
P-0007 S-0007-PL plasma LC-MS metabolomics metabolomics/mzml/S-0007-PL_pos.mzML 2c6e81… M1
P-0007 S-0007-PL plasma LC-MS metabolomics metabolomics/mzml/S-0007-PL_neg.mzML f0a4d9… M1
# EPI-OMICS-2026 Rare epilepsy cohort with paired clinical and multi-omics data. Protocol number, approval, and data use agreement: see project wiki. ## What is here - clinical/phenopackets GA4GH Phenopacket v2, exported from the registry on 2026-06-12 - clinical/omop OMOP CDM v5.4 extract, cut date 2026-06-30, tables in scope only - genomics WGS trios where both parents enrolled, singletons otherwise - transcriptomics fibroblast RNA-seq, all participants with a skin biopsy - proteomics plasma DIA, all participants with a plasma sample - metabolomics plasma LC-MS, same plasma aliquots as proteomics ## Identifiers participant_id (P-xxxx) is the project pseudonym. It appears as phenopacket subject.id and as OMOP person.person_source_value. sample_id (S-xxxx-TYPE) appears as phenopacket biosamples[].id and in every omics file name. Parents are P-xxxxF and P-xxxxM. ## Rules - Raw data folders are read-only. - Write to work/ while analyzing, then copy export candidates to outputs/. - Individual-level records never leave. See the output review policy on the wiki. - Pin pipeline releases and record them in outputs/PROVENANCE.md.
clinical/linkage/participant_sample_map.tsv and pedigrees.ped. The pedigree is in PLINK PED format, which Exomiser and GATK read directly.
participant_id family_id role sample_dna sample_rna sample_plasma omop_person_id P-0007 F-0007 proband S-0007-DNA S-0007-RNA S-0007-PL 100341 P-0007F F-0007 father S-0007F-DNA . . . P-0007M F-0007 mother S-0007M-DNA . . . P-0008 F-0008 proband S-0008-DNA S-0008-RNA S-0008-PL 100342 P-0009 F-0009 proband S-0009-DNA . S-0009-PL 100343 # pedigrees.ped family individual father mother sex affected F-0007 P-0007 P-0007F P-0007M 2 2 F-0007 P-0007F 0 0 1 1 F-0007 P-0007M 0 0 2 1
What the analyst learns from the folder alone
- Which participants have which layers. Not everyone has fibroblast RNA-seq, so any analysis that joins layers starts by counting the overlap from the linkage table rather than assuming it.
- Which participants form trios. De novo variant analysis is only possible where both parents were sequenced.
- Where batch effects will come from. Every manifest row records a batch, and the metabolomics metadata records injection order, because both will go into the models later.
- Where results will be written. Pipelines write to
work/; only reviewed aggregates are copied tooutputs/.
The join key for the whole project
The linkage table is the one place where the participant pseudonym, the OMOP person_id, and every sample identifier sit on one row. Phenopackets are keyed by participant. Omics files are keyed by sample. OMOP tables are keyed by an integer person_id, with the participant pseudonym stored in person_source_value. Every later join goes through this table.
The clinical layer: phenopackets and OMOP
The omics files describe molecules. The clinical layer describes people. The two are joined through the linkage table, and the analyst reads both before defining a single comparison group.
A phenopacket is a JSON document, one per participant, that records phenotypic features as Human Phenotype Ontology (HPO) terms with onset and exclusion flags, the disease diagnosis as a MONDO term, the biosamples taken, and any prior variant interpretations. Some omics tools, Exomiser among them, read a phenopacket directly as input. The OMOP extract holds the longitudinal record: visits, conditions, drug exposures, and lab measurements, coded to standard vocabularies. It answers questions the phenopacket does not, such as which participants were on a given medication when the plasma sample was drawn.
clinical/phenopackets/P-0007.json, abridged. The excluded microcephaly feature is negative evidence, which some tools use.
{
"id": "P-0007",
"subject": {
"id": "P-0007",
"sex": "FEMALE",
"timeAtLastEncounter": { "age": { "iso8601duration": "P4Y2M" } }
},
"phenotypicFeatures": [
{ "type": { "id": "HP:0011097", "label": "Epileptic spasm" },
"onset": { "age": { "iso8601duration": "P5M" } } },
{ "type": { "id": "HP:0001263", "label": "Global developmental delay" } },
{ "type": { "id": "HP:0002376", "label": "Developmental regression" },
"onset": { "age": { "iso8601duration": "P1Y" } } },
{ "type": { "id": "HP:0001250", "label": "Seizure" } },
{ "type": { "id": "HP:0000252", "label": "Microcephaly" }, "excluded": true }
],
"diseases": [
{ "term": { "id": "MONDO:0100062", "label": "developmental and epileptic encephalopathy" } }
],
"biosamples": [
{ "id": "S-0007-DNA", "sampledTissue": { "id": "UBERON:0000178", "label": "blood" } },
{ "id": "S-0007-RNA", "sampledTissue": { "id": "UBERON:0002097", "label": "skin of body" },
"description": "dermal fibroblast culture, passage 4" },
{ "id": "S-0007-PL", "sampledTissue": { "id": "UBERON:0001969", "label": "blood plasma" } }
],
"interpretations": [], ← undiagnosed: no prior molecular finding recorded
"metaData": {
"created": "2026-06-12T09:41:00Z",
"createdBy": "registry-export v3.2",
"resources": [ { "id": "hp", "version": "2026-04-01" }, { "id": "mondo", … }, { "id": "uberon", … } ],
"phenopacketSchemaVersion": "2.0"
}
}
DuckDB reads the Parquet extracts and the mirrored vocabulary without a server. The question here is medication exposure at the time of the plasma draw.
import duckdb
con = duckdb.connect()
q = """
SELECT p.person_source_value AS participant_id,
c.concept_name AS drug,
d.drug_exposure_start_date AS start_date,
d.drug_exposure_end_date AS end_date
FROM read_parquet('clinical/omop/drug_exposure.parquet') d
JOIN read_parquet('clinical/omop/person.parquet') p USING (person_id)
JOIN read_parquet('references/vocab/concept.parquet') c ON c.concept_id = d.drug_concept_id
WHERE p.person_source_value = 'P-0007'
ORDER BY start_date
"""
con.sql(q).df()participant_id drug start_date end_date 0 P-0007 vigabatrin 2022-09-14 2023-03-02 1 P-0007 levetiracetam 2023-01-20 None 2 P-0007 valproic acid 2023-11-08 None
# the plasma sample date comes from the phenopacket biosample, or from the manifest con.sql(""" SELECT p.person_source_value AS participant_id, m.measurement_date, c.concept_name AS lab, m.value_as_number, m.unit_source_value FROM read_parquet('clinical/omop/measurement.parquet') m JOIN read_parquet('clinical/omop/person.parquet') p USING (person_id) JOIN read_parquet('references/vocab/concept.parquet') c ON c.concept_id = m.measurement_concept_id WHERE p.person_source_value = 'P-0007' AND c.concept_name ILIKE '%neurofilament%' """).df()
participant_id measurement_date lab value_as_number unit_source_value 0 P-0007 2026-02-11 Neurofilament light chain [Plasma] 41.7 pg/mL
The comparison group for the whole project is defined once, from the phenopackets, and saved as a table every later stage reads.
import json, glob, pandas as pd rows = [] for f in glob.glob("clinical/phenopackets/P-*.json"): pp = json.load(open(f)) feats = {x["type"]["id"] for x in pp.get("phenotypicFeatures", []) if not x.get("excluded")} rows.append({ "participant_id": pp["subject"]["id"], "sex": pp["subject"]["sex"], "regression": "HP:0002376" in feats, "spasms": "HP:0011097" in feats, "diagnosed": len(pp.get("interpretations", [])) > 0, }) groups = pd.DataFrame(rows).merge( pd.read_csv("clinical/linkage/participant_sample_map.tsv", sep="\t"), on="participant_id") groups.to_csv("work/groups.tsv", sep="\t", index=False) groups.groupby(["regression", "diagnosed"]).size()
regression diagnosed
False False 14
True 21
True False 9
True 16
dtype: int64What each source contributes to the analysis
The phenopacket gives the analysis its phenotype groups (developmental regression present or absent, diagnosed or undiagnosed) and gives Exomiser the HPO terms it needs to rank variants. The OMOP tables give covariates and confounders: medication at sample time, lab values to compare against omics measurements, age at each visit. The analyst writes work/groups.tsv once and every pipeline stage reads it, so a change in how a group is defined happens in one place.
What is inside the OMOP extract
The extract holds only the tables named in the data use agreement, as Parquet files with the standard OMOP CDM v5.4 columns. person holds year of birth, sex, and the participant pseudonym in person_source_value. condition_occurrence holds diagnoses with a standard concept id (usually SNOMED) and start and end dates. drug_exposure holds medication records with an RxNorm concept id, dates, and quantity where known. measurement holds labs with a LOINC concept id, numeric value, and unit. observation holds facts that are neither condition nor measurement, such as developmental milestones recorded at a visit. visit_occurrence holds encounter dates and types. Dates are shifted per participant by a fixed offset inside the enclave, so intervals are preserved and calendar dates are not real.
Genomics: from reads to a ranked candidate
The genomics use case is diagnostic: find a plausible causal variant in an undiagnosed proband by using the trio structure and the phenopacket together.
The pipeline is nf-core/sarek, which takes paired FASTQ reads, aligns them to GRCh38 with BWA-MEM2, marks duplicates, calls variants with DeepVariant and GATK HaplotypeCaller, and joint-calls the trio. Ensembl VEP annotates the calls with gene, consequence, and population frequency. Exomiser then ranks variants by combining variant pathogenicity with how well each gene's known phenotypes match the HPO terms in the phenopacket.
genomics/samplesheet_sarek.csv. One row per lane per sample; the pipeline merges lanes. The status column is 0 for germline.
patient,sex,status,sample,lane,fastq_1,fastq_2 P-0007,XX,0,S-0007-DNA,L001,genomics/fastq/S-0007-DNA_L001_R1.fastq.gz,genomics/fastq/S-0007-DNA_L001_R2.fastq.gz P-0007,XX,0,S-0007-DNA,L002,genomics/fastq/S-0007-DNA_L002_R1.fastq.gz,genomics/fastq/S-0007-DNA_L002_R2.fastq.gz P-0007F,XY,0,S-0007F-DNA,L001,genomics/fastq/S-0007F-DNA_L001_R1.fastq.gz,genomics/fastq/S-0007F-DNA_L001_R2.fastq.gz P-0007M,XX,0,S-0007M-DNA,L001,genomics/fastq/S-0007M-DNA_L001_R1.fastq.gz,genomics/fastq/S-0007M-DNA_L001_R2.fastq.gz # a FASTQ record, four lines per read; a 30x genome is roughly 400 million of these pairs @A01234:88:HXYZ:1:1101:2381:1000 1:N:0:ATCACG GATCTGGCTAGCAACTTAGGTACCTAGGATCTGACCATTAGCATGGATCGATCGATCGGATC... + FFFFFFF:FFFFFFFFFFFFFF:FFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFF...
analyst@enclave-login:~/work$ nextflow run /project/EPI-OMICS-2026/references/nf-core/sarek \ -r <pinned release> -offline -resume \ -profile apptainer,slurm \ --input /project/EPI-OMICS-2026/genomics/samplesheet_sarek.csv \ --genome GATK.GRCh38 --igenomes_base /project/EPI-OMICS-2026/references/igenomes \ --tools deepvariant,haplotypecaller --joint_germline \ --outdir /project/EPI-OMICS-2026/work/sarek N E X T F L O W ~ version 24.10 Launching `references/nf-core/sarek/main.nf` [F-0007-trio] DSL2 executor > slurm (41) [c1/8a2f0e] NFCORE_SAREK:SAREK:FASTQC (S-0007-DNA-L001) [100%] 4 of 4 ✔ [7b/d90c11] NFCORE_SAREK:SAREK:FASTQ_ALIGN_BWAMEM2 (S-0007-DNA) [100%] 4 of 4 ✔ [e4/33b7a9] NFCORE_SAREK:SAREK:GATK4_MARKDUPLICATES (S-0007F-DNA) [100%] 3 of 3 ✔ [12/fc41d8] NFCORE_SAREK:SAREK:BAM_BASERECALIBRATOR (S-0007M-DNA) [100%] 3 of 3 ✔ [9a/0e5b27] NFCORE_SAREK:SAREK:SAMTOOLS_CONVERT_CRAM (S-0007-DNA) [100%] 3 of 3 ✔ [b8/71a3c2] NFCORE_SAREK:SAREK:DEEPVARIANT (S-0007-DNA) [100%] 3 of 3 ✔ [5d/2c9e40] NFCORE_SAREK:SAREK:GATK4_HAPLOTYPECALLER (S-0007-DNA) [100%] 3 of 3 ✔ [03/ab61f7] NFCORE_SAREK:SAREK:GATK4_GENOTYPEGVCFS (F-0007) [100%] 1 of 1 ✔ [cc/4d0e19] NFCORE_SAREK:SAREK:GLNEXUS (F-0007) [100%] 1 of 1 ✔ [6e/9f27b3] NFCORE_SAREK:SAREK:MULTIQC [100%] 1 of 1 ✔ Completed at: 2026-08-19 03:12:44 Duration: 11h 38m CPU hours: 1,206 Succeeded: 41 analyst@enclave-login:~/work$ ls sarek/variant_calling/deepvariant/F-0007 F-0007.deepvariant.joint.vcf.gz F-0007.deepvariant.joint.vcf.gz.tbi analyst@enclave-login:~/work$ vep --offline --cache --dir_cache ../references/vep_cache \ --assembly GRCh38 --everything --vcf \ -i sarek/variant_calling/deepvariant/F-0007/F-0007.deepvariant.joint.vcf.gz \ -o F-0007.vep.vcf.gz --compress_output bgzip
MultiQC gathers per-sample metrics from every step. The analyst checks these before believing any variant.
Sample Reads (M) Mapped % Dup % Mean cov ≥20x % Insert size Contam % Sex check
S-0007-DNA 412.6 99.4 6.1 31.8 95.2 380 0.11 XX ✔
S-0007F-DNA 398.1 99.5 5.8 30.2 94.6 372 0.09 XY ✔
S-0007M-DNA 421.9 99.3 6.4 32.5 95.6 388 0.14 XX ✔
Trio consistency (Mendelian error rate): 0.31 % ✔ within expected range
Relatedness (KING kinship, proband vs father): 0.248 ✔ parent-child
Relatedness (KING kinship, proband vs mother): 0.251 ✔ parent-child
A joint Variant Call Format (VCF) file: a header describing every field, then one row per variant site with one genotype column per sample. Annotation adds the CSQ field.
##fileformat=VCFv4.2 ##reference=GRCh38 ##INFO=<ID=CSQ,Number=.,Type=String,Description="Consequence annotations from Ensembl VEP. Format: Allele|Consequence|IMPACT|SYMBOL|Gene|Feature|HGVSc|HGVSp|gnomADg_AF|CLIN_SIG|..."> ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> ##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for ref and alt"> ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read depth"> ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype quality"> #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT S-0007-DNA S-0007F-DNA S-0007M-DNA chr9 136492817 . C G 58.2 PASS CSQ=G|splice_region_variant&intron_variant|LOW|GENEX1|ENSG…|ENST…|c.1327-3C>G||.|.|... GT:AD:DP:GQ 0/1:16,15:31:99 0/0:29,0:29:87 0/0:33,0:33:93 chr2 165310508 . T C 61.0 PASS CSQ=C|missense_variant|MODERATE|SCN2A|ENSG…|ENST…|c.2941T>C|p.Ser981Pro|0.00012|.|... GT:AD:DP:GQ 0/1:14,17:31:99 0/1:15,16:31:99 0/0:30,0:30:90 # Reading the first row: the proband is heterozygous (0/1) with 15 alternate reads out of 31; # both parents are homozygous reference (0/0) with good depth, so the call is a candidate de novo. # gnomAD frequency is absent (.), meaning the variant was not seen in the population reference. # The second row is inherited from the father and is common enough to be an unlikely cause on its own. analyst@enclave-login:~/work$ bcftools view -H F-0007.vep.vcf.gz | wc -l 4,912,306 # variant sites in the trio; the filtering below reduces this to a handful
analyst@enclave-login:~/work$ exomiser --analysis exomiser_trio.yml \ --sample ../clinical/phenopackets/P-0007.json \ --vcf F-0007.vep.vcf.gz \ --ped ../clinical/linkage/pedigrees.ped \ --output-directory exomiser/P-0007 # exomiser_trio.yml (abridged): frequency ≤ 0.1 %, pathogenicity sources REVEL and MVP, # inheritance modes AD, AR, XD, XR, and de novo via the pedigree; HPO terms come from the phenopacket P-0007.genes.tsv RANK GENE COMBINED PHENO VARIANT MOI HGVS GT(proband,father,mother) gnomAD_AF ClinVar 1 GENEX1 0.9821 0.9412 0.9163 AD c.1327-3C>G (splice reg.) 0/1, 0/0, 0/0 (de novo) . not reported 2 KCNQ2 0.6104 0.7720 0.5011 AD c.1670A>G p.Asn557Ser 0/1, 0/1, 0/0 0.00031 VUS 3 SCN2A 0.4488 0.7115 0.3388 AD c.2941T>C p.Ser981Pro 0/1, 0/1, 0/0 0.00012 not reported 4 CDKL5 0.2917 0.6660 0.2004 XD c.2288C>T p.Ala763Val 0/1, 0/0, 0/1 0.00087 benign … # PHENO is the similarity between the phenopacket HPO terms and each gene's known phenotype profile. # GENEX1 reaches the top through a strong phenotype match and de novo status, even though a splice # region variant has a LOW predicted impact by itself. That is the reason the next stage looks at RNA.
How the phenopacket and pedigree drive the ranking
Exomiser reads the phenopacket file directly, so the HPO terms entered in the registry become the phenotype score without any retyping. The pedigree file, built from the linkage table, tells Exomiser which genotype patterns count as de novo. The excluded microcephaly feature is stored in the phenopacket; whether a given Exomiser release uses excluded features as negative evidence depends on the version, so the analyst checks the release notes rather than assuming.
What the analyst found
A de novo splice region variant in GENEX1 is the top-ranked candidate for P-0007, driven by a strong phenotype match. The variant sits three bases into an intron, where predicted impact is low, so the genomic evidence alone is not enough. The analyst flags it for the transcriptomics stage, where an aberrant splicing test can show whether the variant changes the transcript. Any return of a finding to the participant is a clinical and governance decision outside this workflow; the analyst records the candidate in the project log and continues.
What is inside the genomics output folder
work/sarek/ holds a preprocessing/ folder with one CRAM (compressed alignment) file and index per sample, at roughly 18 GB each; a variant_calling/ folder with per-sample gVCF files and the joint trio VCF; a reports/ folder with FastQC, duplicate metrics, coverage, and the MultiQC HTML; and a pipeline_info/ folder with the execution report, timeline, trace, and the exact software versions, which go into the provenance record at export. Each CRAM can be opened in IGV inside the enclave to look at the reads under a candidate variant.
Transcriptomics: expression outliers and a group comparison
Fibroblast RNA-seq serves two use cases at once. For the undiagnosed proband, it tests whether the candidate variant disrupts splicing. For the cohort, it compares expression between the regression and no-regression groups.
The pipeline is nf-core/rnaseq, run with STAR for alignment and Salmon for transcript quantification. Its main product is a gene-by-sample counts matrix. From there the analyst runs DESeq2 in R for the group comparison, OUTRIDER for single-sample expression outliers, and FRASER for aberrant splicing. The last two are designed for rare disease diagnostics, where the question is "which gene is abnormal in this one person" rather than "which genes differ between groups."
# transcriptomics/samplesheet_rnaseq.csv sample,fastq_1,fastq_2,strandedness S-0007-RNA,transcriptomics/fastq/S-0007-RNA_R1.fastq.gz,transcriptomics/fastq/S-0007-RNA_R2.fastq.gz,reverse S-0008-RNA,transcriptomics/fastq/S-0008-RNA_R1.fastq.gz,transcriptomics/fastq/S-0008-RNA_R2.fastq.gz,reverse … analyst@enclave-login:~/work$ nextflow run /project/EPI-OMICS-2026/references/nf-core/rnaseq \ -r <pinned release> -offline -resume -profile apptainer,slurm \ --input ../transcriptomics/samplesheet_rnaseq.csv \ --fasta ../references/GRCh38/GRCh38.primary_assembly.genome.fa \ --gtf ../references/gencode_v46/gencode.v46.annotation.gtf \ --aligner star_salmon --save_unaligned false \ --outdir rnaseq executor > slurm (188) [a1/9c02d4] NFCORE_RNASEQ:RNASEQ:FASTQ_FASTQC_UMITOOLS_TRIMGALORE (S-0007-RNA) [100%] 44 of 44 ✔ [3f/71be08] NFCORE_RNASEQ:RNASEQ:ALIGN_STAR (S-0007-RNA) [100%] 44 of 44 ✔ [c8/44d1a0] NFCORE_RNASEQ:RNASEQ:QUANTIFY_STAR_SALMON:SALMON_QUANT (S-0007-RNA) [100%] 44 of 44 ✔ [0d/e6f2b7] NFCORE_RNASEQ:RNASEQ:QUANTIFY_STAR_SALMON:SALMON_TX2GENE [100%] 1 of 1 ✔ [77/0a9c3e] NFCORE_RNASEQ:RNASEQ:DESEQ2_QC_STAR_SALMON [100%] 1 of 1 ✔ [b2/5e8d10] NFCORE_RNASEQ:RNASEQ:MULTIQC [100%] 1 of 1 ✔ Completed at: 2026-08-21 14:02:10 Duration: 6h 51m Succeeded: 188 analyst@enclave-login:~/work$ ls rnaseq/star_salmon | head S-0007-RNA/ # per-sample: BAM, salmon quant.sf, bigWig S-0008-RNA/ salmon.merged.gene_counts.tsv # the matrix the R work starts from salmon.merged.gene_tpm.tsv salmon.merged.transcript_counts.tsv tx2gene.tsv
rnaseq/star_salmon/salmon.merged.gene_counts.tsv. One row per gene, one column per sample. Around 60,000 rows including non-coding genes; a typical filter keeps the 15,000 to 20,000 that are expressed in fibroblasts.
gene_id gene_name S-0007-RNA S-0008-RNA S-0011-RNA S-0012-RNA S-0015-RNA … ENSG00000000003.16 TSPAN6 1842.0 1911.5 1766.2 2004.8 1698.3 ENSG00000000419.14 DPM1 3210.4 3402.0 2988.7 3311.2 3120.9 ENSG00000000457.14 SCYL3 688.0 702.1 655.9 741.3 690.4 ENSG00000187608.10 ISG15 412.7 2988.3 398.1 3105.6 455.0 ← high in some samples ENSG00000137959.17 IFI44L 88.2 1220.5 91.0 1301.7 97.4 ENSG00000157601.14 MX1 501.3 4120.9 488.6 4009.2 512.8 ENSG00000XXXXXX.5 GENEX1 1150.2 1204.6 1187.0 1166.3 1173.9 ← placeholder gene … # A per-sample salmon quant.sf, the transcript-level file the matrix is built from Name Length EffectiveLength TPM NumReads ENST00000373020.9 3768 3591.42 12.8043 1701.000 ENST00000494424.1 820 643.55 0.0000 0.000
# groups.tsv came from the phenopackets in the clinical stage; batch from the manifest library(DESeq2); library(readr); library(dplyr) counts <- read_tsv("rnaseq/star_salmon/salmon.merged.gene_counts.tsv") groups <- read_tsv("groups.tsv") %>% filter(!is.na(sample_rna)) manifest <- read_tsv("../MANIFEST.tsv") %>% filter(assay == "RNA-seq") %>% distinct(sample_id, batch) coldata <- groups %>% inner_join(manifest, by = c(sample_rna = "sample_id")) %>% mutate(regression = factor(regression, levels = c(FALSE, TRUE)), batch = factor(batch)) mat <- as.matrix(counts[, coldata$sample_rna]); rownames(mat) <- counts$gene_name dds <- DESeqDataSetFromMatrix(round(mat), coldata, ~ batch + sex + regression) dds <- dds[rowSums(counts(dds) >= 10) >= 8, ] dds <- DESeq(dds) res <- results(dds, name = "regression_TRUE_vs_FALSE") %>% as.data.frame() %>% arrange(padj) %>% head(8) res
baseMean log2FoldChange lfcSE stat pvalue padj
ISG15 1408.2 2.71 0.31 8.74 2.3e-18 3.9e-14
IFI44L 611.5 3.02 0.36 8.39 4.8e-17 4.1e-13
MX1 2140.7 2.44 0.30 8.13 4.4e-16 2.5e-12
IFIT1 902.3 2.58 0.33 7.82 5.3e-15 2.2e-11
OAS1 744.9 2.11 0.29 7.28 3.4e-13 1.2e-9
RSAD2 388.4 2.36 0.34 6.94 3.9e-12 1.1e-8
IFI6 1655.0 1.87 0.28 6.68 2.4e-11 5.8e-8
STAT1 3011.8 1.32 0.21 6.29 3.2e-10 6.8e-7
## Reading this: the regression group shows higher expression of interferon-stimulated genes
## in fibroblasts. The next question is whether medication explains it, which is an OMOP join.library(OUTRIDER); library(FRASER) # OUTRIDER: is any gene an expression outlier in a single sample, after removing hidden confounders? ods <- OutriderDataSet(countData = round(mat)) ods <- filterExpression(ods, minCounts = TRUE) ods <- OUTRIDER(ods) results(ods) %>% filter(sampleID == "S-0007-RNA")
geneID sampleID pValue padjust zScore rawcounts normcounts
1 GENEX1 S-0007-RNA 3.1e-07 0.0041 -4.38 1150 2201 ← about half the expected count
2 COL1A2 S-0007-RNA 8.8e-05 0.31 -3.12 …# FRASER: aberrant splicing, computed from the STAR BAMs. Look at the intron next to the candidate variant. fds <- FraserDataSet(colData = coldata, bamFile = bams, workingDir = "fraser") fds <- countRNAData(fds); fds <- calculatePSIValues(fds); fds <- FRASER(fds) results(fds, sampleIDs = "S-0007-RNA", padjCutoff = 0.05)
seqnames start end sampleID hgncSymbol type pValue padjust deltaPsi 1 chr9 136492521 136492819 S-0007-RNA GENEX1 psi5 2.2e-09 1.9e-05 -0.63 ← the acceptor next to c.1327-3C>G 2 chr9 136492521 136492903 S-0007-RNA GENEX1 psi5 4.0e-08 1.7e-04 +0.58 ← a new acceptor 84 bases downstream ## Reading this: in P-0007 the normal splice acceptor is used far less (deltaPsi -0.63) and a cryptic ## acceptor is used instead, which shifts the reading frame. Reduced expression in OUTRIDER is consistent ## with the abnormal transcript being degraded.
How the clinical layer entered this stage
The DESeq2 design used the regression flag from the phenopackets and the batch from the manifest. The OUTRIDER and FRASER results were filtered to the sample belonging to P-0007, found through the linkage table. The interferon signal now raises a question only OMOP can answer: is it explained by a medication or an infection recorded near the biopsy date? That check happens in the integration stage.
What the analyst found
FRASER shows that the candidate GENEX1 splice region variant from the genomics stage changes splicing in the proband's fibroblasts, with use of a cryptic acceptor and loss of the normal one, and OUTRIDER shows reduced expression of the same gene in the same sample. Two independent layers now point at the same candidate, which is the strongest form of support this project can produce short of functional work. Separately, the regression group shows an interferon-stimulated gene signature that will be tested against covariates and against the proteomics.
Proteomics: plasma proteins between groups
Plasma proteomics asks whether the regression group differs in circulating proteins, and whether any protein-level signal matches the transcript-level signal from fibroblasts or a lab value in the clinical record.
The instrument files arrive already converted from the vendor format to mzML, an open XML format that holds every mass spectrum from a run. The experimental design is described in an SDRF-Proteomics file (sample and data relationship format), which nf-core/quantms reads to know which files belong to which sample and condition. The pipeline runs DIA-NN against a spectral library and produces a protein-by-sample intensity matrix. The analyst then fits linear models with limma in R.
proteomics/sdrf/experiment.sdrf.tsv. One row per raw file. The pipeline reads the factor value column as the condition.
source name characteristics[organism] characteristics[organism part] comment[data file] comment[instrument] comment[label] factor value[phenotype] comment[technical replicate] S-0007-PL Homo sapiens blood plasma S-0007-PL.mzML NT=timsTOF HT;AC=MS:… label free regression 1 S-0008-PL Homo sapiens blood plasma S-0008-PL.mzML NT=timsTOF HT;AC=MS:… label free no regression 1 S-0009-PL Homo sapiens blood plasma S-0009-PL.mzML NT=timsTOF HT;AC=MS:… label free regression 1 … # The factor value column was filled from groups.tsv, so the proteomics design and the RNA design agree.
analyst@enclave-login:~/work$ nextflow run /project/EPI-OMICS-2026/references/nf-core/quantms \ -r <pinned release> -offline -resume -profile apptainer,slurm \ --input ../proteomics/sdrf/experiment.sdrf.tsv \ --database ../references/uniprot_human.fasta \ --acquisition_method dia \ --diann_speclib ../proteomics/spectral_library/plasma_dia.predicted.speclib \ --outdir quantms executor > slurm (66) [4a/1d0e93] NFCORE_QUANTMS:QUANTMS:FILE_PREPARATION:MZMLINDEXING (S-0007-PL) [100%] 58 of 58 ✔ [91/b2c7f0] NFCORE_QUANTMS:QUANTMS:DIA:DIANN_PRELIMINARY_ANALYSIS (S-0007-PL) [100%] 58 of 58 ✔ [e7/30a5d2] NFCORE_QUANTMS:QUANTMS:DIA:ASSEMBLE_EMPIRICAL_LIBRARY [100%] 1 of 1 ✔ [0c/f8e114] NFCORE_QUANTMS:QUANTMS:DIA:DIANN_INDIVIDUAL_FINAL_ANALYSIS [100%] 58 of 58 ✔ [b5/6a92cd] NFCORE_QUANTMS:QUANTMS:DIA:DIANNSUMMARY [100%] 1 of 1 ✔ [3d/c04e77] NFCORE_QUANTMS:QUANTMS:DIA:MSSTATS [100%] 1 of 1 ✔ Completed at: 2026-08-24 22:47:03 Duration: 9h 04m Succeeded: 66 analyst@enclave-login:~/work$ ls quantms/diannsummary diann_report.tsv # one row per precursor per run, the primary output diann_report.pg_matrix.tsv # protein group by sample intensities, MaxLFQ diann_report.pr_matrix.tsv # precursor by sample diann_report.stats.tsv # per-run identification counts and QC
The analyst rarely opens an mzML by hand, but knowing its shape explains why the files are large and why processing takes hours. Around 60,000 spectra per plasma run.
<mzML xmlns="http://psi.hupo.org/ms/mzml" version="1.1.0">
<fileDescription> … source file, instrument, conversion software … </fileDescription>
<run id="S-0007-PL" startTimeStamp="2026-05-14T10:22:41Z">
<spectrumList count="61208">
<spectrum index="0" id="scan=1" defaultArrayLength="4127">
<cvParam accession="MS:1000511" name="ms level" value="1"/>
<scanList><scan><cvParam accession="MS:1000016" name="scan start time" value="0.0041" unitName="minute"/></scan></scanList>
<binaryDataArrayList>
<binaryDataArray> <cvParam name="m/z array"/> <binary>…base64…</binary> </binaryDataArray>
<binaryDataArray> <cvParam name="intensity array"/> <binary>…base64…</binary> </binaryDataArray>
</binaryDataArrayList>
</spectrum>
<spectrum index="1" id="scan=2">
<cvParam accession="MS:1000511" name="ms level" value="2"/>
<precursorList><precursor><isolationWindow> m/z 400.0 to 412.0 </isolationWindow></precursor></precursorList>
…
</spectrum>
…
</spectrumList>
</run>
</mzML>
# In DIA, each MS2 spectrum covers a wide isolation window and mixes fragments from many peptides.
# DIA-NN untangles that mixture by matching against a spectral library, which is why the library
# has to be staged inside the enclave alongside the pipeline.
quantms/diannsummary/diann_report.pg_matrix.tsv. One row per protein group, one column per run, MaxLFQ intensities. Plasma DIA at this depth reports roughly 1,500 to 2,500 protein groups.
Protein.Group Protein.Names Genes First.Protein.Description S-0007-PL S-0008-PL S-0009-PL … P02768 ALBU_HUMAN ALB Albumin 4.21e10 4.35e10 4.08e10 P01023 A2MG_HUMAN A2M Alpha-2-macroglobulin 9.87e9 1.02e10 9.55e9 P07196 NFL_HUMAN NEFL Neurofilament light polypeptide 2.14e5 6.8e4 1.98e5 ← low abundance; near the detection floor P14136 GFAP_HUMAN GFAP Glial fibrillary acidic protein 1.61e5 5.2e4 1.55e5 P05161 ISG15_HUMAN ISG15 Ubiquitin-like protein ISG15 8.9e5 2.3e5 9.4e5 …
library(limma); library(readr); library(dplyr) pg <- read_tsv("quantms/diannsummary/diann_report.pg_matrix.tsv") mat <- log2(as.matrix(pg[, grep("^S-", names(pg))])); rownames(mat) <- pg$Genes mat <- mat[rowMeans(!is.na(mat)) >= 0.7, ] # keep proteins seen in most runs design_df <- groups %>% filter(!is.na(sample_plasma)) %>% arrange(match(sample_plasma, colnames(mat))) design <- model.matrix(~ batch + sex + regression, design_df) fit <- eBayes(lmFit(mat, design)) topTable(fit, coef = "regressionTRUE", n = 6)
logFC AveExpr t P.Value adj.P.Val
NEFL 1.42 17.1 6.02 2.1e-07 3.4e-04
GFAP 1.28 16.9 5.41 1.6e-06 1.3e-03
ISG15 1.63 19.6 5.10 4.9e-06 2.6e-03
CHI3L1 0.94 21.3 4.22 9.8e-05 3.9e-02
B2M 0.51 27.0 3.88 2.8e-04 8.9e-02
C4A -0.37 26.4 -3.41 1.2e-03 3.2e-01
## Reading this: NEFL and GFAP, both markers of neuronal and glial injury, are higher in the regression
## group. ISG15 is higher too, which echoes the interferon-stimulated transcripts in fibroblasts.# cross-check against the clinical lab: OMOP measurement holds a clinical NfL assay for a subset nfl_clin <- con %>% tbl_omop_measurement("Neurofilament light chain [Plasma]") # helper around the DuckDB query from stage 3 cor.test(mat["NEFL", nfl_clin$sample_plasma], log2(nfl_clin$value_as_number))
Pearson's product-moment correlation: r = 0.71, n = 23, p = 1.4e-04
Where the clinical record met the proteome
The SDRF condition column was written from groups.tsv, so proteomics and transcriptomics compared the same groups. The OMOP measurement table then provided an independent clinical assay of one of the top proteins, and the correlation between the research measurement and the clinical one gives the analyst confidence that the low-abundance DIA quantification of NEFL is real rather than noise near the detection floor.
What the analyst found
Plasma NEFL and GFAP are higher in participants with developmental regression, and the research NEFL values agree with clinical NfL assays where both exist. ISG15 is higher at the protein level in the same group that showed interferon-stimulated transcripts in fibroblasts. The pattern is consistent across two tissues and two assays, which makes it a candidate for the integration stage rather than a finding on its own.
Metabolomics: features, annotation, and a confounder
Untargeted plasma metabolomics produces thousands of features, most of them unannotated. The use case is to find metabolic differences between groups and then decide, using the clinical record, which differences are biology and which are treatment.
The instrument files are mzML, the same open format as the proteomics runs, acquired in positive and negative ionization modes. Processing turns raw spectra into a feature table: each feature is a mass-to-charge ratio (m/z) and a retention time with an intensity per sample. The enclave stages nf-core/metaboigniter for this, which wraps peak picking and alignment; many analysts run the same steps directly in R with XCMS. Annotation matches feature m/z and fragmentation spectra against a mirrored copy of the Human Metabolome Database (HMDB). Pooled quality control (QC) injections spread across the run let the analyst correct drift and remove unstable features.
metabolomics/sample_metadata.tsv. Injection order and batch are recorded because signal drifts over a run.
file sample_id sample_type mode batch injection_order acquisition_date BLANK_01_pos.mzML . blank pos M1 1 2026-05-20 QC_pool_01_pos.mzML . qc_pool pos M1 2 2026-05-20 S-0031-PL_pos.mzML S-0031-PL study pos M1 3 2026-05-20 S-0007-PL_pos.mzML S-0007-PL study pos M1 4 2026-05-20 S-0044-PL_pos.mzML S-0044-PL study pos M1 5 2026-05-20 QC_pool_02_pos.mzML . qc_pool pos M1 6 2026-05-20 … # Study samples are injected in randomized order, with a pooled QC every few injections and blanks at the start and end.
analyst@enclave-login:~/work$ nextflow run /project/EPI-OMICS-2026/references/nf-core/metaboigniter \ -r <pinned release> -offline -resume -profile apptainer,slurm \ --input ../metabolomics/sample_metadata.tsv \ --polarity both --ms2_data true \ --outdir metabo executor > slurm (24) [8e/12ab30] NFCORE_METABOIGNITER:METABOIGNITER:CENTROIDING (S-0007-PL_pos) [100%] 132 of 132 ✔ [f2/7c9d4e] NFCORE_METABOIGNITER:METABOIGNITER:PEAK_PICKING (batch M1, pos) [100%] 4 of 4 ✔ [a9/0b3e61] NFCORE_METABOIGNITER:METABOIGNITER:ALIGNMENT_GROUPING [100%] 2 of 2 ✔ [d4/e5f7a8] NFCORE_METABOIGNITER:METABOIGNITER:MS2_MATCHING (HMDB local) [100%] 2 of 2 ✔ [1b/9a4c02] NFCORE_METABOIGNITER:METABOIGNITER:FEATURE_TABLE [100%] 2 of 2 ✔ Completed at: 2026-08-26 08:19:55 Duration: 2h 12m Succeeded: 24 # The equivalent in R, which many analysts prefer for control over parameters: library(xcms) raw <- readMSData(files, mode = "onDisk") peaks <- findChromPeaks(raw, CentWaveParam(ppm = 10, peakwidth = c(5, 30))) peaks <- adjustRtime(peaks, ObiwarpParam()) %>% groupChromPeaks(PeakDensityParam(sampleGroups = md$batch)) feat <- featureValues(fillChromPeaks(peaks), value = "into")
metabo/feature_table_pos.tsv. One row per feature, one column per injection. Several thousand rows in each mode; most stay unannotated.
feature_id mz rt_sec npeaks annotation adduct hmdb_id ms2_match S-0007-PL S-0008-PL S-0009-PL QC_pool_01 … F0001 118.0863 41.2 58 betaine [M+H]+ HMDB0000043 0.92 3.2e7 3.4e7 3.1e7 3.3e7 F0002 162.1125 52.8 58 L-carnitine [M+H]+ HMDB0000062 0.95 8.1e6 2.9e7 7.7e6 2.1e7 ← low in some samples F0003 204.1230 77.4 58 acetyl-L-carnitine [M+H]+ HMDB0000201 0.91 1.4e6 6.2e6 1.3e6 4.4e6 F0004 166.0863 95.1 57 phenylalanine [M+H]+ HMDB0000159 0.96 1.8e7 1.9e7 1.7e7 1.8e7 F0005 205.0972 101.6 58 tryptophan [M+H]+ HMDB0000929 0.94 9.6e6 1.0e7 9.4e6 9.8e6 F0006 144.1019 63.0 44 . . . . 4.2e5 3.9e5 4.4e5 4.1e5 ← unannotated, kept for later F0007 146.1176 258.3 58 . . . . 7.7e5 7.9e5 7.6e5 7.8e5 … # mz is the measured mass-to-charge ratio; rt_sec is the retention time in the chromatography; # npeaks is how many injections the feature was detected in; ms2_match is the fragmentation spectrum similarity to the library.
# keep features that are stable in the pooled QC injections, then correct drift against injection order qc_cv <- apply(feat[, md$sample_type == "qc_pool"], 1, function(x) sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE)) keep <- qc_cv < 0.30 & rowMeans(is.na(feat[, md$sample_type == "study"])) < 0.2 feat <- feat[keep, ] feat_corr <- drift_correct(feat, order = md$injection_order, qc = md$sample_type == "qc_pool") # LOESS on QC, per feature c(before = nrow(feat_raw), after = nrow(feat))
before after
6,412 3,188
## Blank subtraction removed features present in blanks at more than a third of the sample intensity.
## The PCA after correction shows QC pools clustered tightly, which is the check that the correction worked.library(limma) lm_mat <- log2(feat_corr[, study_cols]); colnames(lm_mat) <- md$sample_id[md$sample_type == "study"] design <- model.matrix(~ batch + sex + regression, groups_pl) fit <- eBayes(lmFit(lm_mat, design)) topTable(fit, coef = "regressionTRUE", n = 5) %>% left_join(annot, by = "feature_id")
feature_id annotation logFC adj.P.Val 1 F0002 L-carnitine -1.61 2.2e-05 2 F0003 acetyl-L-carnitine -1.48 6.9e-05 3 F0118 octanoylcarnitine -1.12 8.1e-04 4 F0006 (unannotated) 0.88 3.0e-03 5 F0211 3-hydroxybutyrate 0.74 2.7e-02
# Before calling this biology: which participants were on valproate at the plasma draw? That is an OMOP question. vpa <- omop_exposed_on_date(drug = "valproic acid", date_col = "plasma_draw_date") # DuckDB query on drug_exposure table(groups_pl$regression, vpa$on_valproate)
on_valproate
regression FALSE TRUE
FALSE 26 7
TRUE 6 19
## Valproate is far more common in the regression group. Refit with valproate as a covariate:design2 <- model.matrix(~ batch + sex + on_valproate + regression, cbind(groups_pl, vpa)) fit2 <- eBayes(lmFit(lm_mat, design2)) topTable(fit2, coef = "regressionTRUE", n = 3) %>% left_join(annot, by = "feature_id")
feature_id annotation logFC adj.P.Val
1 F0006 (unannotated) 0.81 1.1e-02
2 F0211 3-hydroxybutyrate 0.69 4.4e-02
3 F0002 L-carnitine -0.21 0.71 ← the carnitine difference was the medicationWhy this stage needed the OMOP tables
The carnitine signal was strong, biologically plausible, and wrong as a group difference. Valproate lowers carnitine, and valproate exposure was unevenly distributed across the phenotype groups. Only the drug_exposure table, joined through the linkage table and restricted to the plasma draw date, could show that. Without the clinical layer inside the enclave, this would have been reported as a finding.
What the analyst found
The largest metabolomic difference between groups was a medication effect and was removed by adjusting for valproate exposure. What remains after adjustment is smaller: a modest increase in 3-hydroxybutyrate and one unannotated feature. The analyst records the unannotated feature's m/z and retention time so that it can be pursued later with an authentic standard, and treats the ketone body signal as exploratory pending a check against dietary records in observation.
Putting the layers together
Each layer produced a table keyed by sample. Integration joins them by participant, fits a model that looks across layers at once, and asks whether the shared structure lines up with the phenotype.
The analyst uses MOFA+ (multi-omics factor analysis), which finds a small number of latent factors that explain variation across several data matrices, each factor with loadings on the features of every layer. A factor that explains variance in both the transcriptome and the proteome, and separates the regression group, is the kind of result that no single layer can produce. Pathway enrichment then names what the loaded features have in common, using gene sets mirrored inside the enclave.
# every layer becomes a features-by-participant matrix; the linkage table maps sample ids to participant ids link <- read_tsv("../clinical/linkage/participant_sample_map.tsv") rna <- vst_mat; colnames(rna) <- link$participant_id[match(colnames(rna), link$sample_rna)] prot <- prot_mat; colnames(prot) <- link$participant_id[match(colnames(prot), link$sample_plasma)] metab <- metab_mat; colnames(metab) <- link$participant_id[match(colnames(metab), link$sample_plasma)] common <- Reduce(intersect, list(colnames(rna), colnames(prot), colnames(metab))) length(common)
[1] 44
## Participants with all three of RNA, proteome, and metabolome. MOFA+ tolerates missing layers,
## so the model below uses everyone and lets the method handle the gaps.library(MOFA2) views <- list(rna = rna[top_var(rna, 3000), ], proteome = prot, metabolome = metab[annotated_rows, ]) mobj <- create_mofa(views) opts <- get_default_model_options(mobj); opts$num_factors <- 8 mobj <- prepare_mofa(mobj, model_options = opts) mobj <- run_mofa(mobj, outfile = "mofa/model.hdf5", use_basilisk = TRUE) # Python env pre-staged in the enclave covs <- groups %>% left_join(vpa) %>% select(participant_id, regression, sex, on_valproate, batch) samples_metadata(mobj) <- covs %>% rename(sample = participant_id) plot_variance_explained(mobj)
Variance explained (%) by factor and view rna proteome metabolome Factor1 18.2 3.1 0.9 ← fibroblast culture and batch Factor2 11.4 9.7 1.2 ← shared across rna and proteome Factor3 1.0 2.4 14.8 ← metabolome only; tracks valproate Factor4 6.3 1.1 0.4 …
correlate_factors_with_covariates(mobj, covariates = c("regression", "sex", "on_valproate", "batch")) plot_top_weights(mobj, view = "rna", factor = 2, nfeatures = 8) plot_top_weights(mobj, view = "proteome", factor = 2, nfeatures = 6)
Factor–covariate association (-log10 p) regression sex on_valproate batch Factor1 0.4 0.2 0.3 6.1 Factor2 5.8 0.1 0.9 0.5 Factor3 1.2 0.3 7.4 0.4 Factor 2, top rna weights: ISG15 IFI44L MX1 IFIT1 OAS1 RSAD2 IFI6 STAT1 Factor 2, top proteome weights: ISG15 NEFL GFAP CHI3L1 B2M CXCL10 ## Reading this: one factor explains variance in both fibroblast transcripts and plasma proteins, ## is associated with the regression phenotype, and is not associated with valproate, sex, or batch. ## Its loadings are interferon-stimulated genes in the RNA and neural injury markers in the proteome.
library(clusterProfiler) reactome <- read.gmt("../references/genesets/reactome_human_2026-03.gmt") # mirrored; no online lookup w <- get_weights(mobj, views = "rna", factors = 2)$rna[, 1] enr <- GSEA(sort(w, decreasing = TRUE), TERM2GENE = reactome, pvalueCutoff = 0.05) head(as.data.frame(enr)[, c("ID", "NES", "p.adjust")], 5)
ID NES p.adjust 1 Interferon alpha/beta signaling 2.61 1.8e-09 2 Interferon Signaling 2.44 4.2e-08 3 Antiviral mechanism by IFN-stimulated genes 2.19 9.7e-06 4 Cytokine Signaling in Immune system 1.83 3.1e-04 5 ISG15 antiviral mechanism 2.02 3.1e-04
Where the clinical layer closed the loop
The regression flag from the phenopackets is the covariate the factor is tested against. The valproate flag from OMOP is what shows Factor 3 to be a medication axis rather than biology, and what clears Factor 2 of the same suspicion. Without both, the analyst would have a factor and no way to say what it means.
What the analyst found
A shared transcript and protein axis, characterized by interferon signaling and neural injury markers, separates participants with developmental regression from those without, independent of medication, sex, and batch. It does not say what causes what: the interferon signature could be a consequence of injury as easily as a cause. What the analyst has is a hypothesis with support from two tissues and two assay types, a defined participant group, and a reproducible pipeline, which is what a grant or a follow-up protocol needs. Separately, the GENEX1 candidate in P-0007 has genomic and transcriptomic support and is recorded for clinical follow-up through the project's governance route.
Leaving the enclave with results
Nothing leaves without review. The analyst copies export candidates to outputs/, writes the provenance file, and submits an airlock request that an output checker reviews against the disclosure policy.
The output review policy in this project allows aggregate tables where every cell meets the enclave's minimum count, summary statistics, figures without individual points that could be traced to a person, model coefficients, pipeline logs, and code. It does not allow individual-level rows, raw variant lines, sample-level intensity matrices, or anything that names a participant pseudonym. The threshold and the exact rules are the enclave's, not the analyst's; the project wiki holds them, and the checker applies them.
/project/EPI-OMICS-2026/outputs/ ├── PROVENANCE.md ├── tables/ │ ├── T1_cohort_by_group.tsv counts by regression × diagnosed × sex, all cells above threshold │ ├── T2_rnaseq_deseq2_regression.tsv gene, log2FC, SE, padj; no per-sample values │ ├── T3_proteome_limma_regression.tsv │ ├── T4_metabolome_limma_adjusted.tsv annotated features only, valproate-adjusted model │ ├── T5_mofa_variance_explained.tsv │ └── T6_pathway_enrichment_factor2.tsv ├── figures/ │ ├── F1_volcano_rnaseq.pdf gene labels only, no sample identifiers │ ├── F2_mofa_factor2_by_group.pdf boxplot with jitter removed; group medians and IQR only │ └── F3_pathway_dotplot.pdf ├── code/ │ ├── notebooks/*.ipynb, *.Rmd outputs cleared │ ├── nextflow_commands.sh │ └── exomiser_trio.yml └── logs/ ├── sarek_pipeline_info/ rnaseq_pipeline_info/ quantms_pipeline_info/ └── software_versions.yml
# PROVENANCE.md
Project: EPI-OMICS-2026
Clinical extract: phenopackets 2026-06-12 (registry-export v3.2); OMOP CDM v5.4 cut 2026-06-30
Reference data: GRCh38 primary assembly; GENCODE v46; Ensembl VEP cache 113; Exomiser data 2506;
UniProt human 2026_02; HMDB mirror 2026-01; Reactome GMT 2026-03
Pipelines (release pinned, run with -offline, Apptainer, Slurm):
nf-core/sarek <release> run 2026-08-18, Nextflow 24.10, log: logs/sarek_pipeline_info
nf-core/rnaseq <release> run 2026-08-21
nf-core/quantms <release> run 2026-08-24
nf-core/metaboigniter <release> run 2026-08-26
Analysis code: code/notebooks (git commit 3f9a2c1 on the enclave's internal repository)
Group definition: work/groups.tsv built from phenopackets, HP:0002376 present and not excluded
Covariates: batch (manifest), sex (phenopacket), valproate exposure at draw date (OMOP drug_exposure)
Software: logs/software_versions.yml (every container and package version)
Reproduce: code/nextflow_commands.sh, then knit code/notebooks in numbered order.
Checker notes on the first submission
Refused figures/F2 (first version) Jittered points allowed a participant with an extreme value
to be identified by anyone who knows the cohort. Resubmitted
with medians and interquartile ranges only.
Refused tables/T1 (first version) One cell (regression × undiagnosed × male) fell below the
minimum count. Resubmitted with sex collapsed for that stratum.
Refused exomiser/P-0007.genes.tsv Individual-level variant data. Not eligible for export under
the data use agreement; recorded in the project log instead.
Approved everything else, second submission, with the provenance file attached.
What follows the analyst out
The exported package holds aggregate results, code, and provenance. What stays behind is everything keyed to a person: the phenopackets, the OMOP extract, the raw and processed omics files, and the per-participant results. Someone reading the manuscript can reproduce the analysis by applying for access to the same enclave and running the same pinned pipelines, which is what the provenance file makes possible.
Tools by use case
A reference table for the whole walkthrough: what the analyst wanted to know, which tool answered it, and what went in and came out.
| Use case | Tool | Input | Output |
|---|---|---|---|
| Run any pipeline reproducibly inside the enclave | Nextflow with nf-core pipelines, Apptainer containers, Slurm | Samplesheet or SDRF, mirrored references, pinned release | Results folder, execution log, software versions |
| Query the clinical record | DuckDB on OMOP Parquet extracts | OMOP CDM tables, mirrored vocabulary | Covariates, exposures at a date, lab values |
| Define phenotype groups | Python or R over phenopacket JSON | HPO terms with onset and exclusion flags | groups.tsv read by every stage |
| Align reads and call germline variants | nf-core/sarek (BWA-MEM2, GATK, DeepVariant, GLnexus) | Paired FASTQ, GRCh38 | CRAM per sample, joint VCF, quality report |
| Annotate variants | Ensembl VEP, offline cache | VCF | VCF with consequence, gene, frequency, ClinVar fields |
| Rank candidate variants by phenotype | Exomiser | Annotated VCF, phenopacket, pedigree | Ranked genes and variants with phenotype and variant scores |
| Look at reads under a variant | IGV inside the enclave | CRAM, GRCh38 | Visual confirmation of the call |
| Quantify gene expression | nf-core/rnaseq (STAR, Salmon, MultiQC) | Paired FASTQ, GENCODE GTF | Gene and transcript count matrices, BAMs |
| Compare expression between groups | DESeq2 (R) | Counts matrix, group and batch design | Fold changes and adjusted p values per gene |
| Find an expression outlier in one participant | OUTRIDER (R) | Counts matrix | Per-sample outlier genes with z scores |
| Find aberrant splicing in one participant | FRASER (R) | STAR BAMs | Per-sample aberrant splice junctions |
| Identify and quantify plasma proteins | nf-core/quantms with DIA-NN | mzML, SDRF, spectral library, UniProt FASTA | Protein group intensity matrix, per-run QC |
| Compare proteins or metabolites between groups | limma (R) | Log intensity matrix, design with covariates | Fold changes and adjusted p values per feature |
| Turn LC-MS spectra into a feature table | nf-core/metaboigniter or XCMS (R) | mzML, sample metadata with injection order | Feature table (m/z, retention time, intensity per sample) |
| Annotate metabolite features | HMDB mirror, MS2 matching | Feature m/z and fragmentation spectra | Putative identities with match scores |
| Correct signal drift and filter unstable features | Pooled QC injections, LOESS in R | Feature table, injection order, QC flags | Corrected, filtered feature table |
| Find structure shared across layers | MOFA+ (R) | Feature matrices per layer keyed by participant | Factors, loadings, variance explained per layer |
| Name the biology behind a signal | clusterProfiler with mirrored gene sets | Ranked genes or loadings | Enriched pathways |
| Take results out | Airlock with output review | outputs/ folder, provenance file, disclosure statement | Approved aggregate package |
Glossary and sources
Terms used on these pages
- Airlock
- The controlled route for moving files into or out of an enclave, with review on the way out.
- CRAM
- A compressed format for aligned sequencing reads, smaller than BAM, that needs the reference genome to read.
- DIA (data-independent acquisition)
- A mass spectrometry mode that fragments everything in wide m/z windows rather than picking individual precursors, so every run measures the same set of peptides.
- FASTQ
- The raw sequencing read format: a name line, the bases, a separator, and a quality string, repeated per read.
- GA4GH
- Global Alliance for Genomics and Health, the body that publishes the phenopacket schema among other standards.
- HMDB
- Human Metabolome Database, the reference the metabolite annotation step matches against.
- HPC
- High-performance computing: the cluster of compute nodes, managed here by the Slurm scheduler, that runs the pipelines.
- HPO
- Human Phenotype Ontology, the controlled vocabulary of phenotypic features used in phenopackets and by Exomiser.
- IRB (Institutional Review Board)
- The committee that reviews and approves research involving human participants. The project's approval and data use agreement set what the enclave may hold and what may leave it.
- LC-MS
- Liquid chromatography mass spectrometry, the platform used here for untargeted metabolomics.
- Linkage table
- The single table mapping participant pseudonym to every sample identifier and to the OMOP person identifier.
- m/z
- Mass-to-charge ratio, the horizontal axis of a mass spectrum.
- MOFA+
- Multi-omics factor analysis, a method that finds latent factors shared across several data matrices.
- mzML
- The open XML format for mass spectrometry data, used for both proteomics and metabolomics runs.
- OMOP CDM
- The Observational Medical Outcomes Partnership Common Data Model, maintained by OHDSI, a standard table structure and vocabulary for observational health data.
- Phenopacket
- A GA4GH standard JSON document describing one person's phenotypes, disease, samples, and interpretations.
- QC pool
- A mixture of all study samples injected repeatedly through a metabolomics run to track and correct drift.
- SDRF-Proteomics
- Sample and data relationship format, the tab-separated file that describes a proteomics experiment's design.
- TRE
- Trusted research environment, also called a data enclave or secure data environment.
- VCF
- Variant Call Format, the text format for genetic variants with one row per site and one column per sample.
- WGS
- Whole genome sequencing.
Sources
The tools and standards shown are documented at the sites below. Versions, parameters, and thresholds on these pages are illustrative; check each project's current documentation before running.
- Nextflow documentation
- nf-core pipelines, including sarek, rnaseq, quantms, and metaboigniter; offline use is described under running offline
- GA4GH Phenopacket schema
- Human Phenotype Ontology
- OMOP Common Data Model and Athena vocabulary browser
- DuckDB documentation
- Exomiser documentation
- Ensembl Variant Effect Predictor
- DESeq2, OUTRIDER, FRASER, limma, xcms, MOFA2, and clusterProfiler on Bioconductor
- DIA-NN
- SDRF-Proteomics specification
- mzML format (HUPO-PSI)
- Human Metabolome Database
- Integrative Genomics Viewer
This tool is part of the Inside the Enclave series at Boyce Data Science. The companion guide, Working Inside a Trusted Research Environment, covers the access, governance, and output review side in more depth.