Issue with load_data dataloader: KeyError 'IR_VJ_1_v_call'

[JiahaoChow]

Hello, thank you for sharing your work on MIST. I am very interested in your project, and I am currently trying to reproduce your algorithm. However, I encountered some issues related to the dataloader.

I created the environment with the following command:

conda create -n mist python=3.10 -y && conda activate mist && pip install git+https://github.com/aapupu/MIST.git

Then I tried running the following code:

from mist.data import load_data

tcr_paths = [
    './data/raw/CD8+Tcells/donor1/vdj_v1_hs_aggregated_donor1_all_contig_annotations.csv', 
    './data/raw/CD8+Tcells/donor2/vdj_v1_hs_aggregated_donor2_all_contig_annotations.csv', 
    './data/raw/CD8+Tcells/donor3/vdj_v1_hs_aggregated_donor3_all_contig_annotations.csv', 
    './data/raw/CD8+Tcells/donor4/vdj_v1_hs_aggregated_donor4_all_contig_annotations.csv'
]

rna_paths = [
    './data/raw/CD8+Tcells/donor1/vdj_v1_hs_aggregated_donor1_filtered_feature_bc_matrix.h5',
    './data/raw/CD8+Tcells/donor2/vdj_v1_hs_aggregated_donor2_filtered_feature_bc_matrix.h5',
    './data/raw/CD8+Tcells/donor3/vdj_v1_hs_aggregated_donor3_filtered_feature_bc_matrix.h5',
    './data/raw/CD8+Tcells/donor4/vdj_v1_hs_aggregated_donor4_filtered_feature_bc_matrix.h5'
]

batches = ['donor1', 'donor2', 'donor3', 'donor4']

adata, *dataloader_tuple = load_data(
    rna_path=rna_paths, 
    tcr_path=tcr_paths, 
    batch=batches,
    rna_data_type='h5',
    tcr_data_type='10X',
)

However, I received the following error:

KeyError: 'IR_VJ_1_v_call'

The traceback points to this line in mist/data.py:

adata_tcr.obs['IR_VJ_1_v_call'] = adata_tcr.obs['IR_VJ_1_v_call'].astype('str')

It seems that the resultant tcr_adata has 0 features and no obs entries.
Do you think this could be due to a version mismatch (e.g., scirpy or related dependencies), or is there something else I might be missing?

Any guidance would be greatly appreciated.

Thank you!

[aapupu]

Hello, Python 3.8 and scirpy 0.12.0 are recommended for MIST. Because the data structure processed by later versions of scirpy has changed. We welcome your further comments and feedback.

[JiahaoChow]

Hello, thank you very much for your timely reply and for clarifying the recommended environment (Python 3.8 and scirpy 0.12.0).

I re-ran the experiment following your suggestion, but unfortunately I still encounter an error. The raw TCR data seems to be empty:

Raw adata_tcr shape: (193656, 0)
Filtering TCR
Traceback (most recent call last):
  ...
    adata_tcr = adata_tcr[adata_tcr.obs['IR_VDJ_1_v_call'].isin(TCR_dict['TRBV'].keys()) &
TypeError: 'NoneType' object is not subscriptable

It looks like adata_tcr ends up with 0 features, so when process_scTCR tries to filter based on IR_VDJ_1_v_call, it fails.

For reproducibility, I obtained the raw data directly from 10x Genomics using the following commands:

curl -O https://cf.10xgenomics.com/samples/cell-vdj/3.0.2/vdj_v1_hs_aggregated_donor1/vdj_v1_hs_aggregated_donor1_all_contig_annotations.csv
curl -O https://cf.10xgenomics.com/samples/cell-vdj/3.0.2/vdj_v1_hs_aggregated_donor2/vdj_v1_hs_aggregated_donor2_all_contig_annotations.csv
curl -O https://cf.10xgenomics.com/samples/cell-vdj/3.0.2/vdj_v1_hs_aggregated_donor3/vdj_v1_hs_aggregated_donor3_all_contig_annotations.csv
curl -O https://cf.10xgenomics.com/samples/cell-vdj/3.0.2/vdj_v1_hs_aggregated_donor4/vdj_v1_hs_aggregated_donor4_all_contig_annotations.csv

curl -O https://cf.10xgenomics.com/samples/cell-vdj/3.0.2/vdj_v1_hs_aggregated_donor1/vdj_v1_hs_aggregated_donor1_filtered_feature_bc_matrix.h5
curl -O https://cf.10xgenomics.com/samples/cell-vdj/3.0.2/vdj_v1_hs_aggregated_donor2/vdj_v1_hs_aggregated_donor2_filtered_feature_bc_matrix.h5
curl -O https://cf.10xgenomics.com/samples/cell-vdj/3.0.2/vdj_v1_hs_aggregated_donor3/vdj_v1_hs_aggregated_donor3_filtered_feature_bc_matrix.h5
curl -O https://cf.10xgenomics.com/samples/cell-vdj/3.0.2/vdj_v1_hs_aggregated_donor4/vdj_v1_hs_aggregated_donor4_filtered_feature_bc_matrix.h5

Do you think this issue might be related to how these Cell Ranger (3.0.2) outputs are parsed, or perhaps differences in the expected column names (e.g., IR_VDJ_1_v_call) between Cell Ranger versions?

Any advice on how to correctly prepare the TCR data or which Cell Ranger/scirpy versions were used in your original experiments would be very helpful.

Thanks again for your guidance!

[aapupu]

adata_tcr without features seems normal, as the information may already be stored in adata_tcr.obs. You can try loading adata_tcr with the Scirpy package to check whether adata_tcr.obs contains the 'IR_VDJ_1_v_call' column. If not, it may be due to a version incompatibility. In addition, please make sure that the TCR file itself is correct.

def load_scTCR(path, batch, type='10X'):
    assert type in ['10X', 'tracer', 'BD', 'h5ad'], 'type of scTCR-seq data must in 10X, tracer, BD, or h5ad.'
    adatas_tcr = {}
    for sample, sample_meta in zip(path, batch):
        if type == '10X':
            adata_tcr = ir.io.read_10x_vdj(sample)
        elif type == 'tracer':
            adata_tcr = ir.io.read_tracer(sample)
        elif type == 'BD':
            adata_tcr = ir.io.read_bd_rhapsody(sample)
        elif type == 'h5ad':
            adata_tcr = ir.io.read_h5ad(sample)
        adata_tcr.obs['batch'] = sample_meta
        adatas_tcr[sample_meta] = adata_tcr

    adata_tcr = anndata.concat(adatas_tcr, index_unique="_")
    return adata_tcr

[JiahaoChow]

Thanks a lot for your timely reply. In the current environment, reading a single TCR file individually works fine. I will further check why it fails specifically inside load_data. Appreciate your help!

[MaybeTommorow]

I encountered the same issue. In scirpy v0.22.5, AIRR data is stored in adata.obsm['airr'] instead of adata.obs . The following function provides a simple fix by transferring the data to adata.obs , which satisfies MIST's input requirements.

"""
Fix TCR h5ad files to add standard AIRR columns expected by MIST.

MIST expects columns like:
  IR_VJ_1_v_call, IR_VJ_1_j_call, IR_VJ_1_junction_aa
  IR_VDJ_1_v_call, IR_VDJ_1_d_call, IR_VDJ_1_j_call, IR_VDJ_1_junction_aa

This script extracts them from scirpy's awkward array (obsm['airr']).
"""

import sys
import scanpy as sc
import awkward as ak
import pandas as pd
import numpy as np

def fix_tcr_h5ad(input_path, output_path=None):

    if output_path is None:
        output_path = input_path

    print(f"Reading {input_path} ...")
    adata = sc.read_h5ad(input_path)

    if 'airr' not in adata.obsm:
        raise ValueError(f"No 'airr' key found in obsm. Is this a scirpy TCR object?")

    airr = adata.obsm['airr']
    n = len(adata)
    results = []

    for i in range(n):
        cell_chains = airr[i]
        num_chains = len(cell_chains)

        tra = None
        trb = None

        for j in range(num_chains):
            locus = str(cell_chains[j]['locus'])
            if locus == 'TRA' and tra is None:
                tra = cell_chains[j]
            elif locus == 'TRB' and trb is None:
                trb = cell_chains[j]

        def get_field(record, field):
            if record is None:
                return None
            val = record[field]
            # Handle awkward None / numpy None
            if val is None or (isinstance(val, (float, np.floating)) and np.isnan(val)):
                return None
            return str(val)

        results.append({
            'IR_VJ_1_v_call': get_field(tra, 'v_call'),
            'IR_VJ_1_j_call': get_field(tra, 'j_call'),
            'IR_VJ_1_junction_aa': get_field(tra, 'junction_aa'),
            'IR_VDJ_1_v_call': get_field(trb, 'v_call'),
            'IR_VDJ_1_d_call': get_field(trb, 'd_call'),
            'IR_VDJ_1_j_call': get_field(trb, 'j_call'),
            'IR_VDJ_1_junction_aa': get_field(trb, 'junction_aa'),
        })

    df = pd.DataFrame(results, index=adata.obs_names)

    # Add columns to obs
    for col in df.columns:
        adata.obs[col] = df[col]

    # Ensure chain_pairing exists (MIST/scirpy uses it)
    if 'chain_pairing' not in adata.obs.columns:
        # Infer from available chains
        def infer_pairing(row):
            has_tra = pd.notna(row['IR_VJ_1_v_call'])
            has_trb = pd.notna(row['IR_VDJ_1_v_call'])
            if has_tra and has_trb:
                return 'single pair'
            elif has_tra or has_trb:
                return 'orphan VJ' if has_tra else 'orphan VDJ'
            else:
                return 'no IR'
        adata.obs['chain_pairing'] = df.apply(infer_pairing, axis=1)

    print(f"Saving to {output_path} ...")
    adata.write_h5ad(output_path)
    print("Done!")
    print("Added columns:", df.columns.tolist())

[aapupu]

Hi MaybeTommorow, I have fixed this issue. Please check it.