Rohit GuptaTidy sample metadata with pandas before any modeling
Multi-platform tumor studies have a lot of moving parts. In a project like our work on gliomas arising in patients with neurofibromatosis type 1 (“Multiplatform molecular analyses refine classification of gliomas arising in patients with neurofibromatosis type 1”, Acta Neuropathologica, 2022), each sample can have sequencing calls, a methylation array, a copy number profile and a row in a clinical spreadsheet. Every one of those comes from a different system, and the only thing tying them together is a sample ID that someone typed.
The models are rarely where things go wrong. The usual failure is quieter: a sample that silently drops out of an inner join, a re-run array counted twice, a grade stored as text so “IV” sorts before “2”, or “no row in the variant table” being read as “wild type” for a sample that was never sequenced. None of these raise an error. They just change your results.
This post is the pandas routine I use to make the metadata boring before any modeling starts. All the data is synthetic and built in code (twelve samples called S001 to S012), so you can paste each block into a notebook and run it top to bottom. The code targets pandas 2.x (it also runs on pandas 3.0), and the Parquet section needs pyarrow.
The synthetic cohort
Four sources, each with the kinds of mess I see in real projects:
- an Illumina-style sample sheet for the methylation arrays, with a header section and a re-run array
- a clinical table exported from a spreadsheet, with inconsistent IDs, grades written three different ways and a free-text date
- a long table of sequencing variant calls, including a sample that is not in the clinical table
- a wide copy number table, one column per gene
from io import StringIO
from pathlib import Path
import pandas as pd
data_dir = Path("demo_data")
data_dir.mkdir(exist_ok=True)
sheet_text = """[Header]
Investigator Name,Demo
Project Name,Synthetic glioma cohort
Date,2020-07-01
[Manifest]
A,MethylationEPIC
[Data]
Sample_Name,Sample_Well,Sample_Plate,Sentrix_ID,Sentrix_Position
S001,A01,PLATE1,203456780001,R01C01
s002,B01,PLATE1,203456780001,R02C01
S-003,C01,PLATE1,203456780001,R03C01
S004 ,D01,PLATE1,203456780001,R04C01
S005,E01,PLATE1,203456780001,R05C01
S006,F01,PLATE1,203456780001,R06C01
S007,G01,PLATE1,203456780001,R07C01
S008,H01,PLATE1,203456780001,R08C01
S009,A02,PLATE1,203456780002,R01C01
S010,B02,PLATE1,203456780002,R02C01
S0011,C02,PLATE1,203456780002,R03C01
S007,D02,PLATE1,203456780002,R04C01
"""
(data_dir / "sample_sheet.csv").write_text(sheet_text)
clinical_raw = pd.DataFrame({
"Sample ID": [" s001", "S002", "S003", "S004", "S005", "S006",
"S007", "S008", "S009", "S010", "S011", "S012"],
"Patient": ["P01", "P01", "P02", "P03", "P04", "P05",
"P06", "P07", "P08", "P09", "P10", "P11"],
"Age at dx": pd.array([8, 8, 34, 51, None, 12, 45, 29, 63, 17, 40, 22], dtype="Int64"),
"Sex": ["F", "F", "M", "m", "F", "M", "F", "M", "F", "unknown", "M", "F"],
"WHO grade": ["1", "1", "II", "grade 4", "3", None,
"4", "2", "IV", "1", "3", "2"],
"Surgery date": ["2015-03-02", "2016-11-20", "2017-01-15", "2017-06-30",
"not recorded", "2018-02-11", "2018-05-09", "2018-09-23",
"2019-01-07", "2019-04-18", "2019-08-30", "2019-12-02"],
"NF1 status": ["germline", "germline", "germline", "sporadic", "germline",
"germline", "sporadic", "germline", "germline", "germline",
"sporadic", "germline"],
})
clinical_raw.to_csv(data_dir / "clinical.csv", index=False)
variants_raw = pd.DataFrame({
"sample": ["S001", "S001", "S002", "S003", "S003", "S004", "S004",
"S006", "S007", "S007", "S009", "s010", "S011", "S013"],
"gene": ["NF1", "FGFR1", "NF1", "NF1", "ATRX", "NF1", "TP53",
"NF1", "NF1", "ATRX", "NF1", "NF1", "TP53", "NF1"],
"variant_class": ["frameshift", "missense", "frameshift", "nonsense",
"frameshift", "splice", "missense", "nonsense",
"frameshift", "nonsense", "missense", "splice",
"missense", "nonsense"],
})
sequenced_raw = pd.Series(["S001", "S002", "S003", "S004", "S006", "S007",
"S008", "S009", "S010", "S011", "S013"])
cn_wide_raw = pd.DataFrame({
"sample_id": ["S001", "S002", "S003", "S004", "S006", "S007", "S008", "S009", "S010"],
"CDKN2A": ["neutral", "neutral", "homdel", "homdel", "neutral",
"loss", "neutral", "homdel", "neutral"],
"PTEN": ["neutral", "neutral", "loss", "neutral", "neutral",
"neutral", "neutral", "loss", "neutral"],
"EGFR": ["neutral", "gain", "neutral", "neutral", "neutral",
"neutral", "amp", "neutral", "neutral"],
})
Every value here is made up. The gene names are there so the tables look familiar, not to say anything about biology.
Reading an Illumina sample sheet
Illumina sample sheets put a few metadata sections ([Header], [Manifest] and so on) above the actual table, which starts after a line reading [Data]. If the header is always exactly the same length you can get away with pd.read_csv(path, skiprows=7). I don’t trust that. Sheets get edited by hand, sections get added, and Excel likes to pad lines with trailing commas ([Data],,,,). Finding the marker is only a few lines more:
def read_illumina_sample_sheet(path):
lines = Path(path).read_text().splitlines()
marker = [i for i, line in enumerate(lines)
if line.strip().strip(",").lower() == "[data]"]
if len(marker) != 1:
raise ValueError(f"{path}: expected one [Data] line, found {len(marker)}")
table = "\n".join(lines[marker[0] + 1:])
return pd.read_csv(StringIO(table), dtype="string")
sheet = read_illumina_sample_sheet(data_dir / "sample_sheet.csv")
sheet.dtypes
Two details matter. First, dtype="string" reads every column as text. Sentrix_ID is a 12-digit barcode, not a number, and the moment pandas parses it as an integer (or worse, Excel shows it in scientific notation) you have a problem waiting for you. Second, the function refuses to guess when there are zero or two [Data] lines. Failing loudly on a malformed file is cheaper than debugging the consequences.
The array identifier that methylation pipelines use (the IDAT “basename”) is the barcode and position joined by an underscore, so I build it right away:
sheet = sheet.assign(array_id=sheet["Sentrix_ID"] + "_" + sheet["Sentrix_Position"])
Reading clinical and variant tables
For CSV exports I tell pandas what types I expect instead of letting it infer them, and I use the nullable dtype backend so missing values don’t turn integer columns into floats:
clinical_in = pd.read_csv(
data_dir / "clinical.csv",
dtype={"Sample ID": "string", "Patient": "string"},
dtype_backend="numpy_nullable",
)
clinical_in.dtypes
Age at dx comes back as Int64 (capital I, the nullable integer) rather than float64, because one value is missing. I’ll come back to nullable types below.
Excel works the same way through pd.read_excel, which needs openpyxl installed for .xlsx files:
clinical_raw.to_excel(data_dir / "clinical.xlsx", sheet_name="Clinical", index=False)
clinical_xl = pd.read_excel(
data_dir / "clinical.xlsx",
sheet_name="Clinical",
dtype={"Sample ID": "string", "Patient": "string"},
)
The pandas side of Excel is fine. The trouble is what Excel did to the file before you got it: leading zeros stripped from IDs, date-like strings converted to dates, and the well-known habit of turning some gene symbols into dates. If you can get a CSV export straight from the source system, take it. If you can’t, read with explicit dtype and check the IDs (next section) before anything else.
Normalizing sample IDs
Look at what we have: " s001", "s002", "S-003", "S004 ", "S0011". Humans read all of these as the same kind of thing. A join does not. I write one function that defines the canonical form and use it on every table:
ID_PATTERN = r"^S[-_ ]?0*(\d+)$"
def normalize_sample_id(raw: pd.Series) -> pd.Series:
cleaned = raw.astype("string").str.strip().str.upper()
number = cleaned.str.extract(ID_PATTERN, expand=False)
return "S" + number.str.zfill(3)
def check_ids(raw: pd.Series, normalized: pd.Series, source: str) -> None:
bad = raw[normalized.isna() & raw.notna()]
if not bad.empty:
raise ValueError(f"{source}: could not parse sample IDs {bad.tolist()}")
ids = normalize_sample_id(sheet["Sample_Name"])
check_ids(sheet["Sample_Name"], ids, "sample sheet")
sheet = sheet.assign(sample_id=ids)
sheet[["Sample_Name", "sample_id"]]
The regex strips an optional separator and any leading zeros, and zfill(3) pads back to three digits, so "S0011" and "S11" both become "S011". Anything that doesn’t match becomes <NA>, and check_ids turns that into an error that names the offending values. The important part is not this particular regex (yours will depend on your ID scheme) but that the rule lives in exactly one place. Ad hoc .str.strip() calls scattered across notebooks are how two tables end up with slightly different ideas of what an ID is.
Validate before you join
Duplicates
The sample sheet has S007 twice: the second array is a re-run. That’s legitimate at the array level and wrong at the sample level. First, look:
sheet[sheet["sample_id"].duplicated(keep=False)]
keep=False marks every copy, not just the second one, which is what you want when inspecting. Then make the rule explicit. Here I keep the last array listed for each sample (the re-run) and print what was dropped, so the decision is visible in the notebook output:
dropped = sheet[sheet["sample_id"].duplicated(keep="last")]
print("Dropping arrays:", dropped["array_id"].tolist())
arrays = sheet.drop_duplicates("sample_id", keep="last")
In a real project the rule might be “keep the array that passed QC” or “keep the higher-purity sample”, and it should come from a QC column, not from row order. Whatever it is, write it down in code.
Missing keys
Before merging, I compare key sets directly. Set differences answer the question “who is going to fall out of this join?” in one line each:
def key_report(left, right, key, names=("left", "right")):
lk, rk = set(left[key].dropna()), set(right[key].dropna())
return {f"only_in_{names[0]}": sorted(lk - rk),
f"only_in_{names[1]}": sorted(rk - lk)}
Merging with validate and indicator
DataFrame.merge has two arguments I now use on almost every join. validate checks the relationship you expect (“one_to_one”, “one_to_many”, “many_to_one”) and raises MergeError if the data disagrees. indicator=True adds a _merge column saying whether each row came from the left table, the right table or both.
Here is validate catching the duplicate array if we had forgotten to resolve it:
clinical_ids = clinical_in.assign(sample_id=normalize_sample_id(clinical_in["Sample ID"]))
try:
clinical_ids.merge(sheet, on="sample_id", validate="one_to_one")
except pd.errors.MergeError as err:
print("MergeError:", err)
Without validate, that merge succeeds and S007 appears twice in everything downstream, including the counts in your Table 1.
Categorical dtypes and ordered grades
WHO grade is ordinal. Stored as text, "IV" sorts before "2" and you can’t ask for “grade 3 or higher”. Stored as numbers, it looks continuous and invites someone to compute a mean grade. An ordered categorical is the right type: it has a fixed set of allowed values and a defined order.
GRADE_MAP = {"1": "1", "I": "1", "2": "2", "II": "2",
"3": "3", "III": "3", "4": "4", "IV": "4"}
GRADE_DTYPE = pd.CategoricalDtype(categories=["1", "2", "3", "4"], ordered=True)
def parse_grade(raw: pd.Series) -> pd.Series:
key = (raw.astype("string").str.upper()
.str.replace(r"^(CNS\s+)?(WHO\s+)?(GRADE\s+)?", "", regex=True)
.str.strip())
unknown = key.notna() & ~key.isin(list(GRADE_MAP))
if unknown.any():
raise ValueError(f"unrecognized grades: {key[unknown].unique().tolist()}")
return key.map(GRADE_MAP).astype(GRADE_DTYPE)
SEX_DTYPE = pd.CategoricalDtype(categories=["F", "M"])
def parse_sex(raw: pd.Series) -> pd.Series:
s = raw.astype("string").str.strip().str.upper()
return s.where(s.isin(["F", "M"])).astype(SEX_DTYPE)
Two choices are baked in here. Unrecognized grades raise, because a new spelling should be looked at by a person. Unrecognized sex values become missing, because “unknown” really is missing for my purposes. You might decide differently; the point is that it’s a decision in code.
Once grade is ordered, comparisons and sorting behave:
grade = parse_grade(clinical_in["WHO grade"])
print(grade.sort_values().tolist())
print((grade >= "3").sum(), "samples are grade 3 or 4")
The categories are the strings "1" to "4" rather than integers. Comparisons follow the category order, not string order, so nothing is lost, and string categories survive the trip to Parquet at the end of this post. Integer categories are less reliable there: with current pyarrow they come back from Parquet as plain numbers.
Categoricals also make value_counts and crosstab show zero counts for categories that exist but don’t appear, which is often exactly the row you need to see.
Parsing dates
Dates in clinical exports arrive as strings, and some of those strings are “not recorded”. I parse with an explicit format and errors="coerce", then check how many values failed:
surgery = pd.to_datetime(clinical_in["Surgery date"], format="%Y-%m-%d", errors="coerce")
failed = clinical_in.loc[surgery.isna() & clinical_in["Surgery date"].notna(), "Surgery date"]
print("Unparsed dates:", failed.tolist())
errors="coerce" on its own is dangerous, since it turns every bad value into NaT without a word. Pairing it with a report of what failed keeps it honest. Passing format matters too: it stops pandas from guessing between day-first and month-first for ambiguous dates.
pd.NA and nullable dtypes
Classic NumPy-backed pandas represents missing values as NaN, which is a float. One missing age turns an integer column into floats, and a boolean column with a gap becomes object. The nullable dtypes ("Int64", "boolean", "string", "Float64") use pd.NA instead and keep their type:
print(pd.Series([8, None, 34]).dtype) # float64
print(pd.Series([8, None, 34], dtype="Int64").dtype) # Int64
pd.NA follows three-valued logic, which is mostly what you want and occasionally surprising:
print(pd.NA | True) # True: the answer is True whatever the missing value is
print(pd.NA & False) # False, for the same reason
print(pd.NA == 1) # <NA>, not False
try:
bool(pd.NA)
except TypeError as err:
print("TypeError:", err)
That last one means if value == something: raises when value is missing, instead of quietly taking the else branch. I consider that a feature. For masks, I fill missing explicitly so the intent is clear: df[(df["age_at_dx"] >= 18).fillna(False)]. DataFrame.convert_dtypes() will convert an existing frame to the nullable types if you didn’t read it that way.
One cleaning function per table, with method chaining
With the pieces defined, each source gets a single function that goes from raw to clean. Method chaining with .assign and .pipe keeps every step in order and avoids the half-modified intermediate frames (and chained-assignment warnings) that come from editing columns in place.
def require_unique(df: pd.DataFrame, key: str) -> pd.DataFrame:
dups = df[df[key].duplicated(keep=False)]
if not dups.empty:
raise ValueError(f"duplicate {key}:\n{dups.sort_values(key)}")
return df
def clean_clinical(raw: pd.DataFrame) -> pd.DataFrame:
out = (
raw
.rename(columns=lambda c: c.strip().lower().replace(" ", "_"))
.rename(columns={"patient": "patient_id"})
.assign(
sample_id=lambda d: normalize_sample_id(d["sample_id"]),
sex=lambda d: parse_sex(d["sex"]),
who_grade=lambda d: parse_grade(d["who_grade"]),
surgery_date=lambda d: pd.to_datetime(
d["surgery_date"], format="%Y-%m-%d", errors="coerce"),
nf1_status=lambda d: d["nf1_status"].astype("category"),
age_at_dx=lambda d: d["age_at_dx"].astype("Int64"),
)
.assign(
days_from_first_surgery=lambda d: (
d["surgery_date"]
- d.groupby("patient_id")["surgery_date"].transform("min")
).dt.days,
)
.pipe(require_unique, "sample_id")
)
check_ids(raw["Sample ID"], out["sample_id"], "clinical")
return out
def clean_variants(raw: pd.DataFrame) -> pd.DataFrame:
return raw.assign(sample_id=lambda d: normalize_sample_id(d["sample"])).drop(columns="sample")
clinical = clean_clinical(clinical_in)
variants = clean_variants(variants_raw)
sequenced = normalize_sample_id(sequenced_raw)
clinical.dtypes
The groupby(...).transform("min") line computes each patient’s first surgery date and broadcasts it back to every row, which gives a per-sample interval since first surgery (P01’s second sample, S002, gets a positive value). Note that this is a patient-level quantity attached to sample rows, which is exactly the kind of thing that gets confused when you have more than one sample per patient. Keep patient_id around even if you think you won’t need it.
Now the joins, each with its expected relationship stated. The methylation classifier output is keyed by array, so I attach it to the deduplicated sheet first:
meth_results = pd.DataFrame({
"array_id": sheet["array_id"],
"meth_class": ["MC1", "MC1", "MC2", "MC3", "MC2", "MC1",
"MC3", "MC1", "MC3", "MC1", "MC2", "MC3"],
"meth_score": [0.97, 0.91, 0.88, 0.95, 0.42, 0.93,
0.61, 0.90, 0.99, 0.86, 0.77, 0.94],
})
arrays_scored = arrays.merge(meth_results, on="array_id", how="left", validate="one_to_one")
print(key_report(clinical, arrays_scored, "sample_id", ("clinical", "arrays")))
cohort = clinical.merge(
arrays_scored[["sample_id", "array_id", "meth_class", "meth_score"]],
on="sample_id", how="outer", validate="one_to_one", indicator=True,
)
print(cohort["_merge"].value_counts())
The key report and the _merge counts tell the same story from two angles: S012 is in the clinical table but has no array. That might be fine (array not run yet) or a typo in the sheet. Either way you now know about it, and you can decide whether downstream analyses should use how="inner" deliberately rather than by accident.
Summaries with groupby, agg and crosstab
Named aggregation in agg gives readable column names and makes the summary self-documenting:
cohort = cohort.assign(meth_class=lambda d: d["meth_class"].astype("category"))
summary = (
cohort
.groupby("meth_class", observed=True)
.agg(
n_samples=("sample_id", "size"),
n_patients=("patient_id", "nunique"),
median_age=("age_at_dx", "median"),
highest_grade=("who_grade", "max"),
low_score=("meth_score", lambda s: (s < 0.9).sum()),
)
)
summary
I pass observed=True explicitly. For categorical group keys, it controls whether categories with no rows show up as empty groups. The default flipped from False to True in pandas 3.0 (pandas 2.x warns when you leave it unset), so spelling it out keeps the output the same on both. Counting samples and patients side by side is a cheap check that catches the “one patient, several samples” issue early.
For two categorical variables, pd.crosstab is the fastest way to see the joint distribution:
pd.crosstab(cohort["meth_class"], cohort["who_grade"], margins=True)
If a cell you expect to be populated is zero, or a margin total doesn’t match the number of samples you think you have, stop and find out why before going further.
Tidy reshaping with melt and pivot_table
Tidy data, in Wickham’s sense, means one variable per column and one observation per row. The copy number table is wide (genes as columns), which is convenient for reading and awkward for filtering or combining with other platforms. melt makes it long:
CN_DTYPE = pd.CategoricalDtype(["homdel", "loss", "neutral", "gain", "amp"], ordered=True)
cn_long = (
cn_wide_raw
.melt(id_vars="sample_id", var_name="gene", value_name="cn_call")
.assign(cn_call=lambda d: d["cn_call"].astype(CN_DTYPE))
)
cn_long[cn_long["cn_call"] <= "loss"]
Long tables from different platforms can then be stacked into a single event table, which is the shape most oncoprint and summary code wants:
events = pd.concat([
variants.assign(platform="sequencing", event=variants["variant_class"])
[["sample_id", "gene", "platform", "event"]],
cn_long[cn_long["cn_call"] != "neutral"]
.assign(platform="copy_number", event=lambda d: d["cn_call"].astype("string"))
[["sample_id", "gene", "platform", "event"]],
], ignore_index=True)
Going the other way, pivot_table turns long calls into a sample by gene matrix for modeling. This is where the “absent row means wild type” trap lives:
mut_matrix = (
variants
.assign(mutated=True)
.pivot_table(index="sample_id", columns="gene", values="mutated",
aggfunc="any", fill_value=False)
.reindex(sequenced, fill_value=False)
.astype("boolean")
.add_prefix("mut_")
.rename_axis(index="sample_id", columns=None)
.reset_index()
)
cohort = cohort.merge(mut_matrix, on="sample_id", how="left", validate="one_to_one")
cohort[["sample_id", "mut_NF1", "mut_TP53"]]
The reindex(sequenced, fill_value=False) step adds rows for samples that were sequenced and had no calls (S008), and fills them with False. Samples that were never sequenced (S005, S012) don’t get a row, so after the left merge they are <NA>. That distinction is the whole point: “we looked and found nothing” and “we didn’t look” are different facts, and a model will treat False as evidence.
Notice also that S013 has a variant call but appears in neither the clinical table nor the arrays. The left merge drops it, which is probably right, but key_report(variants, clinical, "sample_id") is how you find out it existed.
Writing Parquet for reproducibility
CSV forgets everything you just did. Read the cohort back from CSV and the categories, their order, the nullable integers and the dates are all gone. Parquet (through pyarrow) stores the schema with the data:
cohort_out = cohort.drop(columns="_merge")
cohort_out.to_parquet(data_dir / "cohort.parquet", engine="pyarrow", index=False)
back = pd.read_parquet(data_dir / "cohort.parquet")
print(back["who_grade"].dtype == GRADE_DTYPE) # True: categories and order kept
print(back[["age_at_dx", "surgery_date", "mut_NF1"]].dtypes)
I check the columns I care about rather than comparing every dtype exactly. Some details (for example, the dtype of the category labels themselves) can differ after a round trip depending on your pandas and pyarrow versions, without affecting anything downstream.
I treat the Parquet file as the hand-off point between “metadata wrangling” and “analysis”. Analysis notebooks read the Parquet file and never touch the raw spreadsheets, so there is exactly one place where IDs get normalized and duplicates get resolved.
A small data contract
The last step before writing is a function that states what a valid cohort table looks like. It doesn’t need a framework. A function that returns a list of problems is enough:
def validate_cohort(df: pd.DataFrame) -> list[str]:
problems = []
required = {"sample_id", "patient_id", "who_grade", "meth_class", "age_at_dx"}
missing = required - set(df.columns)
if missing:
return [f"missing columns: {sorted(missing)}"]
if df["sample_id"].isna().any():
problems.append("sample_id has missing values")
if df["sample_id"].duplicated().any():
problems.append("sample_id is not unique")
if not df["sample_id"].str.fullmatch(r"S\d{3}").fillna(False).all():
problems.append("sample_id does not match S###")
if df["who_grade"].dtype != GRADE_DTYPE:
problems.append(f"who_grade has dtype {df['who_grade'].dtype}")
if (df["age_at_dx"] < 0).fillna(False).any():
problems.append("negative age_at_dx")
if not df["meth_score"].dropna().between(0, 1).all():
problems.append("meth_score outside [0, 1]")
return problems
problems = validate_cohort(cohort_out)
if problems:
raise ValueError("cohort failed validation:\n" + "\n".join(problems))
I prefer this to a string of bare assert statements for two reasons: it reports every problem at once instead of stopping at the first, and assert is skipped entirely when Python runs with -O. If the contract grows large, libraries such as pandera let you declare schemas instead of writing checks by hand, but the hand-written version covers most of what a single study needs. Run it right before to_parquet and again at the top of every analysis notebook that reads the file.
Checklist
Before any modeling, on every new data drop:
- Read with explicit
dtype(IDs and barcodes as"string"), and usedtype_backend="numpy_nullable"for CSVs. - Locate the
[Data]section in sample sheets rather than hard-codingskiprows. - Normalize sample IDs with one shared function, and raise on anything it can’t parse.
- Inspect duplicates with
duplicated(keep=False)and resolve them with a rule written in code. - Compare key sets before joining, then merge with
validate=andindicator=True. - Convert ordinal and nominal fields to categoricals, with
ordered=Truewhere order means something. - Parse dates with an explicit
format, and report every value thaterrors="coerce"turned intoNaT. - Keep “not tested” (
<NA>) separate from “tested, negative” (False). - Pass
observed=Truetogroupbyon categoricals, and count both samples and patients. - Validate, write Parquet, and have analysis code read only the Parquet file.
The next step after this is the fun part (a classifier, survival models, whatever the study is for), and it goes much better when you trust the table underneath it.
References
- McKinney W. Data Structures for Statistical Computing in Python. Proceedings of the 9th Python in Science Conference (2010).
- Wickham H. Tidy Data. Journal of Statistical Software 59(10) (2014).
- pandas documentation: Merge, join, concatenate and compare, Categorical data, Working with missing data.