<jupyter_start><jupyter_text>Covid Data Dataset containing SARS-CoV-2 T cells described in https://www.medrxiv.org/content/10.1101/2020.12.07.20245274v1Later available under https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE171037<jupyter_code>import pandas as pd
import scanpy as sc
import scipy
import anndata
import scirpy as ir
import numpy as np
import matplotlib.pyplot as pl
from matplotlib import rcParams
import seaborn as sb
sb.set_style('white')<jupyter_output><empty_output><jupyter_text>Single cell data preprocessing Due to limited amount of cells, the data is not filtered in a similar way as the 10x dataset by expression counts, mitochondrial counts and gene counts.Only cells are deleted, that have a non-functional / non-unique receptor, since our models can not deal with that.<jupyter_code>path_base = '../data/Covid/'
path_raw = path_base + 'Covid19_TAs_PBMCs.h5ad'
path_tcr_split = path_base + '03_covidcovidtrain_test_split.h5ad'
path_highlyvar = path_base + '04_covid_highly_var_5000.h5ad'<jupyter_output><empty_output><jupyter_text>Filter cells with productive single IR  Filter cells where there is more or less than one alpha and beta chain present<jupyter_code>adata = sc.read_h5ad(path_raw)
print(f'Number of cells: {adata.n_obs}')
adata = adata[adata.obs['has_tcr']=='True']
print(f'Number of cells with IR: {adata.n_obs}')
print(f'Number of cells: {adata.n_obs}')
adata = adata[adata.obs['multi_chain']=='False']

print(f'Number of cells without >2 chains: {adata.n_obs}')

adata = adata[adata.obs['TRA_2_cdr3']=='None']
adata = adata[adata.obs['TRB_2_cdr3']=='None']

print(f'Number of cells with =1 per chain: {adata.n_obs}')<jupyter_output>Number of cells: 11277
Number of cells without >2 chains: 11147
Number of cells with =1 per chain: 9853<jupyter_text>Filter unpaired Sequences<jupyter_code>print(f'Number of cells: {adata.n_obs}')
adata = adata[adata.obs['TRB_1_cdr3'] != 'None']
print(f'Number of cells with beta chain: {adata.n_obs}')
print(f'Number of cells: {adata.n_obs}')
adata = adata[adata.obs['TRA_1_cdr3'] != 'None']
print(f'Number of cells with alpha chain: {adata.n_obs}')<jupyter_output>Number of cells: 9521
Number of cells with alpha chain: 6713<jupyter_text>Sanity check: only a single CDR3alpha and CDR3beta chain is present<jupyter_code>def check_invalids(sequence):
    allowed_letters = 'GALMFWKQESPVICYHRNDT'
    for letter in sequence:
        if letter not in allowed_letters:
            return False
    return True

for seq_alpha in adata.obs['TRA_1_cdr3']:
    if not check_invalids(seq_alpha):
        print(seq_alpha)
for seq_beta in adata.obs['TRB_1_cdr3']:
    if not check_invalids(seq_beta):
        print(seq_beta)<jupyter_output><empty_output><jupyter_text>Clonotype Annotation Give numerica label to each clonotype<jupyter_code>adata.obs['TRA+TRB'] = adata.obs['TRA_1_cdr3'].astype(str) + '+' + adata.obs['TRB_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'].value_counts()<jupyter_output>Trying to set attribute `.obs` of view, copying.<jupyter_text>Embed TR 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 as tcr
import tcr_embedding.utils_preprocessing as utils
adata.obs['TRA+TRB'] = adata.obs['TRA_1_cdr3'].astype(str) + '+' + adata.obs['TRB_1_cdr3'].astype(str)
pad = adata.obs['TRA+TRB'].str.len().max()
pad = int(pad)
pad
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='TRA+TRB', ohe_col='one_hot', label_col='tcr_seq', length_col='seq_len', pad=pad, aa_to_id=aa_to_id, start_end_symbol=True)

print(adata.obsm['tcr_seq'].shape)
print(adata.obsm['one_hot'].shape)
print(adata.uns['aa_to_id'])
max_length_alpha = adata.obs['TRA_1_cdr3'].str.len().max()
max_length_beta = adata.obs['TRB_1_cdr3'].str.len().max()
max_length = max(max_length_alpha, max_length_beta)
aa_to_id = adata.uns['aa_to_id']
utils.aa_encoding(adata, read_col='TRA_1_cdr3', label_col='alpha_seq', length_col='alpha_len', pad=max_length, aa_to_id=aa_to_id, start_end_symbol=False)
utils.aa_encoding(adata, read_col='TRB_1_cdr3', label_col='beta_seq', length_col='beta_len', pad=max_length, aa_to_id=aa_to_id, start_end_symbol=False)
print(max_length_alpha)
print(max_length_beta)<jupyter_output>20
23<jupyter_text>Train-Val-Test Split Datasets are split s.t. all cells from the same clonotype (here: identical alpha and beta chain) are assigned to the same split. The two samples and the condiation (stimulated vs unstimulated) are distributed evenly over all three sets.<jupyter_code>adata.obs['identifier+condition'] = adata.obs['identifier'].astype(str) + '_' + adata.obs['condition'].astype(str)
adata.obs['identifier+condition'].sample(10)
random_seed = 15
# Create Train-Val and Test set
train_val, test = utils.stratified_group_shuffle_split(adata.obs, stratify_col='identifier+condition', group_col='clonotype', val_split=0.15, random_seed=random_seed)
# Split Train-Val into Train and Val set
train, val = utils.stratified_group_shuffle_split(train_val, stratify_col='identifier+condition', group_col='clonotype', val_split=0.15, random_seed=random_seed)
print('Train Samples: ', len(train))
print('Train Ratio: ', len(train) / len(adata.obs))
print('Val Samples: ', len(val))
print('Val Ratio: ', len(val) / len(adata.obs))
print('Test Samples: ', len(test))
print('Test Ratio: ', len(test) / len(adata.obs))<jupyter_output>Train Samples:  5343
Train Ratio:  0.7959183673469388
Val Samples:  680
Val Ratio:  0.10129599284969462
Test Samples:  690
Test Ratio:  0.10278563980336661<jupyter_text>Visualize Distribution<jupyter_code>pl.figure(figsize=(15,5))
pl.hist([train['identifier'], val['identifier'], test['identifier']], density=True, bins=len(train['identifier'].unique()), label=['train', 'val', 'test'])
pl.xticks(rotation='vertical')
pl.legend(loc='upper right')
pl.show()
pl.figure(figsize=(15,5))
pl.hist([train['condition'], val['condition'], test['condition']], density=True, bins=len(train['condition'].unique()), label=['train', 'val', 'test'])
pl.xticks(rotation='vertical')
pl.legend(loc='upper right')
pl.show()
adata.obs.loc[train.index, 'set'] = 'train'
adata.obs.loc[val.index, 'set'] = 'val'
adata.obs.loc[test.index, 'set'] = 'test'
adata.obs['set'].value_counts()
adata.write_h5ad(path_tcr_split, compression='gzip')<jupyter_output>... storing 'TRA+TRB' as categorical
... storing 'identifier+condition' as categorical
... storing 'set' as categorical<jupyter_text>Highly Variable Genes Reduce the expression data to the 5000 most variable genes.<jupyter_code>adata = sc.read_h5ad(path_tcr_split)
sc.pp.highly_variable_genes(adata, n_top_genes=5000, batch_key='identifier')
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>
{'fraction_non_alphanumeric': 0.15159392529757834, 'fraction_numerical': 0.03187850595156656, 'mean_word_length': 6.971646673936751, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 38, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 1, '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': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '356309', 'license': 'MIT', 'n_tokens_mistral': 2867, 'n_tokens_neox': 2668, 'n_words': 593, 'path': 'preprocessing/covid_preprocessing.ipynb', 'repo_name': 'SchubertLab/mvTCR', 'size': '7309'}