<jupyter_start><jupyter_text>Paths to preprocessed data at different stages.<jupyter_code>path_base = '../data/Borcherding/'
path_merged = path_base + '01_borch_merged.h5ad'
path_filtered = path_base + '02_borch_filtered_tcr.h5ad'
path_highlyvar = path_base + '03_borch_highly_var.h5ad'
path_annotated = path_base + '04_borch_annotated.h5ad'<jupyter_output><empty_output><jupyter_text>Merge Files The dataset of Borcherding contains 1.5 million cells from various joint scRNA and VDJ studies. To convert the dataset to python format refer to the notebook "RDS_2_AnnData.ipynb". In the following they are processed to scanpy / scirpy format.<jupyter_code># Transcriptome data
file_name = path_base + 'rna_combined.h5ad'
adata = sc.read(file_name)
# TCR data
path_anno = path_base + 'tcrs_combined.csv'
anno_tcr = pd.read_csv(path_anno,index_col=0)  
anno_tcr = anno_tcr.set_index('barcode')

for col in ['IR_VJ_2_cdr3', 'IR_VDJ_2_cdr3', 'IR_VDJ_2_c_gene', 'IR_VDJ_2_junction_ins', 'IR_VDJ_2_v_gene',
            'IR_VDJ_2_expr', 'IR_VJ_2_d_gene', 'IR_VDJ_2_j_gene', 'IR_VJ_2_v_gene', 'IR_VJ_2_j_gene', 
            'IR_VJ_2_expr', 'IR_VDJ_2_locus', 'IR_VJ_2_junction_ins', 'IR_VDJ_2_cdr3_nt', 'IR_VJ_2_expr_raw',
            'IR_VDJ_2_expr_raw', 'IR_VJ_2_c_gene', 'IR_VDJ_2_d_gene', 'IR_VDJ_2_cdr3', 'IR_VJ_2_cdr3_nt' 
            'IR_VDJ_1_d_gene', 'IR_VJ_1_junction_ins', 'IR_VDJ_1_expr', 'IR_VJ_1_d_gene', 'IR_VJ_1_expr_raw', 
            'IR_VJ_1_expr', 'IR_VDJ_1_expr_raw', 'IR_VDJ_1_junction_ins', 'IR_VDJ_1_d_gene', 'IR_VJ_2_cdr3_nt']:
    anno_tcr[col] = None

rename_dict = {
#    'TCR1': ,
    'cdr3_aa1': 'IR_VJ_1_cdr3',
    'cdr3_nt1': 'IR_VJ_1_cdr3_nt',
#    'TCR2': ,
    'cdr3_aa2': 'IR_VDJ_1_cdr3', 
    'cdr3_nt2': 'IR_VDJ_1_cdr3_nt',
}
anno_tcr = anno_tcr.rename(rename_dict, axis=1)
anno_tcr['IR_VJ_1_locus'] = 'TRA'
anno_tcr['IR_VJ_2_locus'] = 'TRA'
anno_tcr['IR_VDJ_1_locus'] = 'TRB'
anno_tcr['IR_VDJ_2_locus'] = 'TRB'

anno_tcr['IR_VJ_1_v_gene'] = [el.split('.')[0] if isinstance(el, str) else None for el in anno_tcr['TCR1']]
anno_tcr['IR_VJ_1_c_gene'] = [el.split('.')[2] if isinstance(el, str) else None for el in anno_tcr['TCR1']]
anno_tcr['IR_VJ_1_j_gene'] = [el.split('.')[1] if isinstance(el, str) else None for el in anno_tcr['TCR1']]

anno_tcr['IR_VDJ_1_v_gene'] = [el.split('.')[0] if isinstance(el, str) else None for el in anno_tcr['TCR2']]
anno_tcr['IR_VDJ_1_c_gene'] = [el.split('.')[3] if isinstance(el, str) else None for el in anno_tcr['TCR2']]
anno_tcr['IR_VDJ_1_j_gene'] = [el.split('.')[2] if isinstance(el, str) else None for el in anno_tcr['TCR2']]


anno_tcr['multi_chain'] = False
anno_tcr['has_ir'] = ~anno_tcr['IR_VJ_1_cdr3'].isnull() | ~anno_tcr['IR_VDJ_1_cdr3'].isnull()

anno_tcr.index.name = None
anno_tcr = ann.AnnData(obs=anno_tcr)
anno_tcr.obs_names_make_unique()
anno_tcr
ir.pp.merge_with_ir(adata, anno_tcr)
sc.write(adata=adata, filename=path_merged)
count_irs = sum([1 for x in adata.obs['has_ir'] if x=='True'])
print(f'Total Count: {len(adata)}')
print(f'With IR-Info: {count_irs}')
adata.obs['Tissue'].value_counts()
adata.obs['Type'].value_counts()<jupyter_output><empty_output><jupyter_text>Qualitiy control Filter out cells without alpha and beta chain<jupyter_code>adata = sc.read(path_merged)
print(f'Total cells: {len(adata)}')
adata = adata[adata.obs['has_ir']=='True']
print(f'With IR: {len(adata)}')
adata = adata[adata.obs['IR_VDJ_1_cdr3']!='nan']
print(f'With beta chain: {len(adata)}')
adata = adata[adata.obs['IR_VJ_1_cdr3']!='nan']
print(f'With alpha chain: {len(adata)}')
adata = adata[~adata.obs['IR_VJ_1_cdr3'].str.contains(';')]
print(f'Without dual alpha chains: {len(adata)}')
adata = adata[~adata.obs['IR_VDJ_1_cdr3'].str.contains(';')]
print(f'Without dual beta chains: {len(adata)}')
adata = adata[adata.obs['db.class']=='singlet']
print(f'Without doublets: {len(adata)}')
sc.pp.filter_genes(adata, min_cells=100)
print(f'Number of cells after gene filter: {adata.n_vars}')<jupyter_output>filtered out 37683 genes that are detected in less than 100 cells
Trying to set attribute `.var` of view, copying.<jupyter_text>Load data and calculate the amount of counts, genes and fraction of mitochondrial genes.<jupyter_code>adata.obs['n_counts'] = adata.X.sum(1)
adata.obs['log_counts'] = np.log10(adata.obs['n_counts'])

adata.obs['n_genes'] = (adata.X > 0).sum(1)

mt_gene_mask = [gene.startswith('MT-') for gene in adata.var_names]
mt_gene_index = np.where(mt_gene_mask)[0]
adata.obs['mt_fraction'] = adata.X[:, mt_gene_index].sum(1) / adata.X.sum(1)
rcParams['figure.figsize'] = (15, 5)
sc.pl.violin(adata, ['log_counts'], size=1, log=False, rotation=90)
sc.pl.violin(adata, ['n_genes'], size=1, log=False, rotation=90)
sc.pl.violin(adata, ['mt_fraction'], size=1, log=False, rotation=90)
rcParams['figure.figsize'] = (8, 8)
sc.pl.scatter(adata, y='n_genes', x='n_counts', color='mt_fraction', size=10, show=False)
sc.pl.scatter(adata[np.logical_and(adata.obs['n_genes']<2000, adata.obs['n_counts']<10000)],
             y='n_genes', x='n_counts', color='mt_fraction', size=10, show=False)
sc.write(adata=adata, filename=path_filtered)<jupyter_output><empty_output><jupyter_text>Normalization and highly variable gene selection<jupyter_code>adata = sc.read(path_filtered)
sc.pp.normalize_total(adata, target_sum=1e4)  # target sum? 
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=5000)
print('Shape before: ', adata.shape)
adata = adata[:, adata.var['highly_variable']]
print('Shape after: ', adata.shape)
sc.write(adata=adata, filename=path_highlyvar)<jupyter_output><empty_output><jupyter_text>Model Specific Preproccessing<jupyter_code>adata = sc.read(path_highlyvar)<jupyter_output><empty_output><jupyter_text>Clonotype annotation The clonotype annotation from 10x defines clonotypes not accross all donors. Therefor, public TCRs might leak between training and testing data. Further, cells with the same TCR but multiple IRs are classified as different clonotypes. We therefor assign own clonotype labels based on complete correspondence based on alpha and beta chain.<jupyter_code>adata.obs['TRA+TRB'] = adata.obs['IR_VJ_1_cdr3'].astype(str) + '+' + adata.obs['IR_VDJ_1_cdr3'].astype(str)
clono_dict = {clone: idx for idx, clone in enumerate(adata.obs['TRA+TRB'].unique())}
adata.obs['clonotype'] = adata.obs['TRA+TRB'].map(clono_dict)
print(len(adata.obs['clonotype'].unique()))
adata.obs['clonotype'].head()<jupyter_output>302947<jupyter_text>Embed Proteins Append both chains, add start ('') and seperator ('+') token. For training purpose IR data might be needed as one hot vector.<jupyter_code>import sys
sys.path.append('../')
import tcr_embedding.utils_preprocessing as utils
aa_to_id = {'_': 0, 'A': 1, 'C': 2, 'D': 3, 'E': 4, 'F': 5, 'G': 6, 'H': 7, 'I': 8, 'K': 9, 'L': 10, 'M': 11, 'N': 12, 'P': 13,
            'Q': 14, 'R': 15, 'S': 16, 'T': 17, 'V': 18, 'W': 19, 'Y': 20, '+': 21, '<': 22, '>': 23}

utils.aa_encoding(adata, read_col='IR_VJ_1_cdr3', label_col='alpha_seq', length_col='alpha_len', pad=30, aa_to_id=aa_to_id, start_end_symbol=False)
utils.aa_encoding(adata, read_col='IR_VDJ_1_cdr3', label_col='beta_seq', length_col='beta_len', pad=30, aa_to_id=aa_to_id, start_end_symbol=False)
adata.obsm['beta_seq']
adata.obs['functional.cluster'] = adata.obs['functional.cluster'].cat.add_categories('nan')
adata.obs['functional.cluster'] = adata.obs['functional.cluster'].fillna('nan')
adata.write_h5ad(path_annotated)<jupyter_output>... storing 'TRA+TRB' as categorical
{'fraction_non_alphanumeric': 0.16582177161152614, 'fraction_numerical': 0.029215581643543223, 'mean_word_length': 7.069967707212056, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 32, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '317221', 'license': 'MIT', 'n_tokens_mistral': 3288, 'n_tokens_neox': 3120, 'n_words': 561, 'path': 'preprocessing/Borcherding_preprocessing.ipynb', 'repo_name': 'SchubertLab/mvTCR', 'size': '7496'}