# Filtering and query > **In one line:** boolean masks with `.loc`, or `.query()` for a readable string expression — and the same `&`/`|`/parentheses rules as NumPy. ## Boolean masks ```python df.loc[df["padj"] < 0.05] df.loc[(df["padj"] < 0.05) & (df["log2fc"].abs() > 1)] df.loc[df["chrom"].isin(["chr17", "chr13"])] df.loc[~df["gene"].isna()] df.loc[df["af"].between(0.01, 0.05)] df.loc[df["gene"].str.startswith("BRCA")] ``` **Use `&`, `|`, `~`, never `and`, `or`, `not`**, and parenthesise every comparison. Same reasoning as Boolean Masking: Python's keywords call `bool()` on the whole Series, which raises "truth value is ambiguous". Masks are reusable objects, which is often clearer than one giant expression: ```python is_sig = df["padj"] < 0.05 is_strong = df["log2fc"].abs() > 1 is_coding = df["biotype"] == "protein_coding" df.loc[is_sig & is_strong & is_coding] print(f"sig: {is_sig.sum()}, strong: {is_strong.sum()}, both: {(is_sig & is_strong).sum()}") ``` Naming your masks also gives you the counts for free, which is how you notice that a filter removed 99% of your data. ## query ```python df.query("padj < 0.05 and abs(log2fc) > 1") df.query("chrom in ['chr17', 'chr13']") df.query("gene == @target_gene") # @ references a Python variable df.query("`gene name` == 'TP53'") # backticks for names with spaces df.query("padj < 0.05").query("baseMean > 10") # chainable ``` Inside `query` you use `and`/`or`/`not` (it is a separate mini-language) and you can skip `df[...]` around every column. For a filter with three or more conditions it is substantially more readable. Downsides: no editor autocomplete, typos in column names are runtime errors, and it is marginally slower for small frames (it parses the string). For large frames with `numexpr` installed it can be faster. ## Filtering in a chain ```python result = ( df.query("baseMean > 10") .loc[lambda d: d["padj"] < 0.05] # lambda gets the current frame .sort_values("padj") ) ``` The `lambda` form of `.loc` is what makes filtering work mid-chain — you cannot reference `df` because the intermediate frame has no name. ## Filtering rows by group properties ```python # keep genes seen in >10 samples df.groupby("gene").filter(lambda g: len(g) > 10) # faster for simple counts counts = df["gene"].value_counts() df.loc[df["gene"].isin(counts[counts > 10].index)] # keep the single best row per gene df.loc[df.groupby("gene")["padj"].idxmin()] df.sort_values("padj").drop_duplicates("gene", keep="first") # equivalent, faster ``` That last pattern — sort then `drop_duplicates` — is the idiomatic "best row per group" and much faster than a `groupby().apply()`. ## Top-N ```python df.nlargest(20, "log2fc") df.nsmallest(20, "padj") df.sort_values("padj").head(20) df.nlargest(5, "log2fc", keep="all") # keep ties ``` `nlargest` uses a partial sort — meaningfully faster than a full sort on large frames. ## Deduplication ```python df.duplicated(subset=["chrom", "pos", "ref", "alt"]).sum() df.drop_duplicates(subset=["chrom", "pos", "ref", "alt"], keep="first") df.loc[df.duplicated(subset=["chrom", "pos"], keep=False)] # SHOW all dupes ``` `keep=False` marks *every* member of a duplicate set, which is what you want when investigating rather than removing. ## Common filter recipes for bioinformatics ```python # significant DE genes sig = de.query("padj < 0.05 and abs(log2FoldChange) > 1") # rare, high-impact variants rare = variants.loc[ (variants["gnomad_af"].fillna(0) < 0.001) & (variants["impact"].isin(["HIGH", "MODERATE"])) & (variants["filter"] == "PASS") ] # autosomes only, canonical contigs AUTOSOMES = [f"chr{i}" for i in range(1, 23)] df.loc[df["chrom"].isin(AUTOSOMES)] # expressed genes: >1 CPM in at least 3 samples keep = (cpm > 1).sum(axis=1) >= 3 expr_filtered = expr.loc[keep] # a genomic window df.loc[(df["chrom"] == "chr17") & df["pos"].between(7_660_000, 7_690_000)] # QC pass reads.query("mapq >= 30 and not is_duplicate and not is_secondary") ``` Note `variants["gnomad_af"].fillna(0)` in the second example: a variant absent from gnomAD has NA frequency, and NA fails the `< 0.001` test, so without the fill you would silently *exclude* the rarest variants — the exact opposite of what you want. **Missing values in a filter condition are always excluded**, which is a very easy way to lose the rows you care about most. See Missing Data in pandas. ## Always report what you filtered ```python def report_filter(df, mask, name): kept, total = mask.sum(), len(df) print(f"{name}: kept {kept:,}/{total:,} ({kept/total:.1%})") return df.loc[mask] df = report_filter(df, df["filter"] == "PASS", "PASS variants") df = report_filter(df, df["gnomad_af"].fillna(0) < 0.01, "rare") ``` Two lines that turn "my analysis found nothing" into "oh, step 3 removed 99.8% of rows". Filter accounting belongs in every pipeline and in every methods section. ## Common mistakes - **`and`/`or` instead of `&`/`|`.** - **Missing parentheses** around comparisons. - **NA silently failing every comparison**, dropping rows you wanted. - **Chained assignment after filtering** — `df[mask]["col"] = x` is a silent no-op. See Copy-on-Write. - **Not calling `reset_index(drop=True)`** after filtering, leaving a gappy index. - **`isin` against a list with different formatting** — `chr1` vs `1`, or Ensembl IDs with vs without version suffix. - **Filtering the frame but not a parallel array/list** you are tracking separately. - **Not counting what you removed.** ## See also loc vs iloc · Boolean Masking · Missing Data in pandas · groupby · assign and pipe · Copy-on-Write22:["$","$L1c","filtering-and-query",{"item":{"slug":"filtering-and-query","title":"Filtering and query","track":"pandas","type":"lesson","difficulty":"Beginner","order":6,"contentMarkdown":"$37","starterCode":"","editorSeed":"","testCases":[],"mockFiles":{},"hints":[]},"index":5}] 38:T1732,# Index Objects > **In one line:** the row (and column) labels — an immutable, array-like object that drives alignment, joins and lookups. ## What an Index is ```python idx = pd.Index(["TP53", "BRCA1", "EGFR"], name="gene") idx.name; idx.dtype; idx.is_unique; idx.is_monotonic_increasing idx.tolist(); idx.to_numpy() "TP53" in idx # hash-based membership — O(1) idx.get_loc("BRCA1") # → 1 ``` It is **immutable**. `idx[0] = "X"` raises. To change labels you build a new Index (`df.index = [...]`, `df.rename(...)`). Both axes have one: `df.index` for rows, `df.columns` for columns. They are the same class, which is why `df.columns.str.startswith("sample_")` works. ## The three jobs an Index does **1. Alignment.** Every binary operation between two pandas objects aligns on the index first. This is the source of both pandas' power and its most confusing behaviour. See Series. **2. Fast lookup.** A unique Index builds a hash table, so `df.loc["TP53"]` is O(1) rather than a scan. On a 20,000-gene frame that is the difference between instant and noticeable. **3. Metadata.** The index says what each row *is*. `expr.loc["TP53", "sample_3"]` is self-documenting in a way `expr[4711, 2]` never is. ## Index types | Type | Use | |---|---| | `Index` | general (str, object) | | `RangeIndex` | default 0..n-1; stores only start/stop/step, so it is free | | `DatetimeIndex` | time series; enables `.loc["2026-01"]` partial-string slicing | | `CategoricalIndex` | repeated labels, memory-efficient | | `IntervalIndex` | ranges — genomic intervals! see below | | `MultiIndex` | hierarchical — see MultiIndex | ## Setting and resetting ```python df.set_index("gene") # a column becomes the index df.set_index(["chrom", "pos"]) # → MultiIndex df.reset_index() # index becomes a column again df.reset_index(drop=True) # discard it entirely df.rename_axis("gene_id") # rename the index itself df.index = new_labels # direct assignment ``` **When to set an index:** when you will repeatedly look up by that key, when you will align two frames on it, or when it genuinely identifies the row. A gene ID on an expression matrix, yes. An arbitrary row number, no — leave the RangeIndex. **`reset_index(drop=True)` after filtering** is a habit worth forming. Otherwise your index has gaps (`0, 3, 7, 11...`), which is harmless until something downstream assumes positions and labels agree. ## Uniqueness matters ```python df.index.is_unique df.index.duplicated().sum() df[df.index.duplicated(keep=False)] # show all duplicate rows ``` A non-unique index makes `.loc["TP53"]` return a **DataFrame** instead of a Series, which then breaks the code expecting a scalar. Worse, joins on a non-unique key multiply rows. Check uniqueness after any operation that could introduce duplicates. ## Set operations ```python a.index.intersection(b.index) # genes in both a.index.union(b.index) a.index.difference(b.index) # genes only in a a.index.isin(panel) # boolean mask df.reindex(new_index) # conform to a new index; NaN for missing df.reindex(new_index, fill_value=0) a.align(b, join="inner", axis=0) # align two frames explicitly ``` `intersection` is the clean way to find shared genes across datasets — much better than converting to sets and losing order. `reindex` is subtly different from `.loc[]`: reindex tolerates missing labels (fills NaN), `.loc` with a list of missing labels raises. Use reindex when you want to force a specific ordering, e.g. matching a metadata table's sample order. ## IntervalIndex for genomic ranges ```python bins = pd.IntervalIndex.from_breaks([0, 1e6, 2e6, 3e6], closed="left") binned = pd.cut(variants["pos"], bins) variants.groupby(binned, observed=True).size() ``` `pd.cut` plus an IntervalIndex gives you binning with readable labels. For real interval *overlap* queries (which genes does this variant fall in?), pandas is the wrong tool — use `pyranges`, `bioframe`, or an interval tree. See Applied - GFF and Genomic Intervals. ## Alignment in action ```python expr = pd.read_csv("counts.tsv", index_col=0) # genes × samples meta = pd.read_csv("samples.tsv", index_col=0) # samples × attributes # align samples explicitly rather than assuming the order matches common = expr.columns.intersection(meta.index) expr = expr[common] meta = meta.loc[common] assert (expr.columns == meta.index).all() ``` That assertion is worth writing every time. Sample order mismatch between a count matrix and a metadata table is a classic, silent, result-inverting bug — you compare the wrong groups and get a beautiful, wrong figure. ## Performance ```python df.index.is_monotonic_increasing # sorted? df.sort_index() # sort — enables faster slicing ``` `.loc` on a sorted index can binary-search; on an unsorted non-unique index it must scan. For large frames you slice repeatedly, sorting once pays for itself. ## Common mistakes - **Duplicate index values** silently changing `.loc` return types and multiplying join output. - **Assuming two frames' rows correspond** because they have the same length. Align explicitly. - **Forgetting `reset_index(drop=True)`** after filtering, leaving gaps. - **`df.index = list(...)` of the wrong length** → ValueError, or worse, silently wrong if lengths coincidentally match. - **Setting an index you then need as a column.** `.reset_index()` gets it back, but grouping and merging are often easier with it as a plain column. - **Expecting the index to survive `to_csv`** without `index=True`, or getting a spurious `Unnamed: 0` column when reading back without `index_col=0`. ## See also Series · DataFrame · MultiIndex · loc vs iloc · Merging and Joining · Applied - GFF and Genomic Intervals23:["$","$L1c","index-objects",{"item":{"slug":"index-objects","title":"Index Objects","track":"pandas","type":"lesson","difficulty":"Beginner","order":7,"contentMarkdown":"$38","starterCode":"","editorSeed":"import pandas as pd\n\nidx = pd.Index([\"TP53\", \"BRCA1\", \"EGFR\"], name=\"gene\")\n\nidx.name; idx.dtype; idx.is_unique; idx.is_monotonic_increasing\nidx.tolist(); idx.to_numpy()\n\"TP53\" in idx # hash-based membership — O(1)\nidx.get_loc(\"BRCA1\") # → 1\n","testCases":[],"mockFiles":{},"hints":[]},"index":6}] 39:T196e,# Merging and Joining > **In one line:** `merge` is a SQL join on columns — and if you are not passing `validate=`, you are one duplicate key away from a silently wrong answer. ## The syntax ```python df.merge(other, on="gene", how="left") df.merge(other, left_on="gene_id", right_on="id", how="inner") df.merge(other, on=["chrom", "pos"], how="outer") df.merge(other, left_index=True, right_index=True, how="inner") df.join(other, how="left") # index-based shorthand ``` ## The join types | `how` | Keeps | |---|---| | `"inner"` | keys in **both** (default) | | `"left"` | all of left, NaN where right has no match | | `"right"` | all of right | | `"outer"` | all keys from either side | | `"cross"` | every combination (no key) | For annotation, **left join** is almost always right: keep every one of your variants/genes, and get NaN where the annotation is missing. An inner join silently drops your unannotated rows, which is how a "clean" result set quietly loses 30% of the data. ## validate= — the most important argument ```python df.merge(other, on="gene", how="left", validate="many_to_one") ``` Options: `"one_to_one"`, `"one_to_many"`, `"many_to_one"`, `"many_to_many"`. pandas raises if the actual key cardinality does not match. **Why this matters.** If the right-hand key is not unique, a left join *multiplies* rows: ```python variants # 1,000,000 rows annotations # gene → description, but 50 genes appear twice merged # 1,000,050 rows — and 50 variants are now counted twice ``` No error. No warning. Every downstream count, mean and test is now subtly wrong. `validate="many_to_one"` turns this into an exception at the point of the mistake. Make it a habit: **every merge gets a `validate=`.** If you do not know the cardinality, that itself is the problem. ## indicator= — did it match? ```python merged = df.merge(annot, on="gene", how="left", indicator=True) merged["_merge"].value_counts() # both 950000 # left_only 50000 ← 50k variants with no annotation ``` Combined with `validate`, this is a complete join audit in two arguments. ## Suffixes ```python df.merge(other, on="gene", suffixes=("_tumour", "_normal")) ``` Overlapping non-key column names get suffixed. The default `("_x", "_y")` is uninformative — always set it explicitly, or you will be reading `af_x` in three weeks wondering which file it came from. ## concat vs merge - **`merge`** — combine columns by matching keys (horizontal, key-based). See below. - **`concat`** — stack frames (vertical) or glue columns positionally/by index. See Concat. ## The ID mismatch problem The most common real-world join failure in bioinformatics is not the join mechanics — it is that the keys do not match. ```python # version suffixes "ENSG00000141510.16" vs "ENSG00000141510" df["gene_id"] = df["gene_id"].str.split(".").str[0] # chromosome naming "chr1" (UCSC) vs "1" (Ensembl) df["chrom"] = df["chrom"].str.removeprefix("chr") # gene symbols vs IDs, and symbol aliases "TP53" vs "P53" vs "ENSG00000141510" # case and whitespace df["gene"] = df["gene"].str.strip().str.upper() # Excel damage — genuinely a documented problem in genomics "SEPT9" → "9-Sep", "MARCH1" → "1-Mar" ``` That last one is real: HGNC renamed several genes (SEPT→SEPTIN, MARCH→MARCHF) in 2020 specifically because Excel kept corrupting them. If a collaborator sends you a gene list as .xlsx, check for date-shaped values. **Always check the overlap before you join:** ```python left_keys = set(df["gene_id"]) right_keys = set(annot["gene_id"]) print(f"left only: {len(left_keys - right_keys):,}") print(f"right only: {len(right_keys - left_keys):,}") print(f"shared: {len(left_keys & right_keys):,}") ``` If "shared" is zero or suspiciously small, it is a formatting mismatch, not biology. ## merge_asof — nearest-key joins For "find the closest preceding record", which comes up with genomic positions and time series: ```python pd.merge_asof( variants.sort_values("pos"), genes.sort_values("start"), left_on="pos", right_on="start", by="chrom", # match within chromosome direction="backward", # nearest gene start at or before the variant ) ``` Both frames must be sorted on the key. `by=` restricts matching to within groups, which is exactly what you want for chromosomes. This is *nearest*, not *overlapping*. For true interval overlap ("which gene contains this variant"), pandas is the wrong tool — use `pyranges` or `bioframe`. See Applied - GFF and Genomic Intervals. ## Checking the join afterwards ```python before = len(df) merged = df.merge(annot, on="gene", how="left", validate="many_to_one", indicator=True) assert len(merged) == before, f"row count changed: {before} → {len(merged)}" print(merged["_merge"].value_counts()) ``` Three lines. They will catch more bugs than any amount of careful reading. ## Bioinformatics examples ```python # annotate DE results with gene info de = de.merge(gene_annot, left_index=True, right_on="gene_id", how="left", validate="one_to_one") # add population frequencies to a variant call set variants = variants.merge( gnomad[["chrom", "pos", "ref", "alt", "af", "af_popmax"]], on=["chrom", "pos", "ref", "alt"], how="left", validate="many_to_one", ) print(f"{variants['af'].isna().mean():.1%} not in gnomAD") # attach sample metadata to a long expression table long = long.merge(metadata, on="sample_id", how="left", validate="many_to_one") # tumour vs normal, same gene paired = tumour.merge(normal, on="gene", suffixes=("_tumour", "_normal")) paired["delta"] = paired["expr_tumour"] - paired["expr_normal"] ``` ## Common mistakes - **No `validate=`.** Row multiplication, silently. - **Inner join where you meant left**, dropping unannotated rows. - **Key format mismatch** — version suffixes, `chr` prefix, case, whitespace. - **Not checking `_merge` counts** after a left join. - **Default `("_x", "_y")` suffixes.** - **Merging on a float column.** Floating-point equality is unreliable. - **Merging on a categorical with different categories** — degrades to object or raises. - **Forgetting `merge_asof` needs sorted input.** - **Using `merge_asof` for interval overlap.** It finds the nearest key, not containment. ## See also Concat · groupby · Index Objects · Missing Data in pandas · Applied - Annotating a Variant Table · Applied - GFF and Genomic Intervals24:["$","$L1c","merging-and-joining",{"item":{"slug":"merging-and-joining","title":"Merging and Joining","track":"pandas","type":"lesson","difficulty":"Beginner","order":8,"contentMarkdown":"$39","starterCode":"","editorSeed":"","testCases":[],"mockFiles":{},"hints":[]},"index":7}] 3a:T15c8,# Missing Data in pandas > **In one line:** unlike NumPy, pandas has real missing-value support — but there are two sentinels (`NaN` and `pd.NA`) and knowing which you have matters. ## Detecting ```python df.isna(); df.notna() # element-wise df.isna().sum() # count per column df.isna().sum(axis=1) # count per row df.isna().any().any() # is there any at all? df["af"].isna().mean() # proportion missing df.isna().sum().sort_values(ascending=False).head(10) # worst columns first ``` `isnull`/`notnull` are exact aliases. Use `isna`/`notna` for consistency with `fillna`, `dropna`. ## The two sentinels | | `np.nan` | `pd.NA` | |---|---|---| | Where | NumPy-backed float columns | nullable and Arrow-backed dtypes | | `NA == NA` | `False` | `pd.NA` (propagates!) | | `bool(NA)` | `False` | raises `TypeError` | | In an int column | impossible — promotes to float | fine | ```python pd.Series([1, 2, None]) # float64 with NaN — int promoted pd.Series([1, 2, None], dtype="Int64") # Int64 with — stays integer ``` `pd.NA` propagating through comparisons is the "correct" three-valued logic (like SQL NULL), and it is why `if x == pd.NA` raises instead of quietly being False. In practice you rarely touch it directly — use `.isna()`. **The integer promotion problem** is the one you meet daily: a genomic position column with one missing value becomes `float64`, and `to_csv` writes `7676154.0`. `Int64` fixes it. ## Dropping ```python df.dropna() # any NA in any column df.dropna(how="all") # only fully-empty rows df.dropna(subset=["padj", "log2fc"]) # ← usually what you want df.dropna(axis=1, thresh=len(df) * 0.8) # drop columns >20% missing ``` **Bare `df.dropna()` is almost always too aggressive.** One sparse annotation column will delete most of your rows. Always name the columns you actually require. ## Filling ```python df.fillna(0) df.fillna({"af": 0, "consequence": "unknown"}) # per-column df["x"].fillna(df["x"].median()) df.ffill(); df.bfill() # forward/backward fill df["x"].interpolate(method="linear") # per-group fill — the common real case df["x"] = df["x"].fillna(df.groupby("gene")["x"].transform("mean")) ``` `ffill` on unsorted data propagates values across group boundaries. Sort first, or group: ```python df = df.sort_values(["sample", "pos"]) df["depth"] = df.groupby("sample")["depth"].ffill() ``` ## NA-aware behaviour ```python df["x"].sum() # skipna=True by default → skips NA df["x"].sum(skipna=False) df["x"].mean() # denominator is the non-NA count df["x"].count() # non-NA count (vs len(df) for all rows) df.groupby("gene")["x"].mean() # NA dropped within each group df.groupby("gene", dropna=False).size() # keep NA as its own group key df["gene"].value_counts(dropna=False) # ← always pass this ``` pandas skipping NA by default is the opposite of NumPy, where one NaN poisons everything. Convenient, but it means **your denominator changes silently.** Always report the count alongside the mean: ```python df.groupby("gene")["af"].agg(["count", "mean"]) ``` `value_counts(dropna=False)` deserves to be a reflex. The default hides missing values from your frequency table, which is exactly where you would want to see them. ## Missing does not mean the same thing everywhere In genomics, the *reason* for missingness is usually informative: | Situation | Meaning | Sensible handling | |---|---|---| | `padj` is NA in DESeq2 | Filtered by independent filtering or flagged as an outlier | Do **not** fill with 1 — exclude from the test set | | Genotype `./.` | No call — low coverage | Keep as missing; imputation is a modelling decision | | Zero in scRNA-seq | A real measurement (dropout) | **Not** missing. Do not impute as NA | | Proteomics below LOD | Left-censored, not random | Model it, or impute low, not with the mean | | Empty VCF `INFO` field | Annotation absent | Depends on the annotator | `df.fillna(0)` applied blindly across all of these produces different kinds of wrong. **Decide per column, and write down why.** The DESeq2 case is worth spelling out: `padj = NA` means "this gene was excluded from multiple-testing correction". Filling with 1 makes it look tested-and-not-significant, which changes your denominator and your FDR. Drop those rows instead. ## Interaction with joins ```python merged = df.merge(annot, on="gene", how="left") merged["annotation"].isna().sum() # ← how many failed to match? ``` A left join produces NA for unmatched keys. **Always check that count.** A high number usually means an ID-format mismatch (Ensembl IDs with version suffixes `ENSG00000141510.16` vs without) rather than genuinely absent annotation. See Merging and Joining. ## Common mistakes - **`df.dropna()` with no `subset=`**, deleting most of your data. - **`fillna(0)` on p-values or fold changes.** Zero is a meaningful value there. - **Filling DESeq2 `padj` NAs with 1.** - **Treating scRNA-seq zeros as missing.** - **`value_counts()` without `dropna=False`.** - **Not noticing an integer column became float** because of one NA. - **`ffill` across group boundaries.** - **Reporting a group mean without the count**, so a mean over 2 of 50 samples looks like a mean over 50. - **Not checking post-join NA counts.** ## See also Missing Data in NumPy · pandas dtypes · groupby · Merging and Joining · Reading and Writing Data25:["$","$L1c","missing-data-in-pandas",{"item":{"slug":"missing-data-in-pandas","title":"Missing Data in pandas","track":"pandas","type":"lesson","difficulty":"Beginner","order":9,"contentMarkdown":"$3a","starterCode":"","editorSeed":"import pandas as pd\n\npd.Series([1, 2, None]) # float64 with NaN — int promoted\npd.Series([1, 2, None], dtype=\"Int64\") # Int64 with — stays integer\n","testCases":[],"mockFiles":{},"hints":[]},"index":8}] 3b:T15f1,# MultiIndex > **In one line:** a hierarchical index — several label levels per row — which lets you represent higher-dimensional data in a 2-D frame. ## Creating one ```python df.set_index(["chrom", "pos"]) df.groupby(["gene", "sample"]).mean() # grouping by 2+ keys produces one df.pivot_table(index="gene", columns=["tissue", "condition"]) pd.MultiIndex.from_tuples([("chr1", 100), ("chr1", 200)], names=["chrom", "pos"]) pd.MultiIndex.from_product([["chr1", "chr2"], [100, 200]], names=["chrom", "pos"]) pd.MultiIndex.from_arrays([chroms, positions], names=["chrom", "pos"]) ``` `from_product` is how you build a complete grid — every chromosome × every bin — which is useful for reindexing sparse results into a full matrix. ## Selecting ```python df.loc["chr1"] # whole first level → drops that level df.loc[("chr1", 7676154)] # a full key df.loc[("chr1", 7676154), "ref"] # key + column df.loc["chr1":"chr3"] # slice the outer level (needs sorting) # partial selection on inner levels needs slice(None) or IndexSlice idx = pd.IndexSlice df.loc[idx[:, 7676154], :] # all chromosomes, that position df.loc[idx["chr1", 100:200], :] df.xs("chr1", level="chrom") # cross-section, drops the level df.xs(7676154, level="pos") # select on an inner level directly ``` `.xs` is the readable option for "give me one value of one level". `IndexSlice` is more general but harder to read. ## Sort before you slice ```python df = df.sort_index() df.index.is_monotonic_increasing # must be True for range slicing ``` Slicing an unsorted MultiIndex raises `UnsortedIndexError`, or is quietly slow. **Sort immediately after creating one.** This is the single most common MultiIndex frustration. ## Level manipulation ```python df.index.names # ['chrom', 'pos'] df.index.get_level_values("chrom") # flat array of that level df.index.levels # the unique values per level df.reset_index(level="pos") # one level → column df.reset_index() # all levels → columns df.swaplevel("chrom", "pos").sort_index() df.droplevel("chrom") df.reorder_levels(["pos", "chrom"]) ``` ## Grouping and stacking ```python df.groupby(level="chrom").size() df.groupby(level=[0, 1]).mean() df.sum(level=...) # removed — use groupby(level=...) df.stack() # columns → index level (wider → longer) df.unstack() # index level → columns (longer → wider) df.unstack(level="condition") df.unstack(fill_value=0) ``` `stack`/`unstack` are the MultiIndex-aware cousins of `melt`/`pivot`. See Reshaping with pivot and melt. ## MultiIndex on columns Just as common, and often the more useful direction: ```python expr.columns = pd.MultiIndex.from_arrays( [tissues, conditions, replicates], names=["tissue", "condition", "rep"] ) expr["liver"] # all liver samples expr.loc[:, idx[:, "treated", :]] # all treated, any tissue expr.T.groupby(level=["tissue", "condition"]).mean().T # collapse replicates ``` That last line — averaging technical replicates within condition — is a genuinely elegant use of a column MultiIndex, and awkward without one. ## When it earns its complexity **Worth it:** - Genuinely hierarchical keys: (chromosome, position), (gene, transcript, exon), (patient, timepoint). - The output of a multi-key `groupby` you will immediately `unstack`. - Repeated selection on the outer level of a large frame — the sorted index makes it fast. - Column hierarchies over a wide experimental design. **Not worth it:** - As a habit. Flat columns plus `groupby` handles most work and is far easier to read, merge and debug. - When you will immediately merge with something else. `merge` on MultiIndexes is fiddly; `reset_index()` first. - When collaborators will read your code. MultiIndex selection syntax has a real learning cost. The honest advice: **reach for it when the data is hierarchical and you will exploit that; otherwise `reset_index()` and move on.** ## Bioinformatics examples ```python # a long-format expression table indexed by gene and sample long = counts.melt(id_vars="gene", var_name="sample", value_name="count") mi = long.set_index(["gene", "sample"]).sort_index() mi.loc["TP53"] # all samples for one gene mi.loc[("TP53", "sample_3"), "count"] # one value mi.groupby(level="gene")["count"].mean() # variants keyed by locus v = variants.set_index(["chrom", "pos"]).sort_index() v.loc["chr17"] # all chr17 variants v.loc[("chr17", slice(7_670_000, 7_680_000))] # a genomic window # a summary table that reads well summary = de.groupby(["tissue", "direction"])["gene"].count().unstack(fill_value=0) ``` ## Common mistakes - **Slicing without sorting** → `UnsortedIndexError`. - **Forgetting the tuple.** `df.loc["chr1", 7676154]` is interpreted as (row, column), not as a two-level key. You need `df.loc[("chr1", 7676154)]`. - **Levels retaining unused categories** after filtering — `df.index.levels` still lists them. Use `df.index.remove_unused_levels()`. - **Merging on a MultiIndex** without resetting first. - **`to_csv` round-trip** — a MultiIndex writes fine but needs `index_col=[0, 1]` on read. - **Using it when flat columns would do.** Complexity you pay for on every subsequent line. ## See also Index Objects · groupby · Reshaping with pivot and melt · loc vs iloc · DataFrame26:["$","$L1c","multiindex",{"item":{"slug":"multiindex","title":"MultiIndex","track":"pandas","type":"lesson","difficulty":"Beginner","order":10,"contentMarkdown":"$3b","starterCode":"","editorSeed":"","testCases":[],"mockFiles":{},"hints":[]},"index":9}] 3c:T19ff,# Reading and Writing Data > **In one line:** `read_csv` is the workhorse and it has ~50 parameters — the five that matter are `sep`, `dtype`, `comment`, `na_values` and `usecols`. ## read_csv, properly ```python df = pd.read_csv( "variants.tsv", sep="\t", comment="#", # skip VCF-style header lines dtype={"chrom": "category", "pos": "int32", "sample_id": "str"}, na_values=[".", "NA", "-", ""], # what counts as missing in THIS file usecols=["chrom", "pos", "ref", "alt", "qual"], nrows=1000, # peek at a huge file first low_memory=False, ) ``` **Always pass `dtype=` for genomic files.** Type inference gets these wrong constantly: - A chromosome column containing `1`–`22`, `X`, `Y`, `MT` becomes `str` only because of the letters. Read a chr1-only subset and it silently becomes `int64`, and your merges stop matching. - Sample IDs like `0012` lose their leading zeros if inferred as integers. - A column with one `"."` in it becomes `object` instead of numeric. **`comment="#"`** matters because VCF, GFF and SAM all use `#` header lines. ## Other formats ```python pd.read_csv(...) # also handles .gz, .bz2, .zip, .xz automatically pd.read_parquet("data.parquet") # fast, typed, compressed — the best intermediate pd.read_excel("samples.xlsx", sheet_name="metadata") pd.read_json("results.json", lines=True) pd.read_sql("SELECT * FROM variants", conn) pd.read_hdf("store.h5", key="expr") pd.read_fwf("fixed_width.txt", widths=[10, 8, 6]) pd.read_clipboard() # genuinely useful for one-off pastes ``` `read_csv` detects compression from the extension, so `pd.read_csv("counts.tsv.gz", sep="\t")` just works — no `gzip.open` needed. ## Parquet is the right intermediate format ```python df.to_parquet("intermediate.parquet") # preserves dtypes exactly df = pd.read_parquet("intermediate.parquet") ``` | | CSV/TSV | Parquet | |---|---|---| | Dtypes preserved | no — re-inferred every read | yes | | Size | baseline | ~5–10× smaller | | Read speed | baseline | ~10× faster | | Column subset read | must read all | reads only what you ask | | Human readable | yes | no | | Universal tool support | yes | good and improving | Use TSV at the boundaries (input from collaborators, output for the supplement) and Parquet for everything in between. Re-parsing a 2 GB TSV on every pipeline run is a self-inflicted wound. ## Large files ```python # chunked processing results = [] for chunk in pd.read_csv("huge.tsv", sep="\t", chunksize=100_000): results.append(chunk.query("qual > 30").groupby("chrom").size()) total = pd.concat(results).groupby(level=0).sum() # only the columns you need pd.read_csv("huge.tsv", sep="\t", usecols=["chrom", "pos", "af"]) # check the shape before committing pd.read_csv("huge.tsv", sep="\t", nrows=5) ``` If chunking is not enough, the answer is **polars** or **DuckDB**, both of which query files larger than RAM directly: ```python import duckdb duckdb.sql("SELECT chrom, count(*) FROM 'variants.tsv' GROUP BY chrom").df() ``` ## Writing ```python df.to_csv("out.tsv", sep="\t", index=False) # index=False almost always df.to_csv("out.tsv.gz", sep="\t", index=False) # compression from extension df.to_parquet("out.parquet") df.to_excel("out.xlsx", sheet_name="results", index=False) df.to_markdown() # for a README or a notebook df.to_latex() # for a manuscript table ``` **`index=False` unless the index is real data.** Otherwise you get a spurious unnamed column that reappears as `Unnamed: 0` on the next read, and it accumulates every round-trip. Round-tripping properly: ```python df.to_csv("out.tsv", sep="\t") # writes the index df = pd.read_csv("out.tsv", sep="\t", index_col=0) # reads it back ``` ## Bioinformatics parsing patterns **VCF** — skip the meta lines, find the header: ```python def read_vcf(path): import gzip opener = gzip.open if str(path).endswith(".gz") else open with opener(path, "rt") as fh: for line in fh: if line.startswith("#CHROM"): cols = line.lstrip("#").rstrip().split("\t") break return pd.read_csv(path, sep="\t", comment="#", names=cols, dtype={"CHROM": "category", "POS": "int32"}) ``` See Applied - Reading VCF with pandas for the INFO-field work. **GFF/GTF** — nine fixed columns, no header: ```python GFF_COLS = ["seqid", "source", "type", "start", "end", "score", "strand", "phase", "attributes"] gff = pd.read_csv(path, sep="\t", comment="#", names=GFF_COLS, dtype={"seqid": "category", "type": "category", "start": "int32", "end": "int32"}) ``` **BED** — 0-based, variable column count: ```python BED_COLS = ["chrom", "start", "end", "name", "score", "strand"] bed = pd.read_csv(path, sep="\t", header=None, names=BED_COLS[:n_cols], dtype={"chrom": "category"}) ``` **FASTA/FASTQ** — not tabular. Use SeqIO. ## Validation on load Worth building into any pipeline: ```python def load_variants(path): df = pd.read_csv(path, sep="\t", dtype=SCHEMA, na_values=["."]) assert not df.empty, f"{path} is empty" assert df["pos"].gt(0).all(), "non-positive positions" assert df["chrom"].isin(VALID_CHROMS).all(), \ f"unexpected: {set(df['chrom']) - VALID_CHROMS}" if df.duplicated(["chrom", "pos", "ref", "alt"]).any(): warnings.warn("duplicate variants present") return df ``` Five lines that catch the wrong file, a coordinate-system mix-up, and a chromosome-naming mismatch (`chr1` vs `1`) before they reach your results. That naming mismatch in particular — Ensembl uses `1`, UCSC uses `chr1` — silently produces empty joins. ## Common mistakes - **No `dtype=`** on genomic files. - **No `comment="#"`** on VCF/GFF, so header lines become data rows. - **Wrong `sep`.** The default is `,`; genomics is `\t`. Symptom: one giant column. - **`index=False` forgotten**, giving `Unnamed: 0` columns that multiply. - **Not handling the file's own missing-value convention.** `.` in VCF, `NA` in R output, `-` in some annotation files. - **Re-parsing large TSVs** instead of caching to Parquet. - **`chr1` vs `1`** across files → empty merges with no error. - **Reading a whole 10 GB file** to compute one summary. Chunk, or use DuckDB. ## See also DataFrame · pandas dtypes · Missing Data in pandas · Applied - Reading VCF with pandas · Applied - Genomic File Formats · NumPy File IO27:["$","$L1c","reading-and-writing-data",{"item":{"slug":"reading-and-writing-data","title":"Reading and Writing Data","track":"pandas","type":"lesson","difficulty":"Beginner","order":11,"contentMarkdown":"$3c","starterCode":"","editorSeed":"import pandas as pd\n\ndf = pd.read_csv(\n \"variants.tsv\",\n sep=\"\\t\",\n comment=\"#\", # skip VCF-style header lines\n dtype={\"chrom\": \"category\", \"pos\": \"int32\", \"sample_id\": \"str\"},\n na_values=[\".\", \"NA\", \"-\", \"\"], # what counts as missing in THIS file\n usecols=[\"chrom\", \"pos\", \"ref\", \"alt\", \"qual\"],\n nrows=1000, # peek at a huge file first\n low_memory=False,\n)\n","testCases":[],"mockFiles":{},"hints":[]},"index":10}] 3d:T17ff,# Reshaping with pivot and melt > **In one line:** `melt` goes wide → long (for plotting and grouping), `pivot` goes long → wide (for matrices and human reading). See Tidy Data for why the distinction exists. ## melt: wide → long ```python # wide: gene | sample_1 | sample_2 | sample_3 long = wide.melt( id_vars="gene", # columns to KEEP as identifiers var_name="sample", # name for the old column headers value_name="count", # name for the values ) # long: gene | sample | count ``` ```python wide.melt(id_vars=["gene", "chrom"]) # several id columns wide.melt(id_vars="gene", value_vars=["s1", "s2"]) # only some columns ``` Melt when: the column headers are actually *values* (sample names, timepoints, years), and you are about to plot, group, or run a statistical model. ## pivot: long → wide ```python wide = long.pivot(index="gene", columns="sample", values="count") ``` `pivot` requires the index/columns combination to be **unique**. Duplicates raise `ValueError: Index contains duplicate entries`. ## pivot_table: pivot with aggregation ```python wide = long.pivot_table( index="gene", columns="sample", values="count", aggfunc="mean", # how to combine duplicates fill_value=0, margins=True, # add row/column totals ) ``` `pivot_table` handles duplicates by aggregating them, which is both convenient and a trap: if you did not *expect* duplicates, it silently averages them instead of telling you something is wrong. Use `pivot` when the data should be unique — the error is information. `aggfunc="size"` gives you a contingency table: ```python variants.pivot_table(index="gene", columns="consequence", values="pos", aggfunc="size", fill_value=0) ``` `pd.crosstab` is the shorthand for exactly that: ```python pd.crosstab(variants["gene"], variants["consequence"]) pd.crosstab(variants["chrom"], variants["impact"], normalize="index") # row proportions ``` ## stack and unstack The MultiIndex-aware equivalents: ```python wide.stack() # columns → an index level (wide → long) long_mi.unstack() # innermost index level → columns long_mi.unstack(level="sample") long_mi.unstack(fill_value=0) ``` `stack` on an expression matrix gives you a Series with a (gene, sample) MultiIndex — often more convenient than melt when you want to keep the labels as index rather than columns. See MultiIndex. ## explode: one row per list element Not wide/long exactly, but the third reshaping tool you will need constantly in genomics, because bio formats pack lists into single cells. ```python df = pd.DataFrame({"gene": ["TP53"], "transcripts": [["NM_000546", "NM_001126112"]]}) df.explode("transcripts") # TP53 | NM_000546 # TP53 | NM_001126112 ``` Combined with a string split, this untangles delimited fields: ```python variants.assign(csq=variants["consequence"].str.split("&")).explode("csq") ``` A variant annotated `missense_variant&splice_region_variant` becomes two rows, one per consequence, which is what you need before counting consequence types. See String Accessor. ## The full VCF INFO untangling The canonical bio reshaping problem, using all three tools: ```python info = ( variants["INFO"] .str.split(";") # 'AF=0.3;DP=44' → ['AF=0.3', 'DP=44'] .explode() # one key=value per row .str.split("=", n=1, expand=True) # → two columns .rename(columns={0: "key", 1: "value"}) ) info["variant_idx"] = info.index wide_info = info.pivot_table(index="variant_idx", columns="key", values="value", aggfunc="first") variants = variants.join(wide_info) ``` See Applied - Reading VCF with pandas for the complete treatment including type coercion. ## Bioinformatics examples ```python # expression matrix → long, for seaborn long = expr.reset_index().melt(id_vars="gene", var_name="sample", value_name="count") long = long.merge(metadata, on="sample", how="left", validate="many_to_one") sns.boxplot(data=long, x="condition", y="count", hue="tissue") # long results → a wide matrix for clustering mat = long.pivot(index="gene", columns="sample", values="log_cpm") sns.clustermap(mat, z_score=0, cmap="vlag") # contingency: consequence type by chromosome pd.crosstab(variants["chrom"], variants["impact"], normalize="index") # genotype matrix from long calls gt = calls.pivot(index="variant_id", columns="sample_id", values="genotype") # per-sample-per-metric QC summary qc.pivot_table(index="sample", columns="metric", values="value") ``` ## Choosing a direction | Goal | Shape | |---|---| | seaborn plotting | **long** | | `groupby` on sample/condition | **long** | | Statistical modelling (formula-based) | **long** | | Matrix maths: PCA, correlation, clustering | **wide** | | Heatmap | **wide** | | Human reading / a supplementary table | **wide** | | Storage of a large matrix | **wide** (long repeats every label) | The practical rhythm: **store wide, melt at the point of plotting, discard the long frame.** ## Common mistakes - **`pivot` on non-unique index/column pairs** → ValueError. Use `pivot_table`, but first ask *why* there are duplicates. - **`pivot_table` silently averaging** unexpected duplicates. - **Melting without `id_vars`**, which melts your identifier columns too. - **Forgetting to `reset_index()` before melting** an index-carrying matrix — the gene names are in the index, not a column, so `id_vars="gene"` raises KeyError. - **Long-format memory blowup.** A 20,000 × 500 matrix is 10 million rows long, with the gene name repeated 500 times each. - **`explode` on a string instead of a list.** Split first; `explode` on a plain string does nothing. - **`fill_value=0` in a pivot** where the true meaning of absent is "not measured", not "zero". That distinction matters in proteomics and scRNA-seq. ## See also Tidy Data · groupby · String Accessor · MultiIndex · Applied - Reading VCF with pandas · Applied - Heatmaps and Clustermaps28:["$","$L1c","reshaping-with-pivot-and-melt",{"item":{"slug":"reshaping-with-pivot-and-melt","title":"Reshaping with pivot and melt","track":"pandas","type":"lesson","difficulty":"Beginner","order":12,"contentMarkdown":"$3d","starterCode":"","editorSeed":"import pandas as pd\n\ndf = pd.DataFrame({\"gene\": [\"TP53\"], \"transcripts\": [[\"NM_000546\", \"NM_001126112\"]]})\ndf.explode(\"transcripts\")\n# TP53 | NM_000546\n# TP53 | NM_001126112\n","testCases":[],"mockFiles":{},"hints":[]},"index":11}] 3e:T1559,# Series > **In one line:** a 1-D array plus an index of labels — and the labels participate in every operation, which is the whole difference from a NumPy array. ## Construction ```python s = pd.Series([120, 340, 90], index=["TP53", "BRCA1", "EGFR"], name="counts") s.values # underlying array (NumPy or Arrow) s.to_numpy() # preferred — always gives you a NumPy array s.index # Index(['TP53', 'BRCA1', 'EGFR'], dtype='str') s.dtype # int64 s.name # 'counts' — becomes the column name in a DataFrame ``` `.to_numpy()` over `.values`: `.values` returns whatever the backing store is, which for Arrow-backed or nullable dtypes is not a plain ndarray. ## Access: labels and positions ```python s["TP53"] # by label s[["TP53", "EGFR"]] # multiple labels s.loc["TP53"] # explicit label — use this s.iloc[0] # explicit position — use this s.get("MISSING", 0) # label lookup with a default ``` **Always use `.loc` and `.iloc` explicitly.** Bare `s[0]` is ambiguous: is `0` a label or a position? With an integer index it means label; otherwise position. That ambiguity is a genuine bug source. See loc vs iloc. ## Alignment — the thing that surprises everyone ```python a = pd.Series([1, 2, 3], index=["x", "y", "z"]) b = pd.Series([10, 20, 30], index=["z", "y", "x"]) # reversed order! a + b # x 31 ← matched by LABEL, not position # y 22 # z 13 ``` pandas aligned on the index. NumPy would have added positionally and given you `[11, 22, 33]`. The failure mode: ```python a = pd.Series([1, 2, 3], index=["x", "y", "z"]) b = pd.Series([10, 20, 30], index=[0, 1, 2]) a + b # 0 NaN # 1 NaN # ... all NaN, six rows, no error ``` **"Why is my result all NaN?" is almost always index misalignment.** The fix is to align deliberately (`b.index = a.index`) or drop out of label-land (`a + b.to_numpy()`). ## Vectorized operations Everything from ufuncs works, and the index is carried through: ```python np.log2(s + 1) s * 2 s > 100 s[s > 100] # boolean mask, index preserved s.clip(0, 1000) s.round(2) ``` ## The method catalogue ```python s.sum(); s.mean(); s.median(); s.std() # ddof=1 by default, unlike NumPy! s.min(); s.max(); s.idxmax(); s.idxmin() # idxmax gives the LABEL s.quantile([0.25, 0.5, 0.75]) s.describe() s.value_counts() # frequency table, sorted descending s.value_counts(dropna=False) # ← always use this; NaN counts matter s.value_counts(normalize=True) # proportions s.unique(); s.nunique() s.isna(); s.notna(); s.fillna(0); s.dropna() s.sort_values(); s.sort_index() s.rank(method="dense") s.astype("float32") s.map({"A": "Adenine"}) # elementwise via dict or function s.replace({-9: np.nan}) s.between(0.01, 0.99) s.isin(PANEL_GENES) s.head(); s.tail(); s.sample(5) ``` `s.idxmax()` returning the *label* rather than the position is the useful difference from `np.argmax` — "which gene is highest" answers directly. `s.std()` uses **ddof=1** (sample) while `np.std` uses **ddof=0** (population). Same data, different answer. See Axes and Reductions. ## Accessors Namespaced method groups for particular dtypes: ```python s.str.upper(); s.str.contains("BRCA") # string ops — see [[String Accessor]] s.dt.year; s.dt.dayofweek # datetime — see [[Datetime Handling in pandas]] s.cat.categories; s.cat.codes # categorical — see [[Categorical dtype]] ``` Calling `.str` on a non-string Series raises, which is a useful type check. ## Series vs NumPy array — when to use which | | Series | ndarray | |---|---|---| | Labels | yes | no | | Missing data | first-class | NaN only, floats only | | Mixed types | via object/Arrow | no | | Alignment | automatic | positional | | Speed | small overhead per op | baseline | | n-dimensional | no | yes | For a tight numeric loop over millions of elements, drop to `.to_numpy()`. For anything where losing the labels would be a correctness risk — which in bioinformatics is most things — stay in pandas. ## Bioinformatics examples ```python counts = pd.Series(raw, index=gene_ids, name="counts") cpm = counts / counts.sum() * 1e6 log_cpm = np.log2(cpm + 1) expressed = counts[counts > 10] counts.nlargest(20) # top 20 genes, labels intact counts.rank(ascending=False).loc["TP53"] # where does TP53 rank? df["consequence"].value_counts(dropna=False) df["chrom"].isin([f"chr{i}" for i in range(1, 23)]) # autosomes only ``` The point of `counts.nlargest(20)` over `np.argsort(raw)[-20:]` is that the gene names come along automatically. That is the whole value proposition of pandas in one line. ## Common mistakes - **Index misalignment producing all-NaN.** The single most common pandas confusion. - **Bare `s[0]` on an integer-indexed Series** meaning label, not position. - **Assigning a Series to a DataFrame column and getting NaN** — it aligned on index. Use `.to_numpy()` or `.reset_index(drop=True)` if you meant positional. - **`s.std()` vs `np.std(s)`** disagreeing because of `ddof`. - **`value_counts()` hiding NaN** by default. Pass `dropna=False`. - **`.values` when you wanted `.to_numpy()`.** - **Mutating a Series extracted from a DataFrame** and expecting the frame to change. Under Copy-on-Write it will not. ## See also DataFrame · Index Objects · loc vs iloc · String Accessor · Missing Data in pandas · pandas dtypes · ndarray29:["$","$L1c","series",{"item":{"slug":"series","title":"Series","track":"pandas","type":"lesson","difficulty":"Beginner","order":13,"contentMarkdown":"$3e","starterCode":"","editorSeed":"import pandas as pd\n\ns = pd.Series([120, 340, 90], index=[\"TP53\", \"BRCA1\", \"EGFR\"], name=\"counts\")\n\ns.values # underlying array (NumPy or Arrow)\ns.to_numpy() # preferred — always gives you a NumPy array\ns.index # Index(['TP53', 'BRCA1', 'EGFR'], dtype='str')\ns.dtype # int64\ns.name # 'counts' — becomes the column name in a DataFrame\n","testCases":[],"mockFiles":{},"hints":[]},"index":12}] 3f:T1614,# Sorting and Ranking > **In one line:** `sort_values` for data order, `sort_index` for label order, `rank` for positions within a distribution — and be careful how NaN and ties are handled. ## Sorting ```python df.sort_values("padj") df.sort_values("padj", ascending=False) df.sort_values(["chrom", "pos"]) # multi-key df.sort_values(["chrom", "pos"], ascending=[True, False]) df.sort_values("padj", na_position="first") # default is 'last' df.sort_values("padj", kind="stable") # preserve original order for ties df.sort_index() df.sort_index(axis=1) # sort the columns df.sort_values("gene", key=lambda s: s.str.lower()) # sort by a transformation ``` **`na_position="last"` is the default**, which is usually right for p-values (missing is not significant) and usually wrong for anything where missing should be visible. Decide deliberately. **`kind="stable"`** matters when you sort twice. An unstable sort can reorder tied rows arbitrarily, so results are not reproducible across pandas versions or platforms. If your output includes tied values and you care about determinism, use it — or add a tiebreaker column. ## Ranking ```python s.rank() # average rank for ties (default) s.rank(method="min") # competition ranking: 1, 2, 2, 4 s.rank(method="dense") # no gaps: 1, 2, 2, 3 s.rank(method="first") # break ties by original order s.rank(ascending=False) s.rank(pct=True) # percentile rank in [0, 1] s.rank(na_option="keep") # NaN stays NaN (default) ``` `method="dense"` is the one you usually want for a human-readable ranking. `method="average"` (the default) produces fractional ranks like 2.5, which are correct for rank-based statistics (Spearman, Wilcoxon) and odd-looking in a report. ## Ranking within groups ```python df["rank_in_gene"] = df.groupby("gene")["af"].rank(ascending=False) df["pct_in_sample"] = df.groupby("sample")["count"].rank(pct=True) ``` Per-group ranking is a `groupby().rank()`, not a transform — `rank` already returns one value per row. See groupby. ## Top N without a full sort ```python df.nlargest(20, "log2fc") df.nsmallest(20, "padj") df.nlargest(20, ["log2fc", "baseMean"]) # tiebreak on the second column df.nlargest(5, "x", keep="all") # include all ties ``` `nlargest` uses a partial selection algorithm — meaningfully faster than `sort_values().head(n)` on large frames, same as `np.argpartition` in NumPy. See Fancy Indexing. ## Genomic sort order The recurring problem: **chromosomes do not sort alphabetically.** ```python sorted(["chr1", "chr2", "chr10", "chrX"]) # ['chr1', 'chr10', 'chr2', 'chrX'] ← wrong, and it is wrong in every plot ``` The fix is an ordered categorical, set once: ```python CHROM_ORDER = [f"chr{i}" for i in range(1, 23)] + ["chrX", "chrY", "chrM"] df["chrom"] = pd.Categorical(df["chrom"], categories=CHROM_ORDER, ordered=True) df = df.sort_values(["chrom", "pos"]) ``` This fixes sorting, `groupby` output order, and — importantly — the x-axis order of every seaborn plot, all at once. Do it at load time. See Categorical dtype. Any value not in `CHROM_ORDER` becomes NaN, so decide what to do with scaffolds and alt contigs explicitly: ```python main = df[df["chrom"].isin(CHROM_ORDER)] print(f"dropped {len(df) - len(main):,} rows on non-primary contigs") ``` ## Sorting is not free ```python df.sort_values("pos") # O(n log n), copies the whole frame df.sort_values("pos", ignore_index=True) # skip index bookkeeping df.set_index("pos").sort_index() # pays once, then .loc is fast ``` If you slice a large frame by position repeatedly, sorting once into an index pays for itself. If you sort inside a loop, you have a problem. ## Bioinformatics examples ```python # DE results, most significant first, ties broken deterministically de.sort_values(["padj", "gene"], kind="stable") # canonical genomic order variants = variants.sort_values(["chrom", "pos", "ref", "alt"]) # top expressed genes per sample long.sort_values("count", ascending=False).groupby("sample").head(10) # best transcript per gene (lowest p, then longest) best = (variants .sort_values(["pvalue", "transcript_length"], ascending=[True, False]) .drop_duplicates("gene", keep="first")) # percentile of each gene's expression within its sample long["expr_pct"] = long.groupby("sample")["count"].rank(pct=True) # rank genes for GSEA (signed, by significance) de["gsea_rank"] = np.sign(de["log2FoldChange"]) * -np.log10(de["pvalue"]) ranked = de.sort_values("gsea_rank", ascending=False) ``` The `sort_values().groupby().head(n)` pattern in the third example is the idiomatic "top N per group" and is much faster than `groupby().apply(lambda g: g.nlargest(n, ...))`. ## Common mistakes - **Alphabetical chromosome order.** chr1, chr10, chr11... in a published figure. - **Not setting `kind="stable"`** when tied rows must stay in a fixed order. - **NaN position surprises.** Default is last. - **Sorting inside a loop.** - **Forgetting `reset_index(drop=True)`** after sorting, leaving the index scrambled relative to position. - **`rank()` default returning fractional ranks** in a report. - **Sorting a categorical without `ordered=True`** — you get the category order, which may not be what you set. - **Losing non-primary contigs silently** when applying a categorical chromosome order. ## See also Categorical dtype · groupby · Filtering and query · Index Objects · Fancy Indexing2a:["$","$L1c","sorting-and-ranking",{"item":{"slug":"sorting-and-ranking","title":"Sorting and Ranking","track":"pandas","type":"lesson","difficulty":"Beginner","order":14,"contentMarkdown":"$3f","starterCode":"","editorSeed":"sorted([\"chr1\", \"chr2\", \"chr10\", \"chrX\"])\n# ['chr1', 'chr10', 'chr2', 'chrX'] ← wrong, and it is wrong in every plot\n","testCases":[],"mockFiles":{},"hints":[]},"index":13}] 40:T19eb,# String Accessor > **In one line:** `.str` gives you vectorised string methods on a Series — and it is how you untangle every packed field in every bioinformatics file format. ## The basics ```python s.str.upper(); s.str.lower(); s.str.strip() s.str.len() s.str.startswith("chr"); s.str.endswith("_1") s.str.contains("BRCA") # regex by default s.str.contains("BRCA", regex=False) # literal, faster s.str.replace("chr", "", regex=False) s.str.removeprefix("chr"); s.str.removesuffix(".bam") s.str.split(":") # → Series of lists s.str.split(":", expand=True) # → DataFrame of columns s.str.split(":", n=1, expand=True) # split only on the first occurrence s.str.cat(sep=",") # join all values s.str[0:3] # slice each string s.str.zfill(2) # '1' → '01' s.str.pad(10, side="left") ``` `.str` only works on string dtype. Calling it on a numeric column raises, which is a free type check. Note that in pandas 3.0 string columns are `str` dtype (Arrow-backed) rather than `object`. `.str` works the same way but is considerably faster. See pandas dtypes. ## Extraction with regex ```python s.str.extract(r"(\w+):(\d+)") # first match → DataFrame of groups s.str.extract(r"(?P\w+):(?P\d+)") # named groups → named columns s.str.extractall(r"(\d+)") # ALL matches → MultiIndex s.str.findall(r"\d+") # all matches → Series of lists s.str.count("A") s.str.match(r"^chr\d+$") # anchored boolean test s.str.fullmatch(r"chr\d+") ``` `extract` with named groups is the cleanest way to parse a structured identifier: ```python loc = variants["locus"].str.extract(r"(?Pchr[\dXYM]+):(?P\d+)") loc["pos"] = loc["pos"].astype("int32") ``` **A non-matching row gives NaN, not an error.** Always check: ```python failed = loc["chrom"].isna() if failed.any(): print(f"{failed.sum()} unparseable: {variants.loc[failed, 'locus'].head().tolist()}") ``` ## split + expand: the workhorse ```python df[["chrom", "pos", "ref", "alt"]] = df["variant"].str.split("_", expand=True) ``` `expand=True` turns one column into several in one line. Two cautions: - If rows have *different* numbers of fields, the short ones get None and the column count is set by the longest row. - The results are all strings. Cast afterwards. ```python parts = df["variant"].str.split("_", expand=True) assert parts.shape[1] == 4, f"expected 4 fields, got {parts.shape[1]}" df[["chrom", "pos", "ref", "alt"]] = parts df["pos"] = df["pos"].astype("int32") ``` ## explode: one row per element ```python df.assign(csq=df["consequence"].str.split("&")).explode("csq") ``` A variant annotated `missense_variant&splice_region_variant` becomes two rows. This is how you get to a countable, groupable form. See Reshaping with pivot and melt. ## Bioinformatics parsing recipes ```python # strip Ensembl version suffixes so IDs will actually join df["gene_id"] = df["gene_id"].str.split(".").str[0] # normalise chromosome naming (UCSC ↔ Ensembl) df["chrom"] = df["chrom"].str.removeprefix("chr") df["chrom"] = "chr" + df["chrom"].astype(str) # parse a locus string df[["chrom", "pos"]] = df["locus"].str.split(":", expand=True) # parse a GFF attributes column gff["gene_id"] = gff["attributes"].str.extract(r'gene_id "([^"]+)"') gff["gene_name"] = gff["attributes"].str.extract(r'gene_name "([^"]+)"') # parse a VCF INFO field for one key df["dp"] = df["INFO"].str.extract(r"(?:^|;)DP=(\d+)").astype("Int64") df["af"] = df["INFO"].str.extract(r"(?:^|;)AF=([\d.eE+-]+)").astype("float32") # parse a sample FORMAT field df[["gt", "ad", "dp"]] = df["sample1"].str.split(":", n=2, expand=True) df[["ref_count", "alt_count"]] = df["ad"].str.split(",", expand=True).astype("int32") # HGVS-shaped strings — for FILTERING only, never for interpretation df["accession"] = df["hgvs"].str.split(":").str[0] df["is_coding"] = df["hgvs"].str.contains(":c.") # clean sample names from filenames df["sample"] = df["path"].str.extract(r"([^/]+)\.bam$") # structured sample IDs df[["patient", "tissue", "rep"]] = df["sample_id"].str.split("_", expand=True) ``` The `(?:^|;)` in the INFO extractions matters: without it, `DP=` also matches inside `MQ_DP=` or `ExcessHet;DP=`. Anchoring to a semicolon or the start of the string prevents that. **On the HGVS line:** string manipulation of HGVS is acceptable for coarse filtering (which accession, is it coding). It is **never** acceptable for extracting positions, comparing variants, or deciding equivalence — a regex cannot know that `c.76_78del` and `c.77_79del` are the same event. Use hgvs Parser. See hgvs Pitfalls. ## Missing values ```python s.str.upper() # NaN stays NaN, no error s.fillna("").str.contains("BRCA") # NaN → False s.str.contains("BRCA", na=False) # same, more direct ``` `.str.contains` returns NaN for missing input, and **NaN in a boolean mask raises** when you use it in `.loc`. Always pass `na=False` (or `na=True`) when the result feeds a filter. ## Performance `.str` methods are vectorised at the pandas level but still involve per-element Python string operations for many methods — they are much faster than `.apply(lambda x: ...)` but not as fast as numeric ufuncs. With Arrow-backed strings (the pandas 3.0 default) many operations are now compiled and substantially faster. For very heavy parsing of very large files, consider parsing at read time with a proper format library (`pysam`, `cyvcf2`) rather than regex over a DataFrame. ```python s.str.contains("x", regex=False) # skip the regex engine when you can ``` ## Common mistakes - **NaN propagating into a boolean mask.** Pass `na=False`. - **Forgetting `.str` is regex by default.** `s.str.replace(".", "")` deletes every character. Pass `regex=False` or escape it. - **`split(expand=True)` on ragged rows** producing surprise column counts. - **Not casting after splitting.** Everything comes back as a string, so `pos` sorts as `"1", "10", "100", "2"`. - **Unanchored INFO regexes** matching the wrong key. - **Silent NaN from a failed `extract`.** Always count them. - **Parsing HGVS with regex** for anything semantic. - **Using `.apply(lambda x: x.upper())`** instead of `.str.upper()`. ## See also pandas dtypes · Reshaping with pivot and melt · Applied - Reading VCF with pandas · Applied - GFF and Genomic Intervals · hgvs Parser · Series47:I[18876,["/_next/static/immutable/chunks/1dg34hwq0pupj.js","/_next/static/immutable/chunks/2h8so9g5itegs.js"],"IconMark"] 2b:["$","$L1c","string-accessor",{"item":{"slug":"string-accessor","title":"String Accessor","track":"pandas","type":"lesson","difficulty":"Beginner","order":15,"contentMarkdown":"$40","starterCode":"","editorSeed":"","testCases":[],"mockFiles":{},"hints":[]},"index":14}] 41:T1822,# Window Functions > **In one line:** rolling, expanding and shifted computations — for smoothing coverage, computing running totals, and comparing each row to its neighbours. ## Rolling ```python s.rolling(window=100).mean() s.rolling(window=100, center=True).mean() # centred rather than trailing s.rolling(window=100, min_periods=1).mean() # emit values at the edges s.rolling(window=100).agg(["mean", "std", "min", "max"]) s.rolling(window=100).quantile(0.5) s.rolling(window=100).apply(custom_fn, raw=True) # raw=True passes ndarray — faster ``` **The edge behaviour** is the thing to decide up front. By default the first `window-1` positions are NaN, because there is not a full window yet. `min_periods=1` fills them using whatever data exists, which is usually what you want for a coverage track and usually *not* what you want for a statistic you will interpret. **`center=True`** matters for genomics: an uncentred rolling mean shifts every feature to the right by half a window. On a coverage plot that is a visible, wrong offset. ## Window types ```python s.rolling(100).mean() # flat/rectangular s.rolling(100, win_type="gaussian").mean(std=20) # weighted s.rolling("10min").mean() # time-based (needs DatetimeIndex) s.ewm(span=20).mean() # exponentially weighted s.ewm(halflife=10).mean() ``` `ewm` (exponentially weighted) gives smoothing with no edge NaNs and a smooth response — often nicer than a rectangular rolling mean for visualisation. ## Expanding Cumulative, from the start to each point: ```python s.expanding().mean() # running mean s.expanding().max() # running max s.cumsum(); s.cummax(); s.cumprod() ``` ## Shift and diff ```python s.shift(1) # previous value s.shift(-1) # next value s.diff() # s - s.shift(1) s.diff(2) s.pct_change() df["gap"] = df["pos"] - df["pos"].shift(1) # distance to the previous variant df["is_next_to"] = df["gap"] < 10 ``` `shift` and `diff` are how you compare each row to its neighbours. In genomics that means inter-variant distances, exon gaps, and detecting runs of adjacent features. ## Grouped windows — the crucial bit **Rolling and shifting across group boundaries is a silent bug.** A rolling mean over a position-sorted variant table will happily average the end of chr1 with the start of chr2. ```python # WRONG — smooths across chromosome boundaries df["smooth"] = df["depth"].rolling(100, center=True).mean() # RIGHT df = df.sort_values(["chrom", "pos"]) df["smooth"] = ( df.groupby("chrom")["depth"] .transform(lambda s: s.rolling(100, center=True, min_periods=1).mean()) ) # same for shift/diff df["gap"] = df.groupby("chrom")["pos"].diff() # NaN at each chromosome start ✓ ``` The `groupby().transform()` form is the pattern. See groupby. ## Rolling on a MultiIndex or with a time index ```python df.groupby("sample").rolling(window=10).mean() # produces a MultiIndex df.set_index("time").groupby("sample").rolling("1h").mean() ``` Time-based windows (`"1h"`, `"7D"`) require a sorted DatetimeIndex and are genuinely different from count-based windows when the sampling is irregular. See Datetime Handling in pandas. ## Bioinformatics examples ```python # smooth a per-base coverage track for plotting cov = cov.sort_values(["chrom", "pos"]) cov["smoothed"] = ( cov.groupby("chrom")["depth"] .transform(lambda s: s.rolling(1000, center=True, min_periods=1).mean()) ) # GC content in sliding windows gc = is_gc.rolling(100, center=True).mean() # detect coverage dropouts: sustained low regions low = cov["smoothed"] < 0.2 * cov["smoothed"].median() runs = (low != low.shift()).cumsum() # ← run-length encoding idiom dropouts = cov[low].groupby(runs).agg( start=("pos", "min"), end=("pos", "max"), n=("pos", "size") ).query("n >= 50") # distance to the nearest previous variant, per chromosome variants["dist_prev"] = variants.groupby("chrom")["pos"].diff() # rolling allele-frequency deviation, for LOH detection baf["rolling_dev"] = ( baf.groupby("chrom")["baf"] .transform(lambda s: (s - 0.5).abs().rolling(50, center=True).mean()) ) # running total along a chromosome cov["cumulative"] = cov.groupby("chrom")["depth"].cumsum() ``` The run-length encoding idiom in the dropout example is worth memorising: ```python runs = (flag != flag.shift()).cumsum() ``` Each time the flag changes, the cumulative sum increments, so consecutive identical values share a run id. Grouping by it gives you contiguous blocks — the basis of segment detection, CNV calling, and any "find stretches where X is true" problem. ## Performance ```python s.rolling(100).mean() # compiled — fast s.rolling(100).apply(np.mean) # Python per window — very slow s.rolling(100).apply(f, raw=True) # raw ndarray instead of Series — faster s.rolling(100, engine="numba").mean() # for large data with a custom function ``` Prefer the built-in aggregations. If you need `apply`, always pass `raw=True` — it skips constructing a Series per window. For very large arrays with a simple window, NumPy's `sliding_window_view` is faster still and allocates nothing. See Memory Layout and Strides. ## Common mistakes - **Rolling across chromosome/sample boundaries.** The defining error of this note. - **Not sorting before rolling.** Window functions assume order and will not warn. - **Forgetting `center=True`**, shifting every feature by half a window. - **Leading NaNs** breaking a downstream step. Decide `min_periods` deliberately. - **`apply` without `raw=True`.** - **`shift` inside a group without `groupby`** — the first row of each group picks up the last row of the previous group. - **Interpreting a smoothed track as data.** Smoothing changes the effective resolution; do not call peaks on a smoothed signal without accounting for it. ## See also groupby · Series · Datetime Handling in pandas · Split-Apply-Combine · Applied - Quality Control Plots · Vectorization