Pandas06/21
Starting sandbox…
Beginner·4 graded tests

Filtering and Query

Filtering and query

Core Principle: Data filtering in pandas is executed via boolean series indexing with .loc[] or declarative string expressions via .query(), requiring bitwise logical operators (&, |, ~) and explicit parentheses around conditions.

Boolean masks

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:

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

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

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

# 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

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

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

# 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

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-Write

Exercise

Write report_filter(df, max_af, chroms) for a variant table with columns chrom, pos, gene, af and filter. Apply three filters in order and report each one, so that a step which removes almost everything is visible rather than inferred from the final count.

Apply them in this order: keep rows whose filter equals PASS, then rows whose chrom is in chroms, then rows whose allele frequency is below max_af. A missing af means the variant is absent from the reference population and must be kept by the third filter, not dropped by it. Return the lines PASS <kept>/<entering>, chrom <kept>/<entering> and rare <kept>/<entering>, where entering is the number of rows the step received, followed by genes=<names> listing the surviving genes as a comma-separated, alphabetically sorted list of distinct names, or genes=none when nothing survives.

exerciseEvery submission runs independently in a fresh namespace.

Test cases · 4

#viainputexpected stdout
1stdin0.01 chr17,chr7 chrom pos gene af filter chr17 7676154 TP53 0.002 PASS chr17 43093220 BRCA1 PASS chr7 55191822 EGFR 0.32 PASS chr12 25245350 KRAS 0.001 PASS chr7 140753336 BRAF 0.004 LowQual chr17 7675088 TP53 0.0005 PASS PASS 5/6 chrom 4/5 rare 3/4 genes=BRCA1,TP53
2stdin0.5 chr1,chr2 chrom pos gene af filter chr17 7676154 TP53 0.002 PASS chr7 55191822 EGFR 0.32 LowQual PASS 1/2 chrom 0/1 rare 0/0 genes=none
3stdin0.05 chr17 chrom pos gene af filter chr17 7676154 TP53 PASS chr17 7675088 TP53 PASS chr17 43093220 BRCA1 PASS PASS 3/3 chrom 3/3 rare 3/3 genes=BRCA1,TP53
4stdin0.01 chr17,chr7,chr12 chrom pos gene af filter chr17 7676154 TP53 0.9 PASS PASS 1/1 chrom 1/1 rare 0/1 genes=none

Hints · 2

01Hint

Series.isin(list) builds the chromosome mask; Series.fillna(0) makes an absent frequency compare as rare.

02Hint

Bind each mask to a name before applying it — mask.sum() then gives you the kept count for the report line.