From 10x output to TwINFER input

Convert 10x output (matrix.mtx, features.tsv, barcodes.tsv of one or more samples) and a lineage-barcode / clone table into the input of infer_with_twinfer for real data: one row per cell with

clone_id, cell_id, time_step, {gene}_mRNA, {gene}_mRNA, ...

and normalised expression (log1p of counts per 10,000, with each cell’s total over all genes). The infer_with_twinfer function uses its clones scenario (clones of any size, singletons allowed). Score one time point at a time: infer_with_twinfer(table[table.time_step == t], genes) (see the end of this notebook and the real-data tutorial).

The samples, genes and clone table to use depends on the dataset; the cells to edit are marked edit me. The paths in them are placeholders.

What this notebook does not do: fuzzy / Hamming-distance barcode correction, doublet detection or other clone calling. Clones are assigned by exact barcode match (section 3). If your data need more careful clone calling, do it separately and give this notebook a table with a resolved clone_id column, which section 3 supports directly.

import re
from collections import defaultdict
from pathlib import Path

import numpy as np
import pandas as pd
import scipy.io
import scipy.sparse

from twinfer.utils.paths import get_repo_root

1. Find the samples and give each a time point

The next two cells find one sub-folder per sample and read each sample’s time point from the folder name.

Set samples_root to the folder that holds one sub-folder per sample. A sub-folder counts as a sample when it contains exactly one each of matrix.mtx, features.tsv (or genes.tsv) and barcodes.tsv, gzipped or not; other sub-folders are skipped.
Optional: set sample_name_include to text that every sample folder name contains (for example "LSK") to skip unrelated folders. None keeps all of them.
Set timepoint_regex: a regular expression whose first group is the time point, a whole number. The default _d(\d+)_ reads 2 from LSK_d2_1 and 4 from LSK_d4_2_3. It needs an underscore after the number, so it finds nothing in LSK_d6. Set it to None if the folder names do not contain the time point.
Run the cell and check the table: one row per sample folder found, with the time point read from its name (empty where the regular expression found none).
Correct the table in timepoint_overrides (next cell): list only the samples whose time point is missing or wrong, by folder name, for example {"LSK_d6": 6}. With timepoint_regex = None, list every sample. The cell stops with an error naming any sample that is still without a time point.

# ==================== edit me ====================
samples_root = "/path/to/10x_samples"     # one sub-folder per sample with matrix / features / barcodes

# Only auto-discover folders whose name contains this substring (None = all
# subfolders of samples_root that contain a 10x triplet). Useful when
# samples_root has unrelated folders mixed in.
sample_name_include = None

# Regex whose first capture group is the timepoint, matched against the
# folder name (here: "..._d2_..." -> 2). None = don't try to auto-parse; fill in timepoint_overrides
# for every sample instead.
timepoint_regex = re.compile(r"_d(\d+)_")
# ===================================================


def _find_one(folder: Path, patterns):
    hits = sorted({h for p in patterns for h in folder.glob(p)})
    return hits[0] if len(hits) == 1 else None


def _find_10x_trio(folder: Path):
    mtx = _find_one(folder, ["matrix.mtx", "matrix.mtx.gz"])
    feats = _find_one(folder, ["features.tsv", "features.tsv.gz", "genes.tsv", "genes.tsv.gz"])
    bcs = _find_one(folder, ["barcodes.tsv", "barcodes.tsv.gz"])
    return mtx, feats, bcs


discovered = {}
for sub in sorted(Path(samples_root).iterdir()):
    if not sub.is_dir():
        continue
    if sample_name_include and sample_name_include not in sub.name:
        continue
    mtx, feats, bcs = _find_10x_trio(sub)
    if not (mtx and feats and bcs):
        continue

    auto_time_step = None
    if timepoint_regex is not None:
        m = timepoint_regex.search(sub.name)
        if m:
            auto_time_step = int(m.group(1))

    discovered[sub.name] = {"dir": sub, "time_step": auto_time_step}

print(f"discovered {len(discovered)} sample folder(s) under {samples_root}")
pd.DataFrame(
    {name: {"time_step (auto)": cfg["time_step"]} for name, cfg in discovered.items()}
).T

# ==================== edit me ====================
# Fix anything the regex above got wrong or couldn't parse, or add samples
# by hand. Only samples that need a correction have to be listed here.
timepoint_overrides = {
    # "some_sample_name": 2,
}
# ===================================================

samples = {}
for name, cfg in discovered.items():
    t = timepoint_overrides.get(name, cfg["time_step"])
    samples[name] = {**cfg, "time_step": t}

missing_t = [name for name, cfg in samples.items() if cfg["time_step"] is None]
if missing_t:
    raise ValueError(
        f"No timepoint resolved for: {missing_t}. Add them to timepoint_overrides above."
    )

print(f"{len(samples)} sample(s) ready:")
pd.DataFrame({name: {"time_step": cfg["time_step"]} for name, cfg in samples.items()}).T

Singlet lists (optional, per sample)

If a sample has a doublet-filtering / singlet whitelist (a text file with one cell barcode per line), give its path here. Samples left out of singlet_lists (or explicitly set to None) are not singlet-filtered – all cells in the matrix are kept.

# ==================== edit me ====================
# sample_name -> path to a singlet whitelist, or None / omitted if the
# sample has no such list and should not be singlet-filtered.
singlet_lists = {
    # "sample_A": "/path/to/sample_A_singlets.txt",
}
# ===================================================

unknown_samples = set(singlet_lists) - set(samples)
if unknown_samples:
    raise ValueError(f"singlet_lists has sample(s) not in samples: {sorted(unknown_samples)}")

for name in samples:
    samples[name]["singlet_list"] = singlet_lists.get(name)

pd.DataFrame(
    {name: {"singlet_list": cfg["singlet_list"] or "(none)"} for name, cfg in samples.items()}
).T

2. Genes of interest

A list of genes, or None to keep every gene shared by all samples (the full transcriptome: large). A requested gene missing from any sample stops the notebook (section 4). The library size used for normalisation is always computed over all genes. The example list is a set of mouse haematopoietic regulators; the right list depends on the dataset.

# ==================== edit me ====================
requested_genes = [
    "Gata1", "Gata2", "Gfi1", "Fli1", "Spi1", "Tal1",
    "Cebpa", "Jun", "Egr1", "Nab2", "Klf1", "Zfpm1",
]  # or requested_genes = None to keep every gene common to all samples
# ===================================================

3. Clone assignment (lineage-barcode / clone sheet – required)

TwINFER needs real clones, so a clone sheet is required here – there is no “no info -> treat every cell as its own clone” fallback since that would not have any twin information to analyze.

Two supported sheet shapes, set via clone_sheet_barcode_col / clone_sheet_clone_id_col (exactly one of the two should be set):

Lineage-barcode sheet (clone_sheet_barcode_col set): clones are assigned by exact string match on that column without any additional checks like Hamming-distance correction and subset-compatible merging. If a cell has more than one distinct barcode on record, the most-supported one (by row count) is kept and a count of how often this happened is printed.
Already-resolved clone sheet (clone_sheet_clone_id_col set): used directly, no matching logic at all.

# ==================== edit me ====================
clone_sheet_path = "/path/to/clone_table.csv"   # required -- no clone info, no run

clone_sheet_cell_id_col = "cellID"   # column matching barcodes.tsv values (the 10x cell barcode)
clone_sheet_sample_col = "sample"    # column matching the samples keys from Section 1

clone_sheet_barcode_col = "barcode"  # lineage barcode -> exact-match clone assignment
clone_sheet_clone_id_col = None      # OR set this (and set clone_sheet_barcode_col = None) if the sheet already has resolved clone_id values
# ===================================================


def build_clone_lookup(sheet_path, cell_id_col, sample_col, barcode_col, clone_id_col):
    if sheet_path is None:
        raise ValueError(
            "clone_sheet_path is required -- TwINFER needs real clone/lineage "
            "information, there is no valid 'no sheet' fallback. Set it to a "
            "sheet with a lineage-barcode column (clone_sheet_barcode_col) or "
            "an already-resolved clone_id column (clone_sheet_clone_id_col)."
        )

    if (barcode_col is None) == (clone_id_col is None):
        raise ValueError("Set exactly one of clone_sheet_barcode_col / clone_sheet_clone_id_col")

    id_col = clone_id_col if clone_id_col else barcode_col
    sheet = pd.read_csv(sheet_path, usecols=[sample_col, cell_id_col, id_col])

    unmatched_samples = set(samples) - set(sheet[sample_col].unique())
    if unmatched_samples:
        print(
            f"WARNING: {len(unmatched_samples)} sample(s) don't appear in the clone "
            f"sheet's '{sample_col}' column (no clone_id will be assigned for them): "
            f"{sorted(unmatched_samples)}"
        )

    counts = (
        sheet.groupby([sample_col, cell_id_col, id_col])
        .size()
        .reset_index(name="n_rows")
        .sort_values("n_rows", ascending=False)
    )
    n_ambiguous = counts.groupby([sample_col, cell_id_col]).size().gt(1).sum()
    if n_ambiguous:
        label = "clone_id" if clone_id_col else "lineage barcode"
        print(
            f"NOTE: {n_ambiguous} cell(s) had >1 distinct {label} on record; kept the "
            f"most-supported one (exact match only)."
        )

    dedup = counts.drop_duplicates(subset=[sample_col, cell_id_col], keep="first")
    return {
        (row[sample_col], row[cell_id_col]): row[id_col]
        for _, row in dedup.iterrows()
    }


clone_lookup = build_clone_lookup(
    clone_sheet_path,
    clone_sheet_cell_id_col,
    clone_sheet_sample_col,
    clone_sheet_barcode_col,
    clone_sheet_clone_id_col,
)
print(f"clone_lookup has entries for {len(clone_lookup)} (sample, cell_id) pairs")

4. Load each sample’s 10x matrix, subset to the genes of interest

features.tsv is assumed to be the standard 10x layout (id, symbol, feature_type); gene matching is done on the symbol column.

_sample_cache = {}


def get_sample_raw(sample_name):
    """Load+cache one sample's matrix (cells x genes), gene symbols, cell barcodes."""
    if sample_name in _sample_cache:
        return _sample_cache[sample_name]

    sample_dir = samples[sample_name]["dir"]
    mtx_path, feats_path, bcs_path = _find_10x_trio(sample_dir)

    mat = scipy.io.mmread(mtx_path).tocsr()
    features = pd.read_csv(feats_path, header=None, sep="\t")
    barcodes = pd.read_csv(bcs_path, header=None, sep="\t")

    gene_symbols = (features.iloc[:, 1] if features.shape[1] > 1 else features.iloc[:, 0]).astype(str).to_numpy()
    cell_barcodes = barcodes.iloc[:, 0].astype(str).to_numpy()

    if mat.shape == (len(gene_symbols), len(cell_barcodes)):
        mat = mat.T.tocsr()  # genes x cells -> cells x genes
    elif mat.shape != (len(cell_barcodes), len(gene_symbols)):
        raise ValueError(
            f"[{sample_name}] matrix shape {mat.shape} matches neither "
            f"({len(cell_barcodes)} cells, {len(gene_symbols)} genes) nor its transpose"
        )

    _sample_cache[sample_name] = (mat, gene_symbols, cell_barcodes)
    return _sample_cache[sample_name]


# ---- resolve the gene list against every sample up front ----
sample_gene_sets = {}
for name in samples:
    _, gene_symbols, _ = get_sample_raw(name)
    sample_gene_sets[name] = set(gene_symbols)

genes_in_all_samples = set.intersection(*sample_gene_sets.values())

if requested_genes is None:
    resolved_genes = sorted(genes_in_all_samples)
else:
    missing = defaultdict(list)
    for g in requested_genes:
        if g not in genes_in_all_samples:
            for name, gs in sample_gene_sets.items():
                if g not in gs:
                    missing[g].append(name)
    if missing:
        lines = [f"  {g}: missing from {names}" for g, names in missing.items()]
        raise ValueError(
            "Requested gene(s) missing from at least one sample:\n" + "\n".join(lines)
        )
    resolved_genes = list(requested_genes)

print(f"using {len(resolved_genes)} gene(s): {resolved_genes}")

5. Build the per-cell table (one sample at a time)

Expression is normalised as log1p(counts / total counts of the cell × 10,000), with the total over all genes of the cell (computed before the genes are selected). Set normalise = False only to export raw counts. Cell ids are written as {sample}::{10x barcode}, since 10x barcodes recur across samples and lanes.

normalise = True   # log1p counts per 10,000 (library size over all genes)


def build_sample_table(sample_name, gene_list, clone_lookup):
    cfg = samples[sample_name]
    mat, gene_symbols, cell_barcodes = get_sample_raw(sample_name)

    gene_index = {g: i for i, g in enumerate(gene_symbols)}
    col_idx = [gene_index[g] for g in gene_list]
    sub = mat[:, col_idx]
    if scipy.sparse.issparse(sub):
        sub = sub.toarray()
    sub = np.asarray(sub, dtype=float)
    if normalise:
        totals = np.asarray(mat.sum(axis=1)).ravel()          # library size over all genes, before the gene subset
        sub = np.log1p(sub / np.maximum(totals, 1)[:, None] * 1e4)

    df = pd.DataFrame(sub, columns=[f"{g}_mRNA" for g in gene_list])
    df.insert(0, "time_step", cfg["time_step"])
    df.insert(0, "cell_id", [f"{sample_name}::{bc}" for bc in cell_barcodes])
    df["_raw_barcode"] = cell_barcodes

    if cfg["singlet_list"] is not None:
        singlets = set(pd.read_csv(cfg["singlet_list"], header=None)[0].astype(str))
        n_before = len(df)
        df = df[df["_raw_barcode"].isin(singlets)].copy()
        print(f"[{sample_name}] singlet filter: {n_before} -> {len(df)} cells")

    df["clone_id"] = [clone_lookup.get((sample_name, bc)) for bc in df["_raw_barcode"]]
    n_missing = df["clone_id"].isna().sum()
    if n_missing:
        print(f"[{sample_name}] {n_missing}/{len(df)} cell(s) have no clone_id match -> dropped")
    df = df.dropna(subset=["clone_id"]).copy()

    gene_cols = [f"{g}_mRNA" for g in gene_list]
    return df[["clone_id", "cell_id", "time_step"] + gene_cols]


tables = [build_sample_table(name, resolved_genes, clone_lookup) for name in samples]
combined_df = pd.concat(tables, ignore_index=True)
print(f"\ncombined: {len(combined_df)} cells, {combined_df['clone_id'].nunique()} clones, "
      f"timepoints {sorted(combined_df['time_step'].unique())}")
combined_df.head()

6. Check the table

Required columns, no missing values, no duplicated (cell_id, time_step) rows, and at least one gene column.

def validate_canonical_input(df):
    required = ["clone_id", "cell_id", "time_step"]

    missing = [c for c in required if c not in df.columns]
    if missing:
        raise ValueError(f"missing required column(s): {missing}")

    gene_cols = [c for c in df.columns if c.endswith("_mRNA")]
    if not gene_cols:
        raise ValueError("must contain at least one '*_mRNA' column")

    forbidden = [c for c in ("twin_id", "pair_id", "n_twins", "n_pairs") if c in df.columns]
    if forbidden:
        raise ValueError(f"remove twin-expansion column(s) before saving: {forbidden}")

    check_cols = required + gene_cols
    if df[check_cols].isna().any().any():
        bad = [c for c in check_cols if df[c].isna().any()]
        raise ValueError(f"missing values in: {bad}")

    duplicated = df.duplicated(subset=["cell_id", "time_step"], keep=False)
    if duplicated.any():
        examples = df.loc[duplicated, ["cell_id", "time_step"]].drop_duplicates().head(5).to_dict("records")
        raise ValueError(f"duplicate (cell_id, time_step) rows, e.g.: {examples}")

    print(f"OK: {len(df)} rows, {len(gene_cols)} gene column(s), no duplicates, no missing values")


validate_canonical_input(combined_df)

7. Save

# ==================== edit me ====================
output_path = get_repo_root() / "tutorials" / "tutorial_output" / "convert_10x_to_twinfer_input" / "converted_twinfer_input.csv"   # next to this notebook
# ===================================================

output_path.parent.mkdir(parents=True, exist_ok=True)
combined_df.to_csv(output_path, index=False)
print(f"wrote {output_path}")

8. Score the gene pairs

One time point at a time (for the analysis of the result, see the real-data tutorial):

from twinfer.inference.covariance_decomposition_z_scores import infer_with_twinfer

gene_columns = [c for c in combined_df.columns if c.endswith("_mRNA")]
for t in sorted(combined_df["time_step"].unique()):
    cells = combined_df[combined_df["time_step"] == t]
    result = infer_with_twinfer(cells[["clone_id"] + gene_columns], gene_columns)
    print(f"time point {t}: {len(cells)} cells, {len(result)} gene pairs, scenarios",
          result["inference_scenario"].value_counts().sort_index().to_dict())