Loading and merging¶
Open a cohort¶
from c4.cohort import cohort
demo = cohort("demo")
Cohort names are case-insensitive, and aliases work too (cohort("DEMO")).
Opening a cohort prints everything the catalog knows about it:
Cohort: demo (DEMO) owner: A. Researcher
Catalog: c4/catalog.yaml root: /data/lab
Loadable with .load():
name id params status
──────────────────────── ───────────────── ───────────────────────────────────────────── ──────
id_map - ✓
pheno.main * patient (ID) ✓
pcs array (IID) ✓
hla.hlahd.standardized * exome (seqID) ✓
hla.hibag.raw * array (sample.id) panel: european*, african, asian, multiethnic 3/4 ✗
Genetic data, use .path():
name format build params status
─────────── ────── ───── ───────────────── ───────
dna.imputed vcf b38 chrom: 1..22 (22) 22/22 ✓
dna.exome - - not set
* default: "pheno" -> pheno.main, "hla.hlahd" -> hla.hlahd.standardized, "hla.hibag" -> hla.hibag.raw
How to read it:
- name is what you pass to
.load()or.path(). A*marks a default:"pheno"on its own loadspheno.main(see the footnote). - id is the kind of ID the table uses and the column it is in. This is how the table gets linked to patients (How IDs link).
- params are the values a templated file takes, with the default marked
*. Long lists are shortened (1..22 (22)). - status:
✓all files found ·✗ missing·3/4 ✗some of a templated file's variants are missing ·not setthe catalog hasTODOfor this path.
The table prints again whenever you display the object (demo on its own line
in a notebook, or print(demo)). Pass show=False to cohort() to skip it.
Load a table¶
pheno = demo.load("pheno")
Left out 1 patient(s) marked drop=1 in the id_map (keep_dropped=True to keep them).
patient_id qc_note ID age female site OS_mo OS_event
0 P001 NaN P001 64 0 A 31.2 1
1 P002 NaN P002 58 1 A 44.0 0
2 P003 old exome_id: EX-0003 P003 71 1 B 12.5 1
3 P004 old exome_id: EX-0004 P004 66 0 B 27.9 1
4 P005 NaN P005 49 1 C 50.1 0
5 P007 NaN P007 73 0 A 19.0 1
6 P008 NaN P008 55 0 B 38.6 0
Every loaded table comes back as a pandas DataFrame with patient_id added at
the front, looked up from the table's own ID column through the cohort's ID
map. The original ID column (ID here) is kept. If the ID map has a
qc_note column, that is added too.
The demo has 8 patients but only 7 rows came back: P006 is marked drop=1 in
the ID map, so it is left out of everything you load. See
Dropped patients and QC notes.
Tables keyed by a genetic ID are linked the same way. The PC file uses array IDs:
pcs = demo.load("pcs")
Note: 1 of 7 rows in 'pcs' have no patient_id (array ID not in id_map).
Left out 1 patient(s) marked drop=1 in the id_map (keep_dropped=True to keep them).
patient_id qc_note IID PC1 PC2 PC3
0 P001 NaN ARR-0001 0.012 -0.031 0.004
1 P002 NaN ARR-0002 0.010 -0.028 0.011
2 P003 old exome_id: EX-0003 ARR-0003 -0.044 0.015 -0.002
3 P004 old exome_id: EX-0004 ARR-0004 0.009 -0.030 0.007
4 P007 NaN ARR-0007 -0.051 0.019 0.003
5 NaN NaN CTRL-NA12878 0.011 -0.029 0.005
CTRL-NA12878 is a control sample, not a patient, so it isn't in the ID map.
It is kept with an empty patient_id, and the note tells you it happened.
Templated files¶
Some catalog entries stand for several files, one per value of a parameter. HIBAG calls, for instance, exist for each reference panel. Pass the parameter by name:
hibag = demo.load("hla.hibag", panel="african")
patient_id qc_note sample.id A prob_A
0 P001 NaN ARR-0001 A*02:01/A*24:02 0.97
1 P002 NaN ARR-0002 A*01:01/A*03:01 0.99
2 P003 old exome_id: EX-0003 ARR-0003 A*11:01/A*24:02 0.88
3 P004 old exome_id: EX-0004 ARR-0004 A*02:01/A*02:01 0.95
4 P007 NaN ARR-0007 A*30:01/A*33:03 0.71
Leave the parameter out and its default is used (european here). An invalid
value raises an error that lists the valid ones.
HLA columns are genotype objects
Every table under hla: loads with one column per gene holding
HLAgenotype objects, not text. Raw HIBAG output has allele1_A and
allele2_A; the HIBAG loader combines them into A, and each genotype's
.prob is set from prob_A (also kept as a column). The objects print as
A*02:01/A*24:02 and can do more than text can:
gt = hibag.loc[0, "A"]
gt.prob # 0.97
gt.copies_of("A*24:02") # 1
gt.genotype(format="standard_no_gene") # "02:01/24:02"
Patients missing from an HLA table after a merge get NaN, not a
genotype object.
Merge tables¶
Give .load() several names and the tables are merged on patient_id:
df = demo.load("pheno", "pcs", "hla.hlahd")
Note: 1 of 7 rows in 'pcs' have no patient_id (array ID not in id_map).
Left out 1 patient(s) marked drop=1 in the id_map (keep_dropped=True to keep them).
pheno.main 7 patients
(1 row(s) of pcs without a patient_id dropped before merging)
+ pcs 5 matched (2 without)
+ hla.hlahd.standardized 5 matched (2 without)
= 7 rows (left join on patient_id)
df[["patient_id", "qc_note", "OS_mo", "IID", "PC1", "seqID", "A"]]
patient_id qc_note OS_mo IID PC1 seqID A
0 P001 NaN 31.2 ARR-0001 0.012 EX-0001 A*02:01/A*24:02
1 P002 NaN 44.0 ARR-0002 0.010 EX-0002 A*01:01/A*03:01
2 P003 old exome_id: EX-0003 12.5 ARR-0003 -0.044 EX-0004 A*11:01/A*24:02
3 P004 old exome_id: EX-0004 27.9 ARR-0004 0.009 EX-0003 A*02:01/A*02:01
4 P005 NaN 50.1 NaN NaN EX-0005 A*03:01/A*68:01
5 P007 NaN 19.0 ARR-0007 -0.051 NaN NaN
6 P008 NaN 38.6 NaN NaN NaN NaN
The PCs and HLA calls came from different ID systems (array and exome), but
both landed on the right patients. Notice that P003 got exome EX-0004 and
P004 got EX-0003: those two samples were swapped at sequencing, and the ID
map corrects it (How IDs link). Their qc_note
records the old assignment.
What the merge does¶
- Starts from the first table and adds the others with a left join, so
every row of the first table is kept. Put the table that defines your
study population first (usually
"pheno"). Usehow="inner"to keep only patients present in every table, orhow="outer"for everyone. - Prints a report: how many patients each table matched and how many it
didn't. Read it. It is the quickest way to notice that a table covers far
fewer patients than you expected.
verbose=Falseturns it off. - Drops rows with no
patient_id(like the control sample above) before merging, since they can't be joined, and says how many. - Refuses to merge a table with more than one row per patient. Duplicate patients would silently multiply rows in the result, so instead you get an error naming the patients. Load that table on its own to look at it.
- Prefixes clashing column names. If two tables share a column name (other
than
patient_id), the later table's copy gets a prefix, e.g.pcs.site. The first table keeps its names. The report lists any prefixed columns.
Parameters apply to whichever tables use them:
demo.load("pheno", "hla.hibag", panel="asian")
Dropped patients and QC notes¶
The ID map can have two optional columns about sample QC:
drop:1means "don't use this patient". They are left out of every table.load()returns, and out of survival analyses. Blank or0means keep.qc_note: free text, for the record. It never removes anyone. Use it for things like an old sample ID or why a patient was dropped.
.load() says when it leaves anyone out. To see them, keep them:
demo.load("pheno", keep_dropped=True)
patient_id qc_note drop ID age female site OS_mo OS_event
0 P001 NaN 0 P001 64 0 A 31.2 1
1 P002 NaN 0 P002 58 1 A 44.0 0
2 P003 old exome_id: EX-0003 0 P003 71 1 B 12.5 1
3 P004 old exome_id: EX-0004 0 P004 66 0 B 27.9 1
4 P005 NaN 0 P005 49 1 C 50.1 0
5 P006 sex mismatch (genetic male, recorded female) 1 P006 60 1 C 8.4 1
6 P007 NaN 0 P007 73 0 A 19.0 1
7 P008 NaN 0 P008 55 0 B 38.6 0
With keep_dropped=True the ID map's drop column is added as well, so you
can filter later. .dropped_patients() lists who is marked:
demo.dropped_patients()
{'P006'}
qc_note=False leaves the note column out of what you load.
Genetic data¶
VCF and PLINK files aren't tables you load into pandas; you hand them to
bcftools, plink or the c4 commands. .path() gives you the file:
demo.path("dna.imputed", chrom=6)
/data/lab/demo/dna/chr6.dose.vcf.gz
.paths() gives every variant of a templated entry, keyed by parameter value,
which makes per-chromosome loops easy:
import subprocess
for chrom, vcf in demo.paths("dna.imputed").items():
subprocess.run(["bcftools", "index", "-n", str(vcf)], check=True)
The sample names inside a genetic file are IDs of that data type (imputed
IDs for dna.imputed). To connect them to patients, use the ID map:
ids = demo.id_map() # the full table
to_patient = dict(zip(ids["imputed_id"], ids["patient_id"]))
.samples() gives the sample names to keep, leaving out patients marked
drop=1. Write them to a file to subset a VCF:
samples = demo.samples("imputed") # ["ARR-0001", "ARR-0002", ...]
Path("keep.txt").write_text("\n".join(samples) + "\n")
# bcftools view -S keep.txt chr6.dose.vcf.gz
To pull a few variants out of the VCFs as a table, with patient_id
attached, see Extracting variants.
Data dictionary¶
If a phenotype table has a data dictionary in the catalog, .describe() looks
up columns:
demo.describe("OS_mo")
Overall survival (months) (numeric)
Lookups are case-insensitive, and a near miss suggests the closest match:
demo.describe("os_months")
Column 'os_months' not found. Did you mean: 'OS_mo'?
Several columns print as a table:
demo.describe("age", "female", "OS_event")
+----------+------------------------------------------+
| Column | Definition |
+----------+------------------------------------------+
| age | Age at treatment start (years) (numeric) |
| female | 1 = female, 0 = male (binary) |
| OS_event | 1 = death, 0 = censored (binary) |
+----------+------------------------------------------+
describe uses "pheno" by default; pass table= with any
other phenotype entry that has a data dictionary.
Check the catalog¶
.check() walks every entry and reports anything that doesn't line up:
missing files, paths still marked TODO, ID columns that aren't in the file,
and ID map problems.
problems = demo.check()
demo: 2 problem(s)
✗ dna.exome: path not set (TODO)
✗ hla.hibag.raw: 1 of 4 file(s) missing, e.g. /data/lab/demo/hla/hibag_multiethnic.csv
It also returns the problems as a list (empty when everything is fine), so a script can fail on them. Run it after editing the catalog or moving data.
Older loaders¶
The functions in c4.metadata from before the catalog still work, so existing
notebooks and scripts don't break. They take the same arguments and return
exactly what they always did: the file as stored, without patient_id or
qc_note, and without leaving out drop=1 patients. Only where the files are
found has changed: they now come from the catalog too.
| Older call | New equivalent |
|---|---|
load_pheno("<name>") |
cohort("<cohort>").load("pheno") |
load_pcs("<name>") |
cohort("<cohort>").load("pcs") |
load_hla("<name>", caller="hlahd") |
cohort("<cohort>").load("hla.hlahd") |
load_hla("<name>", panel="african") (HIBAG) |
cohort("<cohort>").load("hla.hibag.raw", panel="african") |
load_pheno("<name>", df=pcs, left_ID=..., right_ID="IID", how="inner") |
cohort("<cohort>").load("pheno", "pcs", how="inner") |
load_qc("<name>") |
cohort("<cohort>").id_map() |
The table in c4/metadata.py lists which old names map to which cohort and
catalog entry.
Use cohort() for new code.