Filtering

We now move on to filtering out BCR contigs (and corresponding cells if necessary) from the BCR data and transcriptome object loaded in scanpy.

[1]:
import os
import scanpy as sc
import pandas as pd
import polars as pl
import dandelion as ddl

# change directory to holding tutorial data
os.chdir("dandelion_tutorial")

Prepare the gene expression data

[2]:
samples = [
    "sc5p_v2_hs_PBMC_1k_b",
    "sc5p_v2_hs_PBMC_10k_b",
    "vdj_v1_hs_pbmc3_b",
    "vdj_nextgem_hs_pbmc3_b",
]
adata_list = []
for sample in samples:
    adata = sc.read_10x_h5(
        sample + "/filtered_feature_bc_matrix.h5", gex_only=True
    )
    adata.obs["sample_id"] = sample
    # rename cells to sample id + barcode
    adata.obs_names = [str(sample) + "_" + str(j) for j in adata.obs_names]
    adata.var_names_make_unique()
    adata_list.append(adata)
adata = sc.concat(adata_list)
# rename the obs_names again, this time cleaving the trailing -#
adata.obs_names = [str(j).split("-")[0] for j in adata.obs_names]
adata
/opt/homebrew/Caskroom/miniforge/base/envs/dandelion/lib/python3.12/site-packages/anndata/_core/anndata.py:1880: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
/opt/homebrew/Caskroom/miniforge/base/envs/dandelion/lib/python3.12/site-packages/anndata/_core/anndata.py:1880: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
/opt/homebrew/Caskroom/miniforge/base/envs/dandelion/lib/python3.12/site-packages/anndata/_core/anndata.py:1880: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
/opt/homebrew/Caskroom/miniforge/base/envs/dandelion/lib/python3.12/site-packages/anndata/_core/anndata.py:1880: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
/opt/homebrew/Caskroom/miniforge/base/envs/dandelion/lib/python3.12/site-packages/anndata/_core/anndata.py:1880: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
/opt/homebrew/Caskroom/miniforge/base/envs/dandelion/lib/python3.12/site-packages/anndata/_core/anndata.py:1880: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
/opt/homebrew/Caskroom/miniforge/base/envs/dandelion/lib/python3.12/site-packages/anndata/_core/anndata.py:1880: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
/opt/homebrew/Caskroom/miniforge/base/envs/dandelion/lib/python3.12/site-packages/anndata/_core/anndata.py:1880: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
[2]:
AnnData object with n_obs × n_vars = 30471 × 31915
    obs: 'sample_id'

I’m using a wrapper called pp.recipe_scanpy_qc to run through a generic scanpy workflow. You can skip this if you already have a pre-processed AnnData object for the subsequent steps.

[3]:
ddl.pp.recipe_scanpy_qc(adata, mito_cutoff=None)  # use a gmm model to decide
# we can continue with those that survive qc
adata = adata[adata.obs["filter_rna"] == "False"].copy()
adata
OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
/opt/homebrew/Caskroom/miniforge/base/envs/dandelion/lib/python3.12/functools.py:912: UserWarning: zero-centering a sparse array/matrix densifies it.
/Users/uqztuong/Documents/GitHub/dandelion/src/dandelion/external/scanpy/__init__.py:142: FutureWarning: In the future, the default backend for leiden will be igraph instead of leidenalg.

 To achieve the future defaults please pass: flavor="igraph" and n_iterations=2.  directed must also be False to work with igraph's implementation.
/opt/homebrew/Caskroom/miniforge/base/envs/dandelion/lib/python3.12/site-packages/numpy/_core/fromnumeric.py:3860: RuntimeWarning: Mean of empty slice.
/opt/homebrew/Caskroom/miniforge/base/envs/dandelion/lib/python3.12/site-packages/numpy/_core/_methods.py:144: RuntimeWarning: invalid value encountered in scalar divide
[3]:
AnnData object with n_obs × n_vars = 25057 × 31915
    obs: 'sample_id', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'gmm_pct_count_clusters_keep', 'scrublet_score', 'is_doublet', 'filter_rna'

Filter cells that are potental doublets and poor quality in both the V(D)J data and transcriptome data

Polars backend update

ddl.pp.check_contigs has been reimplemented on a Polars backend (dandelion.polars). The function signature and behavior remain the same, but the internals have changed significantly:

  • Input: now accepts DandelionPolars | pl.DataFrame | pl.LazyFrame | str.

  • Lazy evaluation: data is converted to a LazyFrame immediately upon loading. Filtering steps (productive status, library type, locus) are composed into a single query plan and only materialized when the actual computation begins, reducing memory overhead on large datasets.

  • Vectorized marking: the core marking logic (mark_ambiguous_contigs_vec) replaces the original row-wise iteration with Polars window functions (over(["cell_id", "locus_type"])), performing UMI/consensus dominance tests, ranking, and flag assignment in bulk. The same applies to duplicate resolution (resolve_duplicates) and chimeric gene detection (check_chimeric_genes_vec), both of which now use pure Polars expressions with no loops or apply. This yields a significant speedup on large contig tables.

  • Output: returns DandelionPolars. To save the contig table, use vdj.data.write_csv("file.tsv", separator="\t").

We use ddl.pp.check_contigs to mark and filter out cells and contigs from both the V(D)J data and transcriptome data in AnnData. The operation will remove bad quality cells based on transcriptome information as well as remove V(D)J doublets (multiplet heavy/long chains, and/or light/short chains) from the V(D)J data. In some situations, a single cell can have multiple heavy/long and light/short chain contigs although they have an identical V(D)J+C alignment; in situations like this, the contigs with lesser UMIs will be dropped and the UMIs transferred to umi_count column. The same procedure is applied to both chains before further checks of the annotation quality, UMI and consensus count distributions.

Cells in the gene expression object without V(D)J information will not be affected which means that the AnnData object can hold non-B/T cells.

[4]:
# first we read in the 4 bcr files
bcr_files = []
for sample in samples:
    file_location = sample + "/dandelion/filtered_contig_dandelion.tsv"
    bcr_files.append(pl.read_csv(file_location, separator="\t"))
bcr = pl.concat(bcr_files)
bcr
[4]:
shape: (7_385, 120)
sequence_idsequencerev_compproductivev_calld_callj_callsequence_alignmentgermline_alignmentjunctionjunction_aav_cigard_cigarj_cigarstop_codonvj_in_framelocusjunction_lengthnp1_lengthnp2_lengthv_sequence_startv_sequence_endv_germline_startv_germline_endd_sequence_startd_sequence_endd_germline_startd_germline_endj_sequence_startj_sequence_endj_germline_startj_germline_endv_scorev_identityv_supportd_scored_identityd_score_blastnd_sequence_alignment_blastnd_germline_alignment_blastnd_call_igblastnd_sourced_support_igblastnd_score_igblastnv_call_genotypedgermline_alignment_d_masksample_idc_callc_sequence_alignmentc_germline_alignmentc_sequence_startc_sequence_endc_scorec_identityc_call_10xjunction_aa_lengthfwr1_aafwr2_aafwr3_aafwr4_aacdr1_aacdr2_aacdr3_aasequence_alignment_aav_sequence_alignment_aad_sequence_alignment_aaj_sequence_alignment_aacomplete_vdjj_call_multimappersj_call_multiplicityj_call_sequence_start_multimappersj_call_sequence_end_multimappersj_call_support_multimappersmu_count
strstrstrstrstrstrstrstrstrstrstrstrstrstrstrstrstri64i64i64i64i64i64i64i64i64i64i64i64i64i64i64f64f64f64f64f64f64strstrstrstrf64f64strstrstrstrstrstri64i64f64f64stri64strstrstrstrstrstrstrstrstrstrstrstrstrf64strstrstrf64
"sc5p_v2_hs_PBMC_1k_b_AACTCCCAG…"ACTGCGGGGGTAAGAGGTTGTGTCCACCAT…"F""T""IGLV5-45*03"null"IGLJ3*02""CAGGCTGTGCTGACTCAGCCGTCTTCC...…"CAGGCTGTGCTGACTCAGCCGTCTTCC...…"TGTATGATTTGGCACAGCAGCGCTTGGGTG…"CMIWHSSAWVF""85S310="null"395S2N36=""F""T""IGL"330null863951334nullnullnullnull396431338484.01.01.2700e-138nullnullnullnullnullnullnullnullnull"IGLV5-45*03""CAGGCTGTGCTGACTCAGCCGTCTTCC...…"sc5p_v2_hs_PBMC_1k_b""IGLC""GTCAGCCCAAGGCTGCCCCCTCGGTCACTC…"GTCAGCCCAAGGCTGCCCCCTCGGTCACTC…432642390.0100.0"IGLC2"11"QAVLTQPSSLSASPGASASLTCTLR""IYWYQQKPGSPPQYLLR""QQGSGVPSRFSGSKDASANAGILLISGLQS…"FGGGTKLTVL""SGINVGTYR""YKSDSDK""MIWHSSAWV""QAVLTQPSSLSASPGASASLTCTLRSGINV…"QAVLTQPSSLSASPGASASLTCTLRSGINV…null"VFGGGTKLTVL"null"["IGLJ3*01"]"1.0"[397]""[431]""[5.29e-13]"0.0
"sc5p_v2_hs_PBMC_1k_b_AACTCCCAG…"ATACTTTCTGAGAGTCCTGGACCTCCTGTG…"F""T""IGHV4-61*02""IGHD3-3*01""IGHJ6*02""CAGGTGCAGCTGCAGGAGTCGGGCCCA...…"CAGGTGCAGCTGCAGGAGTCGGGCCCA...…"TGTGCGAGAGAAAATTACGATTTTTGGAGT…"CARENYDFWSGYYHGADVW""94S299=""395S2N23=""418S11N3=1X5=2X41=""F""T""IGH"57209539313203964183254194701263467.01.01.5300e-13344.91.047.3"ATTACGATTTTTGGAGTGGTTATTA""ATTACGATTTTTGGAGTGGTTATTA""IGHD3-3*01"null1.1100e-844.9"IGHV4-61*02""CAGGTGCAGCTGCAGGAGTCGGGCCCA...…"sc5p_v2_hs_PBMC_1k_b""IGHM""GGAGTGCATCCGCCCCAACCCTTTTCCCCC…"GGAGTGCATCCGCCCCAACCCTTTTCCCCC…471541132.0100.0"IGHM"19"QVQLQESGPGLVKPSQTLSLTCTVS""WSWIRQPAGKGLEWIGR""NYNPSLKSRVTISVDTSKNQFSLKLSSVTA…"WGQGTTVTVSS""GGSISSGSYY""IYTSGST""ARENYDFWSGYYHGADV""QVQLQESGPGLVKPSQTLSLTCTVSGGSIS…"QVQLQESGPGLVKPSQTLSLTCTVSGGSIS…"YDFWSGY""YHGADVWGQGTTVTVSS"null"["IGHJ6*02"]"1.0"[416]""[470]""[3.41e-19]"3.0
"sc5p_v2_hs_PBMC_1k_b_AACTCCCAG…"GGCTGGGGTCTCAGGAGGCAGCGCTCTGGG…"F""F""IGLV2-5*01"null"IGLJ3*02""CAGTCTGCCCTGATTCAGCCTCCCTCC...…"CAGTCTGCCCTGATTCAGCCTCCCTCC...…"TGCTGCTCATATACAAGCAGTGCCACTTTC…"CCSYTSSATFLGV""98S295="null"395S38=""T""F""IGL"402null993931337nullnullnullnull396433138461.01.01.3800e-131nullnullnullnullnullnullnullnullnull"IGLV2-5*01""CAGTCTGCCCTGATTCAGCCTCCCTCC...…"sc5p_v2_hs_PBMC_1k_b""IGLC3""GTCAGCCCAAGGCTGCCCCCTCGGTCACTC…"GTCAGCCCAAGGCTGCCCCCTCGGTCACTC…434644390.0100.0"IGLC3"13"QSALIQPPSVSGSPGQSVTISCTGT""VSWYQQHPGTVPKPMIY""TQPSGVPDRFSGSKSGNTASMTISGLQAED…"RRRDQADRP""SSDVGSYDY""NVN""CSYTSSATFLG""QSALIQPPSVSGSPGQSVTISCTGTSSDVG…"QSALIQPPSVSGSPGQSVTISCTGTSSDVG…null"LGVRRRDQADRP"null"["IGLJ3*02"]"1.0"[396]""[433]""[2.45e-16]"0.0
"sc5p_v2_hs_PBMC_1k_b_AACTCTTGT…"AGCTCTGAGAGAGGAGCCTTAGCCCTGGAT…"F""T""IGHV3-21*01""IGHD3-22*01""IGHJ4*02""GAGGTGCAGCTGGTGGAGTCTGGGGGA...…"GAGGTGCAGCTGGTGGAGTCTGGGGGA...…"TGTGCGAGACGTTACTATGATAGTAGTGGT…"CARRYYDSSGYSANFDYW""136S294=""432S3N23=""461S3N45=""F""T""IGH"54261374301318433455426462506448459.01.03.6200e-13144.91.043.6"TTACTATGATAGTAGTGGTTATT""TTACTATGATAGTAGTGGTTATT""IGHD3-22*01"null1.1900e-844.9"IGHV3-21*01""GAGGTGCAGCTGGTGGAGTCTGGGGGA...…"sc5p_v2_hs_PBMC_1k_b""IGHM""GGAGTGCATCCGCCCCAACCCTTTTCCCCC…"GGAGTGCATCCGCCCCAACCCTTTTCCCCC…507577132.0100.0"IGHM"18"EVQLVESGGGLVKPGGSLRLSCAAS""MNWVRQAPGKGLEWVSS""YYADSVKGRFTISRDNAKNSLYLQMNSLRA…"WGQGTLVTVSS""GFTFSSYS""ISSSSSYI""ARRYYDSSGYSANFDY""EVQLVESGGGLVKPGGSLRLSCAASGFTFS…"EVQLVESGGGLVKPGGSLRLSCAASGFTFS…"YYDSSGY""FDYWGQGTLVTVSS"null"["IGHJ4*02"]"1.0"[462]""[506]""[2.81e-20]"0.0
"sc5p_v2_hs_PBMC_1k_b_AACTCTTGT…"AGAGCTCTGGGGAGTCTGCACCATGGCTTG…"F""T""IGLV4-69*01"null"IGLJ1*01""CAGCTTGTGCTGACTCAATCGCCCTCT...…"CAGCTTGTGCTGACTCAATCGCCCTCT...…"TGTCAGACCTGGGGCACTGGCATTTATGTC…"CQTWGTGIYVF""82S297="null"379S1N37=""F""T""IGL"330null833791333nullnullnullnull380416238464.01.01.5500e-132nullnullnullnullnullnullnullnullnull"IGLV4-69*01""CAGCTTGTGCTGACTCAATCGCCCTCT...…"sc5p_v2_hs_PBMC_1k_b""IGLC1""GTCAGCCCAAGGCCAACCCCACTGTCACTC…"GTCAGCCCAAGGCCAACCCCACTGTCACTC…417627390.0100.0"IGLC1"11"QLVLTQSPSASASLGASVKLTCTLS""IAWHQQQPEKGPRYLMK""SKGDGIPDRFSGSSSGAERYLTISSLQSED…"FGTGTKVTVL""SGHSSYA""LNSDGSH""QTWGTGIYV""QLVLTQSPSASASLGASVKLTCTLSSGHSS…"QLVLTQSPSASASLGASVKLTCTLSSGHSS…null"YVFGTGTKVTVL"null"["IGLJ1*01"]"1.0"[379]""[416]""[2.39e-16]"0.0
"vdj_nextgem_hs_pbmc3_b_TTTGCGC…"ATCACATAACAACCACATTCCTCCTCTAAA…"F""T""IGHV1-69*01,IGHV1-69D*01""IGHD3-22*01""IGHJ4*02""CAGGTGCAGCTGGTGCAGTCTGGGGCT...…"CAGGTGCAGCTGGTGCAGTCTGGGGCT...…"TGTGCGAGGGGGAAGTATTACTATGATAAA…"CARGKYYYDKSGSPPPIYSFDYW""115S29=1X40=1X11=1X6=1X11=1X2=…"414S14=2X5=""456S4N44=""F""T""IGH"696211164081317415435121457500548408.00.945391.2700e-11529.50.9047628.8"GTATTACTATGATAAAAGTGG""GTATTACTATGATAGTAGTGG""IGHD3-22*01"null0.0005529.5"IGHV1-69*01,IGHV1-69D*01""CAGGTGCAGCTGGTGCAGTCTGGGGCT...…"vdj_nextgem_hs_pbmc3_b""IGHG2""CCTCCACCAAGGGCCCATCGGTCTTCCCCC…"CCTCCACCAAGGGCCCATCGGTCTTCCCCC…501672318.0100.0"IGHG2"23"QVQLVQSGAEVKKPGSSVKVSCKVS""ISWVRQAPGQGLEWMAG""NYAQKFQGRVSITADESTTTAYMELSSLRS…"WGQGTLVTVSS""GGIFSSYA""IIPIFGAT""ARGKYYYDKSGSPPPIYSFDY""QVQLVQSGAEVKKPGSSVKVSCKVSGGIFS…"QVQLVQSGAEVKKPGSSVKVSCKVSGGIFS…"YYYDKSG""FDYWGQGTLVTVSS"null"["IGHJ4*02"]"1.0"[457]""[500]""[1.18e-19]"16.0
"vdj_nextgem_hs_pbmc3_b_TTTGGTT…"AGAGCTCTGGAGAAGAGCTGCTCAGTTAGG…"F""T""IGKV3-20*01"null"IGKJ2*01,IGKJ2*02""GAAATTGTGTTGACGCAGTCTCCAGGCACC…"GAAATTGTGTTGACGCAGTCTCCAGGCACC…"TGTCAGCAGTATGATGAGTCACCTCTGACT…"CQQYDESPLTF""104S84=1X3=1X2=1X4=1X1=1X17=1X…null"395S5N34=""F""T""IGK"332null1053931334nullnullnullnull396429639414.00.958481.4100e-117nullnullnullnullnullnullnullnullnull"IGKV3-20*01""GAAATTGTGTTGACGCAGTCTCCAGGCACC…"vdj_nextgem_hs_pbmc3_b""IGKC""CGAACTGTGGCTGCACCATCTGTCTTCATC…"CGAACTGTGGCTGCACCATCTGTCTTCATC…429565254.0100.0"IGKC"11"EIVLTQSPGTLSLSPGERATLSCRAS""LAWYQQKPGQAPRLLIY""SRATGIPDRFSGSGSGTDFTLTISRLVPED…"FGQGTKLEIK""QSLTNSQ""GAS""QQYDESPLT""EIVLTQSPGTLSLSPGERATLSCRASQSLT…"EIVLTQSPGTLSLSPGERATLSCRASQSLT…null"TFGQGTKLEIK"null"["IGKJ2*02"]"1.0"[396]""[429]""[3.59e-14]"11.0
"vdj_nextgem_hs_pbmc3_b_TTTGGTT…"AGCTCTGGGAGAGGAGCCCCAGCCCTGAGA…"F""T""IGHV3-9*01""IGHD5-18*01,IGHD5-5*01""IGHJ6*03""GAAGTGCAGCTGGTGGAGTCTGGGGGA...…"GAAGTGCAGCTGGTGGAGTCTGGGGGA...…"TGTGCAAAAGACGGATACAGCTATCGTTCG…"CAKDGYSYRSSYYFYMDVW""137S59=1X22=1X14=1X45=1X3=1X19…"434S2N12=1X3=""455S8N7=2X46=""F""T""IGH"57151384331320435450318456510963438.00.972971.5500e-12425.70.9375nullnullnull"IGHD5-18*01,IGHD5-5*01"null0.00825.7"IGHV3-9*01""GAAGTGCAGCTGGTGGAGTCTGGGGGA...…"vdj_nextgem_hs_pbmc3_b""IGHA1""CATCCCCGACCAGCCCCAAGGTCTTCCCGC…"CATCCCCGACCAGCCCCAAGGTCTTCCCGC…511670296.0100.0"IGHA1"19"EVQLVESGGGLVQPGRSLRLSCAAS""MHWVRQAPGKGLEWVSS""GYADSVKGRFTISRDNAKNSLYLQMNSLRA…"WGKGTTVTVSS""GFSFDDYV""ISWNSGRT""AKDGYSYRSSYYFYMDV""EVQLVESGGGLVQPGRSLRLSCAASGFSFD…"EVQLVESGGGLVQPGRSLRLSCAASGFSFD…"GYSYR""YYFYMDVWGKGTTVTVSS"null"["IGHJ6*03"]"1.0"[456]""[510]""[1.96e-22]"10.0
"vdj_nextgem_hs_pbmc3_b_TTTGTCA…"AGCTCTGAGAGAGGAGCCCAGCCCTGGGAT…"F""T""IGHV3-23*01,IGHV3-23D*01""IGHD4-17*01""IGHJ4*02""GAGGTGCAGCTGTTGGAGTCTGGGGGA...…"GAGGTGCAGCTGTTGGAGTCTGGGGGA...…"TGTGCGAAAGATTTTAGGTCGCCATACGGT…"CAKDFRSPYGDYYFDYW""136S296=""445S4N10=""455S48=""F""T""IGH"511301374321320446455514456503148463.01.04.1500e-13219.91.0nullnullnull"IGHD4-17*01"null0.3919.9"IGHV3-23*01,IGHV3-23D*01""GAGGTGCAGCTGTTGGAGTCTGGGGGA...…"vdj_nextgem_hs_pbmc3_b""IGHM""GGAGTGCATCCGCCCCAACCCTTTTCCCCC…"GGAGTGCATCCGCCCCAACCCTTTTCCCCC…504574132.0100.0"IGHM"17"EVQLLESGGGLVQPGGSLRLSCAAS""MSWVRQAPGKGLEWVSA""YYADSVKGRFTISRDNSKNTLYLQMNSLRA…"WGQGTLVTVSS""GFTFSSYA""ISGSGGST""AKDFRSPYGDYYFDY""EVQLLESGGGLVQPGGSLRLSCAASGFTFS…"EVQLLESGGGLVQPGGSLRLSCAASGFTFS…"YGD""YFDYWGQGTLVTVSS"null"["IGHJ4*02"]"1.0"[456]""[503]""[6.02e-22]"0.0
"vdj_nextgem_hs_pbmc3_b_TTTGTCA…"GTGGGTCCAGGAGGCAGAACTCTGGGTGTC…"F""T""IGLV3-25*03"null"IGLJ1*01""TCCTATGAGCTGACACAGCCACCCTCG...…"TCCTATGAGCTGACACAGCCACCCTCG...…"TGTCAATCAGCAGACAGCAGTGGTACTTAT…"CQSADSSGTYLYVF""92S289="null"382S38=""F""T""IGL"421null933811340nullnullnullnull383420138452.01.05.9500e-129nullnullnullnullnullnullnullnullnull"IGLV3-25*03""TCCTATGAGCTGACACAGCCACCCTCG...…"vdj_nextgem_hs_pbmc3_b"nullnullnullnullnullnullnull"None"14"SYELTQPPSVSVSPGQTARITCSGD""AYWYQQKPGQAPVLVIY""ERPSGIPERFSGSSSGTTVTLTISGVQAED…"FGTGTKVTVL""ALPKQY""KDS""QSADSSGTYLYV""SYELTQPPSVSVSPGQTARITCSGDALPKQ…"SYELTQPPSVSVSPGQTARITCSGDALPKQ…null"YVFGTGTKVTVL"null"["IGLJ1*01"]"1.0"[383]""[420]""[1.64e-16]"0.0

Annotating V/D/J gene functionality

A new preprocessing function now allows you to annotate the V/D/J gene functionality for each contig in the V(D)J data. The function will parse the v_call, d_call, and j_call columns and extract the gene name and functionality (functional (F), open reading frame (ORF; not considered a confirmed functional gene) and pesudogene (P)) from each call. The resulting annotations will be stored in new columns: <call>_functionality e.g. v_call_functionality. This information maye be useful before downstream analyses, such as filtering out these contigs or identifying specific gene usage patterns.

[5]:
# first, convert bcr to a Dandelion object
bcr = ddl.Dandelion(bcr)
bcr
[5]:
Lazy Dandelion object with n_obs = 3158 and n_contigs = 7385
    data: sequence_id, sequence, rev_comp, productive, v_call, d_call, j_call, sequence_alignment, germline_alignment, junction, junction_aa, v_cigar, d_cigar, j_cigar, stop_codon, vj_in_frame, locus, junction_length, np1_length, np2_length, v_sequence_start, v_sequence_end, v_germline_start, v_germline_end, d_sequence_start, d_sequence_end, d_germline_start, d_germline_end, j_sequence_start, j_sequence_end, j_germline_start, j_germline_end, v_score, v_identity, v_support, d_score, d_identity, d_support, j_score, j_identity, j_support, fwr1, fwr2, fwr3, fwr4, cdr1, cdr2, cdr3, cell_id, consensus_count, umi_count, v_call_10x, d_call_10x, j_call_10x, junction_10x, junction_10x_aa, j_call_blastn, j_identity_blastn, j_alignment_length_blastn, j_number_of_mismatches_blastn, j_number_of_gap_openings_blastn, j_sequence_start_blastn, j_sequence_end_blastn, j_germline_start_blastn, j_germline_end_blastn, j_support_blastn, j_score_blastn, j_sequence_alignment_blastn, j_germline_alignment_blastn, j_call_igblastn, j_source, j_support_igblastn, j_score_igblastn, d_call_blastn, d_identity_blastn, d_alignment_length_blastn, d_number_of_mismatches_blastn, d_number_of_gap_openings_blastn, d_sequence_start_blastn, d_sequence_end_blastn, d_germline_start_blastn, d_germline_end_blastn, d_support_blastn, d_score_blastn, d_sequence_alignment_blastn, d_germline_alignment_blastn, d_call_igblastn, d_source, d_support_igblastn, d_score_igblastn, v_call_genotyped, germline_alignment_d_mask, sample_id, c_call, c_sequence_alignment, c_germline_alignment, c_sequence_start, c_sequence_end, c_score, c_identity, c_call_10x, junction_aa_length, fwr1_aa, fwr2_aa, fwr3_aa, fwr4_aa, cdr1_aa, cdr2_aa, cdr3_aa, sequence_alignment_aa, v_sequence_alignment_aa, d_sequence_alignment_aa, j_sequence_alignment_aa, complete_vdj, j_call_multimappers, j_call_multiplicity, j_call_sequence_start_multimappers, j_call_sequence_end_multimappers, j_call_support_multimappers, mu_count, rearrangement_status
    metadata: cell_id, sample_id, productive_VDJ, productive_VJ, d_call_VDJ, j_call_VDJ, j_call_VJ, junction_VDJ, junction_VJ, junction_aa_VDJ, junction_aa_VJ, locus_VDJ, locus_VJ, v_call_VDJ, v_call_VJ, c_call_VDJ, c_call_VJ, umi_count_VDJ, umi_count_VJ, productive_VDJ_main, productive_VJ_main, d_call_VDJ_main, j_call_VDJ_main, j_call_VJ_main, junction_VDJ_main, junction_VJ_main, junction_aa_VDJ_main, junction_aa_VJ_main, locus_VDJ_main, locus_VJ_main, v_call_genotyped_VDJ_main, v_call_genotyped_VJ_main, c_call_VDJ_main, c_call_VJ_main, umi_count_VDJ_main, umi_count_VJ_main, isotype, isotype_main, isotype_status, locus_status, chain_status, rearrangement_status_VDJ, rearrangement_status_VJ
[6]:
germline_db = (
    "/Users/uqztuong/Documents/GitHub/dandelion/container/database/germlines"
)
ddl.pp.annotate_functionality(bcr, germline=germline_db, org="human", db="imgt")
bcr.data[
    [
        "v_call",
        "v_call_functionality",
        "d_call",
        "d_call_functionality",
        "j_call",
        "j_call_functionality",
    ]
].collect()
[6]:
shape: (7_385, 6)
v_callv_call_functionalityd_calld_call_functionalityj_callj_call_functionality
strstrstrstrstrstr
"IGLV5-45*03""F"nullnull"IGLJ3*02""F"
"IGHV4-61*02""F""IGHD3-3*01""F""IGHJ6*02""F"
"IGLV2-5*01""P"nullnull"IGLJ3*02""F"
"IGHV3-21*01""F""IGHD3-22*01""F""IGHJ4*02""F"
"IGLV4-69*01""F"nullnull"IGLJ1*01""F"
"IGHV1-69*01,IGHV1-69D*01""F""IGHD3-22*01""F""IGHJ4*02""F"
"IGKV3-20*01""F"nullnull"IGKJ2*01,IGKJ2*02""F"
"IGHV3-9*01""F""IGHD5-18*01,IGHD5-5*01""F""IGHJ6*03""F"
"IGHV3-23*01,IGHV3-23D*01""F""IGHD4-17*01""F""IGHJ4*02""F"
"IGLV3-25*03""F"nullnull"IGLJ1*01""F"

Library type

It is recommended to specify the library_type argument as it will remove all contigs that do not belong to the related loci. The rationale is that the choice of the library type should mean that the primers used would most likely amplify those related sequences and if there’s any unexpected loci, they likely represent artifacts and shouldn’t be analysed. The optional argument accepts: ig, tr-ab, tr-gd or None where None means all contigs will be kept.

The main output of this function are two an additional columns in vdj.data, extra and ambiguous, which flags T or F for contigs that were marked accordingly. The rules for marking contigs are as follows:

extra is marked as T if the contig passes the internal QC filters based on umi_count (or consensus_count if there are ties in the umi_count) in a cell. If you are only interested in just the top contig pair, you can set filter_extra=True to remove the extra contigs.

For VDJ chains, the current rule set is to keep the top 1 productive contig with the highest counts and mark the rest as extra (or ambiguous if appropriate). Toggle ntop_vdj to keep the top n (default 1) contigs.

For VJ chains, the current rule set is to keep the top 2 productive contigs with the highest counts and mark the rest as extra (or ambiguous if appropriate). Toggle ntop_vj to keep the top n (default 2) contigs.

ambiguous is marked as T if the contig is of poor quality annotation and would be removed from downstream analysis. Cells with multiple contigs with very low umi_counts and/or consensus_counts are also marked as ambiguous as it is not possible to distinguish which is the most representative contig.

Please note that the default for filter_extra is True. If you want to keep the extra contigs for whatever reasons e.g. interested in T/B-cell development datasets, you need to set filter_extra=False. We are setting this as False in this example because later on we want to visualise these extra contigs.

[7]:
vdj, adata = ddl.pp.check_contigs(
    bcr, adata, library_type="ig", filter_extra=False
)
Filtering contigs...
Marking ambiguous contigs...
Initializing DandelionPolars object...

If you want to filter out contigs that are marked as pseudogenes or ORFs, you can use the filter_pseudo argument in check_contigs. This argument takes a dictionary with keys “v”, “d”, and “j” corresponding to the gene segments, and values as lists of functionality strings to filter out. For example, to only keep functional genes, you would set filter_pseudo={"v": ["P", "ORF"], "d": ["P", "ORF"], "j": ["P", "ORF"]}. Note that this requires that the functionality columns (v_call_functionality, d_call_functionality, j_call_functionality) are present in the AIRR table, which can be populated by running ddl.pp.annotate_functionality first.

Check the DandelionPolars object

[8]:
print(type(vdj))
vdj
<class 'dandelion.polars.core._core.DandelionPolars'>
[8]:
Lazy Dandelion object with n_obs = 2493 and n_contigs = 5785
    data: sequence_id, sequence, rev_comp, productive, v_call, d_call, j_call, sequence_alignment, germline_alignment, junction, junction_aa, v_cigar, d_cigar, j_cigar, stop_codon, vj_in_frame, locus, junction_length, np1_length, np2_length, v_sequence_start, v_sequence_end, v_germline_start, v_germline_end, d_sequence_start, d_sequence_end, d_germline_start, d_germline_end, j_sequence_start, j_sequence_end, j_germline_start, j_germline_end, v_score, v_identity, v_support, d_score, d_identity, d_support, j_score, j_identity, j_support, fwr1, fwr2, fwr3, fwr4, cdr1, cdr2, cdr3, cell_id, consensus_count, umi_count, v_call_10x, d_call_10x, j_call_10x, junction_10x, junction_10x_aa, j_call_blastn, j_identity_blastn, j_alignment_length_blastn, j_number_of_mismatches_blastn, j_number_of_gap_openings_blastn, j_sequence_start_blastn, j_sequence_end_blastn, j_germline_start_blastn, j_germline_end_blastn, j_support_blastn, j_score_blastn, j_sequence_alignment_blastn, j_germline_alignment_blastn, j_call_igblastn, j_source, j_support_igblastn, j_score_igblastn, d_call_blastn, d_identity_blastn, d_alignment_length_blastn, d_number_of_mismatches_blastn, d_number_of_gap_openings_blastn, d_sequence_start_blastn, d_sequence_end_blastn, d_germline_start_blastn, d_germline_end_blastn, d_support_blastn, d_score_blastn, d_sequence_alignment_blastn, d_germline_alignment_blastn, d_call_igblastn, d_source, d_support_igblastn, d_score_igblastn, v_call_genotyped, germline_alignment_d_mask, sample_id, c_call, c_sequence_alignment, c_germline_alignment, c_sequence_start, c_sequence_end, c_score, c_identity, c_call_10x, junction_aa_length, fwr1_aa, fwr2_aa, fwr3_aa, fwr4_aa, cdr1_aa, cdr2_aa, cdr3_aa, sequence_alignment_aa, v_sequence_alignment_aa, d_sequence_alignment_aa, j_sequence_alignment_aa, complete_vdj, j_call_multimappers, j_call_multiplicity, j_call_sequence_start_multimappers, j_call_sequence_end_multimappers, j_call_support_multimappers, mu_count, rearrangement_status, v_call_functionality, d_call_functionality, j_call_functionality, extra, ambiguous
    metadata: cell_id, sample_id, productive_VDJ, productive_VJ, d_call_VDJ, j_call_VDJ, j_call_VJ, junction_VDJ, junction_VJ, junction_aa_VDJ, junction_aa_VJ, locus_VDJ, locus_VJ, v_call_VDJ, v_call_VJ, c_call_VDJ, c_call_VJ, umi_count_VDJ, umi_count_VJ, productive_VDJ_main, productive_VJ_main, d_call_VDJ_main, j_call_VDJ_main, j_call_VJ_main, junction_VDJ_main, junction_VJ_main, junction_aa_VDJ_main, junction_aa_VJ_main, locus_VDJ_main, locus_VJ_main, v_call_genotyped_VDJ_main, v_call_genotyped_VJ_main, c_call_VDJ_main, c_call_VJ_main, umi_count_VDJ_main, umi_count_VJ_main, isotype, isotype_main, isotype_status, locus_status, chain_status, rearrangement_status_VDJ, rearrangement_status_VJ

Check the AnnData object as well

[9]:
adata
[9]:
AnnData object with n_obs × n_vars = 25057 × 31915
    obs: 'sample_id', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'gmm_pct_count_clusters_keep', 'scrublet_score', 'is_doublet', 'filter_rna', 'has_contig', 'productive_VDJ', 'productive_VJ', 'd_call_VDJ', 'j_call_VDJ', 'j_call_VJ', 'junction_VDJ', 'junction_VJ', 'junction_aa_VDJ', 'junction_aa_VJ', 'locus_VDJ', 'locus_VJ', 'v_call_VDJ', 'v_call_VJ', 'c_call_VDJ', 'c_call_VJ', 'umi_count_VDJ', 'umi_count_VJ', 'productive_VDJ_main', 'productive_VJ_main', 'd_call_VDJ_main', 'j_call_VDJ_main', 'j_call_VJ_main', 'junction_VDJ_main', 'junction_VJ_main', 'junction_aa_VDJ_main', 'junction_aa_VJ_main', 'locus_VDJ_main', 'locus_VJ_main', 'v_call_genotyped_VDJ_main', 'v_call_genotyped_VJ_main', 'c_call_VDJ_main', 'c_call_VJ_main', 'umi_count_VDJ_main', 'umi_count_VJ_main', 'isotype', 'isotype_main', 'isotype_status', 'locus_status', 'chain_status', 'rearrangement_status_VDJ', 'rearrangement_status_VJ'

These are the relevant columns for looking at the QC status of the cells and contigs in the .obs slot in the AnnData object (and also .metadata slot in the Dandelion object):

Relevant columns in obs

  • has_contig

  • whether cells have V(D)J chains.

  • locus_status

  • detailed information on chain status pairings (below).

  • chain_status

  • summarised information of the chain locus status pairings (similar to chain_pairing in scirpy).

  • rearrangement_status_VDJ and rearrangement_status_VJ

  • whether or not V(D)J gene usage are standard (i.e. all from the same locus).

So in a standard situation, I would remove cells flagged with Orphan VJ, Orphan VJ-exception, ambiguous in .metadata.chain_status, and also any cell marked as chimeric in the .metadata.rearrangement_status_VDJ and .metadata.rearrangement_status_VJ from downstream cell-level calculations/analysis.

Having said that, you will find that most of Dandelion’s functions will work without the need to requirement to perform additional filtering and filtering can be performed on the final AnnData object (described in the visualisation section).

Let’s take a look at these new columns

[10]:
pd.crosstab(adata.obs["chain_status"], adata.obs["locus_status"])
[10]:
locus_status + Extra VJ Extra VDJ + Extra VJ Extra VDJ + IGK Extra VDJ + IGL IGH + Extra VJ IGH + IGK IGH + IGL Orphan IGH Orphan IGK Orphan IGL
chain_status
Extra pair 14 66 9 5 101 0 0 0 0 0
Orphan VDJ 0 0 0 0 0 0 0 9 0 0
Orphan VJ 0 0 0 0 0 0 0 0 97 46
Single pair 0 0 0 0 0 1248 898 0 0 0

if there are multiple library types, i.e. ddl.pp.check_contigs was run with library_type = None, or if several tcr/bcr Dandelion objects are concatenated, there will be additional columns where the v/d/j/c calls and productive will be split into additional columns to reflect those that belong to a B cell, alpha-beta T cell, or gamma-delta T cell.

Now actually filter the AnnData object and run through a standard workflow starting by filtering genes and normalizing the data

Because the ‘filtered’ AnnData object was returned as a filtered but otherwise unprocessed object, we still need to normalize and run through the usual process here. The following is just a standard scanpy workflow.

[11]:
# standard preprocessing steps
sc.pp.filter_genes(adata, min_cells=3)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
adata.raw = adata
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
adata = adata[:, adata.var.highly_variable].copy()
sc.pp.regress_out(adata, ["total_counts", "pct_counts_mt"])
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, svd_solver="arpack")
sc.pp.neighbors(adata)
sc.tl.umap(adata)
sc.tl.leiden(adata)

Visualizing the clusters and whether or not there’s a corresponding V(D)J receptor

[12]:
sc.pl.umap(adata, color=["sample_id", "leiden", "chain_status"])
../../_images/notebooks_polars_2_dandelion_filtering-10x_data_polars_26_0.png

Visualizing some B cell genes

[13]:
sc.pl.umap(adata, color=["IGHM", "JCHAIN"])
../../_images/notebooks_polars_2_dandelion_filtering-10x_data_polars_28_0.png
[14]:
adata.write_h5ad("adata.h5ad", compression="gzip")

Save Dandelion object

To save the vdj object, we have two options - either save the .data and .metadata slots with polars’ functions:

[15]:
vdj.data.collect().write_csv("filtered_vdj_table.tsv", separator="\t")

Or save the whole DandelionPolars class object using one of the built-in write methods. The default and recommended format is .zipddl (Zarr ZipStore), which .write() now points to. It uses a hybrid storage scheme: Parquet for tabular data (.data and .metadata), Zarr arrays for distance matrices, and embedded HDF5 for graphs, layouts, and germline references. The legacy .h5ddl format (HDF5-based) is still available via .write_h5ddl().

[16]:
vdj.write_zipddl("dandelion_results.zipddl")

Running ddl.pp.check_contigs without AnnData

Finally, ddl.pp.check_contigs can also be run without an AnnData object:

[17]:
vdj3 = ddl.pp.check_contigs(bcr)
vdj3
Filtering contigs...
Marking ambiguous contigs...
Initializing DandelionPolars object...
[17]:
Lazy Dandelion object with n_obs = 3158 and n_contigs = 6384
    data: sequence_id, sequence, rev_comp, productive, v_call, d_call, j_call, sequence_alignment, germline_alignment, junction, junction_aa, v_cigar, d_cigar, j_cigar, stop_codon, vj_in_frame, locus, junction_length, np1_length, np2_length, v_sequence_start, v_sequence_end, v_germline_start, v_germline_end, d_sequence_start, d_sequence_end, d_germline_start, d_germline_end, j_sequence_start, j_sequence_end, j_germline_start, j_germline_end, v_score, v_identity, v_support, d_score, d_identity, d_support, j_score, j_identity, j_support, fwr1, fwr2, fwr3, fwr4, cdr1, cdr2, cdr3, cell_id, consensus_count, umi_count, v_call_10x, d_call_10x, j_call_10x, junction_10x, junction_10x_aa, j_call_blastn, j_identity_blastn, j_alignment_length_blastn, j_number_of_mismatches_blastn, j_number_of_gap_openings_blastn, j_sequence_start_blastn, j_sequence_end_blastn, j_germline_start_blastn, j_germline_end_blastn, j_support_blastn, j_score_blastn, j_sequence_alignment_blastn, j_germline_alignment_blastn, j_call_igblastn, j_source, j_support_igblastn, j_score_igblastn, d_call_blastn, d_identity_blastn, d_alignment_length_blastn, d_number_of_mismatches_blastn, d_number_of_gap_openings_blastn, d_sequence_start_blastn, d_sequence_end_blastn, d_germline_start_blastn, d_germline_end_blastn, d_support_blastn, d_score_blastn, d_sequence_alignment_blastn, d_germline_alignment_blastn, d_call_igblastn, d_source, d_support_igblastn, d_score_igblastn, v_call_genotyped, germline_alignment_d_mask, sample_id, c_call, c_sequence_alignment, c_germline_alignment, c_sequence_start, c_sequence_end, c_score, c_identity, c_call_10x, junction_aa_length, fwr1_aa, fwr2_aa, fwr3_aa, fwr4_aa, cdr1_aa, cdr2_aa, cdr3_aa, sequence_alignment_aa, v_sequence_alignment_aa, d_sequence_alignment_aa, j_sequence_alignment_aa, complete_vdj, j_call_multimappers, j_call_multiplicity, j_call_sequence_start_multimappers, j_call_sequence_end_multimappers, j_call_support_multimappers, mu_count, rearrangement_status, v_call_functionality, d_call_functionality, j_call_functionality, extra, ambiguous
    metadata: cell_id, sample_id, productive_VDJ, productive_VJ, d_call_VDJ, j_call_VDJ, j_call_VJ, junction_VDJ, junction_VJ, junction_aa_VDJ, junction_aa_VJ, locus_VDJ, locus_VJ, v_call_VDJ, v_call_VJ, c_call_VDJ, c_call_VJ, umi_count_VDJ, umi_count_VJ, productive_VDJ_main, productive_VJ_main, d_call_VDJ_main, j_call_VDJ_main, j_call_VJ_main, junction_VDJ_main, junction_VJ_main, junction_aa_VDJ_main, junction_aa_VJ_main, locus_VDJ_main, locus_VJ_main, v_call_genotyped_VDJ_main, v_call_genotyped_VJ_main, c_call_VDJ_main, c_call_VJ_main, umi_count_VDJ_main, umi_count_VJ_main, isotype, isotype_main, isotype_status, locus_status, chain_status, rearrangement_status_VDJ, rearrangement_status_VJ