<jupyter_start><jupyter_text>Compare two repertoires using `tcrdist`, `hierdiff` and `palmotif`This example notebook demonstrates how to: 1. Load/parse the [Dash et al. (2017)](https://pubmed.ncbi.nlm.nih.gov/28636592/) dataset 2. Compute distances between all pairs of unique TCR clones 3. Test for relative enrichment of a TCRs from an experimental condition (e.g. epitope-specific) within neighborhoods around each TCR 4. Test for relative enrichment within clusters defined by nodes of a hierarchical clustering tree 5. Make a sequence logo from CDR3 amino-acid frequencies within a cluster of TCRs 7. Display an interactive hierarhical clustering tree for navigating a TCR dataset<jupyter_code>"""Not required for general use: addition of these paths was helpful for loading dependencies and datasets during development"""
import sys
from os.path import join as opj
from fg_shared import _git
sys.path.append(opj(_git, 'hierdiff'))
sys.path.append(opj(_git, 'tcrdist2'))
sys.path.append(opj(_git, 'palmotif'))
import numpy as np
import pandas as pd
import scipy.spatial
import scipy.cluster.hierarchy as sch
import operator
from os.path import join as opj
import io

from IPython.core.display import SVG, HTML
import IPython
from tcrdist.repertoire import TCRrep
import tcrdist as td

"""hierdiff can be installed from PyPI and is a dependency of tcrdist2
It provides statistical testing and plotting for TCR clusters"""
import hierdiff

"""palmotif can be installed from PyPI and provides functions for computing
and displaying sequence logos as SVG"""
import palmotif<jupyter_output><empty_output><jupyter_text>1. Load and parse the dataThe Dash et al. (2017) dataset is provided in the `datasets` folder and is available for download [here](https://raw.githubusercontent.com/kmayerb/tcrdist2/API2/tcrdist/datasets/dash.csv). Demonstration will use the PB1 and NP mouse epitopes only. V-region and J-region gene usage is inferred from the nucleotide sequences. Unique/deduplicated clones are the input for this analysis.<jupyter_code>df = pd.read_csv(opj(td.__path__[0], 'datasets', 'dash.csv'))
df = df[df.epitope.isin(['PB1', 'NP'])]

#2 create instance of TCRrep class, initializes input as tr.cell_df attribute
tr = TCRrep(cell_df=df, chains=['beta'], organism='mouse')

#3 Infer CDR1,CDR2,CDR2.5 (a.k.a. phmc) from germline v-genes
tr.infer_cdrs_from_v_gene(chain='beta',  imgt_aligned=True)

#4 Define index columns for determining unique clones.
tr.index_cols = ['clone_id', 'subject', 'epitope',
                 'v_b_gene', 'j_b_gene',
                'cdr3_b_aa', 'cdr1_b_aa', 'cdr2_b_aa', 'pmhc_b_aa',
                'cdr3_b_nucseq']

#4 Deduplicate based on index cols, creating tr.clone_df attribute
tr.deduplicate()
tr.clone_df<jupyter_output><empty_output><jupyter_text>2. Compute distances between all pairs of unique TCR beta chainsThe `metric='nw'` defines a TCR distance as the alignment score from a Needleman-Wunsch global alignment between a pair of amino-acid sequences. Distances are computed for the CDR1, CDR2 and CDR3 regions and summed for the total TCR distance. The matrix is stored in the `pmhc_b_aa_pw` attribute of the `TCRRep` object. Note that the number of rows of `tr.clone_df` above and the square dimensions of the pairwise distance matrix are the same and the ordering of TCRs is preserved.<jupyter_code>tr.compute_pairwise_all('beta', metric='nw', processes=3)
pw = tr.pmhc_b_aa_pw
pw.shape<jupyter_output><empty_output><jupyter_text>3. Test for enrichment within TCR distance-defined neighborhoodsWhen TCRs are sampled from two different experimental conditions (e.g. two epitope sorts as they are here), we can look for enrichment of TCRs from a particular condition within neighborhoods around each unique TCR. One can think of this as finding the differences between one or more "repertoires", where each repertoire comes from a different experimental condition (e.g. pre/post treatment or two or more individuals).The neighborhood can be defined using a fixed radius or a fixed number of unique neighbors. The columns containing categorical conditions are specified in `x_cols`. The expansions of each unique TCR can be optionally provided as a `count_col`. Depending on the organization of the data a Chi-squared test, a Cochran–Mantel–Haenszel test or Fisher's Exact test can be applied (see docstring for details).Although this nearest-neighbor (NN) approach tests many nearly identical neighborhoods around similar TCRs, which can reduce overall power after a multiplicity adjustment, it can be helpful to have an enrichment statistic for each TCR that can be plotted relative to other TCR features (e.g. the probability of generation); it also is convenient for overlaying as a color on a scatter plot embedding of TCRs (left for another tutorial!). In the future the NN approach could be implemented for large datasets without ever having to compute and hold all pairwise distances in memory, and could take advantage of existing scalable NN machine learning tools.<jupyter_code>res = td.stats.neighborhood_diff(tr.clone_df,
                          pwmat=tr.pmhc_b_aa_pw,
                          x_cols=['epitope'],
                          count_col='count',
                          test_method='fishers',
                          knn_neighbors=None, knn_radius=5)
res = tr.clone_df.join(res)
res = res.sort_values(by='OR', ascending=False)<jupyter_output><empty_output><jupyter_text>Columns in the results fileThe results `pd.DataFrame` is designed in such a way that the counts could be used in a more sophisticated regression or other model for analysis. The statistical tests we provide are good for scanning repertoires for interesting differences, but it is likely that for real-world experimental designs, a more complex model will be required (e.g. adjusting for more covariates, continuous features, etc.) This format also allows for transporting the data via a CSV file into R where other statistical approaches are available.Each row of the table represents a single cluster of TCRs in a neighborhood around a unique TCR. By default the method tests a neighborhood around every unique TCR, though a subset can be specified. Each row of the table represents an *n x 2* table where *n* is the number of combinations of the `x_cols` values and the 2 columns indicate `MEM+/-` cluster membership. Here is a brief description of the resulting columns: - `epitope`, `v_b_gene`, `j_b_gene`, `cdr3_b_aa` all come from the original `clone_df` which can be joined to provde info about the centroid for each neighborhood - `K_neighbors` number of unique TCRs in the cluster (may be a fixed parameter) - `R_radius` distance to the furthest neighbor in the cluster (may be a fixed parameter) - `levels` contains the ordered levels of each `x_col` and cluster membership - `val_n` contains the values of each level associated with each count; MEM+ indicates cluste membership - `ct_n` contains the counts in each element of the table. associated with each `val_n` - `X|MEM+`, `X|MEM-` contain the  probability of observing `X` conditioned on cluster membership (non-membership), `MEM+` (`MEM-`); these are only provided for 2 x 2 analyses with a single binary varible in `x_cols` - `OR` the odds of being in the cluster for TCRs with the first level of the `x_col` divided by the odds of being in the cluster with the second level of the `x_col`; this is only provided for 2 x 2 analyses with a single binary varible in `x_cols` - `pvalue`, `FWERp` significance determined by the specified test; the family-wise error rate (FWER) adjusted p-value is important as a guide due to the typically large number of clusters/nodes that are tested. A false-discovery rate (FDR) adjusted q-value may also be valuable and can be calculated from the p-values using `td.stats.adjustnonnan` which wraps adjustments available in the `statsmodels` package.<jupyter_code>cols = ['epitope', 'v_b_gene', 'j_b_gene', 'cdr3_b_aa', 'K_neighbors', 'R_radius',
        'levels', 'val_0', 'ct_0', 'val_1', 'ct_1', 'val_2', 'ct_2', 'val_3', 'ct_3',
        'X|MEM+', 'X|MEM-', 'OR', 'pvalue', 'FWERp']
res[cols].drop_duplicates(['cdr3_b_aa']).head(10)<jupyter_output><empty_output><jupyter_text>Above shows clusters enriched with TCRs from the NP epitope sort, while below shows those enriched with TCRs from the PB1 sort.<jupyter_code>res[cols].drop_duplicates(['cdr3_b_aa']).tail(10)<jupyter_output><empty_output><jupyter_text>4. Make frequency-based logo plots of interesting clustersLogo plots are helpful for summarizing a set of similar sequences. Here we use logos to summarize the CDR3 sequences of interesting clusters, based on the neighborhood analysis. Since the CDR3 region has a variable length, we make a logo by aligning each CDR3 to a reference sequence, which here is naturally the centroid of the neighborhood.This first logo is a cluster with OR > 1 indicating enrichment of TCRs identified in the NP tetramer sort. The logo represents the log frequencies of the amino-acids in the cluster CDR3s at each position of the centroid sequence (x-axis), after the alignment. Note that the conserved cysteine-alanine residues at the N-terminus of the CDR3 are emphasized in this plot in addition to other positions.<jupyter_code>ind = 0
motif = palmotif.compute_pal_motif(res.loc[ind, 'cdr3_b_aa'],
                        seqs=res.loc[res.loc[ind, 'neighbors'], 'cdr3_b_aa'])
html = palmotif.svg_logo(motif, return_html=True, svg_height='100px', svg_width='800px')
HTML(data=html)<jupyter_output><empty_output><jupyter_text>The second logo of the same cluster uses symbol heights to represent the log-odds ratio of observing the amino-acid within the cluster vs. outside the cluster; this is helpful for showing which amino-acid residues help differentiate the cluster from all the other sequences. For this analysis, a reference of random TCRs could also be helpful. Note that the conserved residues at the N-terminus are now demphasized because all TCRs within and outside the neighborhood contain those residues. The method also highlights that a cluster could potentially be defined by a depletion of an amino-acid at a particular position (i.e. negative log-OR).<jupyter_code>ind = 0
motif = palmotif.compute_pal_motif(res.loc[ind, 'cdr3_b_aa'],
                        seqs=res.loc[res.loc[ind, 'neighbors'], 'cdr3_b_aa'],
                         refs=res.loc[[i for i in res.index if not i in res.loc[ind, 'neighbors']], 'cdr3_b_aa'])
html = palmotif.svg_logo(motif, return_html=True, svg_height='100px', svg_width='800px')
HTML(data=html)<jupyter_output><empty_output><jupyter_text>We can also look at a cluster with OR < 1 indicating enrichment for PB1-specific TCRs, again with or without using sequences outside the neighborhood as a reference.<jupyter_code>ind = 250
motif = palmotif.compute_pal_motif(res.loc[ind, 'cdr3_b_aa'],
                        seqs=res.loc[res.loc[ind, 'neighbors'], 'cdr3_b_aa'])
html = palmotif.svg_logo(motif, return_html=True, svg_height='100px', svg_width='800px')
HTML(data=html)
ind = 250
motif = palmotif.compute_pal_motif(res.loc[ind, 'cdr3_b_aa'],
                        seqs=res.loc[res.loc[ind, 'neighbors'], 'cdr3_b_aa'],
                         refs=res.loc[[i for i in res.index if not i in res.loc[ind, 'neighbors']], 'cdr3_b_aa'])
html = palmotif.svg_logo(motif, return_html=True, svg_height='100px', svg_width='800px')
HTML(data=html)<jupyter_output><empty_output><jupyter_text>5. Test for enrichment within TCR clusters on a hierarchical clustering treeOne issue with the NN approach is that the nieghborhood size or radius must be pre-specified. A hierarchical clustering approach is able to test both small clusters of similar TCRs near the "leaves" of the tree and larger clusters of less similar TCRs near the "root"; in this way the two approaches are complementary. Here we use `scipy` and the pairwise TCR distances to cluster the TCRs and then test every node of the tree for enrichment/depletion of TCRs associated with an experimental condition(s) in `x_cols`. Otherwise the API is similar to that of `td.stats.neighborhood_diff` and allows a few different simple tests that can help scan two or more repertoires for potentially interesting differences.Note that the rows of the testing results do not match the rows of the input `clone_df` one-to-one, the way they did with the neighborhood approach. Instead we use `td.stats.member_summ` to add columns to the results file that aid interpretation (e.g. top 2 V-genes or most common CDR3).<jupyter_code>hdiff_res, Z = td.stats.hcluster_diff(tr.clone_df,
                          pwmat=tr.cdr3_b_aa_pw,
                          x_cols=['epitope'],
                          count_col='count',
                          test_method='fishers',
                          hclust_method='ward')
"""Adding epitope is superfluous since all the testing is around epitope specificity.
addl_cols could be used to summarize other phenotype info in other columns"""
summ_df = td.stats.member_summ(res_df=hdiff_res, clone_df=tr.clone_df, count_col='count', addl_cols=['epitope'])
hdiff_res = hdiff_res.join(summ_df)
hdiff_res = hdiff_res.sort_values(by='pvalue')<jupyter_output><empty_output><jupyter_text>Columns in the results fileThe results `pd.DataFrame` can be provided directly to `hierdiff` plotting functions, but is also designed in such a way that the counts could be used in a more sophisticated regression or other model for analysis. The statistical tests we provide are good for scanning repertoires for interesting differences, but it is likely that for real-world experimental designs, a more complex model will be required (e.g. adjusting for more covariates, continuous features, etc.) This also allows for transporting the data into R where other statistical approaches are available.Each row of the table represents a single cluster/node of TCRs on the hiearchical clustering tree. The method tests every node of the tree, which can generate many correlated results, but is also helpful for a complete and unbiased look at the data. You may notice that small clusters near the leaves lack sample size and power to detect a significant enrichment. Larger clusters near the root may not be significant if they include too many unrelated TCRs. For epitope specific repertiores we may expect the tight clusters that reflect recognition of the same epitope to be most enriched; however these issues are not evident in this simple comparison of two epitope specific repertoires.Note that there can also be false-discoveries with large clusters. The sample size can be high making the tests overly sensitive to small deviations from the null-hypothesis. One should interpret large clusters with caution because clusters of dissimilar TCRs may not share a specificity or be biologically relevant. Also, TCR data is fundamentally compositional, which means that all frequencies add to one for a given sample and an enrichment of one TCR must mean a depletion of all other TCRs in the sample. Though this can be an issue with any analysis of TCR data, the large nodes of a hierarchical clustering tree are particularly susceptible; if accumulation on one subtree leads to enrichment it may make it look as though there is depletion on the other subtrees. Generally, large datasets (e.g. bulk sequencing) and diverse datasets are less susceptible to issues involving compositionality.Each row of the table represents an `n x 2` table where n is the number of combinations of the `x_cols` values and the 2 columns indicate `MEM+/-` cluster membership. Here is a brief description of the resulting columns: - `epitope`, `v_b_gene`, `j_b_gene`, `cdr3_b_aa` are all summaries of the original `clone_df` generated by `td.stats.member_summ` - `K_neighbors` number of unique TCRs in the cluster - `R_radius` maximum distance between TCRs in the cluster - `levels` contains the ordered levels of each `x_col` and cluster membership - `val_n` contains the values of each level associated with each count; MEM+ indicates cluste membership - `ct_n` contains the counts in each element of the table. associated with each `val_n` - `X|MEM+`, `X|MEM-` contain the  probability of observing `X` conditioned on cluster membership (non-membership), `MEM+` (`MEM-`); these are only provided for 2 x 2 analyses with a single binary varible in `x_cols` - `OR` the odds of being in the cluster for TCRs with the first level of the `x_col` divided by the odds of being in the cluster with the second level of the `x_col`; this is only provided for 2 x 2 analyses with a single binary varible in `x_cols` - `pvalue`, `FWERp` significance determined by the specified test; the family-wise error rate (FWER) adjusted p-value is important as a guide due to the typically large number of clusters/nodes that are tested. A false-discovery rate (FDR) adjusted q-value may also be valuable and can be calculated from the p-values using `td.stats.adjustnonnan` which wraps adjustments available in the `statsmodels` package.<jupyter_code>hdiff_res.columns
cols = ['epitope', 'v_b_gene', 'j_b_gene', 'cdr3_b_aa', 'K_neighbors', 'R_radius',
        'levels', 'val_0', 'ct_0', 'val_1', 'ct_1', 'val_2', 'ct_2', 'val_3', 'ct_3',
        'X|MEM+', 'X|MEM-', 'OR', 'pvalue', 'FWERp']
hdiff_res[cols]
html = hierdiff.plot_hclust_props(Z, title='NP versus PB1 repertoires (mouse)',
                                  res=hdiff_res,
                                  tooltip_cols=['v_b_gene', 'j_b_gene', 'cdr3_b_aa', 'X|MEM+', 'X|MEM-', 'OR', 'FWERp'],
                                  alpha=0.05, alpha_col='FWERp')
# Output HTML can be saved to a file or rendered in the output of a jupyter notebook. For the latter an IFrame is required to isolate the D3/javascript.
# h = IPython.display.HTML(html)
# IPython.display.display_html(h,  metadata={"isolated": True, 'height':700})
"""Save the tree to a local file, which can be read and displayed in an IFrame"""
fn = './PB1_NP_tree.html'
with open(fn, 'w') as fh:
    fh.write(html)
IPython.display.IFrame(src='./PB1_NP_tree.html', height=700, width=1000)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07407815079801872, 'fraction_numerical': 0.0089157952669235, 'mean_word_length': 4.5467032967032965, 'pattern_counts': {'":': 1, '<': 46, '<?xml version=': 0, '>': 46, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '397576', 'license': 'MIT', 'n_tokens_mistral': 5396, 'n_tokens_neox': 4977, 'n_words': 2318, 'path': 'important_notebooks/hierdiff_example.ipynb', 'repo_name': 'agartland/tcrdist2', 'size': '18170'}