From 8ed5a55c8268f46da76a1d1c631e96d051e7491c Mon Sep 17 00:00:00 2001 From: dimalvovs Date: Wed, 17 Jun 2026 19:21:21 -0400 Subject: [PATCH 01/21] init channel-based reports --- modules/local/report/render_notebook.nf | 32 ++++++++++++++++ tests/fixtures/test.qmd | 9 +++++ .../local/report/render_notebook.nf.test | 38 +++++++++++++++++++ 3 files changed, 79 insertions(+) create mode 100644 modules/local/report/render_notebook.nf create mode 100644 tests/fixtures/test.qmd create mode 100644 tests/modules/local/report/render_notebook.nf.test diff --git a/modules/local/report/render_notebook.nf b/modules/local/report/render_notebook.nf new file mode 100644 index 0000000..5b98d4d --- /dev/null +++ b/modules/local/report/render_notebook.nf @@ -0,0 +1,32 @@ +// Generic process to render a Quarto notebook to HTML +process RENDER_NOTEBOOK { + tag "${report_name}" + label 'process_single' + + input: + tuple val(report_name), path(notebook), val(data_dir) + val project_name + val workflow_cmd + + output: + path "${report_name}.html" + + script: + """ + ## copy quarto notebook to working directory + cp $notebook ${report_name}.qmd + + ## render qmd report to html + quarto render ${report_name}.qmd \\ + -P project_name:${project_name} \\ + -P workflow_cmd:'${workflow_cmd}' \\ + -P project_dir:${data_dir} \\ + -P sample_table:${file(params.samplesheet)} \\ + --to html + """ + + stub: + """ + touch ${report_name}.html + """ +} diff --git a/tests/fixtures/test.qmd b/tests/fixtures/test.qmd new file mode 100644 index 0000000..2ab696a --- /dev/null +++ b/tests/fixtures/test.qmd @@ -0,0 +1,9 @@ +--- +title: "test notebook" +format: html +engine: knitr +--- + +```{python} +print("Hello world!") +``` diff --git a/tests/modules/local/report/render_notebook.nf.test b/tests/modules/local/report/render_notebook.nf.test new file mode 100644 index 0000000..6497f1d --- /dev/null +++ b/tests/modules/local/report/render_notebook.nf.test @@ -0,0 +1,38 @@ +nextflow_process { + + name "Test RENDER_NOTEBOOK" + script "modules/local/report/render_notebook.nf" + process "RENDER_NOTEBOOK" + + test("Should render basic notebook") { + // note that on Mac M series, jupyter may fail in container + // so the test notebook uses knitr engine instead + + tag "basic" + tag "notebook" + + when { + params { + samplesheet = "${projectDir}/tests/fixtures/valid_samplesheet.csv" + } + + process { + """ + input[0] = [ + "test_notebook", + file("${projectDir}/tests/fixtures/test.qmd"), + "${projectDir}/tests/fixtures" + ] + input[1] = "TCRtoolkit" + input[2] = "nextflow run main.nf" + """ + } + } + + then { + assert process.success + assert process.out.size() == 1 + } + } + +} From 9e398cf7c57b45bf692d30da7ad86f4d4eded675 Mon Sep 17 00:00:00 2001 From: dimalvovs Date: Wed, 17 Jun 2026 19:31:53 -0400 Subject: [PATCH 02/21] Apply suggestions from code review Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- tests/modules/local/report/render_notebook.nf.test | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tests/modules/local/report/render_notebook.nf.test b/tests/modules/local/report/render_notebook.nf.test index 6497f1d..bbbba00 100644 --- a/tests/modules/local/report/render_notebook.nf.test +++ b/tests/modules/local/report/render_notebook.nf.test @@ -32,6 +32,11 @@ nextflow_process { then { assert process.success assert process.out.size() == 1 + + def html = path(process.out.get(0)) + assert html.getFileName().toString() == 'test_notebook.html' + assert html.exists() + assert html.text.contains('Hello world!') } } From b23a86ec1786bae86145b25ed254f09112806ec2 Mon Sep 17 00:00:00 2001 From: dimalvovs Date: Wed, 17 Jun 2026 21:20:08 -0400 Subject: [PATCH 03/21] Revert "Apply suggestions from code review" This reverts commit 9e398cf7c57b45bf692d30da7ad86f4d4eded675. --- tests/modules/local/report/render_notebook.nf.test | 5 ----- 1 file changed, 5 deletions(-) diff --git a/tests/modules/local/report/render_notebook.nf.test b/tests/modules/local/report/render_notebook.nf.test index bbbba00..6497f1d 100644 --- a/tests/modules/local/report/render_notebook.nf.test +++ b/tests/modules/local/report/render_notebook.nf.test @@ -32,11 +32,6 @@ nextflow_process { then { assert process.success assert process.out.size() == 1 - - def html = path(process.out.get(0)) - assert html.getFileName().toString() == 'test_notebook.html' - assert html.exists() - assert html.text.contains('Hello world!') } } From 291f56ab364bbe4d916bc72820e4772a77a1cfc3 Mon Sep 17 00:00:00 2001 From: dimalvovs Date: Thu, 18 Jun 2026 16:59:45 -0400 Subject: [PATCH 04/21] clean up notebooks --- notebooks/compare_stats_template.qmd | 190 --- notebooks/gliph2_report_template.qmd | 56 - notebooks/sample_stats_template.qmd | 316 ---- notebooks/template_qc.qmd | 2328 ++++++++++++++++++++++++++ 4 files changed, 2328 insertions(+), 562 deletions(-) delete mode 100644 notebooks/compare_stats_template.qmd delete mode 100644 notebooks/gliph2_report_template.qmd delete mode 100644 notebooks/sample_stats_template.qmd create mode 100644 notebooks/template_qc.qmd diff --git a/notebooks/compare_stats_template.qmd b/notebooks/compare_stats_template.qmd deleted file mode 100644 index 01b2c66..0000000 --- a/notebooks/compare_stats_template.qmd +++ /dev/null @@ -1,190 +0,0 @@ ---- -title: "Comparative T Cell Repertoire statistics" -format: - html: - theme: flatly - toc: true - toc_depth: 3 - code-fold: show - embed-resources: true - number-sections: true - smooth-scroll: true - grid: - body-width: 1000px - margin-width: 300px - -jupyter: python3 ---- - -Thank you for using TCRtoolkit! This report is generated from sample data and metadata you provided. The report is divided into two sections: - -@sec-heatmap : Heatmap of sample to sample repertoire similarity using Jaccard, Sorensen, and Morisita indices (v1.0) - -# Report Setup - -```{python, echo=false} -#| tags: [parameters] -#| echo: false - -## 2. Pipeline Parameters -#Default inputs are overwritten at the command line in `modules/local/plot_sample.nf` -workflow_cmd='' -project_name='path/to/project_name' -jaccard_mat='path/to/jaccard_mat.csv' -sorensen_mat='path/to/sorensen_mat.csv' -morisita_mat='path/to/morisita_mat.csv' -# jsd_mat='path/to/jsd_mat.csv' -sample_utf8='path/to/sample_utf8.csv' -``` - -```{python} -#| tags: [setup] -#| warning: false - -# 1. Load Packages -import os -import shutil -import datetime -import sys -import numpy as np -import pandas as pd -import seaborn as sns -# from matplotlib.colors import LinearSegmentedColormap -# import scipy.cluster.hierarchy as sch - -# 2. Print Pipeline Information -print('Pipeline information and parameters:' + '\n') -print('Project Name: ' + project_name) -print('Workflow command: ' + workflow_cmd) -print('Date and time: ' + str(datetime.datetime.now())) - -# 3. Importing similarity data -## 3a. jaccard similarity matrix -jaccard_df = pd.read_csv(jaccard_mat, sep=',', header=0, index_col=0) - -## 3b. sorensen similarity matrix -sorensen_df = pd.read_csv(sorensen_mat, sep=',', header=0, index_col=0) - -## 3c. morisita similarity matrix -morisita_df = pd.read_csv(morisita_mat, sep=',', header=0, index_col=0) - -``` - -# Analysis - -## Overall Repertoire Similarity {#sec-heatmap} - -Similarity metrics such as Jaccard, Sorensen, and Morisita are often used to compare the similarity between two samples. Here, we compare the similarity of TCR repertoires between samples using these three metrics. Details on how each metric is calculated can be found below the figure. - -```{python} -import plotly.express as px -import plotly.graph_objects as go -from plotly.figure_factory import create_dendrogram -from plotly.subplots import make_subplots -import scipy.spatial.distance as ssd -import matplotlib.pyplot as plt - -# preprocess the data prior to clustering -jaccard_numeric = jaccard_df.apply(pd.to_numeric, errors='coerce') -sorensen_numeric = sorensen_df.apply(pd.to_numeric, errors='coerce') -morisita_numeric = morisita_df.apply(pd.to_numeric, errors='coerce') - -# Assuming jaccard_numeric, sorensen_numeric, and morisita_numeric are your DataFrames -# and sns_cluster_jaccard, sns_cluster_sorensen, and sns_cluster_morisita are the corresponding clustermaps -sns_cluster_jaccard = sns.clustermap(jaccard_numeric) -plt.close() -sns_cluster_sorensen = sns.clustermap(sorensen_numeric) -plt.close() -sns_cluster_morisita = sns.clustermap(morisita_numeric) -plt.close() -# sns_cluster_jsd = sns.clustermap(jsd_numeric) -# plt.close() - -# Create a subplot with 3 rows -fig = make_subplots(rows=3, cols=1) - -# Reindex the dataframes to match the clustering -jaccard_clustered = jaccard_numeric.iloc[sns_cluster_jaccard.dendrogram_row.reordered_ind, sns_cluster_jaccard.dendrogram_col.reordered_ind] -sorensen_clustered = sorensen_numeric.iloc[sns_cluster_sorensen.dendrogram_row.reordered_ind, sns_cluster_sorensen.dendrogram_col.reordered_ind] -morisita_clustered = morisita_numeric.iloc[sns_cluster_morisita.dendrogram_row.reordered_ind, sns_cluster_morisita.dendrogram_col.reordered_ind] -# jsd_clustered = jsd_numeric.iloc[sns_cluster_jsd.dendrogram_row.reordered_ind, sns_cluster_jsd.dendrogram_col.reordered_ind] - -# Create individual heatmaps -heatmap_jaccard = go.Heatmap( - z=jaccard_clustered, - x=jaccard_clustered.columns, - y=jaccard_clustered.index, - coloraxis="coloraxis", - visible=False) -heatmap_sorensen = go.Heatmap( - z=sorensen_clustered, - x=sorensen_clustered.columns, - y=sorensen_clustered.index, - coloraxis="coloraxis", - visible=False) -heatmap_morisita = go.Heatmap( - z=morisita_clustered, - x=morisita_clustered.columns, - y=morisita_clustered.index, - coloraxis="coloraxis", - visible=True) - -# Add the heatmaps to the figure -fig = go.Figure(data=[heatmap_jaccard, heatmap_sorensen, heatmap_morisita]) - -# Create buttons to switch between the heatmaps -buttons = [ - dict(label="Jaccard", method="update", - args=[{"visible": [True, False, False]}, {"title": "Jaccard"}]), - dict(label="Sorensen", method="update", - args=[{"visible": [False, True, False]}, {"title": "Sorensen"}]), - dict(label="Morisita", method="update", - args=[{"visible": [False, False, True]}, {"title": "Morisita"}]) -] - -# Update the layout of the figure -fig.update_layout( - updatemenus=[dict(type="buttons", showactive=True, buttons=buttons)], - title='Similarity Matrices', - xaxis_title='Sample ID', - yaxis_title='Sample ID', - autosize=False, - width=950, - height=950, - coloraxis=dict(colorscale='Viridis', colorbar=dict(title="Log Scale")) -) - -fig.show() -``` - -**Jaccard Index**: - -$$ -J(A,B) = \frac{|A \cap B|}{|A \cup B|} -$$ - -Where $A$ and $B$ are the sets of unique CDR3 amino acid sequences (TCRs) in samples A and B, respectively. The Jaccard Index is defined as the ratio of the number of common elements between two sets to the total number of distinct elements in the two sets. - -**Sorensen Index**: - -$$ -S(A,B) = \frac{2|A \cap B|}{|A| + |B|} -$$ - -Where $A$ and $B$ are the sets of unique CDR3 amino acid sequences (TCRs) in samples A and B, respectively. The difference between the sorensen index and the jaccard index is that the sorensen index takes into account the size of the two sets being compared, while the jaccard index only considers the number of common elements between the two sets. - -**Morisita-Horn Index:** - -$$ -M(A,B)=\frac{2\sum_{i=1}^{S}a_{i}b_{i}}{(D_{a}+D_{b})AB} ; D_{a}=\frac{\sum_{i=1}^{S}a_{i}^{2}}{A^2}, D_{b}=\frac{\sum_{i=1}^{S}b_{i}^{2}}{B^2} -$$ - -Where: - -- $A$ and $B$ are the sets of unique CDR3 amino acid sequences (TCRs) in samples A and B, - -- $a_{i}$ ($b_{i}$) is the number of times TCR $i$ is represented in the total $A$ ($B$) from one sample. - -- $S$ is the total number of unique TCRs in the two samples. - -- $D_{a}$ and $D_{b}$ are the Simpson Index values for samples A and B, respectively. diff --git a/notebooks/gliph2_report_template.qmd b/notebooks/gliph2_report_template.qmd deleted file mode 100644 index e39d248..0000000 --- a/notebooks/gliph2_report_template.qmd +++ /dev/null @@ -1,56 +0,0 @@ ---- -title: "Comparative T Cell Repertoire statistics" -format: - html: - theme: flatly - toc: true - toc_depth: 3 - code-fold: show - embed-resources: true - number-sections: true - smooth-scroll: true - grid: - body-width: 1000px - margin-width: 300px - -jupyter: python3 ---- - -Thank you for using TCRtoolkit! This report is generated from sample data and metadata you provided. The report is divided into two sections: - -@sec-heatmap : Heatmap of sample to sample repertoire similarity using Jaccard, Sorensen, and Morisita indices (v1.0) - -# Report Setup - -```{python, echo=false} -#| tags: [parameters] -#| echo: false - -## 2. Pipeline Parameters -#Default inputs are overwritten at the command line in `modules/local/plot_gliph2.nf` -workflow_cmd='' -project_name='path/to/project_name' -clusters='path/to/{project_name}_cluster.csv' -cluster_stats='path/to/{project_name}_cluster.txt' -``` - -```{python} -#| tags: [setup] -#| warning: false - -# 1. Load Packages -import os -import shutil -import datetime -import sys -import numpy as np -import pandas as pd -import seaborn as sns -# from matplotlib.colors import LinearSegmentedColormap -# import scipy.cluster.hierarchy as sch - -# 2. Print Pipeline Information -print('Pipeline information and parameters:' + '\n') -print('Project Name: ' + project_name) -print('Workflow command: ' + workflow_cmd) -print('Date and time: ' + str(datetime.datetime.now())) diff --git a/notebooks/sample_stats_template.qmd b/notebooks/sample_stats_template.qmd deleted file mode 100644 index a85d7ef..0000000 --- a/notebooks/sample_stats_template.qmd +++ /dev/null @@ -1,316 +0,0 @@ ---- -title: "Sample Level T Cell Repertoire statistics" -format: - html: - theme: flatly - toc: true - toc_depth: 3 - code-fold: show - embed-resources: true - number-sections: true - smooth-scroll: true - grid: - body-width: 1000px - margin-width: 300px - -jupyter: python3 ---- - -Thank you for using TCRtoolkit! This report is generated from sample data and metadata you provided. The report is divided into three sections: - -@sec-report-setup : Code to setup the report. This section includes the parameters you used to run the pipeline, loading necessary packages, data, etc. - -@sec-sample-level-stats : Typical sample level T cell repertoire statistics. Each plot has a description about the statistic shown and formulas used to calculate the metric. - -@sec-gene-family-usage : TCR V gene family usage. The x-axis shows the timepoint collected for each individual, and the y-axis shows the proportion of TCRs that use each V gene family. - -# Report Setup {#sec-report-setup} - -```{python, echo=false} -#| tags: [parameters] -#| echo: false - -workflow_cmd='' -project_name='' -sample_table='' -sample_stats_csv='' -v_family_csv='' - -samplechart_x_col='' -samplechart_color_col='' -vgene_subject_col='' -vgene_x_cols='' -``` - -```{python} -#| code-fold: true - -# 1. Load Packages -from IPython.display import Image -import os -import datetime -import sys -import pandas as pd -import math -import matplotlib.pyplot as plt -import seaborn as sns -from matplotlib.colors import LinearSegmentedColormap -import plotly.express as px -import plotly.graph_objects as go - -import warnings -warnings.filterwarnings( - "ignore", - category=FutureWarning, - module="plotly" -) - -# 2. Print pipeline parameters - -print('Project Name: ' + project_name) -print('Workflow command: ' + workflow_cmd) -print('Date and time: ' + str(datetime.datetime.now())) - -# 3. Loading data -## reading sample metadata -meta = pd.read_csv(sample_table, sep=',') -meta_cols = meta.columns.tolist() - -## reading combined repertoire statistics -df = pd.read_csv(sample_stats_csv, sep=',') -df = pd.merge(df, meta, on='sample', how='left') -df = df[meta_cols + [c for c in df.columns if c not in meta_cols]] - -## reading V gene family usage -v_family = pd.read_csv(v_family_csv, sep=',') -v_family = pd.merge(v_family, meta, on='sample', how='left') -v_family = v_family[meta_cols + [c for c in v_family.columns if c not in meta_cols]] -v_family = v_family.sort_values(by=[vgene_subject_col]) -``` - -# Sample level statistics {#sec-sample-level-stats} - -Below are plots showing basic T cell repertoire statistics. Each plot has a description about the statistic shown and formulas used to calculate the metric when applicable. Specific biological interpretation of each plot is left to the user. - -Version 3 of these plots features plotly express interactive plots. This version is exploratory and may be updated in the future. - -```{python} -#| code-fold: true - - -x_category = df[samplechart_x_col].unique().tolist() -x_category.sort() - -def samplechart(df, y_col): - fig = px.box( - df, - x=samplechart_x_col, - y=y_col, - color=samplechart_color_col, - points='all', - hover_data='sample', - category_orders={f'{samplechart_x_col}': x_category} - ) - - fig.update_layout( - margin=dict(b=80) - ) - - fig.show() -``` - -## Number of clones - -```{python} -#| code-fold: true - -samplechart(df, y_col='num_clones') -``` - -**Figure 1. Number of clones across sample groupings.** A clone is defined as a T cell with a unique CDR3 amino acid sequence. The number of clones is shown on the y-axis. The x-axis represents sample groupings defined by user-specified metadata fields (eg. origin, timepoint). - -## Clonality - -```{python} -#| code-fold: true - -samplechart(df, y_col='clonality') -``` - -**Figure 2. The clonality of samples across sample groupings.** Clonality is a measure of T cell clonal expansion and reflects the degree to which the sample is dominated by 1 or more T cell clones. Clonality is calculated via: $$Clonality = \frac {1-H} {\log_{2} N} \quad\text{,}\quad H = -\sum\limits_{i=1}^N p_i \log_{2}{p_i}$$ where $H$ is the Shannon entropy of a given sample, $N$ is the number of unique TCRs in the sample, and $p_i$ is the frequency of the $i$ th unique TCR in the sample. - -## Simpson Index - -```{python} -#| code-fold: true - -samplechart(df, y_col='simpson_index_corrected') -``` - -**Figure 3. Corrected Simpson Index.** The Simpson Index is a measure of diversity that takes into account the number of clones and the relative abundance of each clone in a sample. The corrected Simpson Index, $D$, is calculated as: - -$$D = \sum\limits_{i=1}^N \frac {p_i(p_i - 1)} {T(T - 1)} \quad\text{,}\quad T = \sum\limits_{i=1}^N p_i$$ - -Where $N$ is the number of unique TCRs in the sample, $p_i$ is the frequency of the $i$ th unique TCR in the sample, and $T$ is the total number of T Cells counted in the sample. - -## Percent of productive TCRs - -```{python} -#| code-fold: true - -samplechart(df, y_col='pct_prod') -``` - -**Figure 4. Percent of productive TCRs.** A productive TCR is a DNA/RNA sequence that can be translated into a protein sequence, i.e. it does not contain a premature stop codon or an out of frame rearrangement. The percent of productive TCRs is calculated as: - -$$ Percent \text{ } productive \text{ } TCRs = \frac P N $$ - -where $P$ is the number of productive TCRs and $N$ is the total number of TCRs in a given sample. - -## Average productive CDR3 Length - -```{python} -#| code-fold: true - -samplechart(df, y_col='productive_cdr3_avg_len') -``` - -**Figure 5. Average Productive CDR3 Length** The average length of the CDR3 region of the TCR for productive clones. The CDR3 region is the most variable region of the TCR and is the region that determines antigen specificity. - -## TCR Convergence - -```{python} -#| code-fold: true - -samplechart(df, y_col='ratio_convergent') -``` - -**Figure 6. TCR Convergence** The ratio of convergent TCRs to total TCRs. A convergent TCR is a TCR that is generated via 2 or more unique nucleotide sequences via codon degeneracy. - -# Gene Family Usage {#sec-gene-family-usage} - -## V gene family usage - -The V gene family usage of the TCRs in each sample is shown in the plots below. The x-axis shows the timepoint collected for each individual, and the y-axis shows the proportion of TCRs that use each V gene family. - -The V gene usage proportion, $V_k$, is calculated via: - -$$ -V_k = \frac {N_{k}} {T} \quad\text{,}\quad T = \sum\limits_{i=1}^N p_i -$$ - -where $N_{k}$ is the number of TCRs that use the $k$ th V gene, and T is the total number of TCRs in the sample. - -```{python} -#| code-fold: true - -## code adapted from https://www.moritzkoerber.com/posts/plotly-grouped-stacked-bar-chart/ -colors = ["#fafa70","#fdef6b","#ffe566","#ffda63","#ffd061","#ffc660","#ffbb5f","#fdb15f","#fba860","#f79e61","#f39562","#ef8c63","#e98365","#e37b66","#dd7367","#d66b68","#ce6469","#c65e6a","#bd576b","#b4526b","#ab4c6b","#a1476a","#974369","#8c3e68","#823a66","#773764","#6d3361","#62305e","#572c5a","#4d2956"] -vgene_x_cols = vgene_x_cols.split(',') - -## calculate proportions and add to v_family_long -v_family_long = pd.melt(v_family, - id_vars=meta_cols, - value_vars=[c for c in v_family.columns if c.startswith('TRBV')], - var_name='v_gene', - value_name='count') -v_family_long['proportion'] = v_family_long.groupby(meta_cols)['count'].transform(lambda x: x / x.sum()) - -## add in the total number of v genes for each sample -total_v_genes = v_family_long.groupby(meta_cols)['count'].sum().reset_index().rename(columns={'count': 'total_v_genes'}) -v_family_long = v_family_long.merge(total_v_genes, on=meta_cols) - -subjects = ( - v_family_long[vgene_subject_col] - .dropna() - .unique() -) - -for subject in subjects: - current = v_family_long[v_family_long[vgene_subject_col] == subject] - fig = go.Figure() - fig.update_layout( - template="simple_white", - title_text=f"{vgene_subject_col}: {subject}", - xaxis=dict(title_text=f"{', '.join(vgene_x_cols)}"), - yaxis=dict(title_text="proportion"), - barmode="stack", - margin=dict(b=80, r=115), - legend=dict( - orientation="v", - x=1.02, - y=1, - yanchor="top", - itemsizing="constant", - traceorder="reversed" - ) - ) - - fig.update_xaxes( - automargin=True, - title_standoff=30 - ) - - if len(vgene_x_cols) == 1: - x_vals = current[vgene_x_cols[0]] - else: - x_vals = current[vgene_x_cols].astype(str).agg('_'.join, axis=1) - summary = ( - current - .assign(x=x_vals) - .drop_duplicates(subset=meta_cols) - .groupby('x', as_index=False) - .agg( - total_v_genes=('total_v_genes', 'sum'), - ) - ) - y_df = ( - current - .assign(x=x_vals) - .groupby('x', as_index=False) - .agg( - y_top=('proportion', 'sum') - ) - ) - summary = y_df.merge(summary, on='x') - - for g, c in zip(current.v_gene.unique(), colors): - plot_df = current[current.v_gene == g] - if len(vgene_x_cols) == 1: - x = plot_df[vgene_x_cols[0]] - else: - x = [plot_df[col] for col in vgene_x_cols] - fig.add_trace( - go.Bar( - x=x, - y=plot_df['proportion'], - name=g, - marker_color=c, - legendgroup=g, - - customdata=plot_df[['sample']], - hovertemplate=( - "V gene: %{fullData.name}
" - "Sample: %{customdata[0]}
" - "Proportion: %{y}
" - "" - ), - - text=plot_df['total_v_genes'] if ((g == 'TRBV30') and (len(vgene_x_cols) != 1)) else None, - textposition='outside', - ) - ) - - if len(vgene_x_cols) == 1: - for _, row in summary.iterrows(): - fig.add_annotation( - x=row['x'], - y=row['y_top'] + 0.02, - text=str(row['total_v_genes']), - showarrow=False, - yanchor='bottom', - font=dict(size=12) - ) - - fig.show() -``` diff --git a/notebooks/template_qc.qmd b/notebooks/template_qc.qmd new file mode 100644 index 0000000..7928457 --- /dev/null +++ b/notebooks/template_qc.qmd @@ -0,0 +1,2328 @@ +--- +title: "Quality Control" +format: + html: + theme: flatly + toc: true + toc_depth: 3 + code-fold: true + embed-resources: true + number-sections: true + smooth-scroll: true + grid: + body-width: 1000px + margin-width: 300px + +jupyter: python3 +--- + + + +Thank you for using TCRtoolkit! This report is generated from the data provided. + +:::{.callout-note collapse="true"} +## Document Information +**Current Version:** 1.0-beta +**Last Updated:** March 2026 +**Maintainer:** BTC Data Science Team +**Notes:** +::: + +```{python} +#| tags: [parameters] +#| include: false + +# --------------------------------------------------------- +# BASE PARAMETERS +# --------------------------------------------------------- +workflow_cmd = '' +project_name='' +project_dir='' +sample_table='' + +timepoint_col = 'timepoint' +timepoint_order_col = 'timepoint_order' +alias_col = 'alias' +subject_col = 'subject_id' + +``` + +```{python} +#| include: false + +# --------------------------------------------------------- +# DERIVED PATHS +# --------------------------------------------------------- + +# Define files +project_dir=f"{project_dir}/{project_name}" +sample_stats_csv = f"{project_dir}/sample/sample_stats.csv" +concat_file = f"{project_dir}/annotate/concatenated_cdr3_sorted.tsv" + + +``` + +```{python} +#| code-fold: true + +from IPython.display import Image +import os +import datetime +import sys +import pandas as pd +import math +import matplotlib.pyplot as plt +import seaborn as sns +from matplotlib.colors import LinearSegmentedColormap +import plotly.express as px +import plotly.graph_objects as go +import glob +import itertools +import h5py +import igraph as ig +import matplotlib.ticker as ticker +import numpy as np +import scipy.cluster.hierarchy as sch +from IPython.display import display, Markdown +from scipy.sparse import csr_matrix +from scipy.stats import gaussian_kde +from scipy.stats import entropy +from scipy.stats import skew, wasserstein_distance +from scipy.stats import pearsonr +from sklearn.decomposition import PCA +from sklearn.preprocessing import StandardScaler +from scipy import sparse +from sklearn.preprocessing import normalize +from scipy.stats import wilcoxon +from scipy.stats import mannwhitneyu +from itertools import combinations +from scipy.spatial.distance import pdist +from scipy.cluster.hierarchy import linkage, leaves_list +import plotly.figure_factory as ff + +import warnings + +meta = pd.read_csv(sample_table, sep=',') +meta_cols = meta.columns.tolist() + +df = pd.read_csv(sample_stats_csv, sep=',') +df = pd.merge(df, meta, on='sample', how='left') +df = df[meta_cols + [c for c in df.columns if c not in meta_cols]] + +concat_df = pd.read_csv(concat_file, sep='\t') + +``` + +# Technical QC & Biological Relevance 🚩/🦠 {#sec-qc-bio} +## Diversity metrics {#sec-div-metrics} + +Analyzing the architecture of a T-cell repertoire requires moving beyond simply counting the total number of unique sequences (Richness). A healthy immune system is highly diverse and relatively evenly distributed. However, upon antigen encounter, specific T-cells rapidly proliferate, causing the repertoire to become highly skewed. +To quantify this structural shift, we rely on ecological diversity metrics and inequality indices. + +- **Shannon Entropy ($H$)** +Answers the question: "How hard is it to predict the identity of a randomly drawn T-cell?" It quantifies the uncertainty of the clonotype distribution. + - ***High Entropy (High Uncertainty):*** The repertoire is very diverse. There are many different clones with comparable frequencies, so you have no idea which one you will pick next. This signifies a healthy, polyclonal repertoire. + - ***Low Entropy (Low Uncertainty):*** The repertoire is dominated by a few large clones. You can easily guess that the next T-cell will likely belong to the dominant clone. This signifies oligoclonality or clonality. + +It is calculated as: +$H = - \sum_{i=1}^{R} p_i \ln(p_i)$ +Where $R$ is the number of unique clonotypes (Richness), and $p_i$ is the frequency of the $i$-th clonotype. + +- **Inverse Simpson Index ($1/D$)** +Answers the question: "How many TCR clones actually matter in this sample?" It ignores the "long tail" of rare, single-read clones (which might just be sequencing errors or debris) and focuses on the dominant, expanded clones that actually drive the immune response. +While the standard Simpson Index ($D$) measures the probability of collision (picking the same clone twice), the Inverse Simpson converts that probability into an "effective number". + - **High Value:** Indicates a high effective number of TCR clones. The sample is diverse and even in clone size. + - **Low Value (approaching 1):** Indicates the sample behaves as if it only contains 1 clone. It is highly dominated by a single expansion. +It is calculated as: +$\frac{1}{D} = \frac{1}{\sum_{i=1}^{R} p_i^2}$ +Where $p_i$ is the frequency of the $i$-th clonotype. + +::: {.callout-tip title="Pro Tip"} +Compare with sample Richness (# of unique clones) + +If Richness = 5,000 but Inverse Simpson = 1.2: Your sample behaves as if it only has 1.2 clones. It is massively dominated by a single expansion. + +If Richness = 5,000 and Inverse Simpson = 4,500: Your sample behaves as if it has 4,500 equally sized clones. It is extremely even and polyclonal. +::: + + +- **Gini Coefficient** +Answers the question: "Has the immune system picked a winner yet?" Measures inequality. Gini is your best metric for detecting Clonal Expansions independent of how many clones you actually sequenced. + - ***Value close to 0 (Perfect Equality):*** Every clone has the exact same frequency. The immune system is "resting" or "naive." It hasn't been triggered by a specific threat, so no single clone has started to divide rapidly. + - ***Medium value (0.3 - 0.7):*** Reactive. The immune system is fighting something. A few clones have expanded to fight a virus or tumor, but the background diversity is still there. + - ***High value (> 0.8):*** Monoclonal / Oligoclonal. A few clones have taken over completely. If this is a tumor sample, it might indicate a TIL (Tumor Infiltrating Lymphocyte) expansion. +It is calculated as: +$G = \frac{\sum_{i=1}^{R} (2i - R - 1) p_i}{R \sum_{i=1}^{R} p_i}$ + +Where the frequencies $p_i$ are sorted in ascending order, and $i$ is the rank. + +::: {.callout-tip title="Pro Tip"} +Why use this instead of Shannon? +Shannon Entropy is hard to compare if your library sizes (depth) are wildly different. Gini is normalized (always 0 to 1). +::: + +- **Hill Numbers ($^qD$)** +Answer the question: "What is the effective number of species in the sample when we ignore (penalize to degree $q$) rare clones?" This is a unified metric that combines richness, entropy, and dominance into a single scale (counts of effective species). + - ***$q=0$ (Richness):*** Counts every unique clone equally, regardless of frequency. Sensitive to sequencing depth and errors. + - ***$q=1$ (Exponential Shannon):*** Weighs clones by their frequency. Represents the number of "common" clones. + - ***$q=2$ (Inverse Simpson):*** Heavily weighs dominant clones. Represents the number of "very dominant" clones. +The general formula for Hill numbers of order $q$ is: +$^qD = \left( \sum_{i=1}^{R} p_i^q \right)^{1/(1-q)}$ +For $q=0$, this simplifies to $R$ (Total unique clonotypes). For $q=1$, the limit is undefined, so we use $\exp(H)$. For $q=2$, this simplifies to $1 / \sum p_i^2$ (Inverse Simpson). + + +```{python} +#| label: diversity-nested +#| output: asis + +# --- 1. CALCULATION FUNCTION --- +def calculate_diversity_metrics(df): + results = [] + # Check if 'counts' exists, if not try 'clone_count' or similar + count_col = 'counts' if 'counts' in df.columns else df.columns[0] # Fallback + + grouped = df.groupby('sample') + for sample, data in grouped: + counts = data[count_col].values + if counts.sum() == 0: continue + + p = counts / counts.sum() + shannon = -np.sum(p * np.log(p)) + inv_simpson = 1.0 / np.sum(p**2) + sorted_p = np.sort(p) + n = len(p) + index = np.arange(1, n + 1) + gini = ((2 * index - n - 1) * sorted_p).sum() / (n * sorted_p.sum()) + q0 = len(p) + q1 = np.exp(shannon) + q2 = inv_simpson + results.append({ + 'sample': sample, + 'shannon_entropy': shannon, + 'inverse_simpson': inv_simpson, + 'gini_coefficient': gini, + 'hill_q0': q0, + 'hill_q1': q1, + 'hill_q2': q2 + }) + return pd.DataFrame(results).set_index('sample') + +# --- 2. DATA PREPARATION --- +metrics = calculate_diversity_metrics(concat_df) +metrics = metrics.reset_index() + +# Ensure numeric columns are float +num_cols = metrics.select_dtypes(include=[np.number]).columns +metrics[num_cols] = metrics[num_cols].astype(float) + +# Merge with metadata +plot_df = pd.merge(metrics, meta, on='sample', how='inner') + +metrics_to_plot = [ + 'shannon_entropy', 'inverse_simpson', 'gini_coefficient', + 'hill_q0', 'hill_q1', 'hill_q2' +] + +# --- 3. SETUP GROUPS & ORDERING --- + +# Create Timepoint Mapping +if timepoint_col in plot_df.columns and timepoint_order_col in plot_df.columns: + time_order_map = plot_df.set_index(timepoint_col)[timepoint_order_col].to_dict() +else: + time_order_map = {} + +# Define Exclusions +exclude_cols = ['sample', 'file', 'total_counts'] +exclude_from_all_samples = [timepoint_col, timepoint_order_col, 'protocol_day', alias_col] + +# Groups for "All Samples" tab +group_opts_all = [ + col for col in meta.columns + if col not in exclude_cols + and col not in exclude_from_all_samples + and meta[col].nunique() < 35 +] + +# Groups for "By Patient" tab +group_opts_patient = [timepoint_col] if timepoint_col in plot_df.columns else [] + +unique_patients = sorted(plot_df[subject_col].dropna().unique().tolist()) + +# --- 4. PLOTTING FUNCTION --- + +def create_diversity_plot(data, x_col, y_metric, custom_order=None, plot_type="box"): + + if data.empty or data[x_col].dropna().empty: + return + + # Determine order + if custom_order: + unique_cats = [x for x in custom_order if x in data[x_col].unique()] + else: + unique_cats = sorted(data[x_col].dropna().unique().tolist()) + + # --- BASE PLOT --- + if plot_type == "box": + fig = px.box( + data, + x=x_col, + y=y_metric, + color=x_col, + points="all", + hover_data=['alias'], + category_orders={x_col: unique_cats}, + template="simple_white", + width=700, + height=600, + title=f"{y_metric} by {x_col}" + ) + + fig.update_traces(width=0.5, marker=dict(size=10, opacity=0.7)) + + for trace in fig.data: + if isinstance(trace, go.Box): + trace.pointpos = 0 + trace.jitter = 0.2 + + elif plot_type == "line": + # Sort for proper line connection + data = data.copy() + data[x_col] = pd.Categorical(data[x_col], categories=unique_cats, ordered=True) + data = data.sort_values(x_col) + + fig = px.scatter( + data, + x=x_col, + y=y_metric, + color_discrete_sequence=["black"], + hover_data=['alias'], + template="simple_white", + width=700, + height=500, + title=f"{y_metric} by {x_col}" + ) + + # Add dotted line + fig.add_trace(go.Scatter( + x=data[x_col], + y=data[y_metric], + mode='lines+markers', + line=dict(dash='dot', width=2), + marker=dict(size=18, color='black'), + showlegend=False + )) + + # --- STATISTICS (ONLY FOR BOXPLOTS) --- + stack_counter = 0 + + if plot_type == "box" and len(unique_cats) >= 2: + + pairs = list(combinations(unique_cats, 2)) + + y_max = data[y_metric].max() + y_min = data[y_metric].min() + y_range = y_max - y_min + if y_range == 0: + y_range = 1 + + base_offset = y_range * 0.08 # was 0.2 → MUCH closer to boxes + step_size = y_range * 0.06 # was 0.15 → tighter stacking + text_offset = y_range * 0.015 # was 0.03 → less gap above line + + for t1, t2 in pairs: + + g1 = data[data[x_col] == t1][y_metric].dropna() + g2 = data[data[x_col] == t2][y_metric].dropna() + + if len(g1) < 2 or len(g2) < 2: + continue + + try: + stat, p_value = mannwhitneyu(g1, g2, alternative='two-sided') + except ValueError: + continue + + if p_value >= 0.05: + continue + + if p_value < 0.001: + symbol = '***' + elif p_value < 0.01: + symbol = '**' + else: + symbol = '*' + + y_bracket = y_max + base_offset + (stack_counter * step_size) + y_text = y_bracket + text_offset + + fig.add_shape( + type="line", xref="x", yref="y", + x0=t1, y0=y_bracket, x1=t2, y1=y_bracket, + line=dict(color="black", width=1.5) + ) + + tick_len = y_range * 0.0015 + + fig.add_shape(type="line", xref="x", yref="y", + x0=t1, y0=y_bracket, x1=t1, y1=y_bracket - tick_len) + + fig.add_shape(type="line", xref="x", yref="y", + x0=t2, y0=y_bracket, x1=t2, y1=y_bracket - tick_len) + + try: + x_center = (unique_cats.index(t1) + unique_cats.index(t2)) / 2 + except: + x_center = 0 + + fig.add_annotation( + x=x_center, + y=y_text, + text=f"{symbol}
p={p_value:.3f}", + showarrow=False, + font=dict(size=10) + ) + + stack_counter += 1 + + extra_space = y_range * (0.15 + 0.07 * stack_counter) + fig.update_yaxes(range=[y_min, y_max + extra_space]) + + # --- FINAL LAYOUT --- + dynamic_height = 600 + (stack_counter * 40) + + fig.update_layout( + showlegend=False, + yaxis_title=y_metric, + xaxis_title=x_col, + xaxis=dict(tickfont=dict(size=15)), + height=dynamic_height, + margin=dict(t=80), + plot_bgcolor='rgba(0,0,0,0)' + ) + + fig.update_yaxes(showgrid=True, gridcolor='lightgrey') + + fig.show() + +# --- 5. EXECUTION LOOP --- +print(":::::: {.panel-tabset}\n") # Level 1 Start + +for metric in metrics_to_plot: + print(f"## {metric}\n") + + if len(unique_patients) <= 10: + print("::::: {.panel-tabset}\n") # Level 2 Start + + # --- TAB A: ALL SAMPLES --- + print("### All Samples\n") + print(":::: {.panel-tabset}\n") # Level 3 Start + for group in group_opts_all: + print(f"#### {group}\n") + create_diversity_plot(plot_df, group, metric) + print("\n") + print("::::\n") # Level 3 End + + # --- TAB B: BY PATIENT --- + print("### By Patient\n") + print(":::: {.panel-tabset}\n") # Level 3 Start + + for pat in unique_patients: + print(f"#### {pat}\n") + pat_df = plot_df[plot_df[subject_col] == pat] + + has_multi_origin = ( + 'origin' in pat_df.columns and + pat_df['origin'].nunique() > 1 + ) + + if has_multi_origin: + print("::: {.panel-tabset}\n") # Level 4 Start + origins = sorted(pat_df['origin'].dropna().unique()) + + for origin in origins: + print(f"##### {origin}\n") + origin_df = pat_df[pat_df['origin'] == origin] + + for group in group_opts_patient: + if origin_df[group].nunique() > 0: + custom_order = None + if group == timepoint_col and time_order_map: + pat_tps = origin_df[group].dropna().unique().tolist() + custom_order = sorted( + pat_tps, + key=lambda x: time_order_map.get(x, 999) + ) + create_diversity_plot( + origin_df, group, metric, + custom_order, plot_type="line" + ) + else: + print(f"No data for {group} in {origin}, patient {pat}.") + print("\n") + + print(":::\n") # Level 4 End ← was missing in some branches + + else: + for group in group_opts_patient: + if pat_df[group].nunique() > 0: + custom_order = None + if group == timepoint_col and time_order_map: + pat_tps = pat_df[group].dropna().unique().tolist() + custom_order = sorted( + pat_tps, + key=lambda x: time_order_map.get(x, 999) + ) + create_diversity_plot( + pat_df, group, metric, + custom_order, plot_type="line" + ) + else: + print(f"No data for {group} in patient {pat}.") + print("\n") + + print("::::\n") # Level 3 End + print(":::::\n") # Level 2 End ← moved OUTSIDE patient loop + + else: + # Standard View (>10 patients) + print("::::: {.panel-tabset}\n") # Level 2 Start + for group in group_opts_all: + print(f"### {group}\n") + create_diversity_plot(plot_df, group, metric) + print("\n") + print(":::::\n") # Level 2 End + +print("::::::\n") # Level 1 End + +``` +**Figure 1. Comparative Analysis of TCR Repertoire Diversity Metrics.** Boxplots display the distribution of metrics across samples. The interactive panel allows toggling between diversity indices and grouping variables. This visualization facilitates the assessment of repertoire diversity and evenness across different experimental conditions and biological replicates. + +**Use the tabs to toggle between metadata groupings:** + +- The distributions of **Batches** L0, L1, etc., should largely overlap. If L0 is consistently higher than L1 across all sample types, you have a Batch Effect. You cannot compare samples across batches without computational correction (we do not provide such correction). + +- If one **patient** (e.g. Patient01) has drastically lower diversity than all others across all timepoints and tissues, it is likely a biological anomaly or a collection issue (degraded sample). Treat this patient as a separate cohort or exclude them. + +- In longitudinal studies, diversity often drops slightly over time due to treatment. Diversity metrics rarely jump from 0.8 to 0.1 (or vice versa) in a short window **between timepoints**. If so, that sample is suspect. There might be some causes that could explain it: 1) Sample swap of sample timepoint labels for that patient. 2) The patient contracted a viral infection at T2, causing a massive expansion of non-tumor clones. + +## Hill diversity Profile + +It allows to see how "TCR diversity" changes depending on how much you care about rare clones vs. dominant clones. You use Hill Profiles to catch artifacts that single metrics hide. +$q = 0$ (Richness): Ignore frequency, every unique CDR3 counts as 1. +$q = 1$ (Shannon Entropy): Clones are weighted by their frequency. +$q = 2$ (Simpson): Rare clones are mathematically ignored. The top expanded clones dictate the score. + +Two patients might share the exact same Shannon Entropy score, yet one has a healthy, broad repertoire while the other is battling a massive leukemia clone, for example. The Hill Diversity Profile reveals this difference by visualizing the trade-off between Richness (total number of clones) and Evenness (how equally those clones are distributed) in a single curve. + +**How to Read the Curve?** +Think of the **X-axis** ($q$) as a **"sensitivity dial" for dominance**. At the far left ($q=0$), the metric cares only about presence; it counts every unique clone equally, regardless of whether it appears once or a million times. This point represents your total Richness. **As you move right toward $q=2$, the mathematical weight shifts drastically toward the most abundant clones**, effectively ignoring the rare ones. **The slope of the line tells you the story**: a flat line indicates a perfectly even population where everyone is equal, while **a steep**, crashing slope **reveals a population dominated by a few massive "bully" clones.** + +**Here are 2 possible scenarios:** + +- **Scenario A:** The "Depth Trap" +Imagine you are **comparing Sample A and Sample B**, and you immediately notice that **Sample A has a much higher Hill $q=0$ value**, which represents richness. You might be **tempted to conclude that Sample A is biologically "more diverse"** and has a healthier T cell repertoire. However, when you look at the **Hill Profile as it moves toward $q=1$ and $q=2$**, you see the **lines rapidly converge or even cross, showing no difference in the dominant clones.** This reveals that the **"diversity" in Sample A was an illusion caused by sequencing depth.** You simply spent more money sequencing Sample A, allowing the machine to pick up more rare, single-copy clones that were missed in Sample B. The structure of the immune system was identical in both; one was just measured with a magnifying glass while the other was measured with the naked eye. + +- **Scenario B:** Real Biology +Consider a case where you are analyzing a **"Responder" patient versus a "Non-Responder" in a cancer trial**. If you only looked at a single metric like **Shannon Entropy ($q=1$), the two patients might look identical, leading you to think the treatment had no effect on the repertoire structure.**. But the Hill Profile tells a different story. At $q=0$, the Non-Responder is much higher, indicating a vast, unfocused collection of rare T cells that aren't doing much. As you move to **$q=2$, the lines cross, and the Responder suddenly shoots up**. This crossover reveals that **while the Non-Responder has more "types" of cells, the Responder has successfully expanded a specific "army" of effector clones** to fight the tumor**. The profile proves that their immune systems are structurally opposite, a critical biological insight that a single summary statistic would have completely hidden. + +```{python} +#| output: asis +#| fig-width: 6 +#| fig-height: 4.5 +#| label: hill + +# --- 1. Setup Data & Grouping --- +hill_cols = {'hill_q0': 0, 'hill_q1': 1, 'hill_q2': 2} + +exclude_cols = ['sample', 'file', 'total_counts', 'filename', 'sample_id', 'alias', + 'shannon_entropy', 'gini_coefficient', 'inverse_simpson', + 'hill_q0', 'hill_q1', 'hill_q2', 'mean_len', 'std_len'] + +group_opts = [ + col for col in plot_df.columns + if col not in exclude_cols and plot_df[col].nunique() < 35 +] + +sns.set_theme(style="whitegrid", context="paper", font_scale=1.1) + +# --- 2. Start Quarto Tabs --- +print("::: {.panel-tabset}") + +for col in group_opts: + print(f"## {col}") + + # Data Prep + subset = plot_df[[col, 'alias', 'hill_q0', 'hill_q1', 'hill_q2']].copy() + + melted = subset.melt( + id_vars=['alias', col], + value_vars=['hill_q0', 'hill_q1', 'hill_q2'], + var_name='q_metric', + value_name='Diversity' + ) + melted['q'] = melted['q_metric'].map(hill_cols) + melted[col] = melted[col].astype(str) + + # --- Y-AXIS SYNCHRONIZATION --- + # Calculate global min/max for this tab to lock axes + y_min = melted['Diversity'].min() + y_max = melted['Diversity'].max() + + # Add 5% log-padding so points aren't cut off + log_min = np.log10(y_min) + log_max = np.log10(y_max) + pad = (log_max - log_min) * 0.05 + + # Matplotlib limits (Linear values) + mpl_ylim = (10**(log_min - pad), 10**(log_max + pad)) + + # Plotly limits (Log10 values) + plotly_yrange = [log_min - pad, log_max + pad] + + # --- COLOR CONSISTENCY --- + unique_cats = sorted(melted[col].unique()) + palette_tuples = sns.color_palette("viridis", n_colors=len(unique_cats)) + hex_colors = [f'#{int(r*255):02x}{int(g*255):02x}{int(b*255):02x}' for r,g,b in palette_tuples] + color_map = dict(zip(unique_cats, hex_colors)) + + # =========================== + # PLOT 1: SEABORN (Summary) + # =========================== + plt.figure(figsize=(6, 4.5)) # Standard size + + sns.lineplot( + data=melted, + x='q', y='Diversity', + hue=col, style=col, + markers=True, dashes=False, + linewidth=2.5, + palette=color_map, + err_style='band', errorbar=('ci', 95) + ) + + plt.yscale('log') + plt.ylim(mpl_ylim) + plt.xticks([0, 1, 2], ['q=0\n(Richness)', 'q=1\n(Shannon)', 'q=2\n(Simpson)']) + plt.title(f'Average Profile by {col} (with 95% CI)') + plt.ylabel('Effective Clones (Log Scale)') + plt.xlabel('') + + plt.legend(bbox_to_anchor=(1.02, 1), loc='upper left', title=col, + fontsize='x-small', title_fontsize='small', frameon=False) + plt.tight_layout() + plt.show() + + # =========================== + # PLOT 2: PLOTLY (Interactive) + # =========================== + + fig = px.line( + melted, + x='q', + y='Diversity', + color=col, + line_group='alias', # Keep sample here to ensure lines are grouped by the unique identifier + hover_name='alias', + + hover_data={col: True, 'q': False, 'Diversity': ':.2f', 'alias': True}, + log_y=True, + title=f'Individual Sample Profiles (Hover to Identify)', + color_discrete_map=color_map + ) + + fig.update_traces( + mode="lines+markers", + line=dict(width=3), + opacity=0.6, + marker=dict(size=6) + ) + + fig.update_layout( + template="simple_white", + width=600, + height=450, + xaxis=dict( + tickmode='array', + tickvals=[0, 1, 2], + ticktext=['q=0 (Richness)', 'q=1 (Shannon)', 'q=2 (Simpson)'], + title='Sensitivity (q)' + ), + yaxis=dict( + title='Effective Clones (Log Scale)', + range=plotly_yrange + ), + legend=dict(title=col) + ) + + fig.show() + + print("\n") + +print(":::") + +``` +**Figure 2. Hill Diversity Profiles of TCR Repertoires.** Diversity profiles display the effective number of clones (log scale) across sensitivity orders $q=0$ (Richness), $q=1$ (Shannon), and $q=2$ (Simpson). The top panel illustrates the mean repertoire structure for each subject with 95% confidence intervals, highlighting differences in evenness and dominance. The bottom panel visualizes individual sample trajectories, revealing heterogeneity and outlier profiles within each cohort. + +**Using the above metrics to identify technical failures 🚩** +Before attempting biological interpretation, you must validate that your metrics reflect immune biology and not library preparation artifacts. Use these metrics to flag and exclude compromised samples. + +**1. Amplification Bias** +**Symptom:** A sample with moderate Richness ($q=0$) but extreme Dominance ($q=2$). +**Metric Signature:** + +- Shannon ($H$): Disproportionately low compared to the cohort baseline. +- Gini Coefficient: Approaches 1.0 (>0.9 is suspicious in non-tumor samples). +- Hill Profile: The curve starts high at $q=0$ (indicating presence of unique sequences) but crashes vertically as $q \to 2$. + +**Diagnosis:** A "Jackpot" event occurred where a random RNA molecule was preferentially amplified during PCR. The "richness" is likely sequencing noise/errors, while the reads are consumed by a single artifact. +**Action:** Discard sample. + +**2. Sequencing Bias** +**Symptom:** Two samples appear to have different diversities, but the difference disappears when you ignore rare clones. +**Metric Signature:** + +- Hill Profile: Sample A has a much higher $q=0$ (Richness) than Sample B, but the lines converge at $q=1$ and $q=2$. + +**Diagnosis:** This is not a biological difference. Sample A was simply sequenced deeper, detecting more single-copy clones (noise). The structural effective diversity is identical. +**Action:** Do not claim Sample A is "more diverse." Rely on Shannon ($q=1$) or Simpson ($q=2$) for comparison, or downsample (rarefy) the libraries to equal depth. + +**3. Library Failure** +**Symptom**: A sample that looks "empty" compared to the cohort. +**Metric Signature:** + +- Richness ($R$): Drastically lower than the cohort mean. +- Hill Profile: A flat line hovering near the bottom of the Y-axis. +**Diagnosis:** Poor RNA extraction, degradation, or a biopsy that captured mostly fat/stroma (low T-cell content). +**Action:** Exclude from analysis. + +**4. Batch Effects** +**Symptom:** Diversity metrics cluster by processing date rather than biological group. +**Metric Signature:**Visual Check: Boxplots of Shannon Entropy grouped by Batch_ID. If Batch 0 is consistently higher than 1 across all sample types (tumor, blood, healthy), you have a technical batch effect. +**Action:** You cannot compare raw metrics across these batches. + + +**Biological Interpretation of the above metrics 🦠** +Once QC is cleared, these metrics describe the shape of the immune repertoire. + +**1. Naive vs. Reactive Profiles** +**Naive / Resting State:** + +- Signature: High Richness, High Shannon, Low Gini (< 0.3). +- Interpretation: The army is standing by. There is no dominant clone because no specific threat has triggered an expansion. The Hill profile will show a "gentle decline." +**Reactive / Tumor Infiltrating:** + +- Signature: Low Shannon, High Inverse Simpson, High Gini (> 0.6). +- Interpretation: Clonal Expansion. The immune system has "picked a winner." A few specific clones (e.g., tumor-reactive or viral-specific) have expanded to occupy a large fraction of the repertoire (20–50%). + +**2. Inverse Simpson** +The Inverse Simpson Index ($1/D$) converts the probability of coincidence into a concrete "count" of effective clones. It filters out the long tail of rare sequences to quantify the number of clones driving the response. +**High Effective Number** ($1/D \gg 100$): + +- Signature: $1/D$ is close to Richness ($R$; i.e. hill_q0). +- Interpretation: Polyclonal/Diverse. The sample behaves as if it is composed of hundreds or thousands of equally abundant clones. Resources are distributed broadly. + +**Low Effective Number ($1/D < 20$):** + +- Signature: $1/D$ is a small fraction of Richness ($R$). +- Interpretation: Oligoclonal/Focused. Despite potentially containing thousands of unique sequences ($R$), the sample behaves biologically as if it only contains ~10-20 active clones. This confirms a highly focused effector response. + +**3. Hill Profile** +This is the most powerful visualization for distinguishing Responder vs. Non-Responder dynamics in longitudinal data. +**The Scenario:** A Non-Responder may have higher Richness ($q=0$) due to a chaotic, unfocused repertoire. A Responder might have lower Richness but higher Dominance. +**The Signature:** The Hill curves will cross.At $q=0$ (Richness): Non-Responder is higher.At $q=2$ (Dominance): Responder shoots up. +**Interpretation**: The crossover proves the Responder has successfully funneled resources into a specific "army" of effector clones, structurally reorganizing the repertoire to fight the tumor. + +## Combining metrics {#sec-combine-metrics} + +### Evenness Check (Gini vs Shannon) + +Under normal biological conditions, the Gini coefficient and Shannon entropy should exhibit a strong negative correlation. As the immune system focuses on a specific threat, a few clones expand to dominate the repertoire, causing inequality to rise (High Gini) while overall diversity naturally falls (Low Shannon). + +```{python} +#| output: asis +#| fig-width: 6 +#| fig-height: 5 + +# 1. Setup Grouping Options +exclude_cols = ['sample', 'file', 'total_counts', 'timepoint_order', 'filename', 'sample_id', + 'shannon_entropy', 'gini_coefficient', 'inverse_simpson', + 'hill_q0', 'hill_q1', 'hill_q2', 'mean_len', 'std_len'] + +group_opts = [ + col for col in plot_df.columns + if col not in exclude_cols and plot_df[col].nunique() < 35 +] + +# 2. Start Quarto Tabset +print("::: {.panel-tabset}") + +for col in group_opts: + print(f"## {col}") + + fig = px.scatter( + plot_df, + x='shannon_entropy', + y='gini_coefficient', + color=col, + symbol=col, + template="simple_white", + width=600, + height=500, + title=f'Evenness by {col}', + hover_data={ + 'sample': True, + 'alias': True if 'alias' in plot_df.columns else False, + 'shannon_entropy': ':.3f', + 'gini_coefficient': ':.3f', + col: True + } + ) + + # Make points look like your seaborn version + fig.update_traces( + marker=dict(size=10, line=dict(width=1, color='black')), + opacity=0.8 + ) + + fig.update_layout( + xaxis_title='Shannon (Diversity)', + yaxis_title='Gini (Inequality)', + legend_title=col + ) + + fig.show() + print("\n") + +print(":::") + +``` + +**Figure 3. Evaluation of Repertoire Structure via Gini-Shannon Correlation.** This scatter plot contrasts the Gini coefficient (inequality) against Shannon entropy (diversity) for individual samples, colored by metadata variables. The visualization assesses the expected inverse relationship between the two metrics, where antigen-driven clonal expansion typically results in higher inequality and reduced diversity. + +Healthy, **naive, or polyclonal samples** will naturally cluster in the **bottom-right** of the plot. These repertoires are **characterized by high diversity** ($H > 8$) and low inequality ($G < 0.2$), reflecting a deep pool of clones with relatively even frequencies. In contrast, **samples undergoing an active immune response**—such as those from tumors, acute infections, or leukemia—will drift toward the **top-left**. In this "Expansion Zone," the biological dominance of a few clones drives the Gini coefficient up ($> 0.6$) and suppresses the entropy, a signature that is expected in disease states but signals potential PCR bias ("jackpotting") if observed in negative controls. + +The most critical **QC artifact to watch** for is the "Ghost Library" in the **bottom-left corner**. These samples show **low diversity but paradoxically low inequality**. This almost always indicates severe **undersampling**: if a library only contains 50 reads and every read is unique, the sample will appear "perfectly equal" (Gini $\approx$ 0) simply because it lacks the depth to reveal the true distribution. These samples are statistical phantoms and should be **discarded**. + +Finally, samples falling in the **top-right (High Diversity + High Inequality)** represent a **mathematical contradiction for T-cell data**. It is impossible to have a massive, diverse "tail" of clones while simultaneously having a single clone occupy the majority of reads. Samples in this region often indicate **bioinformatic errors**, such as merging incompatible library files or calculating metrics on raw counts rather than frequencies. + +### Richness vs 1/D + +The Inverse Simpson ($1/D$) is used to visualize the magnitude of the response. It ignores the long tail of rare clones. It is a good practice to compare $1/D$ to Richness ($R$). +- If $R = 5000$ and $1/D = 4500$: The sample is polyclonal (many clones of equal size). +- If $R = 5000$ and $1/D = 2$: The sample is oligoclonal. Despite having 5000 unique sequences, it behaves biologically as if it only has 2 clones. + +```{python} +#| output: asis +#| fig-width: 6 +#| fig-height: 5 + +# 1. Setup Grouping Options +exclude_cols = ['sample', 'file', 'total_counts', 'timepoint_order', 'filename', 'sample_id', + 'shannon_entropy', 'gini_coefficient', 'inverse_simpson', + 'hill_q0', 'hill_q1', 'hill_q2', 'mean_len', 'std_len'] + +# Automatically detect categorical columns with reasonable cardinality +group_opts = [ + col for col in plot_df.columns + if col not in exclude_cols and plot_df[col].nunique() < 35 +] + +# ============================================================================== +# PLOT SET 1: RICHNESS vs INVERSE SIMPSON (The "Depth Trap" Check) +# ============================================================================== +print("::: {.panel-tabset}") + +for col in group_opts: + print(f"## {col}") + + fig = px.scatter( + plot_df, + x='inverse_simpson', + y='hill_q0', + color=col, + symbol=col, + template="simple_white", + width=600, + height=500, + title=f'Richness vs. Effective Clones by {col}', + hover_data={ + 'sample': True, + 'alias': True if 'alias' in plot_df.columns else False, + 'inverse_simpson': ':.3f', + 'hill_q0': ':.3f', + col: True + } + ) + + # Match seaborn styling + fig.update_traces( + marker=dict(size=10, line=dict(width=1, color='black')), + opacity=0.8 + ) + + # --- Add diagonal line (y = x) --- + max_val = max( + plot_df['hill_q0'].max(), + plot_df['inverse_simpson'].max() + ) + + fig.add_trace(go.Scatter( + x=[0, max_val], + y=[0, max_val], + mode='lines', + line=dict(dash='dash', color='gray'), + name='Perfect Evenness (y=x)', + hoverinfo='skip' # keeps hover clean + )) + + # Layout + fig.update_layout( + xaxis_title='Inverse Simpson (Effective # Clones)', + yaxis_title='Hill q=0 (Total Richness)', + legend_title=col + ) + + fig.show() + print("\n") + +print(":::") +print("\n") + +``` + +**Figure 4. Analysis of Clonal Dominance vs. Effective Clone Count.** This scatter plot compares the total clonal richness against the effective number of clones, colored by metadata variables. The diagonal line ($y=x$) represents a theoretical state of perfect evenness, where every clone has equal frequency. Deviations from this line indicate the degree of clonal dominance. + +### Gini Coefficient vs. Effective Clones + +This plot maps the trade-off between unevenness and effective size. Typically, these metrics are inversely correlated: as inequality (Gini) rises, the number of effective clones (Inverse Simpson) drops. Deviations from this curve reveal samples with unusual structural properties + +```{python} +#| output: asis +#| fig-width: 6 +#| fig-height: 5 + +print("::: {.panel-tabset}") + +for col in group_opts: + print(f"## {col}") + + # Create Plotly scatter + fig = px.scatter( + plot_df, + x='inverse_simpson', + y='gini_coefficient', + color=col, + symbol=col, + template="simple_white", + width=600, + height=500, + title=f'Inequality vs. Diversity by {col}', + hover_data={ + 'sample': True, + 'alias': True if 'alias' in plot_df.columns else False, + 'inverse_simpson': ':.3f', + 'gini_coefficient': ':.3f', + col: True + }, + color_continuous_scale='magma' if plot_df[col].dtype != object else None + ) + + # Match marker style from seaborn + fig.update_traces( + marker=dict(size=10, line=dict(width=1, color='black')), + opacity=0.8 + ) + + # Fix y-axis for Gini (bounded [0,1]) + fig.update_yaxes(range=[0, 1.05], title='Gini Coefficient (Inequality)') + + # x-axis label + fig.update_xaxes(title='Inverse Simpson (Effective # Clones)') + + # Layout adjustments + fig.update_layout( + legend_title=col, + title=dict(font=dict(size=14)), + width=600, + height=500 + ) + + fig.show() + print("\n") + +print(":::") + +``` + +**Figure 5. Relationship Between Repertoire Inequality and Effective Diversity**. The visualization demonstrates the inverse correlation typically observed in immune repertoires, where high inequality (dominance by a few clones) corresponds to a low effective population size. Deviations from this structural trend can identify samples with aberrant clonal distributions. + +**Combined metrics: Quality Control 🚩** +Use these 2D plots to distinguish between sequencing artifacts and true biological signal. + +**1. Richness vs. Effective Clones** +This plot exposes the "useless tail" of your library—the sequences that appear only once or twice and do not contribute to the effective immune response. + +- The Diagonal Line ($y=x$): A sample on this line has no rare clones; every unique sequence is equally abundant. +- Vertical Separation (The Warning): If Sample A is far above Sample B on the Y-axis (Richness) but they align perfectly on the X-axis (Effective Clones), Sample A is not more diverse. It was simply sequenced deeper. The extra richness is noise (singletons). +- Action: Trust the X-axis (Inverse Simpson) for biological comparison. Ignore the Y-axis offset. Samples clustered in the bottom-left corner ($<100$ on both axes) often indicate failed library preparation or low-input samples (e.g., fibrous tissue with no T-cells). + +**2. Gini vs. Effective Clones Plot** +This plot validates the structural integrity of the library. + +- **The "Impossible" Quadrant:** You should rarely see High Gini (>0.8) combined with High Effective Clones (>1000). Mathematically, you cannot have extreme inequality while simultaneously maintaining thousands of effectively equal clones. If you see this, check for chimeric reads or alignment errors. + +- **The "False Dominance":** If a Baseline/PBMC sample appears in the Top-Left corner (High Gini, Low Effective Count), it is a red flag. Unless the patient has a known blood malignancy, a resting control sample should never look like a tumor. This indicates a PCR Jackpot artifact. + + +**Combined metrics: Biological Interpretation 🦠 ** +Once QC is cleared, the position of a sample in these 2D spaces defines its immunological state. + +**1. Gini vs. Effective Clones (1/D)** +**The Naive/Polyclonal State (Bottom-Right):** +The immune system is in surveillance mode. No specific antigen has triggered a response. Resources are distributed broadly across thousands of clones. +**The Reactive/focused State (Top-Left):** +Clonal Expansion. The immune system has identified a threat (tumor/virus) and "picked a winner." The repertoire is dominated by a focused army of effector clones, physically crowding out the diversity. + +**2. Tracking Response** +In a successful immunotherapy response, track the patient's movement across the plot over time: +**The "Focusing" Shift:** A responder should move diagonally Up and Left. They start with a broad repertoire (Bottom-Right) and shift toward dominance (Top-Left) as tumor-reactive clones expand. +**The "Ignorance" Stasis:** Non-responders often remain stuck in the Bottom-Right quadrant (high diversity, low inequality), indicating the treatment failed to trigger T-cell expansion. + +## TCR Overlap {#sec-tcr-overlap} + +We use two distinct metrics to compare repertoires. One measures membership (who is there?), and the other measures structure (who is dominant?). + +- **Jaccard Index:** Measures the overlap of unique sequences, ignoring frequency. It treats a singleton and a top clone equally. +**Question:** "What fraction of the unique CDR3s are found in both samples?" +**Range:** $0$ (No shared sequences) $\to$ $1$ (Identical unique sequence list). +**Formula:** +$$J(A,B) = \frac{|A \cap B|}{|A \cup B|}$$ +Where $|A \cap B|$ is the count of unique clonotypes shared by both samples, and $|A \cup B|$ is the total count of unique clonotypes in the union of both samples. + +```{python} +#| output: asis +#| fig-width: 7 +#| fig-height: 7 + +# --- 1. CALCULATION FUNCTION (Only Jaccard) --- +def calculate_overlaps(df): + samples = sorted(df['sample'].unique()) + clones = sorted(df['CDR3b'].unique()) + + sample_map = {s: i for i, s in enumerate(samples)} + clone_map = {c: i for i, c in enumerate(clones)} + + row_idx = df['sample'].map(sample_map).values + col_idx = df['CDR3b'].map(clone_map).values + values = df['counts'].values + + n_samples = len(samples) + n_clones = len(clones) + + mat = sparse.coo_matrix((values, (row_idx, col_idx)), shape=(n_samples, n_clones)).tocsr() + + # Jaccard (Binary) + mat_bin = mat.copy() + mat_bin.data[:] = 1 + intersection = mat_bin.dot(mat_bin.T).toarray() + row_sums = mat_bin.sum(axis=1).A1 + union_matrix = row_sums[:, None] + row_sums[None, :] - intersection + + jaccard = np.divide(intersection, union_matrix, + out=np.zeros_like(intersection, dtype=float), + where=union_matrix!=0) + + return pd.DataFrame(jaccard, index=samples, columns=samples) + + +# --- 2. EXECUTE CALCULATION --- +jaccard_df = calculate_overlaps(concat_df) + +# --- 3. RENAME TO ALIAS --- +if alias_col in meta.columns: + sample_to_alias = meta.set_index('sample')[alias_col].to_dict() + jaccard_df = jaccard_df.rename(index=sample_to_alias, columns=sample_to_alias) + +# --- 4. CLUSTERING HELPER --- +def get_clustered_order(matrix): + """Return hierarchical clustering order of samples.""" + dist = pdist(matrix.values, metric='cityblock') + link = linkage(dist, method='average') + leaves = leaves_list(link) + return matrix.index[leaves].tolist() + +# --- 5. INTERACTIVE PLOTTING FUNCTION --- +def plot_interactive_heatmap(matrix, title): + try: + ordered_labels = get_clustered_order(matrix) + matrix = matrix.loc[ordered_labels, ordered_labels] + except Exception as e: + print(f"Clustering failed, using default order: {e}") + + vals = matrix.values + mask = ~np.eye(vals.shape[0], dtype=bool) + max_val = vals[mask].max() if mask.any() else 1.0 + if max_val < 0.05: max_val = 0.1 + + fig = px.imshow( + matrix, + labels=dict(x="Sample", y="Sample", color="Jaccard"), + x=matrix.columns, + y=matrix.index, + color_continuous_scale="Reds", + range_color=[0, max_val], + title=f"{title}
(Max Overlap: {max_val:.3f})" + ) + + fig.update_layout( + width=700, + height=700, + xaxis=dict(tickangle=90, tickfont=dict(size=10)), + yaxis=dict(tickfont=dict(size=10)), + plot_bgcolor='black' + ) + fig.update_traces(xgap=1, ygap=1) + fig.show() + +# --- 6. OUTPUT TABSET (Only Jaccard) --- +print("::: {.panel-tabset}\n") +print("## Jaccard\n") +plot_interactive_heatmap(jaccard_df, "Jaccard Index") +print("\n") +print(":::\n") + +``` +**Figure 6. Assessment of Clonal Overlap and Cross-Contamination.** +Heatmaps display the pairwise similarity between TCR repertoires using the Jaccard (overlap) index. Hierarchical clustering groups samples based on their shared clonal content, facilitating the identification of biological replicates or potential cross-contamination events. The color intensity reflects the degree of similarity, where darker red indicates a higher proportion of shared clonotypes between sample pairs. + +::: {.callout-tip title="Note"} +When calculating Jaccard for this specific purpose (contamination), we are calculating it on the nucleotide sequence, not the amino acid sequence, **whenever the TCRseq was done using DNA as input**. Convergent recombination (different DNA making the same Amino Acid) is common in high-depth samples and will give you a false positive "contamination" signal if you stick to Amino Acids. Contamination should be measured at the physical level (DNA). +::: + +**Technical QC 🚩** +Use these overlap metrics to detect the two most dangerous errors in sequencing: **Sample Swapping and Physical Contamination.** + +- **Contamination Detection: Jaccard Index** + - Red Flag: A block of red squares connecting unrelated samples (Overlap $> 0.05 - 0.1$ Pink/Red). + - Diagnosis: Cross-Contamination. Because TCRs are hypervariable, distinct individuals should almost never share rare sequences (Jaccard should be near 0). + - Scenario A (Splash): Patient_A and Patient_B share 20% of their unique sequences. Physical liquid likely splashed between wells during PCR setup. + - Scenario B (Barcode Hopping): If an entire batch shows faint overlap (~5-10%) with a high-concentration sample, it indicates index hopping during sequencing. + - Action: Flag samples with unexpected Jaccard scores $>0.05$. + + +**Biological Interpretation 🦠** +Once technical artifacts are ruled out, overlap quantifies stability and shared biology. + +- **Public Clones: Jaccard** + - Low but Real Overlap (< 0.01):While distinct humans rarely share repertoires, they may share identical "Public Clones" (often viral-specific, e.g., CMV or EBV). + - Interpretation: A very faint Jaccard signal (0.001 - 0.01) across a cohort can indicate a shared antigen exposure (e.g., a common viral infection in the population) or convergent recombination events. + + +# Technical QC Relevance (🚩) {#sec-qc} +The following section presents metrics designed to assess the processing quality of the samples. These technical QC checks are essential for differentiating between genuine biological signals and technical artifacts, including library failures, undersampling, batch effects, and outliers. However, we acknowledge that some of these metrics can also carry biological meaning, depending on the specific biological context of the analysis. + +## Percent of productive TCRs {#sec-cdr3-prod} + +The percent of productive TCRs is one of the first metrics you should check to validate your data's reliability before proceeding to any biological analysis. + +As a Quality Filter: You should set a QC threshold (e.g., >75% productive reads). Samples that fail to meet this threshold should be flagged for review or potentially excluded from the analysis. Drawing conclusions about T-cell diversity or clonality from a sample with low productivity is unreliable, as the data is likely noisy and not a true representation of the functional T-cell repertoire. + +To Troubleshoot Experiments: If you find that an entire batch of samples has a low productive percentage, it points to a systematic issue in your experimental protocol, most commonly an ineffective gDNA removal step or problems with RNA integrity. + +To Ensure Confidence in Results: By confirming that your samples have a high percentage of productive TCRs, you can be confident that the clonotypes you identify and analyze represent real, functional T-cells that are actively participating in the immune response. + +```{python} +#| output: asis +#| echo: false + +import plotly.express as px +import plotly.graph_objects as go +from scipy.stats import mannwhitneyu +from itertools import combinations +import pandas as pd + +# Create a mapping for timepoint order using the 'meta' dataframe +if 'meta' in locals() and timepoint_col in meta.columns and timepoint_order_col in meta.columns: + time_order_map = meta.set_index(timepoint_col)[timepoint_order_col].to_dict() +else: + # Fallback: if 'meta' isn't loaded, try to use columns in the main df + if timepoint_order_col in df.columns: + time_order_map = df.set_index(timepoint_col)[timepoint_order_col].to_dict() + else: + time_order_map = {} # No sorting map available + +# 1. Define General Exclusions +exclude_cols = ['sample', 'file', 'clonality', 'counts', 'total_counts', + 'num_clones', 'num_TCRs', 'simpson_index', 'simpson_index_corrected', + 'num_prod', 'num_nonprod', 'pct_prod', 'pct_nonprod', + 'productive_cdr3_avg_len', 'num_convergent','ratio_convergent'] + +# 2. Define Groups for "All Samples" Tab (Exclude timepoint/order) +exclude_from_all_samples = [timepoint_col, timepoint_order_col, 'protocol_day', alias_col] + +group_opts_all = [ + col for col in df.columns + if col not in exclude_cols + and col not in exclude_from_all_samples + and df[col].nunique() < 35 +] + +# 3. Define Groups for "By Patient" Tab (Only timepoint) +group_opts_patient = [timepoint_col] if timepoint_col in df.columns else [] + +unique_patients = sorted(df[subject_col].dropna().unique().tolist()) + +# --- 1. DEFINE PLOTTING FUNCTION --- +def create_pct_prod_plot(data, group_col, custom_order=None, plot_type="box"): + """ + custom_order: Optional list of x-axis values in the desired order. + plot_type: "box" for standard comparisons, "line" for patient timelines. + """ + if data.empty or data[group_col].dropna().empty: + return + + # Determine order: Use custom list if provided, otherwise sort alphabetically + if custom_order: + unique_cats = [x for x in custom_order if x in data[group_col].unique()] + else: + unique_cats = sorted(data[group_col].dropna().unique().tolist()) + + # --- BASE PLOT --- + if plot_type == "box": + fig = px.box( + data, + x=group_col, + y='pct_prod', + color=group_col, + points='all', + hover_data=['alias'], + category_orders={group_col: unique_cats}, + template="simple_white", + title=f"Productive TCRs (%) by {group_col}" + ) + + fig.update_traces(width=0.5, marker=dict(size=10, opacity=0.7)) + + for trace in fig.data: + if isinstance(trace, go.Box): + trace.pointpos = 0 + trace.jitter = 0.2 + + elif plot_type == "line": + # Sort for proper line connection + data = data.copy() + data[group_col] = pd.Categorical(data[group_col], categories=unique_cats, ordered=True) + data = data.sort_values(group_col) + + fig = px.scatter( + data, + x=group_col, + y='pct_prod', + color_discrete_sequence=["black"], + hover_data=['alias'], + template="simple_white", + title=f"Productive TCRs (%) by {group_col}" + ) + + # Add dotted line + fig.add_trace(go.Scatter( + x=data[group_col], + y=data['pct_prod'], + mode='lines+markers', + line=dict(dash='dot', width=2), + marker=dict(size=18, color='black'), + showlegend=False + )) + + # --- STATISTICS (ONLY FOR BOXPLOTS) --- + stack_counter = 0 + + if plot_type == "box" and len(unique_cats) >= 2: + pairs = list(combinations(unique_cats, 2)) + + y_max = data['pct_prod'].max() + y_min = data['pct_prod'].min() + y_range = y_max - y_min + if y_range == 0: + y_range = 1 + + # Updated to tighter spacing format + base_offset = y_range * 0.08 + step_size = y_range * 0.06 + text_offset = y_range * 0.015 + + for t1, t2 in pairs: + group1 = data[data[group_col] == t1]['pct_prod'].dropna() + group2 = data[data[group_col] == t2]['pct_prod'].dropna() + + if len(group1) < 2 or len(group2) < 2: + continue + + try: + stat, p_value = mannwhitneyu(group1, group2, alternative='two-sided') + except ValueError: + continue + + if p_value >= 0.05: + continue + + if p_value < 0.001: symbol = '***' + elif p_value < 0.01: symbol = '**' + else: symbol = '*' + + y_bracket = y_max + base_offset + (stack_counter * step_size) + y_text = y_bracket + text_offset + + # Draw Bracket Line + fig.add_shape( + type="line", xref="x", yref="y", + x0=t1, y0=y_bracket, x1=t2, y1=y_bracket, + line=dict(color="black", width=1.5) + ) + + # Draw Bracket Ticks + tick_len = y_range * 0.015 + fig.add_shape(type="line", xref="x", yref="y", + x0=t1, y0=y_bracket, x1=t1, y1=y_bracket - tick_len, + line=dict(color="black", width=1.5)) + fig.add_shape(type="line", xref="x", yref="y", + x0=t2, y0=y_bracket, x1=t2, y1=y_bracket - tick_len, + line=dict(color="black", width=1.5)) + + try: + x_center = (unique_cats.index(t1) + unique_cats.index(t2)) / 2 + except: + x_center = 0 + + # Add P-Value Text + fig.add_annotation( + x=x_center, y=y_text, text=f"{symbol}
p={p_value:.3f}", + showarrow=False, font=dict(size=10, color="black") + ) + + stack_counter += 1 + + extra_space = y_range * (0.15 + 0.07 * stack_counter) + fig.update_yaxes(range=[y_min, y_max + extra_space]) + + # --- FINAL LAYOUT --- + dynamic_height = 600 + (stack_counter * 40) + + fig.update_layout( + xaxis_title=group_col, + yaxis_title="Productive TCRs (%)", + xaxis=dict(tickfont=dict(size=15)), + margin=dict(t=80), + showlegend=False, + width=700, + height=dynamic_height, + plot_bgcolor='rgba(0,0,0,0)' + ) + + fig.update_yaxes(showgrid=True, gridcolor='lightgrey') + fig.show() + +# --- 2. GENERATE TABS --- + +if len(unique_patients) <= 10: + print("::::: {.panel-tabset}\n") + + # --- TAB A: ALL SAMPLES --- + print("## All Samples\n") + print(":::: {.panel-tabset}\n") + for group in group_opts_all: + print(f"### {group}\n") + create_pct_prod_plot(df, group, plot_type="box") + print("\n") + print("::::\n") + + # --- TAB B: BY PATIENT --- + print("## By Patient\n") + print(":::: {.panel-tabset}\n") + + for pat in unique_patients: + print(f"### {pat}\n") + + pat_df = df[df[subject_col] == pat] + + if 'origin' in pat_df.columns and pat_df['origin'].nunique() > 1: + print(":::: {.panel-tabset}\n") # Level 4 Start (Origins) + + origins = sorted(pat_df['origin'].dropna().unique()) + + for origin in origins: + print(f"#### {origin}\n") + + origin_df = pat_df[pat_df['origin'] == origin] + + for group in group_opts_patient: + if origin_df[group].nunique() > 0: + custom_order = None + if group == timepoint_col and time_order_map: + pat_tps = origin_df[group].dropna().unique().tolist() + custom_order = sorted( + pat_tps, + key=lambda x: time_order_map.get(x, 999) + ) + + create_pct_prod_plot( + origin_df, + group, + custom_order, + plot_type="line" + ) + else: + print(f"No data for {group} in {origin}, patient {pat}.") + print("\n") + + print("::::\n") # Level 4 End + + else: + # Fallback if no origin variation + for group in group_opts_patient: + if pat_df[group].nunique() > 0: + + # --- APPLY CUSTOM SORTING IF AVAILABLE --- + custom_order = None + if group == timepoint_col and time_order_map: + pat_tps = pat_df[group].dropna().unique().tolist() + custom_order = sorted(pat_tps, key=lambda x: time_order_map.get(x, 999)) + + create_pct_prod_plot( + pat_df, + group, + custom_order, + plot_type="line" + ) + else: + print(f"No data for {group} in patient {pat}.") + print("\n") + + print("::::\n") + print(":::::\n") + +else: + # --- STANDARD VIEW (>10 Patients) --- + print("::: {.panel-tabset}\n") + for group in group_opts_all: + print(f"## {group}\n") + create_pct_prod_plot(df, group, plot_type="box") + print("\n") + print(":::\n") + +``` + +**Figure 7. Percent of productive TCRs.** Distribution of productive TCRs across sample timepoints. A productive TCR is a DNA/RNA sequence that can be translated into a protein sequence, i.e. it does not contain a premature stop codon or an out of frame rearrangement. The percent of productive TCRs is calculated as: + +$$ Percent \text{ } productive \text{ } TCRs = \frac P N $$ + +where $P$ is the number of productive TCRs and $N$ is the total number of TCRs in a given sample. + + +## Average productive CDR3 Length {#sec-cdr3-len} + +The incredible diversity is generated during the development of T-cells through a process of genetic recombination, where different gene segments (V, D, and J) are pieced together. The CDR3 region spans the junction of these segments and includes random nucleotide additions and deletions, creating a unique sequence for almost every T-cell. **The length of this CDR3 region is a direct consequence of this recombination process.** + +The CDR3 length distribution serves as a primary metric for assessing library complexity and detecting technical biases introduced during library preparation or sequencing. While biological clonality impacts this metric, deviations in non-expanded samples often signal workflow failures rather than biological reality. + +```{python} +#| output: asis +#| echo: false + +# Create a mapping for timepoint order using the 'meta' dataframe +if 'meta' in locals() and timepoint_col in meta.columns and timepoint_order_col in meta.columns: + time_order_map = meta.set_index(timepoint_col)[timepoint_order_col].to_dict() +else: + # Fallback: if 'meta' isn't loaded, try to use columns in the main df + if timepoint_order_col in df.columns: + time_order_map = df.set_index(timepoint_col)[timepoint_order_col].to_dict() + else: + time_order_map = {} # No sorting map available + +# 1. Define General Exclusions +exclude_cols = ['sample', 'file', 'clonality', 'counts', 'total_counts', + 'num_clones', 'num_TCRs', 'simpson_index', 'simpson_index_corrected', + 'num_prod', 'num_nonprod', 'pct_prod', 'pct_nonprod', + 'productive_cdr3_avg_len', 'num_convergent','ratio_convergent'] + +# 2. Define Groups for "All Samples" Tab (Exclude timepoint/order) +exclude_from_all_samples = [timepoint_col, timepoint_order_col, 'protocol_day', alias_col] + +group_opts_all = [ + col for col in df.columns + if col not in exclude_cols + and col not in exclude_from_all_samples + and df[col].nunique() < 35 +] + +# 3. Define Groups for "By Patient" Tab (Only timepoint) +group_opts_patient = [timepoint_col] if timepoint_col in df.columns else [] + +unique_patients = sorted(df[subject_col].dropna().unique().tolist()) + +# --- 1. DEFINE PLOTTING FUNCTION --- +def create_cdr3_len_plot(data, group_col, custom_order=None, plot_type="box"): + """ + custom_order: Optional list of x-axis values in the desired order. + plot_type: "box" for standard comparisons, "line" for patient timelines. + """ + if data.empty or data[group_col].dropna().empty: + return + + # Determine order: Use custom list if provided, otherwise sort alphabetically + if custom_order: + unique_cats = [x for x in custom_order if x in data[group_col].unique()] + else: + unique_cats = sorted(data[group_col].dropna().unique().tolist()) + + # --- BASE PLOT --- + if plot_type == "box": + fig = px.box( + data, + x=group_col, + y='productive_cdr3_avg_len', + color=group_col, + points='all', + hover_data=['alias'], + category_orders={group_col: unique_cats}, + template="simple_white", + title=f"CDR3 Length by {group_col}" + ) + + # --- BIGGER DOTS (size=10) --- + fig.update_traces(width=0.5, marker=dict(size=10, opacity=0.7)) + + for trace in fig.data: + if isinstance(trace, go.Box): + trace.pointpos = 0 + trace.jitter = 0.2 + + elif plot_type == "line": + # Sort for proper line connection + data = data.copy() + data[group_col] = pd.Categorical(data[group_col], categories=unique_cats, ordered=True) + data = data.sort_values(group_col) + + fig = px.scatter( + data, + x=group_col, + y='productive_cdr3_avg_len', + color_discrete_sequence=["black"], + hover_data=['alias'], + template="simple_white", + title=f"CDR3 Length by {group_col}" + ) + + # Add dotted line + fig.add_trace(go.Scatter( + x=data[group_col], + y=data['productive_cdr3_avg_len'], + mode='lines+markers', + line=dict(dash='dot', width=2), + marker=dict(size=18, color='black'), + showlegend=False + )) + + # --- STATISTICS (ONLY FOR BOXPLOTS) --- + stack_counter = 0 + + if plot_type == "box" and len(unique_cats) >= 2: + pairs = list(combinations(unique_cats, 2)) + + y_max = data['productive_cdr3_avg_len'].max() + y_min = data['productive_cdr3_avg_len'].min() + y_range = y_max - y_min + if y_range == 0: + y_range = 1 + + # Tighter spacing format + base_offset = y_range * 0.08 + step_size = y_range * 0.06 + text_offset = y_range * 0.015 + + for t1, t2 in pairs: + group1 = data[data[group_col] == t1]['productive_cdr3_avg_len'].dropna() + group2 = data[data[group_col] == t2]['productive_cdr3_avg_len'].dropna() + + if len(group1) < 2 or len(group2) < 2: + continue + + try: + stat, p_value = mannwhitneyu(group1, group2, alternative='two-sided') + except ValueError: + continue + + if p_value >= 0.05: + continue + + if p_value < 0.001: symbol = '***' + elif p_value < 0.01: symbol = '**' + else: symbol = '*' + + y_bracket = y_max + base_offset + (stack_counter * step_size) + y_text = y_bracket + text_offset + + # Draw Bracket Line + fig.add_shape( + type="line", xref="x", yref="y", + x0=t1, y0=y_bracket, x1=t2, y1=y_bracket, + line=dict(color="black", width=1.5) + ) + + # Draw Bracket Ticks + tick_len = y_range * 0.015 + fig.add_shape(type="line", xref="x", yref="y", + x0=t1, y0=y_bracket, x1=t1, y1=y_bracket - tick_len, + line=dict(color="black", width=1.5)) + fig.add_shape(type="line", xref="x", yref="y", + x0=t2, y0=y_bracket, x1=t2, y1=y_bracket - tick_len, + line=dict(color="black", width=1.5)) + + try: + x_center = (unique_cats.index(t1) + unique_cats.index(t2)) / 2 + except: + x_center = 0 + + # Add P-Value Text + fig.add_annotation( + x=x_center, y=y_text, text=f"{symbol}
p={p_value:.3f}", + showarrow=False, font=dict(size=10, color="black") + ) + + stack_counter += 1 + + extra_space = y_range * (0.15 + 0.07 * stack_counter) + fig.update_yaxes(range=[y_min, y_max + extra_space]) + + # --- FINAL LAYOUT --- + dynamic_height = 600 + (stack_counter * 40) + + fig.update_layout( + xaxis_title=group_col, + yaxis_title="Avg CDR3 Length (ntds)", + xaxis=dict(tickfont=dict(size=15)), + margin=dict(t=80), + showlegend=False, + width=700, + height=dynamic_height, + plot_bgcolor='rgba(0,0,0,0)' + ) + fig.update_yaxes(showgrid=True, gridcolor='lightgrey') + fig.show() + +# --- 2. GENERATE TABS --- + +if len(unique_patients) <= 10: + print("::::: {.panel-tabset}\n") + + # --- TAB A: ALL SAMPLES --- + print("## All Samples\n") + print(":::: {.panel-tabset}\n") + for group in group_opts_all: + print(f"### {group}\n") + create_cdr3_len_plot(df, group, plot_type="box") + print("\n") + print("::::\n") + + # --- TAB B: BY PATIENT --- + print("## By Patient\n") + print(":::: {.panel-tabset}\n") + + for pat in unique_patients: + print(f"### {pat}\n") + + pat_df = df[df[subject_col] == pat] + + # --- Origin tab layer --- + if 'origin' in pat_df.columns and pat_df['origin'].nunique() > 1: + print(":::: {.panel-tabset}\n") # Level 4 Start (Origins) + + origins = sorted(pat_df['origin'].dropna().unique()) + + for origin in origins: + print(f"#### {origin}\n") + + origin_df = pat_df[pat_df['origin'] == origin] + + for group in group_opts_patient: + if origin_df[group].nunique() > 0: + custom_order = None + if group == timepoint_col and time_order_map: + pat_tps = origin_df[group].dropna().unique().tolist() + custom_order = sorted( + pat_tps, + key=lambda x: time_order_map.get(x, 999) + ) + + create_cdr3_len_plot( + origin_df, + group, + custom_order, + plot_type="line" + ) + else: + print(f"No data for {group} in {origin}, patient {pat}.") + print("\n") + + print("::::\n") # Level 4 End + + else: + # Fallback if no origin variation + for group in group_opts_patient: + if pat_df[group].nunique() > 0: + + # --- APPLY CUSTOM SORTING IF AVAILABLE --- + custom_order = None + if group == timepoint_col and time_order_map: + pat_tps = pat_df[group].dropna().unique().tolist() + custom_order = sorted(pat_tps, key=lambda x: time_order_map.get(x, 999)) + + create_cdr3_len_plot( + pat_df, + group, + custom_order, + plot_type="line" + ) + else: + print(f"No data for {group} in patient {pat}.") + print("\n") + + print("::::\n") + print(":::::\n") + +else: + # --- STANDARD VIEW (>10 Patients) --- + print("::: {.panel-tabset}\n") + for group in group_opts_all: + print(f"## {group}\n") + create_cdr3_len_plot(df, group, plot_type="box") + print("\n") + print(":::\n") + +``` +**Figure 8. Average Productive CDR3 Length** Visualizes the frequency of CDR3 nucleotide lengths. This plot is used to verify the expected Gaussian distribution typical of polyclonal repertoires and to identify skewness indicative of technical artifacts or clonal dominance + +**Quality Control Interpretations** + +**Validation of Library Complexity (Gaussian Fit):** A high-quality, polyclonal library must exhibit a near-Gaussian (bell-shaped) distribution of CDR3 lengths. + +- Significant deviations (e.g., flat, multimodal, or ragged distributions) in control or baseline samples indicate poor library complexity, RNA degradation, or sampling bottlenecks (insufficient template starting material). + +**Distinguishing Expansion from Bias:** While clonal expansion biologically causes skew, technical artifacts can mimic this signal. + +- A sharp peak at a single length in a sample expected to be naive or healthy suggests PCR bias or amplicon contamination rather than true biological expansion. +- An abundance of unusually short (<15 ntds or <5 AA) or long (>90 ntds or >30 AA) CDR3s often indicates misalignment, primer dimers, or non-specific amplification products that survived upstream filtering. + +## Number of unique TCRs per sample (Richness) {#sec-richness} + + +```{python} +#| output: asis + + +# --- 1. Define General Exclusions --- +# Exclude continuous metrics, identifiers, and derived stats from being used as tabs +exclude_cols = [ + 'sample', 'file', 'total_counts', 'filename', 'sample_id', 'alias', + 'shannon_entropy', 'gini_coefficient', 'inverse_simpson', + 'hill_q0', 'hill_q1', 'hill_q2', 'mean_len', 'std_len', + 'clonality', 'counts', 'num_clones', 'num_TCRs', 'simpson_index', + 'simpson_index_corrected', 'num_prod', 'num_nonprod', 'pct_prod', 'pct_nonprod', + 'productive_cdr3_avg_len', 'num_convergent', 'ratio_convergent' +] + +# Add timepoint_order_col if it exists in the environment +if 'timepoint_order_col' in locals() or 'timepoint_order_col' in globals(): + exclude_cols.append(timepoint_order_col) + +# --- 2. Setup Grouping Options --- +# Automatically grab any categorical columns with reasonable cardinality (< 35 unique values) +group_opts = [ + col for col in df.columns + if col not in exclude_cols and df[col].nunique() < 35 +] + +# --- 3. DEFINE PLOTTING FUNCTION --- +def create_richness_plot(data, group_col, custom_order=None): + """ + Plots the number of TCRs (richness) as a boxplot with individual sample dots. + """ + if data.empty or data[group_col].dropna().empty: + return + + # Determine order: Use custom list if provided, otherwise sort alphabetically + if custom_order: + unique_cats = [x for x in custom_order if x in data[group_col].unique()] + else: + unique_cats = sorted(data[group_col].dropna().unique().tolist()) + + # Safely handle hover data depending on available columns + hover_cols = ['sample'] + if 'alias' in data.columns: + hover_cols.append('alias') + + # --- BASE PLOT --- + fig = px.box( + data, + x=group_col, + y='num_TCRs', + color=group_col, + points='all', + hover_data=hover_cols, + category_orders={group_col: unique_cats}, + template="simple_white", + title=f"TCR Richness by {group_col}" + ) + + # --- DOT STYLING --- + fig.update_traces(width=0.5, marker=dict(size=10, opacity=0.7)) + + # Center the dots and apply jitter so they don't overlap entirely + for trace in fig.data: + if isinstance(trace, go.Box): + trace.pointpos = 0 + trace.jitter = 0.2 + + # --- FINAL LAYOUT --- + fig.update_layout( + xaxis_title=group_col, + yaxis_title="Number of Unique TCRs", + xaxis=dict(tickfont=dict(size=15)), + margin=dict(t=60), + showlegend=False, + width=700, + height=600, + plot_bgcolor='rgba(0,0,0,0)' + ) + + fig.update_yaxes(showgrid=True, gridcolor='lightgrey') + + fig.show() + +# --- 4. EXECUTION LOOP --- +print("::: {.panel-tabset}\n\n") + +for group in group_opts: + if df[group].dropna().empty: + continue + + print(f"## {group}\n\n") + create_richness_plot(df, group) + print("\n") + +print(":::\n\n") + + +``` +**Figure 9. TCR Richness** Visualizes the number of unique CDR3 sequence. Richness is the simplest measure of diversity, counting how many different clones are present in a sample. It is highly sensitive to sequencing depth and undersampling, which can lead to an underestimation of true diversity if not properly accounted for. + +## Sequencing Depth vs Richness {#sec-depth-richness} + +**Does a sample appear "diverse" because it is biologically complex, or simply because it was sequenced more deeply than others?** +In an ideal experiment, **we would sequence until we reach saturation** (the point where sequencing more reads no longer reveals new clonotypes). On this plot (log-log scale), a saturated cohort appears as a horizontal plateau: increasing Total Counts (x-axis) does not significantly increase Richness (y-axis). This indicates that the library has captured the true biological ceiling of the sample, allowing for valid comparisons between patients. + +```{python} +#| output: asis +#| fig-width: 6 +#| fig-height: 5 + +# --- 1. Merge the missing counts data --- +total_counts = concat_df.groupby('sample')['counts'].sum().to_frame(name='total_counts') +bias_plot_df = plot_df.merge(total_counts, on='sample', how='inner') + +# --- 2. Setup Grouping Options --- +exclude_cols = ['sample', 'file', 'total_counts', 'filename', 'sample_id', + 'shannon_entropy', 'gini_coefficient', 'inverse_simpson', + 'hill_q0', 'hill_q1', 'hill_q2', 'mean_len', 'std_len', timepoint_order_col] + +group_opts = [ + col for col in bias_plot_df.columns + if col not in exclude_cols and bias_plot_df[col].nunique() < 35 +] + +# --- 3. CALCULATE UNIFIED AXIS LIMITS --- +# Find the absolute min and max across BOTH axes to ensure a perfect 1:1 square +all_vals = pd.concat([bias_plot_df['total_counts'], bias_plot_df['hill_q0']]) +global_min = all_vals[all_vals > 0].min() +global_max = all_vals.max() + +# Plotly expects log bounds in log10 format +log_min = np.log10(global_min) - 0.15 +log_max = np.log10(global_max) + 0.15 + +# --- 4. Start Quarto Tabset --- +print("::: {.panel-tabset}\n") + +for col in group_opts: + print(f"## {col}\n") + + # Calculate Correlation (Log-Log) + x_val = bias_plot_df['total_counts'] + y_val = bias_plot_df['hill_q0'] # Richness + + # Handle zeros safely for log calculation + x_log = np.log1p(x_val) + y_log = np.log1p(y_val) + + if len(x_log) > 1: + corr, _ = pearsonr(x_log, y_log) + status = "⚠️ BIAS" if corr > 0.8 else "Pass" + title_str = f'Depth Bias by {col}
R = {corr:.2f} ({status})' + else: + corr = 0 + title_str = f'Depth Bias by {col}
(Not enough data)' + + # Force Categorical Coloring + bias_plot_df[col] = bias_plot_df[col].astype(str) + + # Determine hover columns safely + hover_cols = ['sample'] + if 'alias' in bias_plot_df.columns: + hover_cols.append('alias') + + # --- PLOTLY CONVERSION --- + fig = px.scatter( + bias_plot_df, + x='total_counts', + y='hill_q0', + color=col, + symbol=col, + hover_data=hover_cols, + log_x=True, + log_y=True, + template="simple_white", + title=title_str, + labels={ + 'total_counts': 'Sequencing Depth (Total Counts)', + 'hill_q0': 'Richness (Unique Clones)' + } + ) + + # Marker styling + fig.update_traces(marker=dict(size=10, opacity=0.8, line=dict(width=1, color='DarkSlateGrey'))) + + # Add diagonal warning if correlation is high + if corr > 0.8: + fig.add_annotation( + x=0.05, + y=0.95, + xref="paper", + yref="paper", + text="Strong Depth Bias", + showarrow=False, + font=dict(color="red", size=14) + ) + + # --- ENFORCE 1:1 SCALE & CLEAN LOG GRIDS --- + fig.update_layout( + width=700, + height=700, + margin=dict(t=80), + plot_bgcolor='rgba(0,0,0,0)' + ) + + fig.update_xaxes( + showgrid=True, + gridcolor='lightgrey', + range=[log_min, log_max], + dtick=1, # Forces major ticks to exact powers of 10 + minor=dict(showgrid=False) # Hides the clustered minor grid lines + ) + + fig.update_yaxes( + showgrid=True, + gridcolor='lightgrey', + range=[log_min, log_max], + scaleanchor="x", + scaleratio=1, + dtick=1, # Forces major ticks to exact powers of 10 + minor=dict(showgrid=False) # Hides the clustered minor grid lines + ) + + fig.show() + + print("\n") + +print(":::\n") + +``` + +**Figure 10. Evaluation of Sequencing Depth Bias.** This scatter plot correlates sequencing depth (Total Counts, x-axis) with observed richness (Unique Clones, y-axis) on a log-log scale. + +A **strong positive correlation** (Pearson’s $R > 0.8$) **typically signals a technical failure in experimental design**. If your samples form a strict **diagonal line from the bottom-left to the top-right**, your **"diversity" metric is merely a proxy for library size**. In this scenario, **samples with higher read counts will artificially appear more diverse, not because of biology, but because the sampling depth was insufficient to capture the rare tail of the repertoire**. Diversity **comparisons** in such a cohort **are statistically invalid without subsampling (rarefaction)**. + +Samples with **both low sequencing depth and low richness** represent **technical dropouts**. These libraries likely suffered from **low RNA input or poor amplification efficiency**. They often lack sufficient statistical power to characterize the repertoire and **should be excluded**. + +Samples that exhibit **high sequencing depth but paradoxically low richness** are of particular interest. This discordant position—deep sampling yielding few unique clones—indicates that the repertoire is dominated by a **massive clonal expansion** (e.g., leukemia or an acute viral response) or suffers from **technical PCR "jackpotting"**. +Finding a sample with **low read counts but extreme richness** is mathematically suspicious. It implies that almost every read is a unique clonotype, a pattern often caused by **sequencing errors** (phantoms) or **contamination** with high-diversity amplicon debris. + +### Rarefaction Curve + +You **cannot fairly compare the diversity of a sample with 1 million reads to one with 10,000 reads**. The **deeper sample will always look richer simply because you looked harder**. Rarefaction simulates what your deep samples would look like if you had stopped sequencing earlier. + +```{python} + +# --------------------------------------------------------- +# Compute Basic Metrics: Depth & Richness +# --------------------------------------------------------- + +qc_stats = concat_df.groupby('sample').agg( + sequencing_depth=('counts', 'sum'), + clonotype_richness=('CDR3b', 'nunique') +).reset_index() + +# We extract just the metadata columns we need from 'df' +meta_safe = meta[['sample', 'alias', 'origin', timepoint_col]].drop_duplicates() + +# Merge metadata into qc_stats +qc_stats = qc_stats.merge(meta_safe, left_on='sample', right_on='sample', how='left') + +# --------------------------------------------------------- +# Rarefaction Curve Calculation +# --------------------------------------------------------- +# Rarefaction answers: "If we sequenced everyone to the same depth, how many clones would we find?" + + +def get_rarefaction_curve(sample_df, sample, steps=50): + """ + Simulates rarefaction by expanding counts and shuffling. + """ + # Flat array where each element is a clone ID, repeated by its count + # Use integer IDs for memory efficiency instead of strings + counts = sample_df['counts'].values + n_reads = counts.sum() + + if n_reads == 0: + return pd.DataFrame() + + # Generate clone IDs (0 to N-1) + ids = np.arange(len(counts)) + + # e.g., if clone 0 has 2 counts -> [0, 0] + # WARNING: High memory usage for very deep sequencing + flat_repertoire = np.repeat(ids, counts) + + # Shuffle to simulate random sampling + np.random.shuffle(flat_repertoire) + + # Define depths to check (from 0 to total reads) + depths = np.linspace(0, n_reads, num=steps, dtype=int) + depths = np.unique(depths)[1:] # Remove 0 and duplicates + + results = [] + for d in depths: + subsample = flat_repertoire[:d] + n_unique = len(np.unique(subsample)) + results.append({'sample': sample, 'depth': d, 'richness': n_unique}) + + return pd.DataFrame(results) + +rarefaction_data = [] + +unique_samples = concat_df['sample'].unique() +for s in unique_samples: + s_df = concat_df[concat_df['sample'] == s] + r_df = get_rarefaction_curve(s_df, s) + rarefaction_data.append(r_df) + +df_rarefaction = pd.concat(rarefaction_data) + +# Merge metadata into rarefaction data using the safe source +df_rarefaction = df_rarefaction.merge(meta_safe, on='sample', how='left') + +# ----- Rarefaction Curves + +# --------------------------------------------------------- +# Dynamic Plotting with Adaptive Axes +# --------------------------------------------------------- + +# 1. Calculate dynamic limits based on the actual data distribution +# Instead of a hard 90th percentile, we can use the max depth plus a small buffer +absolute_max_depth = df_rarefaction['depth'].max() +absolute_max_richness = df_rarefaction['richness'].max() + +# 2. Define the plot +fig_rare = px.line(df_rarefaction, + x='depth', + y='richness', + color='alias', + title='Rarefaction Curves (Diversity vs Depth)', + labels={'depth': 'Sequencing Depth (Total Counts)', 'richness': 'Identified Clonotypes'}, + hover_data=['origin'], + template='plotly_white') + +# 3. Add Vertical Line for the "Saturation Point" (Minimum Sample Depth) +min_depth = qc_stats['sequencing_depth'].min() +fig_rare.add_vline(x=min_depth, line_dash="dash", line_color="gray", + annotation_text=f"Min Depth: {min_depth}", + annotation_position="top left") + +# 4. Update Axes Dynamically +# Setting range to None (default) allows Plotly to auto-scale, +# but if you want to cap it to the 95th percentile to avoid outliers stretching the plot: +suggested_max_x = qc_stats['sequencing_depth'].quantile(0.95) + +fig_rare.update_layout( + xaxis=dict( + # Use 'range' if you want a specific cutoff, otherwise omit for auto-scaling + range=[0, suggested_max_x * 1.05], # 5% buffer + autorange=False, + rangemode="tozero" + ), + yaxis=dict( + autorange=True, # Automatically adjust Y based on the visible X range + rangemode="tozero" + ) +) + +fig_rare.show() + +``` +**Figure 11. Rarefaction Analysis of TCR Repertoire Diversity.** Rarefaction curves depict the accumulation of unique clonotypes (richness) as a function of sequencing depth (Total Counts) for each sample. The dotted vertical line marks the depth of the smallest library. + +Follow each sample's line from left to right to track how new clonotypes are discovered as sequencing depth increases. In an **ideal scenario**, you want to see a **saturation Plateau**: the curve shoots up initially and then flattens out into a horizontal line. This indicates success; you have sequenced the sample deeply enough to exhaust the available diversity, and spending more money on reads would yield no new information. + +You will often see ***two distinct deviations from this ideal***: +- **Lazy Climber:** A curve that **rises very slowly**. This signals **clonal dominance**; the **sequencer is reading millions of molecules, but they are all the same few sequences**. This is a biological reality, not a technical failure, so these samples should be kept. +- **Unsatisfied Climber:** A curve that **shoots diagonally upward** even at maximum depth, never bending toward horizontal. This indicates **undersampling**; the library is far more diverse than your sequencing budget allowed. Because the "Richness" value here is just a lower bound of an unknown total, these samples often cannot be accurately compared and may need to be discarded. + +The **vertical dashed line** on your plot marks the **Lowest Sequencing Depth in your entire cohort**. To perform a valid statistical comparison, you must mathematically "downsample" every patient to this threshold. This forces a difficult trade-off: if your lowest sample has an unusable depth (e.g., < 1,000 reads), keeping it would require you to throw away 90% of the data from your high-quality samples just to make them comparable. In such cases, it is better to **delete the single bad sample rather than degrade the resolution of the entire cohort**. + +**Re-sequencing?** +This visualization also serves as the financial evidence for deciding whether to re-sequence a dropout sample or not: +- If the curve is **rising steeply**, the sample was **undersampled**. **Re-sequencing is possible** and will yield missing data. +- If the **curve has plateaued**, the **sample is naturally low-diversity** (or the library prep failed). **Re-sequencing would be a waste of money**, as deeper reads will not recover new clonotypes. + +**Expected Richness** +Always verify the magnitude of your Y-axis ("Identified Clonotypes"). For a typical bulk human blood sample, you should expect to see between $10^3$ and $10^5$ unique clonotypes. If your plot shows a maximum of only ~90 clones (as seen in some synthetic or filtered datasets), this is a major red flag for real biological data. Such a low count in primary tissue usually implies a library prep failure, an oligoclonal cell line, or an artifact of synthetic simulation. + + +::: {.callout-caution title="Warning"} +Not all diversity metrics are equally sensitive to library size: + +**High Sensitivity (Depth Matters!) $\rightarrow$ Use Rarefied Data** +Metrics that rely on counting the **number of unique clones** are heavily biased by sequencing depth. +* **Richness ($q=0$):** A direct count of unique clonotypes. This is the most sensitive metric. +* **Shannon Entropy ($q=1$):** While it considers frequency, it is still significantly influenced by the "long tail" of rare clones. +* **Clonality:** Because it is calculated as $1 - \text{Normalized Entropy}$ (which relies on Richness), it inherits this bias. A deep sample will often appear "less clonal" than a shallow one simply due to the math. + +**Low Sensitivity (Depth is Irrelevant) $\rightarrow$ Use Full Data** +Metrics that measure **dominance** or **inequality** focus on the most abundant clones. Because rare clones (singletons) contribute almost zero to these calculations, adding more reads does not significantly change the result. +- **Simpson Index ($q=2$):** Measures the probability of picking two identical clones. It is mathematically weighted toward the top clones, making it extremely robust to depth differences. +- **Gini Coefficient:** Measures the inequality of the distribution (Lorenz curve). It stabilizes quickly and does not require downsampling. +- **Gene Usage (V/J):** Based on percentages (frequency), not raw counts. +::: + +## V-Gene Usage PCA {#sec-pca} + +```{python} +#| output: asis +#| fig-width: 6 +#| fig-height: 5 + +# --- 1. PCA Calculation Function (Unchanged) --- +def analyze_v_usage(df): + """ + Expects df with columns: ['sample', 'TRBV', 'counts'] + """ + matrix = df.groupby(['sample', 'TRBV'])['counts'].sum().unstack(fill_value=0) + freq_matrix = matrix.div(matrix.sum(axis=1), axis=0) + + scaler = StandardScaler() + scaled_data = scaler.fit_transform(freq_matrix) + + pca = PCA(n_components=2) + components = pca.fit_transform(scaled_data) + + pca_df = pd.DataFrame(data=components, columns=['PC1', 'PC2'], index=freq_matrix.index) + return freq_matrix, pca_df + +# --- 2. Run Analysis & Prepare Data --- +_, pca_df = analyze_v_usage(concat_df) + +if pca_df.index.name == 'sample' or 'sample' not in pca_df.columns: + pca_df = pca_df.reset_index() + if 'index' in pca_df.columns and 'sample' not in pca_df.columns: + pca_df = pca_df.rename(columns={'index': 'sample'}) + +plot_data = pd.merge(pca_df, meta, on='sample', how='inner') + +# --- 3. Determine Grouping Columns --- +exclude_cols = ['sample', 'file', 'total_counts', 'filename', 'sample_id'] +group_opts = [ + col for col in meta.columns + if col not in exclude_cols and meta[col].nunique() < 35 +] + +# --- 4. Generate Tabs & Plotly Charts --- +print("::: {.panel-tabset}\n") + +for col in group_opts: + print(f"## {col}\n") + + # Force categorical coloring + plot_data[col] = plot_data[col].astype(str) + + # Safely determine hover columns + hover_cols = ['sample'] + if 'alias' in plot_data.columns: + hover_cols.append('alias') + + # --- PLOTLY SCATTER PLOT --- + fig = px.scatter( + plot_data, + x='PC1', + y='PC2', + color=col, + symbol=col, + hover_data=hover_cols, # <-- Adds the hover over info + template='simple_white', + title=f'TRBV Gene Usage PCA (Colored by {col})' + ) + + # Match the Seaborn styling (large dots, transparency, black outlines) + fig.update_traces(marker=dict(size=12, opacity=0.8, line=dict(width=1, color='DarkSlateGrey'))) + + # --- OUTLIER LABELING LOGIC --- + outlier_mask = (np.abs(plot_data['PC1']) > 3 * plot_data['PC1'].std()) | \ + (np.abs(plot_data['PC2']) > 3 * plot_data['PC2'].std()) + + outliers = plot_data[outlier_mask] + + # Add permanent text labels next to outlier points + for _, row in outliers.iterrows(): + fig.add_annotation( + x=row['PC1'], + y=row['PC2'], + text=f"{row['sample']}", + showarrow=False, + xshift=15, # Shifts text to the right + yshift=15, # Shifts text up + font=dict(size=10, color="black") + ) + + # --- FINAL LAYOUT --- + fig.update_layout( + xaxis_title="PC1", + yaxis_title="PC2", + width=700, + height=600, + margin=dict(t=60), + plot_bgcolor='rgba(0,0,0,0)' + ) + + # Add gridlines and emphasize the 0,0 origin axes + fig.update_xaxes(showgrid=True, gridcolor='lightgrey', zeroline=True, zerolinecolor='gray', zerolinewidth=1.5) + fig.update_yaxes(showgrid=True, gridcolor='lightgrey', zeroline=True, zerolinecolor='gray', zerolinewidth=1.5) + + fig.show() + + print("\n") + +print(":::\n") + +``` +**Figure 12. Dimensionality Reduction of TRBV Gene Usage.** Principal Component Analysis (PCA) plot visualizing the variation in T-cell Receptor Beta Variable (TRBV) gene usage frequencies across samples, colored by metadata variables. + +This plot compresses the high-dimensional complexity of T-cell receptor usage (e.g., TRBV1, TRBV2... TRBV30) into a 2D map. Each dot represents a sample, and the **distance between dots represents the similarity of their "V-gene fingerprints."** **If two samples are close together, they are using the V-genes in almost identical proportions**. If they are far apart, their repertoire structures are radically different. + +**Technical Failures 🚩** +One **purpose** of this visualization is to **detect technical artifacts**. When you toggle the tabs to color by **Sequencing Batch for same condition**, you should ideally see a randomly **mixed cloud**. This indicates that the V-gene usage is **consistent** across your **experiments**. However, **if you see distinct, non-overlapping islands** corresponding to different batches (e.g., all "Batch A" samples cluster on the left, and "Batch B" samples on the right), **you have a severe Batch Effect**. This usually implies a change in the wet-lab protocol, such as a new set of multiplex primers that amplify genes with different efficiencies. If observed, **you cannot compare samples across these batches without statistical correction** (e.g., ComBat). + +**Biological Interpretation 🦠** +In the **context of cancer treatment** (especially Checkpoint Inhibitors like PD-1 blockade), **we are looking for Remodeling**. We want to **know if the drug "woke up" the immune system and caused it to expand new T-cell armies (clones) to fight the tumor. For example, if "Responders" form a distinct island separate from "Non-Responders", this suggests that having a specific bias toward a certain V-gene, TRBV-19 for example, pre-disposes a patient to fight the cancer effectively. + +An outlier on this plot can have a **specific biological cause**: + +- ***Monoclonal Expansion***. If a patient has a **massive expansion of a single clone** (e.g., a clone using TRBV7-2 takes up 50% of the repertoire), that sample will be pulled violently toward the "TRBV7-2 direction" in PCA space. Finding such an **outlier in a "healthy" control often indicates a PCR jackpotting error or contamination**. + +## Outlier Detection {#sec-outlier-flag} + +This aggregates everything. We calculate Z-scores for every metric. If a sample is $>3$ standard deviations away from the mean, we flag it. +```{python} + +# CDR3 Length Distribution & Distances +def analyze_cdr3_lengths(df): + """ + Expects df with columns: ['sample', 'CDR3b'] + """ + # Calculate lengths + df['cdr3_len'] = df['CDR3b'].apply(len) + + # 1. Basic Stats per sample + stats = df.groupby('sample')['cdr3_len'].agg( + mean_len='mean', + std_len='std', + skew_len=lambda x: skew(x) + ) + + # 2. Cohort "Ideal" Distribution + # We aggregate ALL samples to create a reference distribution + cohort_distribution = df['cdr3_len'].values + + # 3. Distance Calculation + # How far is each sample from the cohort average? + distances = [] + for sample in stats.index: + sample_dist = df[df['sample'] == sample]['cdr3_len'].values + + # Earth Mover's Distance (Wasserstein) + # Robust metric for histogram similarity + wd = wasserstein_distance(sample_dist, cohort_distribution) + distances.append(wd) + + stats['cdr3_len_dist_to_cohort'] = distances + return stats + +def run_qc_pipeline(df): + div_df = calculate_diversity_metrics(df) + len_df = analyze_cdr3_lengths(df) + _, pca_df = analyze_v_usage(df) + + # Merge all metrics into one master QC table + qc_table = div_df.join(len_df).join(pca_df) + + # Outlier Detection Logic (Z-Scores) + # We define which columns we care about for outliers + metrics_to_check = ['shannon_entropy', 'inverse_simpson', 'mean_len', 'cdr3_len_dist_to_cohort', 'PC1', 'PC2'] + + flags = [] + + for metric in metrics_to_check: + mu = qc_table[metric].mean() + sigma = qc_table[metric].std() + + # Calculate Z-score + z_score = (qc_table[metric] - mu) / sigma + + # Find outliers (|z| > 3 is standard, maybe use 2.5 for stricter QC) + outliers = qc_table[np.abs(z_score) > 3].index.tolist() + + for out in outliers: + flags.append({ + 'sample_id': out, + 'metric': metric, + 'value': qc_table.loc[out, metric], + 'z_score': z_score[out], + 'flag': f"Extreme {metric}" + }) + + return qc_table, pd.DataFrame(flags) + +qc_table,flags = run_qc_pipeline(concat_df) +print("\nMetrics values for all samples\n") +print(qc_table) +print("\nThe following samples had a metric that deviates from the rest (|z-score|>3):\n") +print(flags) +``` From d00be437aeb97f65466c583390c94c6734a9083c Mon Sep 17 00:00:00 2001 From: dimalvovs Date: Mon, 22 Jun 2026 14:06:06 -0400 Subject: [PATCH 05/21] revert back to enable tests --- notebooks/compare_stats_template.qmd | 190 +++ notebooks/gliph2_report_template.qmd | 56 + notebooks/sample_stats_template.qmd | 316 ++++ notebooks/template_qc.qmd | 2328 -------------------------- 4 files changed, 562 insertions(+), 2328 deletions(-) create mode 100644 notebooks/compare_stats_template.qmd create mode 100644 notebooks/gliph2_report_template.qmd create mode 100644 notebooks/sample_stats_template.qmd delete mode 100644 notebooks/template_qc.qmd diff --git a/notebooks/compare_stats_template.qmd b/notebooks/compare_stats_template.qmd new file mode 100644 index 0000000..01b2c66 --- /dev/null +++ b/notebooks/compare_stats_template.qmd @@ -0,0 +1,190 @@ +--- +title: "Comparative T Cell Repertoire statistics" +format: + html: + theme: flatly + toc: true + toc_depth: 3 + code-fold: show + embed-resources: true + number-sections: true + smooth-scroll: true + grid: + body-width: 1000px + margin-width: 300px + +jupyter: python3 +--- + +Thank you for using TCRtoolkit! This report is generated from sample data and metadata you provided. The report is divided into two sections: + +@sec-heatmap : Heatmap of sample to sample repertoire similarity using Jaccard, Sorensen, and Morisita indices (v1.0) + +# Report Setup + +```{python, echo=false} +#| tags: [parameters] +#| echo: false + +## 2. Pipeline Parameters +#Default inputs are overwritten at the command line in `modules/local/plot_sample.nf` +workflow_cmd='' +project_name='path/to/project_name' +jaccard_mat='path/to/jaccard_mat.csv' +sorensen_mat='path/to/sorensen_mat.csv' +morisita_mat='path/to/morisita_mat.csv' +# jsd_mat='path/to/jsd_mat.csv' +sample_utf8='path/to/sample_utf8.csv' +``` + +```{python} +#| tags: [setup] +#| warning: false + +# 1. Load Packages +import os +import shutil +import datetime +import sys +import numpy as np +import pandas as pd +import seaborn as sns +# from matplotlib.colors import LinearSegmentedColormap +# import scipy.cluster.hierarchy as sch + +# 2. Print Pipeline Information +print('Pipeline information and parameters:' + '\n') +print('Project Name: ' + project_name) +print('Workflow command: ' + workflow_cmd) +print('Date and time: ' + str(datetime.datetime.now())) + +# 3. Importing similarity data +## 3a. jaccard similarity matrix +jaccard_df = pd.read_csv(jaccard_mat, sep=',', header=0, index_col=0) + +## 3b. sorensen similarity matrix +sorensen_df = pd.read_csv(sorensen_mat, sep=',', header=0, index_col=0) + +## 3c. morisita similarity matrix +morisita_df = pd.read_csv(morisita_mat, sep=',', header=0, index_col=0) + +``` + +# Analysis + +## Overall Repertoire Similarity {#sec-heatmap} + +Similarity metrics such as Jaccard, Sorensen, and Morisita are often used to compare the similarity between two samples. Here, we compare the similarity of TCR repertoires between samples using these three metrics. Details on how each metric is calculated can be found below the figure. + +```{python} +import plotly.express as px +import plotly.graph_objects as go +from plotly.figure_factory import create_dendrogram +from plotly.subplots import make_subplots +import scipy.spatial.distance as ssd +import matplotlib.pyplot as plt + +# preprocess the data prior to clustering +jaccard_numeric = jaccard_df.apply(pd.to_numeric, errors='coerce') +sorensen_numeric = sorensen_df.apply(pd.to_numeric, errors='coerce') +morisita_numeric = morisita_df.apply(pd.to_numeric, errors='coerce') + +# Assuming jaccard_numeric, sorensen_numeric, and morisita_numeric are your DataFrames +# and sns_cluster_jaccard, sns_cluster_sorensen, and sns_cluster_morisita are the corresponding clustermaps +sns_cluster_jaccard = sns.clustermap(jaccard_numeric) +plt.close() +sns_cluster_sorensen = sns.clustermap(sorensen_numeric) +plt.close() +sns_cluster_morisita = sns.clustermap(morisita_numeric) +plt.close() +# sns_cluster_jsd = sns.clustermap(jsd_numeric) +# plt.close() + +# Create a subplot with 3 rows +fig = make_subplots(rows=3, cols=1) + +# Reindex the dataframes to match the clustering +jaccard_clustered = jaccard_numeric.iloc[sns_cluster_jaccard.dendrogram_row.reordered_ind, sns_cluster_jaccard.dendrogram_col.reordered_ind] +sorensen_clustered = sorensen_numeric.iloc[sns_cluster_sorensen.dendrogram_row.reordered_ind, sns_cluster_sorensen.dendrogram_col.reordered_ind] +morisita_clustered = morisita_numeric.iloc[sns_cluster_morisita.dendrogram_row.reordered_ind, sns_cluster_morisita.dendrogram_col.reordered_ind] +# jsd_clustered = jsd_numeric.iloc[sns_cluster_jsd.dendrogram_row.reordered_ind, sns_cluster_jsd.dendrogram_col.reordered_ind] + +# Create individual heatmaps +heatmap_jaccard = go.Heatmap( + z=jaccard_clustered, + x=jaccard_clustered.columns, + y=jaccard_clustered.index, + coloraxis="coloraxis", + visible=False) +heatmap_sorensen = go.Heatmap( + z=sorensen_clustered, + x=sorensen_clustered.columns, + y=sorensen_clustered.index, + coloraxis="coloraxis", + visible=False) +heatmap_morisita = go.Heatmap( + z=morisita_clustered, + x=morisita_clustered.columns, + y=morisita_clustered.index, + coloraxis="coloraxis", + visible=True) + +# Add the heatmaps to the figure +fig = go.Figure(data=[heatmap_jaccard, heatmap_sorensen, heatmap_morisita]) + +# Create buttons to switch between the heatmaps +buttons = [ + dict(label="Jaccard", method="update", + args=[{"visible": [True, False, False]}, {"title": "Jaccard"}]), + dict(label="Sorensen", method="update", + args=[{"visible": [False, True, False]}, {"title": "Sorensen"}]), + dict(label="Morisita", method="update", + args=[{"visible": [False, False, True]}, {"title": "Morisita"}]) +] + +# Update the layout of the figure +fig.update_layout( + updatemenus=[dict(type="buttons", showactive=True, buttons=buttons)], + title='Similarity Matrices', + xaxis_title='Sample ID', + yaxis_title='Sample ID', + autosize=False, + width=950, + height=950, + coloraxis=dict(colorscale='Viridis', colorbar=dict(title="Log Scale")) +) + +fig.show() +``` + +**Jaccard Index**: + +$$ +J(A,B) = \frac{|A \cap B|}{|A \cup B|} +$$ + +Where $A$ and $B$ are the sets of unique CDR3 amino acid sequences (TCRs) in samples A and B, respectively. The Jaccard Index is defined as the ratio of the number of common elements between two sets to the total number of distinct elements in the two sets. + +**Sorensen Index**: + +$$ +S(A,B) = \frac{2|A \cap B|}{|A| + |B|} +$$ + +Where $A$ and $B$ are the sets of unique CDR3 amino acid sequences (TCRs) in samples A and B, respectively. The difference between the sorensen index and the jaccard index is that the sorensen index takes into account the size of the two sets being compared, while the jaccard index only considers the number of common elements between the two sets. + +**Morisita-Horn Index:** + +$$ +M(A,B)=\frac{2\sum_{i=1}^{S}a_{i}b_{i}}{(D_{a}+D_{b})AB} ; D_{a}=\frac{\sum_{i=1}^{S}a_{i}^{2}}{A^2}, D_{b}=\frac{\sum_{i=1}^{S}b_{i}^{2}}{B^2} +$$ + +Where: + +- $A$ and $B$ are the sets of unique CDR3 amino acid sequences (TCRs) in samples A and B, + +- $a_{i}$ ($b_{i}$) is the number of times TCR $i$ is represented in the total $A$ ($B$) from one sample. + +- $S$ is the total number of unique TCRs in the two samples. + +- $D_{a}$ and $D_{b}$ are the Simpson Index values for samples A and B, respectively. diff --git a/notebooks/gliph2_report_template.qmd b/notebooks/gliph2_report_template.qmd new file mode 100644 index 0000000..e39d248 --- /dev/null +++ b/notebooks/gliph2_report_template.qmd @@ -0,0 +1,56 @@ +--- +title: "Comparative T Cell Repertoire statistics" +format: + html: + theme: flatly + toc: true + toc_depth: 3 + code-fold: show + embed-resources: true + number-sections: true + smooth-scroll: true + grid: + body-width: 1000px + margin-width: 300px + +jupyter: python3 +--- + +Thank you for using TCRtoolkit! This report is generated from sample data and metadata you provided. The report is divided into two sections: + +@sec-heatmap : Heatmap of sample to sample repertoire similarity using Jaccard, Sorensen, and Morisita indices (v1.0) + +# Report Setup + +```{python, echo=false} +#| tags: [parameters] +#| echo: false + +## 2. Pipeline Parameters +#Default inputs are overwritten at the command line in `modules/local/plot_gliph2.nf` +workflow_cmd='' +project_name='path/to/project_name' +clusters='path/to/{project_name}_cluster.csv' +cluster_stats='path/to/{project_name}_cluster.txt' +``` + +```{python} +#| tags: [setup] +#| warning: false + +# 1. Load Packages +import os +import shutil +import datetime +import sys +import numpy as np +import pandas as pd +import seaborn as sns +# from matplotlib.colors import LinearSegmentedColormap +# import scipy.cluster.hierarchy as sch + +# 2. Print Pipeline Information +print('Pipeline information and parameters:' + '\n') +print('Project Name: ' + project_name) +print('Workflow command: ' + workflow_cmd) +print('Date and time: ' + str(datetime.datetime.now())) diff --git a/notebooks/sample_stats_template.qmd b/notebooks/sample_stats_template.qmd new file mode 100644 index 0000000..a85d7ef --- /dev/null +++ b/notebooks/sample_stats_template.qmd @@ -0,0 +1,316 @@ +--- +title: "Sample Level T Cell Repertoire statistics" +format: + html: + theme: flatly + toc: true + toc_depth: 3 + code-fold: show + embed-resources: true + number-sections: true + smooth-scroll: true + grid: + body-width: 1000px + margin-width: 300px + +jupyter: python3 +--- + +Thank you for using TCRtoolkit! This report is generated from sample data and metadata you provided. The report is divided into three sections: + +@sec-report-setup : Code to setup the report. This section includes the parameters you used to run the pipeline, loading necessary packages, data, etc. + +@sec-sample-level-stats : Typical sample level T cell repertoire statistics. Each plot has a description about the statistic shown and formulas used to calculate the metric. + +@sec-gene-family-usage : TCR V gene family usage. The x-axis shows the timepoint collected for each individual, and the y-axis shows the proportion of TCRs that use each V gene family. + +# Report Setup {#sec-report-setup} + +```{python, echo=false} +#| tags: [parameters] +#| echo: false + +workflow_cmd='' +project_name='' +sample_table='' +sample_stats_csv='' +v_family_csv='' + +samplechart_x_col='' +samplechart_color_col='' +vgene_subject_col='' +vgene_x_cols='' +``` + +```{python} +#| code-fold: true + +# 1. Load Packages +from IPython.display import Image +import os +import datetime +import sys +import pandas as pd +import math +import matplotlib.pyplot as plt +import seaborn as sns +from matplotlib.colors import LinearSegmentedColormap +import plotly.express as px +import plotly.graph_objects as go + +import warnings +warnings.filterwarnings( + "ignore", + category=FutureWarning, + module="plotly" +) + +# 2. Print pipeline parameters + +print('Project Name: ' + project_name) +print('Workflow command: ' + workflow_cmd) +print('Date and time: ' + str(datetime.datetime.now())) + +# 3. Loading data +## reading sample metadata +meta = pd.read_csv(sample_table, sep=',') +meta_cols = meta.columns.tolist() + +## reading combined repertoire statistics +df = pd.read_csv(sample_stats_csv, sep=',') +df = pd.merge(df, meta, on='sample', how='left') +df = df[meta_cols + [c for c in df.columns if c not in meta_cols]] + +## reading V gene family usage +v_family = pd.read_csv(v_family_csv, sep=',') +v_family = pd.merge(v_family, meta, on='sample', how='left') +v_family = v_family[meta_cols + [c for c in v_family.columns if c not in meta_cols]] +v_family = v_family.sort_values(by=[vgene_subject_col]) +``` + +# Sample level statistics {#sec-sample-level-stats} + +Below are plots showing basic T cell repertoire statistics. Each plot has a description about the statistic shown and formulas used to calculate the metric when applicable. Specific biological interpretation of each plot is left to the user. + +Version 3 of these plots features plotly express interactive plots. This version is exploratory and may be updated in the future. + +```{python} +#| code-fold: true + + +x_category = df[samplechart_x_col].unique().tolist() +x_category.sort() + +def samplechart(df, y_col): + fig = px.box( + df, + x=samplechart_x_col, + y=y_col, + color=samplechart_color_col, + points='all', + hover_data='sample', + category_orders={f'{samplechart_x_col}': x_category} + ) + + fig.update_layout( + margin=dict(b=80) + ) + + fig.show() +``` + +## Number of clones + +```{python} +#| code-fold: true + +samplechart(df, y_col='num_clones') +``` + +**Figure 1. Number of clones across sample groupings.** A clone is defined as a T cell with a unique CDR3 amino acid sequence. The number of clones is shown on the y-axis. The x-axis represents sample groupings defined by user-specified metadata fields (eg. origin, timepoint). + +## Clonality + +```{python} +#| code-fold: true + +samplechart(df, y_col='clonality') +``` + +**Figure 2. The clonality of samples across sample groupings.** Clonality is a measure of T cell clonal expansion and reflects the degree to which the sample is dominated by 1 or more T cell clones. Clonality is calculated via: $$Clonality = \frac {1-H} {\log_{2} N} \quad\text{,}\quad H = -\sum\limits_{i=1}^N p_i \log_{2}{p_i}$$ where $H$ is the Shannon entropy of a given sample, $N$ is the number of unique TCRs in the sample, and $p_i$ is the frequency of the $i$ th unique TCR in the sample. + +## Simpson Index + +```{python} +#| code-fold: true + +samplechart(df, y_col='simpson_index_corrected') +``` + +**Figure 3. Corrected Simpson Index.** The Simpson Index is a measure of diversity that takes into account the number of clones and the relative abundance of each clone in a sample. The corrected Simpson Index, $D$, is calculated as: + +$$D = \sum\limits_{i=1}^N \frac {p_i(p_i - 1)} {T(T - 1)} \quad\text{,}\quad T = \sum\limits_{i=1}^N p_i$$ + +Where $N$ is the number of unique TCRs in the sample, $p_i$ is the frequency of the $i$ th unique TCR in the sample, and $T$ is the total number of T Cells counted in the sample. + +## Percent of productive TCRs + +```{python} +#| code-fold: true + +samplechart(df, y_col='pct_prod') +``` + +**Figure 4. Percent of productive TCRs.** A productive TCR is a DNA/RNA sequence that can be translated into a protein sequence, i.e. it does not contain a premature stop codon or an out of frame rearrangement. The percent of productive TCRs is calculated as: + +$$ Percent \text{ } productive \text{ } TCRs = \frac P N $$ + +where $P$ is the number of productive TCRs and $N$ is the total number of TCRs in a given sample. + +## Average productive CDR3 Length + +```{python} +#| code-fold: true + +samplechart(df, y_col='productive_cdr3_avg_len') +``` + +**Figure 5. Average Productive CDR3 Length** The average length of the CDR3 region of the TCR for productive clones. The CDR3 region is the most variable region of the TCR and is the region that determines antigen specificity. + +## TCR Convergence + +```{python} +#| code-fold: true + +samplechart(df, y_col='ratio_convergent') +``` + +**Figure 6. TCR Convergence** The ratio of convergent TCRs to total TCRs. A convergent TCR is a TCR that is generated via 2 or more unique nucleotide sequences via codon degeneracy. + +# Gene Family Usage {#sec-gene-family-usage} + +## V gene family usage + +The V gene family usage of the TCRs in each sample is shown in the plots below. The x-axis shows the timepoint collected for each individual, and the y-axis shows the proportion of TCRs that use each V gene family. + +The V gene usage proportion, $V_k$, is calculated via: + +$$ +V_k = \frac {N_{k}} {T} \quad\text{,}\quad T = \sum\limits_{i=1}^N p_i +$$ + +where $N_{k}$ is the number of TCRs that use the $k$ th V gene, and T is the total number of TCRs in the sample. + +```{python} +#| code-fold: true + +## code adapted from https://www.moritzkoerber.com/posts/plotly-grouped-stacked-bar-chart/ +colors = ["#fafa70","#fdef6b","#ffe566","#ffda63","#ffd061","#ffc660","#ffbb5f","#fdb15f","#fba860","#f79e61","#f39562","#ef8c63","#e98365","#e37b66","#dd7367","#d66b68","#ce6469","#c65e6a","#bd576b","#b4526b","#ab4c6b","#a1476a","#974369","#8c3e68","#823a66","#773764","#6d3361","#62305e","#572c5a","#4d2956"] +vgene_x_cols = vgene_x_cols.split(',') + +## calculate proportions and add to v_family_long +v_family_long = pd.melt(v_family, + id_vars=meta_cols, + value_vars=[c for c in v_family.columns if c.startswith('TRBV')], + var_name='v_gene', + value_name='count') +v_family_long['proportion'] = v_family_long.groupby(meta_cols)['count'].transform(lambda x: x / x.sum()) + +## add in the total number of v genes for each sample +total_v_genes = v_family_long.groupby(meta_cols)['count'].sum().reset_index().rename(columns={'count': 'total_v_genes'}) +v_family_long = v_family_long.merge(total_v_genes, on=meta_cols) + +subjects = ( + v_family_long[vgene_subject_col] + .dropna() + .unique() +) + +for subject in subjects: + current = v_family_long[v_family_long[vgene_subject_col] == subject] + fig = go.Figure() + fig.update_layout( + template="simple_white", + title_text=f"{vgene_subject_col}: {subject}", + xaxis=dict(title_text=f"{', '.join(vgene_x_cols)}"), + yaxis=dict(title_text="proportion"), + barmode="stack", + margin=dict(b=80, r=115), + legend=dict( + orientation="v", + x=1.02, + y=1, + yanchor="top", + itemsizing="constant", + traceorder="reversed" + ) + ) + + fig.update_xaxes( + automargin=True, + title_standoff=30 + ) + + if len(vgene_x_cols) == 1: + x_vals = current[vgene_x_cols[0]] + else: + x_vals = current[vgene_x_cols].astype(str).agg('_'.join, axis=1) + summary = ( + current + .assign(x=x_vals) + .drop_duplicates(subset=meta_cols) + .groupby('x', as_index=False) + .agg( + total_v_genes=('total_v_genes', 'sum'), + ) + ) + y_df = ( + current + .assign(x=x_vals) + .groupby('x', as_index=False) + .agg( + y_top=('proportion', 'sum') + ) + ) + summary = y_df.merge(summary, on='x') + + for g, c in zip(current.v_gene.unique(), colors): + plot_df = current[current.v_gene == g] + if len(vgene_x_cols) == 1: + x = plot_df[vgene_x_cols[0]] + else: + x = [plot_df[col] for col in vgene_x_cols] + fig.add_trace( + go.Bar( + x=x, + y=plot_df['proportion'], + name=g, + marker_color=c, + legendgroup=g, + + customdata=plot_df[['sample']], + hovertemplate=( + "V gene: %{fullData.name}
" + "Sample: %{customdata[0]}
" + "Proportion: %{y}
" + "" + ), + + text=plot_df['total_v_genes'] if ((g == 'TRBV30') and (len(vgene_x_cols) != 1)) else None, + textposition='outside', + ) + ) + + if len(vgene_x_cols) == 1: + for _, row in summary.iterrows(): + fig.add_annotation( + x=row['x'], + y=row['y_top'] + 0.02, + text=str(row['total_v_genes']), + showarrow=False, + yanchor='bottom', + font=dict(size=12) + ) + + fig.show() +``` diff --git a/notebooks/template_qc.qmd b/notebooks/template_qc.qmd deleted file mode 100644 index 7928457..0000000 --- a/notebooks/template_qc.qmd +++ /dev/null @@ -1,2328 +0,0 @@ ---- -title: "Quality Control" -format: - html: - theme: flatly - toc: true - toc_depth: 3 - code-fold: true - embed-resources: true - number-sections: true - smooth-scroll: true - grid: - body-width: 1000px - margin-width: 300px - -jupyter: python3 ---- - - - -Thank you for using TCRtoolkit! This report is generated from the data provided. - -:::{.callout-note collapse="true"} -## Document Information -**Current Version:** 1.0-beta -**Last Updated:** March 2026 -**Maintainer:** BTC Data Science Team -**Notes:** -::: - -```{python} -#| tags: [parameters] -#| include: false - -# --------------------------------------------------------- -# BASE PARAMETERS -# --------------------------------------------------------- -workflow_cmd = '' -project_name='' -project_dir='' -sample_table='' - -timepoint_col = 'timepoint' -timepoint_order_col = 'timepoint_order' -alias_col = 'alias' -subject_col = 'subject_id' - -``` - -```{python} -#| include: false - -# --------------------------------------------------------- -# DERIVED PATHS -# --------------------------------------------------------- - -# Define files -project_dir=f"{project_dir}/{project_name}" -sample_stats_csv = f"{project_dir}/sample/sample_stats.csv" -concat_file = f"{project_dir}/annotate/concatenated_cdr3_sorted.tsv" - - -``` - -```{python} -#| code-fold: true - -from IPython.display import Image -import os -import datetime -import sys -import pandas as pd -import math -import matplotlib.pyplot as plt -import seaborn as sns -from matplotlib.colors import LinearSegmentedColormap -import plotly.express as px -import plotly.graph_objects as go -import glob -import itertools -import h5py -import igraph as ig -import matplotlib.ticker as ticker -import numpy as np -import scipy.cluster.hierarchy as sch -from IPython.display import display, Markdown -from scipy.sparse import csr_matrix -from scipy.stats import gaussian_kde -from scipy.stats import entropy -from scipy.stats import skew, wasserstein_distance -from scipy.stats import pearsonr -from sklearn.decomposition import PCA -from sklearn.preprocessing import StandardScaler -from scipy import sparse -from sklearn.preprocessing import normalize -from scipy.stats import wilcoxon -from scipy.stats import mannwhitneyu -from itertools import combinations -from scipy.spatial.distance import pdist -from scipy.cluster.hierarchy import linkage, leaves_list -import plotly.figure_factory as ff - -import warnings - -meta = pd.read_csv(sample_table, sep=',') -meta_cols = meta.columns.tolist() - -df = pd.read_csv(sample_stats_csv, sep=',') -df = pd.merge(df, meta, on='sample', how='left') -df = df[meta_cols + [c for c in df.columns if c not in meta_cols]] - -concat_df = pd.read_csv(concat_file, sep='\t') - -``` - -# Technical QC & Biological Relevance 🚩/🦠 {#sec-qc-bio} -## Diversity metrics {#sec-div-metrics} - -Analyzing the architecture of a T-cell repertoire requires moving beyond simply counting the total number of unique sequences (Richness). A healthy immune system is highly diverse and relatively evenly distributed. However, upon antigen encounter, specific T-cells rapidly proliferate, causing the repertoire to become highly skewed. -To quantify this structural shift, we rely on ecological diversity metrics and inequality indices. - -- **Shannon Entropy ($H$)** -Answers the question: "How hard is it to predict the identity of a randomly drawn T-cell?" It quantifies the uncertainty of the clonotype distribution. - - ***High Entropy (High Uncertainty):*** The repertoire is very diverse. There are many different clones with comparable frequencies, so you have no idea which one you will pick next. This signifies a healthy, polyclonal repertoire. - - ***Low Entropy (Low Uncertainty):*** The repertoire is dominated by a few large clones. You can easily guess that the next T-cell will likely belong to the dominant clone. This signifies oligoclonality or clonality. - -It is calculated as: -$H = - \sum_{i=1}^{R} p_i \ln(p_i)$ -Where $R$ is the number of unique clonotypes (Richness), and $p_i$ is the frequency of the $i$-th clonotype. - -- **Inverse Simpson Index ($1/D$)** -Answers the question: "How many TCR clones actually matter in this sample?" It ignores the "long tail" of rare, single-read clones (which might just be sequencing errors or debris) and focuses on the dominant, expanded clones that actually drive the immune response. -While the standard Simpson Index ($D$) measures the probability of collision (picking the same clone twice), the Inverse Simpson converts that probability into an "effective number". - - **High Value:** Indicates a high effective number of TCR clones. The sample is diverse and even in clone size. - - **Low Value (approaching 1):** Indicates the sample behaves as if it only contains 1 clone. It is highly dominated by a single expansion. -It is calculated as: -$\frac{1}{D} = \frac{1}{\sum_{i=1}^{R} p_i^2}$ -Where $p_i$ is the frequency of the $i$-th clonotype. - -::: {.callout-tip title="Pro Tip"} -Compare with sample Richness (# of unique clones) - -If Richness = 5,000 but Inverse Simpson = 1.2: Your sample behaves as if it only has 1.2 clones. It is massively dominated by a single expansion. - -If Richness = 5,000 and Inverse Simpson = 4,500: Your sample behaves as if it has 4,500 equally sized clones. It is extremely even and polyclonal. -::: - - -- **Gini Coefficient** -Answers the question: "Has the immune system picked a winner yet?" Measures inequality. Gini is your best metric for detecting Clonal Expansions independent of how many clones you actually sequenced. - - ***Value close to 0 (Perfect Equality):*** Every clone has the exact same frequency. The immune system is "resting" or "naive." It hasn't been triggered by a specific threat, so no single clone has started to divide rapidly. - - ***Medium value (0.3 - 0.7):*** Reactive. The immune system is fighting something. A few clones have expanded to fight a virus or tumor, but the background diversity is still there. - - ***High value (> 0.8):*** Monoclonal / Oligoclonal. A few clones have taken over completely. If this is a tumor sample, it might indicate a TIL (Tumor Infiltrating Lymphocyte) expansion. -It is calculated as: -$G = \frac{\sum_{i=1}^{R} (2i - R - 1) p_i}{R \sum_{i=1}^{R} p_i}$ - -Where the frequencies $p_i$ are sorted in ascending order, and $i$ is the rank. - -::: {.callout-tip title="Pro Tip"} -Why use this instead of Shannon? -Shannon Entropy is hard to compare if your library sizes (depth) are wildly different. Gini is normalized (always 0 to 1). -::: - -- **Hill Numbers ($^qD$)** -Answer the question: "What is the effective number of species in the sample when we ignore (penalize to degree $q$) rare clones?" This is a unified metric that combines richness, entropy, and dominance into a single scale (counts of effective species). - - ***$q=0$ (Richness):*** Counts every unique clone equally, regardless of frequency. Sensitive to sequencing depth and errors. - - ***$q=1$ (Exponential Shannon):*** Weighs clones by their frequency. Represents the number of "common" clones. - - ***$q=2$ (Inverse Simpson):*** Heavily weighs dominant clones. Represents the number of "very dominant" clones. -The general formula for Hill numbers of order $q$ is: -$^qD = \left( \sum_{i=1}^{R} p_i^q \right)^{1/(1-q)}$ -For $q=0$, this simplifies to $R$ (Total unique clonotypes). For $q=1$, the limit is undefined, so we use $\exp(H)$. For $q=2$, this simplifies to $1 / \sum p_i^2$ (Inverse Simpson). - - -```{python} -#| label: diversity-nested -#| output: asis - -# --- 1. CALCULATION FUNCTION --- -def calculate_diversity_metrics(df): - results = [] - # Check if 'counts' exists, if not try 'clone_count' or similar - count_col = 'counts' if 'counts' in df.columns else df.columns[0] # Fallback - - grouped = df.groupby('sample') - for sample, data in grouped: - counts = data[count_col].values - if counts.sum() == 0: continue - - p = counts / counts.sum() - shannon = -np.sum(p * np.log(p)) - inv_simpson = 1.0 / np.sum(p**2) - sorted_p = np.sort(p) - n = len(p) - index = np.arange(1, n + 1) - gini = ((2 * index - n - 1) * sorted_p).sum() / (n * sorted_p.sum()) - q0 = len(p) - q1 = np.exp(shannon) - q2 = inv_simpson - results.append({ - 'sample': sample, - 'shannon_entropy': shannon, - 'inverse_simpson': inv_simpson, - 'gini_coefficient': gini, - 'hill_q0': q0, - 'hill_q1': q1, - 'hill_q2': q2 - }) - return pd.DataFrame(results).set_index('sample') - -# --- 2. DATA PREPARATION --- -metrics = calculate_diversity_metrics(concat_df) -metrics = metrics.reset_index() - -# Ensure numeric columns are float -num_cols = metrics.select_dtypes(include=[np.number]).columns -metrics[num_cols] = metrics[num_cols].astype(float) - -# Merge with metadata -plot_df = pd.merge(metrics, meta, on='sample', how='inner') - -metrics_to_plot = [ - 'shannon_entropy', 'inverse_simpson', 'gini_coefficient', - 'hill_q0', 'hill_q1', 'hill_q2' -] - -# --- 3. SETUP GROUPS & ORDERING --- - -# Create Timepoint Mapping -if timepoint_col in plot_df.columns and timepoint_order_col in plot_df.columns: - time_order_map = plot_df.set_index(timepoint_col)[timepoint_order_col].to_dict() -else: - time_order_map = {} - -# Define Exclusions -exclude_cols = ['sample', 'file', 'total_counts'] -exclude_from_all_samples = [timepoint_col, timepoint_order_col, 'protocol_day', alias_col] - -# Groups for "All Samples" tab -group_opts_all = [ - col for col in meta.columns - if col not in exclude_cols - and col not in exclude_from_all_samples - and meta[col].nunique() < 35 -] - -# Groups for "By Patient" tab -group_opts_patient = [timepoint_col] if timepoint_col in plot_df.columns else [] - -unique_patients = sorted(plot_df[subject_col].dropna().unique().tolist()) - -# --- 4. PLOTTING FUNCTION --- - -def create_diversity_plot(data, x_col, y_metric, custom_order=None, plot_type="box"): - - if data.empty or data[x_col].dropna().empty: - return - - # Determine order - if custom_order: - unique_cats = [x for x in custom_order if x in data[x_col].unique()] - else: - unique_cats = sorted(data[x_col].dropna().unique().tolist()) - - # --- BASE PLOT --- - if plot_type == "box": - fig = px.box( - data, - x=x_col, - y=y_metric, - color=x_col, - points="all", - hover_data=['alias'], - category_orders={x_col: unique_cats}, - template="simple_white", - width=700, - height=600, - title=f"{y_metric} by {x_col}" - ) - - fig.update_traces(width=0.5, marker=dict(size=10, opacity=0.7)) - - for trace in fig.data: - if isinstance(trace, go.Box): - trace.pointpos = 0 - trace.jitter = 0.2 - - elif plot_type == "line": - # Sort for proper line connection - data = data.copy() - data[x_col] = pd.Categorical(data[x_col], categories=unique_cats, ordered=True) - data = data.sort_values(x_col) - - fig = px.scatter( - data, - x=x_col, - y=y_metric, - color_discrete_sequence=["black"], - hover_data=['alias'], - template="simple_white", - width=700, - height=500, - title=f"{y_metric} by {x_col}" - ) - - # Add dotted line - fig.add_trace(go.Scatter( - x=data[x_col], - y=data[y_metric], - mode='lines+markers', - line=dict(dash='dot', width=2), - marker=dict(size=18, color='black'), - showlegend=False - )) - - # --- STATISTICS (ONLY FOR BOXPLOTS) --- - stack_counter = 0 - - if plot_type == "box" and len(unique_cats) >= 2: - - pairs = list(combinations(unique_cats, 2)) - - y_max = data[y_metric].max() - y_min = data[y_metric].min() - y_range = y_max - y_min - if y_range == 0: - y_range = 1 - - base_offset = y_range * 0.08 # was 0.2 → MUCH closer to boxes - step_size = y_range * 0.06 # was 0.15 → tighter stacking - text_offset = y_range * 0.015 # was 0.03 → less gap above line - - for t1, t2 in pairs: - - g1 = data[data[x_col] == t1][y_metric].dropna() - g2 = data[data[x_col] == t2][y_metric].dropna() - - if len(g1) < 2 or len(g2) < 2: - continue - - try: - stat, p_value = mannwhitneyu(g1, g2, alternative='two-sided') - except ValueError: - continue - - if p_value >= 0.05: - continue - - if p_value < 0.001: - symbol = '***' - elif p_value < 0.01: - symbol = '**' - else: - symbol = '*' - - y_bracket = y_max + base_offset + (stack_counter * step_size) - y_text = y_bracket + text_offset - - fig.add_shape( - type="line", xref="x", yref="y", - x0=t1, y0=y_bracket, x1=t2, y1=y_bracket, - line=dict(color="black", width=1.5) - ) - - tick_len = y_range * 0.0015 - - fig.add_shape(type="line", xref="x", yref="y", - x0=t1, y0=y_bracket, x1=t1, y1=y_bracket - tick_len) - - fig.add_shape(type="line", xref="x", yref="y", - x0=t2, y0=y_bracket, x1=t2, y1=y_bracket - tick_len) - - try: - x_center = (unique_cats.index(t1) + unique_cats.index(t2)) / 2 - except: - x_center = 0 - - fig.add_annotation( - x=x_center, - y=y_text, - text=f"{symbol}
p={p_value:.3f}", - showarrow=False, - font=dict(size=10) - ) - - stack_counter += 1 - - extra_space = y_range * (0.15 + 0.07 * stack_counter) - fig.update_yaxes(range=[y_min, y_max + extra_space]) - - # --- FINAL LAYOUT --- - dynamic_height = 600 + (stack_counter * 40) - - fig.update_layout( - showlegend=False, - yaxis_title=y_metric, - xaxis_title=x_col, - xaxis=dict(tickfont=dict(size=15)), - height=dynamic_height, - margin=dict(t=80), - plot_bgcolor='rgba(0,0,0,0)' - ) - - fig.update_yaxes(showgrid=True, gridcolor='lightgrey') - - fig.show() - -# --- 5. EXECUTION LOOP --- -print(":::::: {.panel-tabset}\n") # Level 1 Start - -for metric in metrics_to_plot: - print(f"## {metric}\n") - - if len(unique_patients) <= 10: - print("::::: {.panel-tabset}\n") # Level 2 Start - - # --- TAB A: ALL SAMPLES --- - print("### All Samples\n") - print(":::: {.panel-tabset}\n") # Level 3 Start - for group in group_opts_all: - print(f"#### {group}\n") - create_diversity_plot(plot_df, group, metric) - print("\n") - print("::::\n") # Level 3 End - - # --- TAB B: BY PATIENT --- - print("### By Patient\n") - print(":::: {.panel-tabset}\n") # Level 3 Start - - for pat in unique_patients: - print(f"#### {pat}\n") - pat_df = plot_df[plot_df[subject_col] == pat] - - has_multi_origin = ( - 'origin' in pat_df.columns and - pat_df['origin'].nunique() > 1 - ) - - if has_multi_origin: - print("::: {.panel-tabset}\n") # Level 4 Start - origins = sorted(pat_df['origin'].dropna().unique()) - - for origin in origins: - print(f"##### {origin}\n") - origin_df = pat_df[pat_df['origin'] == origin] - - for group in group_opts_patient: - if origin_df[group].nunique() > 0: - custom_order = None - if group == timepoint_col and time_order_map: - pat_tps = origin_df[group].dropna().unique().tolist() - custom_order = sorted( - pat_tps, - key=lambda x: time_order_map.get(x, 999) - ) - create_diversity_plot( - origin_df, group, metric, - custom_order, plot_type="line" - ) - else: - print(f"No data for {group} in {origin}, patient {pat}.") - print("\n") - - print(":::\n") # Level 4 End ← was missing in some branches - - else: - for group in group_opts_patient: - if pat_df[group].nunique() > 0: - custom_order = None - if group == timepoint_col and time_order_map: - pat_tps = pat_df[group].dropna().unique().tolist() - custom_order = sorted( - pat_tps, - key=lambda x: time_order_map.get(x, 999) - ) - create_diversity_plot( - pat_df, group, metric, - custom_order, plot_type="line" - ) - else: - print(f"No data for {group} in patient {pat}.") - print("\n") - - print("::::\n") # Level 3 End - print(":::::\n") # Level 2 End ← moved OUTSIDE patient loop - - else: - # Standard View (>10 patients) - print("::::: {.panel-tabset}\n") # Level 2 Start - for group in group_opts_all: - print(f"### {group}\n") - create_diversity_plot(plot_df, group, metric) - print("\n") - print(":::::\n") # Level 2 End - -print("::::::\n") # Level 1 End - -``` -**Figure 1. Comparative Analysis of TCR Repertoire Diversity Metrics.** Boxplots display the distribution of metrics across samples. The interactive panel allows toggling between diversity indices and grouping variables. This visualization facilitates the assessment of repertoire diversity and evenness across different experimental conditions and biological replicates. - -**Use the tabs to toggle between metadata groupings:** - -- The distributions of **Batches** L0, L1, etc., should largely overlap. If L0 is consistently higher than L1 across all sample types, you have a Batch Effect. You cannot compare samples across batches without computational correction (we do not provide such correction). - -- If one **patient** (e.g. Patient01) has drastically lower diversity than all others across all timepoints and tissues, it is likely a biological anomaly or a collection issue (degraded sample). Treat this patient as a separate cohort or exclude them. - -- In longitudinal studies, diversity often drops slightly over time due to treatment. Diversity metrics rarely jump from 0.8 to 0.1 (or vice versa) in a short window **between timepoints**. If so, that sample is suspect. There might be some causes that could explain it: 1) Sample swap of sample timepoint labels for that patient. 2) The patient contracted a viral infection at T2, causing a massive expansion of non-tumor clones. - -## Hill diversity Profile - -It allows to see how "TCR diversity" changes depending on how much you care about rare clones vs. dominant clones. You use Hill Profiles to catch artifacts that single metrics hide. -$q = 0$ (Richness): Ignore frequency, every unique CDR3 counts as 1. -$q = 1$ (Shannon Entropy): Clones are weighted by their frequency. -$q = 2$ (Simpson): Rare clones are mathematically ignored. The top expanded clones dictate the score. - -Two patients might share the exact same Shannon Entropy score, yet one has a healthy, broad repertoire while the other is battling a massive leukemia clone, for example. The Hill Diversity Profile reveals this difference by visualizing the trade-off between Richness (total number of clones) and Evenness (how equally those clones are distributed) in a single curve. - -**How to Read the Curve?** -Think of the **X-axis** ($q$) as a **"sensitivity dial" for dominance**. At the far left ($q=0$), the metric cares only about presence; it counts every unique clone equally, regardless of whether it appears once or a million times. This point represents your total Richness. **As you move right toward $q=2$, the mathematical weight shifts drastically toward the most abundant clones**, effectively ignoring the rare ones. **The slope of the line tells you the story**: a flat line indicates a perfectly even population where everyone is equal, while **a steep**, crashing slope **reveals a population dominated by a few massive "bully" clones.** - -**Here are 2 possible scenarios:** - -- **Scenario A:** The "Depth Trap" -Imagine you are **comparing Sample A and Sample B**, and you immediately notice that **Sample A has a much higher Hill $q=0$ value**, which represents richness. You might be **tempted to conclude that Sample A is biologically "more diverse"** and has a healthier T cell repertoire. However, when you look at the **Hill Profile as it moves toward $q=1$ and $q=2$**, you see the **lines rapidly converge or even cross, showing no difference in the dominant clones.** This reveals that the **"diversity" in Sample A was an illusion caused by sequencing depth.** You simply spent more money sequencing Sample A, allowing the machine to pick up more rare, single-copy clones that were missed in Sample B. The structure of the immune system was identical in both; one was just measured with a magnifying glass while the other was measured with the naked eye. - -- **Scenario B:** Real Biology -Consider a case where you are analyzing a **"Responder" patient versus a "Non-Responder" in a cancer trial**. If you only looked at a single metric like **Shannon Entropy ($q=1$), the two patients might look identical, leading you to think the treatment had no effect on the repertoire structure.**. But the Hill Profile tells a different story. At $q=0$, the Non-Responder is much higher, indicating a vast, unfocused collection of rare T cells that aren't doing much. As you move to **$q=2$, the lines cross, and the Responder suddenly shoots up**. This crossover reveals that **while the Non-Responder has more "types" of cells, the Responder has successfully expanded a specific "army" of effector clones** to fight the tumor**. The profile proves that their immune systems are structurally opposite, a critical biological insight that a single summary statistic would have completely hidden. - -```{python} -#| output: asis -#| fig-width: 6 -#| fig-height: 4.5 -#| label: hill - -# --- 1. Setup Data & Grouping --- -hill_cols = {'hill_q0': 0, 'hill_q1': 1, 'hill_q2': 2} - -exclude_cols = ['sample', 'file', 'total_counts', 'filename', 'sample_id', 'alias', - 'shannon_entropy', 'gini_coefficient', 'inverse_simpson', - 'hill_q0', 'hill_q1', 'hill_q2', 'mean_len', 'std_len'] - -group_opts = [ - col for col in plot_df.columns - if col not in exclude_cols and plot_df[col].nunique() < 35 -] - -sns.set_theme(style="whitegrid", context="paper", font_scale=1.1) - -# --- 2. Start Quarto Tabs --- -print("::: {.panel-tabset}") - -for col in group_opts: - print(f"## {col}") - - # Data Prep - subset = plot_df[[col, 'alias', 'hill_q0', 'hill_q1', 'hill_q2']].copy() - - melted = subset.melt( - id_vars=['alias', col], - value_vars=['hill_q0', 'hill_q1', 'hill_q2'], - var_name='q_metric', - value_name='Diversity' - ) - melted['q'] = melted['q_metric'].map(hill_cols) - melted[col] = melted[col].astype(str) - - # --- Y-AXIS SYNCHRONIZATION --- - # Calculate global min/max for this tab to lock axes - y_min = melted['Diversity'].min() - y_max = melted['Diversity'].max() - - # Add 5% log-padding so points aren't cut off - log_min = np.log10(y_min) - log_max = np.log10(y_max) - pad = (log_max - log_min) * 0.05 - - # Matplotlib limits (Linear values) - mpl_ylim = (10**(log_min - pad), 10**(log_max + pad)) - - # Plotly limits (Log10 values) - plotly_yrange = [log_min - pad, log_max + pad] - - # --- COLOR CONSISTENCY --- - unique_cats = sorted(melted[col].unique()) - palette_tuples = sns.color_palette("viridis", n_colors=len(unique_cats)) - hex_colors = [f'#{int(r*255):02x}{int(g*255):02x}{int(b*255):02x}' for r,g,b in palette_tuples] - color_map = dict(zip(unique_cats, hex_colors)) - - # =========================== - # PLOT 1: SEABORN (Summary) - # =========================== - plt.figure(figsize=(6, 4.5)) # Standard size - - sns.lineplot( - data=melted, - x='q', y='Diversity', - hue=col, style=col, - markers=True, dashes=False, - linewidth=2.5, - palette=color_map, - err_style='band', errorbar=('ci', 95) - ) - - plt.yscale('log') - plt.ylim(mpl_ylim) - plt.xticks([0, 1, 2], ['q=0\n(Richness)', 'q=1\n(Shannon)', 'q=2\n(Simpson)']) - plt.title(f'Average Profile by {col} (with 95% CI)') - plt.ylabel('Effective Clones (Log Scale)') - plt.xlabel('') - - plt.legend(bbox_to_anchor=(1.02, 1), loc='upper left', title=col, - fontsize='x-small', title_fontsize='small', frameon=False) - plt.tight_layout() - plt.show() - - # =========================== - # PLOT 2: PLOTLY (Interactive) - # =========================== - - fig = px.line( - melted, - x='q', - y='Diversity', - color=col, - line_group='alias', # Keep sample here to ensure lines are grouped by the unique identifier - hover_name='alias', - - hover_data={col: True, 'q': False, 'Diversity': ':.2f', 'alias': True}, - log_y=True, - title=f'Individual Sample Profiles (Hover to Identify)', - color_discrete_map=color_map - ) - - fig.update_traces( - mode="lines+markers", - line=dict(width=3), - opacity=0.6, - marker=dict(size=6) - ) - - fig.update_layout( - template="simple_white", - width=600, - height=450, - xaxis=dict( - tickmode='array', - tickvals=[0, 1, 2], - ticktext=['q=0 (Richness)', 'q=1 (Shannon)', 'q=2 (Simpson)'], - title='Sensitivity (q)' - ), - yaxis=dict( - title='Effective Clones (Log Scale)', - range=plotly_yrange - ), - legend=dict(title=col) - ) - - fig.show() - - print("\n") - -print(":::") - -``` -**Figure 2. Hill Diversity Profiles of TCR Repertoires.** Diversity profiles display the effective number of clones (log scale) across sensitivity orders $q=0$ (Richness), $q=1$ (Shannon), and $q=2$ (Simpson). The top panel illustrates the mean repertoire structure for each subject with 95% confidence intervals, highlighting differences in evenness and dominance. The bottom panel visualizes individual sample trajectories, revealing heterogeneity and outlier profiles within each cohort. - -**Using the above metrics to identify technical failures 🚩** -Before attempting biological interpretation, you must validate that your metrics reflect immune biology and not library preparation artifacts. Use these metrics to flag and exclude compromised samples. - -**1. Amplification Bias** -**Symptom:** A sample with moderate Richness ($q=0$) but extreme Dominance ($q=2$). -**Metric Signature:** - -- Shannon ($H$): Disproportionately low compared to the cohort baseline. -- Gini Coefficient: Approaches 1.0 (>0.9 is suspicious in non-tumor samples). -- Hill Profile: The curve starts high at $q=0$ (indicating presence of unique sequences) but crashes vertically as $q \to 2$. - -**Diagnosis:** A "Jackpot" event occurred where a random RNA molecule was preferentially amplified during PCR. The "richness" is likely sequencing noise/errors, while the reads are consumed by a single artifact. -**Action:** Discard sample. - -**2. Sequencing Bias** -**Symptom:** Two samples appear to have different diversities, but the difference disappears when you ignore rare clones. -**Metric Signature:** - -- Hill Profile: Sample A has a much higher $q=0$ (Richness) than Sample B, but the lines converge at $q=1$ and $q=2$. - -**Diagnosis:** This is not a biological difference. Sample A was simply sequenced deeper, detecting more single-copy clones (noise). The structural effective diversity is identical. -**Action:** Do not claim Sample A is "more diverse." Rely on Shannon ($q=1$) or Simpson ($q=2$) for comparison, or downsample (rarefy) the libraries to equal depth. - -**3. Library Failure** -**Symptom**: A sample that looks "empty" compared to the cohort. -**Metric Signature:** - -- Richness ($R$): Drastically lower than the cohort mean. -- Hill Profile: A flat line hovering near the bottom of the Y-axis. -**Diagnosis:** Poor RNA extraction, degradation, or a biopsy that captured mostly fat/stroma (low T-cell content). -**Action:** Exclude from analysis. - -**4. Batch Effects** -**Symptom:** Diversity metrics cluster by processing date rather than biological group. -**Metric Signature:**Visual Check: Boxplots of Shannon Entropy grouped by Batch_ID. If Batch 0 is consistently higher than 1 across all sample types (tumor, blood, healthy), you have a technical batch effect. -**Action:** You cannot compare raw metrics across these batches. - - -**Biological Interpretation of the above metrics 🦠** -Once QC is cleared, these metrics describe the shape of the immune repertoire. - -**1. Naive vs. Reactive Profiles** -**Naive / Resting State:** - -- Signature: High Richness, High Shannon, Low Gini (< 0.3). -- Interpretation: The army is standing by. There is no dominant clone because no specific threat has triggered an expansion. The Hill profile will show a "gentle decline." -**Reactive / Tumor Infiltrating:** - -- Signature: Low Shannon, High Inverse Simpson, High Gini (> 0.6). -- Interpretation: Clonal Expansion. The immune system has "picked a winner." A few specific clones (e.g., tumor-reactive or viral-specific) have expanded to occupy a large fraction of the repertoire (20–50%). - -**2. Inverse Simpson** -The Inverse Simpson Index ($1/D$) converts the probability of coincidence into a concrete "count" of effective clones. It filters out the long tail of rare sequences to quantify the number of clones driving the response. -**High Effective Number** ($1/D \gg 100$): - -- Signature: $1/D$ is close to Richness ($R$; i.e. hill_q0). -- Interpretation: Polyclonal/Diverse. The sample behaves as if it is composed of hundreds or thousands of equally abundant clones. Resources are distributed broadly. - -**Low Effective Number ($1/D < 20$):** - -- Signature: $1/D$ is a small fraction of Richness ($R$). -- Interpretation: Oligoclonal/Focused. Despite potentially containing thousands of unique sequences ($R$), the sample behaves biologically as if it only contains ~10-20 active clones. This confirms a highly focused effector response. - -**3. Hill Profile** -This is the most powerful visualization for distinguishing Responder vs. Non-Responder dynamics in longitudinal data. -**The Scenario:** A Non-Responder may have higher Richness ($q=0$) due to a chaotic, unfocused repertoire. A Responder might have lower Richness but higher Dominance. -**The Signature:** The Hill curves will cross.At $q=0$ (Richness): Non-Responder is higher.At $q=2$ (Dominance): Responder shoots up. -**Interpretation**: The crossover proves the Responder has successfully funneled resources into a specific "army" of effector clones, structurally reorganizing the repertoire to fight the tumor. - -## Combining metrics {#sec-combine-metrics} - -### Evenness Check (Gini vs Shannon) - -Under normal biological conditions, the Gini coefficient and Shannon entropy should exhibit a strong negative correlation. As the immune system focuses on a specific threat, a few clones expand to dominate the repertoire, causing inequality to rise (High Gini) while overall diversity naturally falls (Low Shannon). - -```{python} -#| output: asis -#| fig-width: 6 -#| fig-height: 5 - -# 1. Setup Grouping Options -exclude_cols = ['sample', 'file', 'total_counts', 'timepoint_order', 'filename', 'sample_id', - 'shannon_entropy', 'gini_coefficient', 'inverse_simpson', - 'hill_q0', 'hill_q1', 'hill_q2', 'mean_len', 'std_len'] - -group_opts = [ - col for col in plot_df.columns - if col not in exclude_cols and plot_df[col].nunique() < 35 -] - -# 2. Start Quarto Tabset -print("::: {.panel-tabset}") - -for col in group_opts: - print(f"## {col}") - - fig = px.scatter( - plot_df, - x='shannon_entropy', - y='gini_coefficient', - color=col, - symbol=col, - template="simple_white", - width=600, - height=500, - title=f'Evenness by {col}', - hover_data={ - 'sample': True, - 'alias': True if 'alias' in plot_df.columns else False, - 'shannon_entropy': ':.3f', - 'gini_coefficient': ':.3f', - col: True - } - ) - - # Make points look like your seaborn version - fig.update_traces( - marker=dict(size=10, line=dict(width=1, color='black')), - opacity=0.8 - ) - - fig.update_layout( - xaxis_title='Shannon (Diversity)', - yaxis_title='Gini (Inequality)', - legend_title=col - ) - - fig.show() - print("\n") - -print(":::") - -``` - -**Figure 3. Evaluation of Repertoire Structure via Gini-Shannon Correlation.** This scatter plot contrasts the Gini coefficient (inequality) against Shannon entropy (diversity) for individual samples, colored by metadata variables. The visualization assesses the expected inverse relationship between the two metrics, where antigen-driven clonal expansion typically results in higher inequality and reduced diversity. - -Healthy, **naive, or polyclonal samples** will naturally cluster in the **bottom-right** of the plot. These repertoires are **characterized by high diversity** ($H > 8$) and low inequality ($G < 0.2$), reflecting a deep pool of clones with relatively even frequencies. In contrast, **samples undergoing an active immune response**—such as those from tumors, acute infections, or leukemia—will drift toward the **top-left**. In this "Expansion Zone," the biological dominance of a few clones drives the Gini coefficient up ($> 0.6$) and suppresses the entropy, a signature that is expected in disease states but signals potential PCR bias ("jackpotting") if observed in negative controls. - -The most critical **QC artifact to watch** for is the "Ghost Library" in the **bottom-left corner**. These samples show **low diversity but paradoxically low inequality**. This almost always indicates severe **undersampling**: if a library only contains 50 reads and every read is unique, the sample will appear "perfectly equal" (Gini $\approx$ 0) simply because it lacks the depth to reveal the true distribution. These samples are statistical phantoms and should be **discarded**. - -Finally, samples falling in the **top-right (High Diversity + High Inequality)** represent a **mathematical contradiction for T-cell data**. It is impossible to have a massive, diverse "tail" of clones while simultaneously having a single clone occupy the majority of reads. Samples in this region often indicate **bioinformatic errors**, such as merging incompatible library files or calculating metrics on raw counts rather than frequencies. - -### Richness vs 1/D - -The Inverse Simpson ($1/D$) is used to visualize the magnitude of the response. It ignores the long tail of rare clones. It is a good practice to compare $1/D$ to Richness ($R$). -- If $R = 5000$ and $1/D = 4500$: The sample is polyclonal (many clones of equal size). -- If $R = 5000$ and $1/D = 2$: The sample is oligoclonal. Despite having 5000 unique sequences, it behaves biologically as if it only has 2 clones. - -```{python} -#| output: asis -#| fig-width: 6 -#| fig-height: 5 - -# 1. Setup Grouping Options -exclude_cols = ['sample', 'file', 'total_counts', 'timepoint_order', 'filename', 'sample_id', - 'shannon_entropy', 'gini_coefficient', 'inverse_simpson', - 'hill_q0', 'hill_q1', 'hill_q2', 'mean_len', 'std_len'] - -# Automatically detect categorical columns with reasonable cardinality -group_opts = [ - col for col in plot_df.columns - if col not in exclude_cols and plot_df[col].nunique() < 35 -] - -# ============================================================================== -# PLOT SET 1: RICHNESS vs INVERSE SIMPSON (The "Depth Trap" Check) -# ============================================================================== -print("::: {.panel-tabset}") - -for col in group_opts: - print(f"## {col}") - - fig = px.scatter( - plot_df, - x='inverse_simpson', - y='hill_q0', - color=col, - symbol=col, - template="simple_white", - width=600, - height=500, - title=f'Richness vs. Effective Clones by {col}', - hover_data={ - 'sample': True, - 'alias': True if 'alias' in plot_df.columns else False, - 'inverse_simpson': ':.3f', - 'hill_q0': ':.3f', - col: True - } - ) - - # Match seaborn styling - fig.update_traces( - marker=dict(size=10, line=dict(width=1, color='black')), - opacity=0.8 - ) - - # --- Add diagonal line (y = x) --- - max_val = max( - plot_df['hill_q0'].max(), - plot_df['inverse_simpson'].max() - ) - - fig.add_trace(go.Scatter( - x=[0, max_val], - y=[0, max_val], - mode='lines', - line=dict(dash='dash', color='gray'), - name='Perfect Evenness (y=x)', - hoverinfo='skip' # keeps hover clean - )) - - # Layout - fig.update_layout( - xaxis_title='Inverse Simpson (Effective # Clones)', - yaxis_title='Hill q=0 (Total Richness)', - legend_title=col - ) - - fig.show() - print("\n") - -print(":::") -print("\n") - -``` - -**Figure 4. Analysis of Clonal Dominance vs. Effective Clone Count.** This scatter plot compares the total clonal richness against the effective number of clones, colored by metadata variables. The diagonal line ($y=x$) represents a theoretical state of perfect evenness, where every clone has equal frequency. Deviations from this line indicate the degree of clonal dominance. - -### Gini Coefficient vs. Effective Clones - -This plot maps the trade-off between unevenness and effective size. Typically, these metrics are inversely correlated: as inequality (Gini) rises, the number of effective clones (Inverse Simpson) drops. Deviations from this curve reveal samples with unusual structural properties - -```{python} -#| output: asis -#| fig-width: 6 -#| fig-height: 5 - -print("::: {.panel-tabset}") - -for col in group_opts: - print(f"## {col}") - - # Create Plotly scatter - fig = px.scatter( - plot_df, - x='inverse_simpson', - y='gini_coefficient', - color=col, - symbol=col, - template="simple_white", - width=600, - height=500, - title=f'Inequality vs. Diversity by {col}', - hover_data={ - 'sample': True, - 'alias': True if 'alias' in plot_df.columns else False, - 'inverse_simpson': ':.3f', - 'gini_coefficient': ':.3f', - col: True - }, - color_continuous_scale='magma' if plot_df[col].dtype != object else None - ) - - # Match marker style from seaborn - fig.update_traces( - marker=dict(size=10, line=dict(width=1, color='black')), - opacity=0.8 - ) - - # Fix y-axis for Gini (bounded [0,1]) - fig.update_yaxes(range=[0, 1.05], title='Gini Coefficient (Inequality)') - - # x-axis label - fig.update_xaxes(title='Inverse Simpson (Effective # Clones)') - - # Layout adjustments - fig.update_layout( - legend_title=col, - title=dict(font=dict(size=14)), - width=600, - height=500 - ) - - fig.show() - print("\n") - -print(":::") - -``` - -**Figure 5. Relationship Between Repertoire Inequality and Effective Diversity**. The visualization demonstrates the inverse correlation typically observed in immune repertoires, where high inequality (dominance by a few clones) corresponds to a low effective population size. Deviations from this structural trend can identify samples with aberrant clonal distributions. - -**Combined metrics: Quality Control 🚩** -Use these 2D plots to distinguish between sequencing artifacts and true biological signal. - -**1. Richness vs. Effective Clones** -This plot exposes the "useless tail" of your library—the sequences that appear only once or twice and do not contribute to the effective immune response. - -- The Diagonal Line ($y=x$): A sample on this line has no rare clones; every unique sequence is equally abundant. -- Vertical Separation (The Warning): If Sample A is far above Sample B on the Y-axis (Richness) but they align perfectly on the X-axis (Effective Clones), Sample A is not more diverse. It was simply sequenced deeper. The extra richness is noise (singletons). -- Action: Trust the X-axis (Inverse Simpson) for biological comparison. Ignore the Y-axis offset. Samples clustered in the bottom-left corner ($<100$ on both axes) often indicate failed library preparation or low-input samples (e.g., fibrous tissue with no T-cells). - -**2. Gini vs. Effective Clones Plot** -This plot validates the structural integrity of the library. - -- **The "Impossible" Quadrant:** You should rarely see High Gini (>0.8) combined with High Effective Clones (>1000). Mathematically, you cannot have extreme inequality while simultaneously maintaining thousands of effectively equal clones. If you see this, check for chimeric reads or alignment errors. - -- **The "False Dominance":** If a Baseline/PBMC sample appears in the Top-Left corner (High Gini, Low Effective Count), it is a red flag. Unless the patient has a known blood malignancy, a resting control sample should never look like a tumor. This indicates a PCR Jackpot artifact. - - -**Combined metrics: Biological Interpretation 🦠 ** -Once QC is cleared, the position of a sample in these 2D spaces defines its immunological state. - -**1. Gini vs. Effective Clones (1/D)** -**The Naive/Polyclonal State (Bottom-Right):** -The immune system is in surveillance mode. No specific antigen has triggered a response. Resources are distributed broadly across thousands of clones. -**The Reactive/focused State (Top-Left):** -Clonal Expansion. The immune system has identified a threat (tumor/virus) and "picked a winner." The repertoire is dominated by a focused army of effector clones, physically crowding out the diversity. - -**2. Tracking Response** -In a successful immunotherapy response, track the patient's movement across the plot over time: -**The "Focusing" Shift:** A responder should move diagonally Up and Left. They start with a broad repertoire (Bottom-Right) and shift toward dominance (Top-Left) as tumor-reactive clones expand. -**The "Ignorance" Stasis:** Non-responders often remain stuck in the Bottom-Right quadrant (high diversity, low inequality), indicating the treatment failed to trigger T-cell expansion. - -## TCR Overlap {#sec-tcr-overlap} - -We use two distinct metrics to compare repertoires. One measures membership (who is there?), and the other measures structure (who is dominant?). - -- **Jaccard Index:** Measures the overlap of unique sequences, ignoring frequency. It treats a singleton and a top clone equally. -**Question:** "What fraction of the unique CDR3s are found in both samples?" -**Range:** $0$ (No shared sequences) $\to$ $1$ (Identical unique sequence list). -**Formula:** -$$J(A,B) = \frac{|A \cap B|}{|A \cup B|}$$ -Where $|A \cap B|$ is the count of unique clonotypes shared by both samples, and $|A \cup B|$ is the total count of unique clonotypes in the union of both samples. - -```{python} -#| output: asis -#| fig-width: 7 -#| fig-height: 7 - -# --- 1. CALCULATION FUNCTION (Only Jaccard) --- -def calculate_overlaps(df): - samples = sorted(df['sample'].unique()) - clones = sorted(df['CDR3b'].unique()) - - sample_map = {s: i for i, s in enumerate(samples)} - clone_map = {c: i for i, c in enumerate(clones)} - - row_idx = df['sample'].map(sample_map).values - col_idx = df['CDR3b'].map(clone_map).values - values = df['counts'].values - - n_samples = len(samples) - n_clones = len(clones) - - mat = sparse.coo_matrix((values, (row_idx, col_idx)), shape=(n_samples, n_clones)).tocsr() - - # Jaccard (Binary) - mat_bin = mat.copy() - mat_bin.data[:] = 1 - intersection = mat_bin.dot(mat_bin.T).toarray() - row_sums = mat_bin.sum(axis=1).A1 - union_matrix = row_sums[:, None] + row_sums[None, :] - intersection - - jaccard = np.divide(intersection, union_matrix, - out=np.zeros_like(intersection, dtype=float), - where=union_matrix!=0) - - return pd.DataFrame(jaccard, index=samples, columns=samples) - - -# --- 2. EXECUTE CALCULATION --- -jaccard_df = calculate_overlaps(concat_df) - -# --- 3. RENAME TO ALIAS --- -if alias_col in meta.columns: - sample_to_alias = meta.set_index('sample')[alias_col].to_dict() - jaccard_df = jaccard_df.rename(index=sample_to_alias, columns=sample_to_alias) - -# --- 4. CLUSTERING HELPER --- -def get_clustered_order(matrix): - """Return hierarchical clustering order of samples.""" - dist = pdist(matrix.values, metric='cityblock') - link = linkage(dist, method='average') - leaves = leaves_list(link) - return matrix.index[leaves].tolist() - -# --- 5. INTERACTIVE PLOTTING FUNCTION --- -def plot_interactive_heatmap(matrix, title): - try: - ordered_labels = get_clustered_order(matrix) - matrix = matrix.loc[ordered_labels, ordered_labels] - except Exception as e: - print(f"Clustering failed, using default order: {e}") - - vals = matrix.values - mask = ~np.eye(vals.shape[0], dtype=bool) - max_val = vals[mask].max() if mask.any() else 1.0 - if max_val < 0.05: max_val = 0.1 - - fig = px.imshow( - matrix, - labels=dict(x="Sample", y="Sample", color="Jaccard"), - x=matrix.columns, - y=matrix.index, - color_continuous_scale="Reds", - range_color=[0, max_val], - title=f"{title}
(Max Overlap: {max_val:.3f})" - ) - - fig.update_layout( - width=700, - height=700, - xaxis=dict(tickangle=90, tickfont=dict(size=10)), - yaxis=dict(tickfont=dict(size=10)), - plot_bgcolor='black' - ) - fig.update_traces(xgap=1, ygap=1) - fig.show() - -# --- 6. OUTPUT TABSET (Only Jaccard) --- -print("::: {.panel-tabset}\n") -print("## Jaccard\n") -plot_interactive_heatmap(jaccard_df, "Jaccard Index") -print("\n") -print(":::\n") - -``` -**Figure 6. Assessment of Clonal Overlap and Cross-Contamination.** -Heatmaps display the pairwise similarity between TCR repertoires using the Jaccard (overlap) index. Hierarchical clustering groups samples based on their shared clonal content, facilitating the identification of biological replicates or potential cross-contamination events. The color intensity reflects the degree of similarity, where darker red indicates a higher proportion of shared clonotypes between sample pairs. - -::: {.callout-tip title="Note"} -When calculating Jaccard for this specific purpose (contamination), we are calculating it on the nucleotide sequence, not the amino acid sequence, **whenever the TCRseq was done using DNA as input**. Convergent recombination (different DNA making the same Amino Acid) is common in high-depth samples and will give you a false positive "contamination" signal if you stick to Amino Acids. Contamination should be measured at the physical level (DNA). -::: - -**Technical QC 🚩** -Use these overlap metrics to detect the two most dangerous errors in sequencing: **Sample Swapping and Physical Contamination.** - -- **Contamination Detection: Jaccard Index** - - Red Flag: A block of red squares connecting unrelated samples (Overlap $> 0.05 - 0.1$ Pink/Red). - - Diagnosis: Cross-Contamination. Because TCRs are hypervariable, distinct individuals should almost never share rare sequences (Jaccard should be near 0). - - Scenario A (Splash): Patient_A and Patient_B share 20% of their unique sequences. Physical liquid likely splashed between wells during PCR setup. - - Scenario B (Barcode Hopping): If an entire batch shows faint overlap (~5-10%) with a high-concentration sample, it indicates index hopping during sequencing. - - Action: Flag samples with unexpected Jaccard scores $>0.05$. - - -**Biological Interpretation 🦠** -Once technical artifacts are ruled out, overlap quantifies stability and shared biology. - -- **Public Clones: Jaccard** - - Low but Real Overlap (< 0.01):While distinct humans rarely share repertoires, they may share identical "Public Clones" (often viral-specific, e.g., CMV or EBV). - - Interpretation: A very faint Jaccard signal (0.001 - 0.01) across a cohort can indicate a shared antigen exposure (e.g., a common viral infection in the population) or convergent recombination events. - - -# Technical QC Relevance (🚩) {#sec-qc} -The following section presents metrics designed to assess the processing quality of the samples. These technical QC checks are essential for differentiating between genuine biological signals and technical artifacts, including library failures, undersampling, batch effects, and outliers. However, we acknowledge that some of these metrics can also carry biological meaning, depending on the specific biological context of the analysis. - -## Percent of productive TCRs {#sec-cdr3-prod} - -The percent of productive TCRs is one of the first metrics you should check to validate your data's reliability before proceeding to any biological analysis. - -As a Quality Filter: You should set a QC threshold (e.g., >75% productive reads). Samples that fail to meet this threshold should be flagged for review or potentially excluded from the analysis. Drawing conclusions about T-cell diversity or clonality from a sample with low productivity is unreliable, as the data is likely noisy and not a true representation of the functional T-cell repertoire. - -To Troubleshoot Experiments: If you find that an entire batch of samples has a low productive percentage, it points to a systematic issue in your experimental protocol, most commonly an ineffective gDNA removal step or problems with RNA integrity. - -To Ensure Confidence in Results: By confirming that your samples have a high percentage of productive TCRs, you can be confident that the clonotypes you identify and analyze represent real, functional T-cells that are actively participating in the immune response. - -```{python} -#| output: asis -#| echo: false - -import plotly.express as px -import plotly.graph_objects as go -from scipy.stats import mannwhitneyu -from itertools import combinations -import pandas as pd - -# Create a mapping for timepoint order using the 'meta' dataframe -if 'meta' in locals() and timepoint_col in meta.columns and timepoint_order_col in meta.columns: - time_order_map = meta.set_index(timepoint_col)[timepoint_order_col].to_dict() -else: - # Fallback: if 'meta' isn't loaded, try to use columns in the main df - if timepoint_order_col in df.columns: - time_order_map = df.set_index(timepoint_col)[timepoint_order_col].to_dict() - else: - time_order_map = {} # No sorting map available - -# 1. Define General Exclusions -exclude_cols = ['sample', 'file', 'clonality', 'counts', 'total_counts', - 'num_clones', 'num_TCRs', 'simpson_index', 'simpson_index_corrected', - 'num_prod', 'num_nonprod', 'pct_prod', 'pct_nonprod', - 'productive_cdr3_avg_len', 'num_convergent','ratio_convergent'] - -# 2. Define Groups for "All Samples" Tab (Exclude timepoint/order) -exclude_from_all_samples = [timepoint_col, timepoint_order_col, 'protocol_day', alias_col] - -group_opts_all = [ - col for col in df.columns - if col not in exclude_cols - and col not in exclude_from_all_samples - and df[col].nunique() < 35 -] - -# 3. Define Groups for "By Patient" Tab (Only timepoint) -group_opts_patient = [timepoint_col] if timepoint_col in df.columns else [] - -unique_patients = sorted(df[subject_col].dropna().unique().tolist()) - -# --- 1. DEFINE PLOTTING FUNCTION --- -def create_pct_prod_plot(data, group_col, custom_order=None, plot_type="box"): - """ - custom_order: Optional list of x-axis values in the desired order. - plot_type: "box" for standard comparisons, "line" for patient timelines. - """ - if data.empty or data[group_col].dropna().empty: - return - - # Determine order: Use custom list if provided, otherwise sort alphabetically - if custom_order: - unique_cats = [x for x in custom_order if x in data[group_col].unique()] - else: - unique_cats = sorted(data[group_col].dropna().unique().tolist()) - - # --- BASE PLOT --- - if plot_type == "box": - fig = px.box( - data, - x=group_col, - y='pct_prod', - color=group_col, - points='all', - hover_data=['alias'], - category_orders={group_col: unique_cats}, - template="simple_white", - title=f"Productive TCRs (%) by {group_col}" - ) - - fig.update_traces(width=0.5, marker=dict(size=10, opacity=0.7)) - - for trace in fig.data: - if isinstance(trace, go.Box): - trace.pointpos = 0 - trace.jitter = 0.2 - - elif plot_type == "line": - # Sort for proper line connection - data = data.copy() - data[group_col] = pd.Categorical(data[group_col], categories=unique_cats, ordered=True) - data = data.sort_values(group_col) - - fig = px.scatter( - data, - x=group_col, - y='pct_prod', - color_discrete_sequence=["black"], - hover_data=['alias'], - template="simple_white", - title=f"Productive TCRs (%) by {group_col}" - ) - - # Add dotted line - fig.add_trace(go.Scatter( - x=data[group_col], - y=data['pct_prod'], - mode='lines+markers', - line=dict(dash='dot', width=2), - marker=dict(size=18, color='black'), - showlegend=False - )) - - # --- STATISTICS (ONLY FOR BOXPLOTS) --- - stack_counter = 0 - - if plot_type == "box" and len(unique_cats) >= 2: - pairs = list(combinations(unique_cats, 2)) - - y_max = data['pct_prod'].max() - y_min = data['pct_prod'].min() - y_range = y_max - y_min - if y_range == 0: - y_range = 1 - - # Updated to tighter spacing format - base_offset = y_range * 0.08 - step_size = y_range * 0.06 - text_offset = y_range * 0.015 - - for t1, t2 in pairs: - group1 = data[data[group_col] == t1]['pct_prod'].dropna() - group2 = data[data[group_col] == t2]['pct_prod'].dropna() - - if len(group1) < 2 or len(group2) < 2: - continue - - try: - stat, p_value = mannwhitneyu(group1, group2, alternative='two-sided') - except ValueError: - continue - - if p_value >= 0.05: - continue - - if p_value < 0.001: symbol = '***' - elif p_value < 0.01: symbol = '**' - else: symbol = '*' - - y_bracket = y_max + base_offset + (stack_counter * step_size) - y_text = y_bracket + text_offset - - # Draw Bracket Line - fig.add_shape( - type="line", xref="x", yref="y", - x0=t1, y0=y_bracket, x1=t2, y1=y_bracket, - line=dict(color="black", width=1.5) - ) - - # Draw Bracket Ticks - tick_len = y_range * 0.015 - fig.add_shape(type="line", xref="x", yref="y", - x0=t1, y0=y_bracket, x1=t1, y1=y_bracket - tick_len, - line=dict(color="black", width=1.5)) - fig.add_shape(type="line", xref="x", yref="y", - x0=t2, y0=y_bracket, x1=t2, y1=y_bracket - tick_len, - line=dict(color="black", width=1.5)) - - try: - x_center = (unique_cats.index(t1) + unique_cats.index(t2)) / 2 - except: - x_center = 0 - - # Add P-Value Text - fig.add_annotation( - x=x_center, y=y_text, text=f"{symbol}
p={p_value:.3f}", - showarrow=False, font=dict(size=10, color="black") - ) - - stack_counter += 1 - - extra_space = y_range * (0.15 + 0.07 * stack_counter) - fig.update_yaxes(range=[y_min, y_max + extra_space]) - - # --- FINAL LAYOUT --- - dynamic_height = 600 + (stack_counter * 40) - - fig.update_layout( - xaxis_title=group_col, - yaxis_title="Productive TCRs (%)", - xaxis=dict(tickfont=dict(size=15)), - margin=dict(t=80), - showlegend=False, - width=700, - height=dynamic_height, - plot_bgcolor='rgba(0,0,0,0)' - ) - - fig.update_yaxes(showgrid=True, gridcolor='lightgrey') - fig.show() - -# --- 2. GENERATE TABS --- - -if len(unique_patients) <= 10: - print("::::: {.panel-tabset}\n") - - # --- TAB A: ALL SAMPLES --- - print("## All Samples\n") - print(":::: {.panel-tabset}\n") - for group in group_opts_all: - print(f"### {group}\n") - create_pct_prod_plot(df, group, plot_type="box") - print("\n") - print("::::\n") - - # --- TAB B: BY PATIENT --- - print("## By Patient\n") - print(":::: {.panel-tabset}\n") - - for pat in unique_patients: - print(f"### {pat}\n") - - pat_df = df[df[subject_col] == pat] - - if 'origin' in pat_df.columns and pat_df['origin'].nunique() > 1: - print(":::: {.panel-tabset}\n") # Level 4 Start (Origins) - - origins = sorted(pat_df['origin'].dropna().unique()) - - for origin in origins: - print(f"#### {origin}\n") - - origin_df = pat_df[pat_df['origin'] == origin] - - for group in group_opts_patient: - if origin_df[group].nunique() > 0: - custom_order = None - if group == timepoint_col and time_order_map: - pat_tps = origin_df[group].dropna().unique().tolist() - custom_order = sorted( - pat_tps, - key=lambda x: time_order_map.get(x, 999) - ) - - create_pct_prod_plot( - origin_df, - group, - custom_order, - plot_type="line" - ) - else: - print(f"No data for {group} in {origin}, patient {pat}.") - print("\n") - - print("::::\n") # Level 4 End - - else: - # Fallback if no origin variation - for group in group_opts_patient: - if pat_df[group].nunique() > 0: - - # --- APPLY CUSTOM SORTING IF AVAILABLE --- - custom_order = None - if group == timepoint_col and time_order_map: - pat_tps = pat_df[group].dropna().unique().tolist() - custom_order = sorted(pat_tps, key=lambda x: time_order_map.get(x, 999)) - - create_pct_prod_plot( - pat_df, - group, - custom_order, - plot_type="line" - ) - else: - print(f"No data for {group} in patient {pat}.") - print("\n") - - print("::::\n") - print(":::::\n") - -else: - # --- STANDARD VIEW (>10 Patients) --- - print("::: {.panel-tabset}\n") - for group in group_opts_all: - print(f"## {group}\n") - create_pct_prod_plot(df, group, plot_type="box") - print("\n") - print(":::\n") - -``` - -**Figure 7. Percent of productive TCRs.** Distribution of productive TCRs across sample timepoints. A productive TCR is a DNA/RNA sequence that can be translated into a protein sequence, i.e. it does not contain a premature stop codon or an out of frame rearrangement. The percent of productive TCRs is calculated as: - -$$ Percent \text{ } productive \text{ } TCRs = \frac P N $$ - -where $P$ is the number of productive TCRs and $N$ is the total number of TCRs in a given sample. - - -## Average productive CDR3 Length {#sec-cdr3-len} - -The incredible diversity is generated during the development of T-cells through a process of genetic recombination, where different gene segments (V, D, and J) are pieced together. The CDR3 region spans the junction of these segments and includes random nucleotide additions and deletions, creating a unique sequence for almost every T-cell. **The length of this CDR3 region is a direct consequence of this recombination process.** - -The CDR3 length distribution serves as a primary metric for assessing library complexity and detecting technical biases introduced during library preparation or sequencing. While biological clonality impacts this metric, deviations in non-expanded samples often signal workflow failures rather than biological reality. - -```{python} -#| output: asis -#| echo: false - -# Create a mapping for timepoint order using the 'meta' dataframe -if 'meta' in locals() and timepoint_col in meta.columns and timepoint_order_col in meta.columns: - time_order_map = meta.set_index(timepoint_col)[timepoint_order_col].to_dict() -else: - # Fallback: if 'meta' isn't loaded, try to use columns in the main df - if timepoint_order_col in df.columns: - time_order_map = df.set_index(timepoint_col)[timepoint_order_col].to_dict() - else: - time_order_map = {} # No sorting map available - -# 1. Define General Exclusions -exclude_cols = ['sample', 'file', 'clonality', 'counts', 'total_counts', - 'num_clones', 'num_TCRs', 'simpson_index', 'simpson_index_corrected', - 'num_prod', 'num_nonprod', 'pct_prod', 'pct_nonprod', - 'productive_cdr3_avg_len', 'num_convergent','ratio_convergent'] - -# 2. Define Groups for "All Samples" Tab (Exclude timepoint/order) -exclude_from_all_samples = [timepoint_col, timepoint_order_col, 'protocol_day', alias_col] - -group_opts_all = [ - col for col in df.columns - if col not in exclude_cols - and col not in exclude_from_all_samples - and df[col].nunique() < 35 -] - -# 3. Define Groups for "By Patient" Tab (Only timepoint) -group_opts_patient = [timepoint_col] if timepoint_col in df.columns else [] - -unique_patients = sorted(df[subject_col].dropna().unique().tolist()) - -# --- 1. DEFINE PLOTTING FUNCTION --- -def create_cdr3_len_plot(data, group_col, custom_order=None, plot_type="box"): - """ - custom_order: Optional list of x-axis values in the desired order. - plot_type: "box" for standard comparisons, "line" for patient timelines. - """ - if data.empty or data[group_col].dropna().empty: - return - - # Determine order: Use custom list if provided, otherwise sort alphabetically - if custom_order: - unique_cats = [x for x in custom_order if x in data[group_col].unique()] - else: - unique_cats = sorted(data[group_col].dropna().unique().tolist()) - - # --- BASE PLOT --- - if plot_type == "box": - fig = px.box( - data, - x=group_col, - y='productive_cdr3_avg_len', - color=group_col, - points='all', - hover_data=['alias'], - category_orders={group_col: unique_cats}, - template="simple_white", - title=f"CDR3 Length by {group_col}" - ) - - # --- BIGGER DOTS (size=10) --- - fig.update_traces(width=0.5, marker=dict(size=10, opacity=0.7)) - - for trace in fig.data: - if isinstance(trace, go.Box): - trace.pointpos = 0 - trace.jitter = 0.2 - - elif plot_type == "line": - # Sort for proper line connection - data = data.copy() - data[group_col] = pd.Categorical(data[group_col], categories=unique_cats, ordered=True) - data = data.sort_values(group_col) - - fig = px.scatter( - data, - x=group_col, - y='productive_cdr3_avg_len', - color_discrete_sequence=["black"], - hover_data=['alias'], - template="simple_white", - title=f"CDR3 Length by {group_col}" - ) - - # Add dotted line - fig.add_trace(go.Scatter( - x=data[group_col], - y=data['productive_cdr3_avg_len'], - mode='lines+markers', - line=dict(dash='dot', width=2), - marker=dict(size=18, color='black'), - showlegend=False - )) - - # --- STATISTICS (ONLY FOR BOXPLOTS) --- - stack_counter = 0 - - if plot_type == "box" and len(unique_cats) >= 2: - pairs = list(combinations(unique_cats, 2)) - - y_max = data['productive_cdr3_avg_len'].max() - y_min = data['productive_cdr3_avg_len'].min() - y_range = y_max - y_min - if y_range == 0: - y_range = 1 - - # Tighter spacing format - base_offset = y_range * 0.08 - step_size = y_range * 0.06 - text_offset = y_range * 0.015 - - for t1, t2 in pairs: - group1 = data[data[group_col] == t1]['productive_cdr3_avg_len'].dropna() - group2 = data[data[group_col] == t2]['productive_cdr3_avg_len'].dropna() - - if len(group1) < 2 or len(group2) < 2: - continue - - try: - stat, p_value = mannwhitneyu(group1, group2, alternative='two-sided') - except ValueError: - continue - - if p_value >= 0.05: - continue - - if p_value < 0.001: symbol = '***' - elif p_value < 0.01: symbol = '**' - else: symbol = '*' - - y_bracket = y_max + base_offset + (stack_counter * step_size) - y_text = y_bracket + text_offset - - # Draw Bracket Line - fig.add_shape( - type="line", xref="x", yref="y", - x0=t1, y0=y_bracket, x1=t2, y1=y_bracket, - line=dict(color="black", width=1.5) - ) - - # Draw Bracket Ticks - tick_len = y_range * 0.015 - fig.add_shape(type="line", xref="x", yref="y", - x0=t1, y0=y_bracket, x1=t1, y1=y_bracket - tick_len, - line=dict(color="black", width=1.5)) - fig.add_shape(type="line", xref="x", yref="y", - x0=t2, y0=y_bracket, x1=t2, y1=y_bracket - tick_len, - line=dict(color="black", width=1.5)) - - try: - x_center = (unique_cats.index(t1) + unique_cats.index(t2)) / 2 - except: - x_center = 0 - - # Add P-Value Text - fig.add_annotation( - x=x_center, y=y_text, text=f"{symbol}
p={p_value:.3f}", - showarrow=False, font=dict(size=10, color="black") - ) - - stack_counter += 1 - - extra_space = y_range * (0.15 + 0.07 * stack_counter) - fig.update_yaxes(range=[y_min, y_max + extra_space]) - - # --- FINAL LAYOUT --- - dynamic_height = 600 + (stack_counter * 40) - - fig.update_layout( - xaxis_title=group_col, - yaxis_title="Avg CDR3 Length (ntds)", - xaxis=dict(tickfont=dict(size=15)), - margin=dict(t=80), - showlegend=False, - width=700, - height=dynamic_height, - plot_bgcolor='rgba(0,0,0,0)' - ) - fig.update_yaxes(showgrid=True, gridcolor='lightgrey') - fig.show() - -# --- 2. GENERATE TABS --- - -if len(unique_patients) <= 10: - print("::::: {.panel-tabset}\n") - - # --- TAB A: ALL SAMPLES --- - print("## All Samples\n") - print(":::: {.panel-tabset}\n") - for group in group_opts_all: - print(f"### {group}\n") - create_cdr3_len_plot(df, group, plot_type="box") - print("\n") - print("::::\n") - - # --- TAB B: BY PATIENT --- - print("## By Patient\n") - print(":::: {.panel-tabset}\n") - - for pat in unique_patients: - print(f"### {pat}\n") - - pat_df = df[df[subject_col] == pat] - - # --- Origin tab layer --- - if 'origin' in pat_df.columns and pat_df['origin'].nunique() > 1: - print(":::: {.panel-tabset}\n") # Level 4 Start (Origins) - - origins = sorted(pat_df['origin'].dropna().unique()) - - for origin in origins: - print(f"#### {origin}\n") - - origin_df = pat_df[pat_df['origin'] == origin] - - for group in group_opts_patient: - if origin_df[group].nunique() > 0: - custom_order = None - if group == timepoint_col and time_order_map: - pat_tps = origin_df[group].dropna().unique().tolist() - custom_order = sorted( - pat_tps, - key=lambda x: time_order_map.get(x, 999) - ) - - create_cdr3_len_plot( - origin_df, - group, - custom_order, - plot_type="line" - ) - else: - print(f"No data for {group} in {origin}, patient {pat}.") - print("\n") - - print("::::\n") # Level 4 End - - else: - # Fallback if no origin variation - for group in group_opts_patient: - if pat_df[group].nunique() > 0: - - # --- APPLY CUSTOM SORTING IF AVAILABLE --- - custom_order = None - if group == timepoint_col and time_order_map: - pat_tps = pat_df[group].dropna().unique().tolist() - custom_order = sorted(pat_tps, key=lambda x: time_order_map.get(x, 999)) - - create_cdr3_len_plot( - pat_df, - group, - custom_order, - plot_type="line" - ) - else: - print(f"No data for {group} in patient {pat}.") - print("\n") - - print("::::\n") - print(":::::\n") - -else: - # --- STANDARD VIEW (>10 Patients) --- - print("::: {.panel-tabset}\n") - for group in group_opts_all: - print(f"## {group}\n") - create_cdr3_len_plot(df, group, plot_type="box") - print("\n") - print(":::\n") - -``` -**Figure 8. Average Productive CDR3 Length** Visualizes the frequency of CDR3 nucleotide lengths. This plot is used to verify the expected Gaussian distribution typical of polyclonal repertoires and to identify skewness indicative of technical artifacts or clonal dominance - -**Quality Control Interpretations** - -**Validation of Library Complexity (Gaussian Fit):** A high-quality, polyclonal library must exhibit a near-Gaussian (bell-shaped) distribution of CDR3 lengths. - -- Significant deviations (e.g., flat, multimodal, or ragged distributions) in control or baseline samples indicate poor library complexity, RNA degradation, or sampling bottlenecks (insufficient template starting material). - -**Distinguishing Expansion from Bias:** While clonal expansion biologically causes skew, technical artifacts can mimic this signal. - -- A sharp peak at a single length in a sample expected to be naive or healthy suggests PCR bias or amplicon contamination rather than true biological expansion. -- An abundance of unusually short (<15 ntds or <5 AA) or long (>90 ntds or >30 AA) CDR3s often indicates misalignment, primer dimers, or non-specific amplification products that survived upstream filtering. - -## Number of unique TCRs per sample (Richness) {#sec-richness} - - -```{python} -#| output: asis - - -# --- 1. Define General Exclusions --- -# Exclude continuous metrics, identifiers, and derived stats from being used as tabs -exclude_cols = [ - 'sample', 'file', 'total_counts', 'filename', 'sample_id', 'alias', - 'shannon_entropy', 'gini_coefficient', 'inverse_simpson', - 'hill_q0', 'hill_q1', 'hill_q2', 'mean_len', 'std_len', - 'clonality', 'counts', 'num_clones', 'num_TCRs', 'simpson_index', - 'simpson_index_corrected', 'num_prod', 'num_nonprod', 'pct_prod', 'pct_nonprod', - 'productive_cdr3_avg_len', 'num_convergent', 'ratio_convergent' -] - -# Add timepoint_order_col if it exists in the environment -if 'timepoint_order_col' in locals() or 'timepoint_order_col' in globals(): - exclude_cols.append(timepoint_order_col) - -# --- 2. Setup Grouping Options --- -# Automatically grab any categorical columns with reasonable cardinality (< 35 unique values) -group_opts = [ - col for col in df.columns - if col not in exclude_cols and df[col].nunique() < 35 -] - -# --- 3. DEFINE PLOTTING FUNCTION --- -def create_richness_plot(data, group_col, custom_order=None): - """ - Plots the number of TCRs (richness) as a boxplot with individual sample dots. - """ - if data.empty or data[group_col].dropna().empty: - return - - # Determine order: Use custom list if provided, otherwise sort alphabetically - if custom_order: - unique_cats = [x for x in custom_order if x in data[group_col].unique()] - else: - unique_cats = sorted(data[group_col].dropna().unique().tolist()) - - # Safely handle hover data depending on available columns - hover_cols = ['sample'] - if 'alias' in data.columns: - hover_cols.append('alias') - - # --- BASE PLOT --- - fig = px.box( - data, - x=group_col, - y='num_TCRs', - color=group_col, - points='all', - hover_data=hover_cols, - category_orders={group_col: unique_cats}, - template="simple_white", - title=f"TCR Richness by {group_col}" - ) - - # --- DOT STYLING --- - fig.update_traces(width=0.5, marker=dict(size=10, opacity=0.7)) - - # Center the dots and apply jitter so they don't overlap entirely - for trace in fig.data: - if isinstance(trace, go.Box): - trace.pointpos = 0 - trace.jitter = 0.2 - - # --- FINAL LAYOUT --- - fig.update_layout( - xaxis_title=group_col, - yaxis_title="Number of Unique TCRs", - xaxis=dict(tickfont=dict(size=15)), - margin=dict(t=60), - showlegend=False, - width=700, - height=600, - plot_bgcolor='rgba(0,0,0,0)' - ) - - fig.update_yaxes(showgrid=True, gridcolor='lightgrey') - - fig.show() - -# --- 4. EXECUTION LOOP --- -print("::: {.panel-tabset}\n\n") - -for group in group_opts: - if df[group].dropna().empty: - continue - - print(f"## {group}\n\n") - create_richness_plot(df, group) - print("\n") - -print(":::\n\n") - - -``` -**Figure 9. TCR Richness** Visualizes the number of unique CDR3 sequence. Richness is the simplest measure of diversity, counting how many different clones are present in a sample. It is highly sensitive to sequencing depth and undersampling, which can lead to an underestimation of true diversity if not properly accounted for. - -## Sequencing Depth vs Richness {#sec-depth-richness} - -**Does a sample appear "diverse" because it is biologically complex, or simply because it was sequenced more deeply than others?** -In an ideal experiment, **we would sequence until we reach saturation** (the point where sequencing more reads no longer reveals new clonotypes). On this plot (log-log scale), a saturated cohort appears as a horizontal plateau: increasing Total Counts (x-axis) does not significantly increase Richness (y-axis). This indicates that the library has captured the true biological ceiling of the sample, allowing for valid comparisons between patients. - -```{python} -#| output: asis -#| fig-width: 6 -#| fig-height: 5 - -# --- 1. Merge the missing counts data --- -total_counts = concat_df.groupby('sample')['counts'].sum().to_frame(name='total_counts') -bias_plot_df = plot_df.merge(total_counts, on='sample', how='inner') - -# --- 2. Setup Grouping Options --- -exclude_cols = ['sample', 'file', 'total_counts', 'filename', 'sample_id', - 'shannon_entropy', 'gini_coefficient', 'inverse_simpson', - 'hill_q0', 'hill_q1', 'hill_q2', 'mean_len', 'std_len', timepoint_order_col] - -group_opts = [ - col for col in bias_plot_df.columns - if col not in exclude_cols and bias_plot_df[col].nunique() < 35 -] - -# --- 3. CALCULATE UNIFIED AXIS LIMITS --- -# Find the absolute min and max across BOTH axes to ensure a perfect 1:1 square -all_vals = pd.concat([bias_plot_df['total_counts'], bias_plot_df['hill_q0']]) -global_min = all_vals[all_vals > 0].min() -global_max = all_vals.max() - -# Plotly expects log bounds in log10 format -log_min = np.log10(global_min) - 0.15 -log_max = np.log10(global_max) + 0.15 - -# --- 4. Start Quarto Tabset --- -print("::: {.panel-tabset}\n") - -for col in group_opts: - print(f"## {col}\n") - - # Calculate Correlation (Log-Log) - x_val = bias_plot_df['total_counts'] - y_val = bias_plot_df['hill_q0'] # Richness - - # Handle zeros safely for log calculation - x_log = np.log1p(x_val) - y_log = np.log1p(y_val) - - if len(x_log) > 1: - corr, _ = pearsonr(x_log, y_log) - status = "⚠️ BIAS" if corr > 0.8 else "Pass" - title_str = f'Depth Bias by {col}
R = {corr:.2f} ({status})' - else: - corr = 0 - title_str = f'Depth Bias by {col}
(Not enough data)' - - # Force Categorical Coloring - bias_plot_df[col] = bias_plot_df[col].astype(str) - - # Determine hover columns safely - hover_cols = ['sample'] - if 'alias' in bias_plot_df.columns: - hover_cols.append('alias') - - # --- PLOTLY CONVERSION --- - fig = px.scatter( - bias_plot_df, - x='total_counts', - y='hill_q0', - color=col, - symbol=col, - hover_data=hover_cols, - log_x=True, - log_y=True, - template="simple_white", - title=title_str, - labels={ - 'total_counts': 'Sequencing Depth (Total Counts)', - 'hill_q0': 'Richness (Unique Clones)' - } - ) - - # Marker styling - fig.update_traces(marker=dict(size=10, opacity=0.8, line=dict(width=1, color='DarkSlateGrey'))) - - # Add diagonal warning if correlation is high - if corr > 0.8: - fig.add_annotation( - x=0.05, - y=0.95, - xref="paper", - yref="paper", - text="Strong Depth Bias", - showarrow=False, - font=dict(color="red", size=14) - ) - - # --- ENFORCE 1:1 SCALE & CLEAN LOG GRIDS --- - fig.update_layout( - width=700, - height=700, - margin=dict(t=80), - plot_bgcolor='rgba(0,0,0,0)' - ) - - fig.update_xaxes( - showgrid=True, - gridcolor='lightgrey', - range=[log_min, log_max], - dtick=1, # Forces major ticks to exact powers of 10 - minor=dict(showgrid=False) # Hides the clustered minor grid lines - ) - - fig.update_yaxes( - showgrid=True, - gridcolor='lightgrey', - range=[log_min, log_max], - scaleanchor="x", - scaleratio=1, - dtick=1, # Forces major ticks to exact powers of 10 - minor=dict(showgrid=False) # Hides the clustered minor grid lines - ) - - fig.show() - - print("\n") - -print(":::\n") - -``` - -**Figure 10. Evaluation of Sequencing Depth Bias.** This scatter plot correlates sequencing depth (Total Counts, x-axis) with observed richness (Unique Clones, y-axis) on a log-log scale. - -A **strong positive correlation** (Pearson’s $R > 0.8$) **typically signals a technical failure in experimental design**. If your samples form a strict **diagonal line from the bottom-left to the top-right**, your **"diversity" metric is merely a proxy for library size**. In this scenario, **samples with higher read counts will artificially appear more diverse, not because of biology, but because the sampling depth was insufficient to capture the rare tail of the repertoire**. Diversity **comparisons** in such a cohort **are statistically invalid without subsampling (rarefaction)**. - -Samples with **both low sequencing depth and low richness** represent **technical dropouts**. These libraries likely suffered from **low RNA input or poor amplification efficiency**. They often lack sufficient statistical power to characterize the repertoire and **should be excluded**. - -Samples that exhibit **high sequencing depth but paradoxically low richness** are of particular interest. This discordant position—deep sampling yielding few unique clones—indicates that the repertoire is dominated by a **massive clonal expansion** (e.g., leukemia or an acute viral response) or suffers from **technical PCR "jackpotting"**. -Finding a sample with **low read counts but extreme richness** is mathematically suspicious. It implies that almost every read is a unique clonotype, a pattern often caused by **sequencing errors** (phantoms) or **contamination** with high-diversity amplicon debris. - -### Rarefaction Curve - -You **cannot fairly compare the diversity of a sample with 1 million reads to one with 10,000 reads**. The **deeper sample will always look richer simply because you looked harder**. Rarefaction simulates what your deep samples would look like if you had stopped sequencing earlier. - -```{python} - -# --------------------------------------------------------- -# Compute Basic Metrics: Depth & Richness -# --------------------------------------------------------- - -qc_stats = concat_df.groupby('sample').agg( - sequencing_depth=('counts', 'sum'), - clonotype_richness=('CDR3b', 'nunique') -).reset_index() - -# We extract just the metadata columns we need from 'df' -meta_safe = meta[['sample', 'alias', 'origin', timepoint_col]].drop_duplicates() - -# Merge metadata into qc_stats -qc_stats = qc_stats.merge(meta_safe, left_on='sample', right_on='sample', how='left') - -# --------------------------------------------------------- -# Rarefaction Curve Calculation -# --------------------------------------------------------- -# Rarefaction answers: "If we sequenced everyone to the same depth, how many clones would we find?" - - -def get_rarefaction_curve(sample_df, sample, steps=50): - """ - Simulates rarefaction by expanding counts and shuffling. - """ - # Flat array where each element is a clone ID, repeated by its count - # Use integer IDs for memory efficiency instead of strings - counts = sample_df['counts'].values - n_reads = counts.sum() - - if n_reads == 0: - return pd.DataFrame() - - # Generate clone IDs (0 to N-1) - ids = np.arange(len(counts)) - - # e.g., if clone 0 has 2 counts -> [0, 0] - # WARNING: High memory usage for very deep sequencing - flat_repertoire = np.repeat(ids, counts) - - # Shuffle to simulate random sampling - np.random.shuffle(flat_repertoire) - - # Define depths to check (from 0 to total reads) - depths = np.linspace(0, n_reads, num=steps, dtype=int) - depths = np.unique(depths)[1:] # Remove 0 and duplicates - - results = [] - for d in depths: - subsample = flat_repertoire[:d] - n_unique = len(np.unique(subsample)) - results.append({'sample': sample, 'depth': d, 'richness': n_unique}) - - return pd.DataFrame(results) - -rarefaction_data = [] - -unique_samples = concat_df['sample'].unique() -for s in unique_samples: - s_df = concat_df[concat_df['sample'] == s] - r_df = get_rarefaction_curve(s_df, s) - rarefaction_data.append(r_df) - -df_rarefaction = pd.concat(rarefaction_data) - -# Merge metadata into rarefaction data using the safe source -df_rarefaction = df_rarefaction.merge(meta_safe, on='sample', how='left') - -# ----- Rarefaction Curves - -# --------------------------------------------------------- -# Dynamic Plotting with Adaptive Axes -# --------------------------------------------------------- - -# 1. Calculate dynamic limits based on the actual data distribution -# Instead of a hard 90th percentile, we can use the max depth plus a small buffer -absolute_max_depth = df_rarefaction['depth'].max() -absolute_max_richness = df_rarefaction['richness'].max() - -# 2. Define the plot -fig_rare = px.line(df_rarefaction, - x='depth', - y='richness', - color='alias', - title='Rarefaction Curves (Diversity vs Depth)', - labels={'depth': 'Sequencing Depth (Total Counts)', 'richness': 'Identified Clonotypes'}, - hover_data=['origin'], - template='plotly_white') - -# 3. Add Vertical Line for the "Saturation Point" (Minimum Sample Depth) -min_depth = qc_stats['sequencing_depth'].min() -fig_rare.add_vline(x=min_depth, line_dash="dash", line_color="gray", - annotation_text=f"Min Depth: {min_depth}", - annotation_position="top left") - -# 4. Update Axes Dynamically -# Setting range to None (default) allows Plotly to auto-scale, -# but if you want to cap it to the 95th percentile to avoid outliers stretching the plot: -suggested_max_x = qc_stats['sequencing_depth'].quantile(0.95) - -fig_rare.update_layout( - xaxis=dict( - # Use 'range' if you want a specific cutoff, otherwise omit for auto-scaling - range=[0, suggested_max_x * 1.05], # 5% buffer - autorange=False, - rangemode="tozero" - ), - yaxis=dict( - autorange=True, # Automatically adjust Y based on the visible X range - rangemode="tozero" - ) -) - -fig_rare.show() - -``` -**Figure 11. Rarefaction Analysis of TCR Repertoire Diversity.** Rarefaction curves depict the accumulation of unique clonotypes (richness) as a function of sequencing depth (Total Counts) for each sample. The dotted vertical line marks the depth of the smallest library. - -Follow each sample's line from left to right to track how new clonotypes are discovered as sequencing depth increases. In an **ideal scenario**, you want to see a **saturation Plateau**: the curve shoots up initially and then flattens out into a horizontal line. This indicates success; you have sequenced the sample deeply enough to exhaust the available diversity, and spending more money on reads would yield no new information. - -You will often see ***two distinct deviations from this ideal***: -- **Lazy Climber:** A curve that **rises very slowly**. This signals **clonal dominance**; the **sequencer is reading millions of molecules, but they are all the same few sequences**. This is a biological reality, not a technical failure, so these samples should be kept. -- **Unsatisfied Climber:** A curve that **shoots diagonally upward** even at maximum depth, never bending toward horizontal. This indicates **undersampling**; the library is far more diverse than your sequencing budget allowed. Because the "Richness" value here is just a lower bound of an unknown total, these samples often cannot be accurately compared and may need to be discarded. - -The **vertical dashed line** on your plot marks the **Lowest Sequencing Depth in your entire cohort**. To perform a valid statistical comparison, you must mathematically "downsample" every patient to this threshold. This forces a difficult trade-off: if your lowest sample has an unusable depth (e.g., < 1,000 reads), keeping it would require you to throw away 90% of the data from your high-quality samples just to make them comparable. In such cases, it is better to **delete the single bad sample rather than degrade the resolution of the entire cohort**. - -**Re-sequencing?** -This visualization also serves as the financial evidence for deciding whether to re-sequence a dropout sample or not: -- If the curve is **rising steeply**, the sample was **undersampled**. **Re-sequencing is possible** and will yield missing data. -- If the **curve has plateaued**, the **sample is naturally low-diversity** (or the library prep failed). **Re-sequencing would be a waste of money**, as deeper reads will not recover new clonotypes. - -**Expected Richness** -Always verify the magnitude of your Y-axis ("Identified Clonotypes"). For a typical bulk human blood sample, you should expect to see between $10^3$ and $10^5$ unique clonotypes. If your plot shows a maximum of only ~90 clones (as seen in some synthetic or filtered datasets), this is a major red flag for real biological data. Such a low count in primary tissue usually implies a library prep failure, an oligoclonal cell line, or an artifact of synthetic simulation. - - -::: {.callout-caution title="Warning"} -Not all diversity metrics are equally sensitive to library size: - -**High Sensitivity (Depth Matters!) $\rightarrow$ Use Rarefied Data** -Metrics that rely on counting the **number of unique clones** are heavily biased by sequencing depth. -* **Richness ($q=0$):** A direct count of unique clonotypes. This is the most sensitive metric. -* **Shannon Entropy ($q=1$):** While it considers frequency, it is still significantly influenced by the "long tail" of rare clones. -* **Clonality:** Because it is calculated as $1 - \text{Normalized Entropy}$ (which relies on Richness), it inherits this bias. A deep sample will often appear "less clonal" than a shallow one simply due to the math. - -**Low Sensitivity (Depth is Irrelevant) $\rightarrow$ Use Full Data** -Metrics that measure **dominance** or **inequality** focus on the most abundant clones. Because rare clones (singletons) contribute almost zero to these calculations, adding more reads does not significantly change the result. -- **Simpson Index ($q=2$):** Measures the probability of picking two identical clones. It is mathematically weighted toward the top clones, making it extremely robust to depth differences. -- **Gini Coefficient:** Measures the inequality of the distribution (Lorenz curve). It stabilizes quickly and does not require downsampling. -- **Gene Usage (V/J):** Based on percentages (frequency), not raw counts. -::: - -## V-Gene Usage PCA {#sec-pca} - -```{python} -#| output: asis -#| fig-width: 6 -#| fig-height: 5 - -# --- 1. PCA Calculation Function (Unchanged) --- -def analyze_v_usage(df): - """ - Expects df with columns: ['sample', 'TRBV', 'counts'] - """ - matrix = df.groupby(['sample', 'TRBV'])['counts'].sum().unstack(fill_value=0) - freq_matrix = matrix.div(matrix.sum(axis=1), axis=0) - - scaler = StandardScaler() - scaled_data = scaler.fit_transform(freq_matrix) - - pca = PCA(n_components=2) - components = pca.fit_transform(scaled_data) - - pca_df = pd.DataFrame(data=components, columns=['PC1', 'PC2'], index=freq_matrix.index) - return freq_matrix, pca_df - -# --- 2. Run Analysis & Prepare Data --- -_, pca_df = analyze_v_usage(concat_df) - -if pca_df.index.name == 'sample' or 'sample' not in pca_df.columns: - pca_df = pca_df.reset_index() - if 'index' in pca_df.columns and 'sample' not in pca_df.columns: - pca_df = pca_df.rename(columns={'index': 'sample'}) - -plot_data = pd.merge(pca_df, meta, on='sample', how='inner') - -# --- 3. Determine Grouping Columns --- -exclude_cols = ['sample', 'file', 'total_counts', 'filename', 'sample_id'] -group_opts = [ - col for col in meta.columns - if col not in exclude_cols and meta[col].nunique() < 35 -] - -# --- 4. Generate Tabs & Plotly Charts --- -print("::: {.panel-tabset}\n") - -for col in group_opts: - print(f"## {col}\n") - - # Force categorical coloring - plot_data[col] = plot_data[col].astype(str) - - # Safely determine hover columns - hover_cols = ['sample'] - if 'alias' in plot_data.columns: - hover_cols.append('alias') - - # --- PLOTLY SCATTER PLOT --- - fig = px.scatter( - plot_data, - x='PC1', - y='PC2', - color=col, - symbol=col, - hover_data=hover_cols, # <-- Adds the hover over info - template='simple_white', - title=f'TRBV Gene Usage PCA (Colored by {col})' - ) - - # Match the Seaborn styling (large dots, transparency, black outlines) - fig.update_traces(marker=dict(size=12, opacity=0.8, line=dict(width=1, color='DarkSlateGrey'))) - - # --- OUTLIER LABELING LOGIC --- - outlier_mask = (np.abs(plot_data['PC1']) > 3 * plot_data['PC1'].std()) | \ - (np.abs(plot_data['PC2']) > 3 * plot_data['PC2'].std()) - - outliers = plot_data[outlier_mask] - - # Add permanent text labels next to outlier points - for _, row in outliers.iterrows(): - fig.add_annotation( - x=row['PC1'], - y=row['PC2'], - text=f"{row['sample']}", - showarrow=False, - xshift=15, # Shifts text to the right - yshift=15, # Shifts text up - font=dict(size=10, color="black") - ) - - # --- FINAL LAYOUT --- - fig.update_layout( - xaxis_title="PC1", - yaxis_title="PC2", - width=700, - height=600, - margin=dict(t=60), - plot_bgcolor='rgba(0,0,0,0)' - ) - - # Add gridlines and emphasize the 0,0 origin axes - fig.update_xaxes(showgrid=True, gridcolor='lightgrey', zeroline=True, zerolinecolor='gray', zerolinewidth=1.5) - fig.update_yaxes(showgrid=True, gridcolor='lightgrey', zeroline=True, zerolinecolor='gray', zerolinewidth=1.5) - - fig.show() - - print("\n") - -print(":::\n") - -``` -**Figure 12. Dimensionality Reduction of TRBV Gene Usage.** Principal Component Analysis (PCA) plot visualizing the variation in T-cell Receptor Beta Variable (TRBV) gene usage frequencies across samples, colored by metadata variables. - -This plot compresses the high-dimensional complexity of T-cell receptor usage (e.g., TRBV1, TRBV2... TRBV30) into a 2D map. Each dot represents a sample, and the **distance between dots represents the similarity of their "V-gene fingerprints."** **If two samples are close together, they are using the V-genes in almost identical proportions**. If they are far apart, their repertoire structures are radically different. - -**Technical Failures 🚩** -One **purpose** of this visualization is to **detect technical artifacts**. When you toggle the tabs to color by **Sequencing Batch for same condition**, you should ideally see a randomly **mixed cloud**. This indicates that the V-gene usage is **consistent** across your **experiments**. However, **if you see distinct, non-overlapping islands** corresponding to different batches (e.g., all "Batch A" samples cluster on the left, and "Batch B" samples on the right), **you have a severe Batch Effect**. This usually implies a change in the wet-lab protocol, such as a new set of multiplex primers that amplify genes with different efficiencies. If observed, **you cannot compare samples across these batches without statistical correction** (e.g., ComBat). - -**Biological Interpretation 🦠** -In the **context of cancer treatment** (especially Checkpoint Inhibitors like PD-1 blockade), **we are looking for Remodeling**. We want to **know if the drug "woke up" the immune system and caused it to expand new T-cell armies (clones) to fight the tumor. For example, if "Responders" form a distinct island separate from "Non-Responders", this suggests that having a specific bias toward a certain V-gene, TRBV-19 for example, pre-disposes a patient to fight the cancer effectively. - -An outlier on this plot can have a **specific biological cause**: - -- ***Monoclonal Expansion***. If a patient has a **massive expansion of a single clone** (e.g., a clone using TRBV7-2 takes up 50% of the repertoire), that sample will be pulled violently toward the "TRBV7-2 direction" in PCA space. Finding such an **outlier in a "healthy" control often indicates a PCR jackpotting error or contamination**. - -## Outlier Detection {#sec-outlier-flag} - -This aggregates everything. We calculate Z-scores for every metric. If a sample is $>3$ standard deviations away from the mean, we flag it. -```{python} - -# CDR3 Length Distribution & Distances -def analyze_cdr3_lengths(df): - """ - Expects df with columns: ['sample', 'CDR3b'] - """ - # Calculate lengths - df['cdr3_len'] = df['CDR3b'].apply(len) - - # 1. Basic Stats per sample - stats = df.groupby('sample')['cdr3_len'].agg( - mean_len='mean', - std_len='std', - skew_len=lambda x: skew(x) - ) - - # 2. Cohort "Ideal" Distribution - # We aggregate ALL samples to create a reference distribution - cohort_distribution = df['cdr3_len'].values - - # 3. Distance Calculation - # How far is each sample from the cohort average? - distances = [] - for sample in stats.index: - sample_dist = df[df['sample'] == sample]['cdr3_len'].values - - # Earth Mover's Distance (Wasserstein) - # Robust metric for histogram similarity - wd = wasserstein_distance(sample_dist, cohort_distribution) - distances.append(wd) - - stats['cdr3_len_dist_to_cohort'] = distances - return stats - -def run_qc_pipeline(df): - div_df = calculate_diversity_metrics(df) - len_df = analyze_cdr3_lengths(df) - _, pca_df = analyze_v_usage(df) - - # Merge all metrics into one master QC table - qc_table = div_df.join(len_df).join(pca_df) - - # Outlier Detection Logic (Z-Scores) - # We define which columns we care about for outliers - metrics_to_check = ['shannon_entropy', 'inverse_simpson', 'mean_len', 'cdr3_len_dist_to_cohort', 'PC1', 'PC2'] - - flags = [] - - for metric in metrics_to_check: - mu = qc_table[metric].mean() - sigma = qc_table[metric].std() - - # Calculate Z-score - z_score = (qc_table[metric] - mu) / sigma - - # Find outliers (|z| > 3 is standard, maybe use 2.5 for stricter QC) - outliers = qc_table[np.abs(z_score) > 3].index.tolist() - - for out in outliers: - flags.append({ - 'sample_id': out, - 'metric': metric, - 'value': qc_table.loc[out, metric], - 'z_score': z_score[out], - 'flag': f"Extreme {metric}" - }) - - return qc_table, pd.DataFrame(flags) - -qc_table,flags = run_qc_pipeline(concat_df) -print("\nMetrics values for all samples\n") -print(qc_table) -print("\nThe following samples had a metric that deviates from the rest (|z-score|>3):\n") -print(flags) -``` From 749a0732befb55e28ea6bc150660d8150f94e3d5 Mon Sep 17 00:00:00 2001 From: dimalvovs Date: Mon, 22 Jun 2026 14:08:46 -0400 Subject: [PATCH 06/21] fix quarto jupyter on macs --- Dockerfile | 29 +++++++++++++++++++---------- conf/base.config | 2 +- env.yml | 8 ++++---- nextflow.config | 3 +++ 4 files changed, 27 insertions(+), 15 deletions(-) diff --git a/Dockerfile b/Dockerfile index 50f3cb3..c8c5a0b 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,4 +1,4 @@ -FROM mambaorg/micromamba:1.5.8 +FROM mambaorg/micromamba:2 # Ensure we run as root for apt USER root @@ -50,14 +50,22 @@ RUN git init /opt/GIANA && \ chmod +x /opt/GIANA/GIANA4.1.py && \ ln -s /opt/GIANA/GIANA4.1.py /usr/local/bin/GIANA -# Install quarto -RUN mkdir -p /opt/quarto/1.6.42 \ - && curl -o /tmp/quarto.tar.gz -L \ - "https://github.com/quarto-dev/quarto-cli/releases/download/v1.6.42/quarto-1.6.42-linux-amd64.tar.gz" \ - && tar -zxvf /tmp/quarto.tar.gz \ - -C "/opt/quarto/1.6.42" \ - --strip-components=1 \ - && rm /tmp/quarto.tar.gz +# Install quarto (match the container architecture) +RUN set -eux; \ + arch="$(dpkg --print-architecture)"; \ + case "$arch" in \ + amd64) quarto_arch="amd64" ;; \ + arm64) quarto_arch="arm64" ;; \ + *) echo "Unsupported architecture for Quarto: $arch"; exit 1 ;; \ + esac; \ + mkdir -p /opt/quarto/1.9.38; \ + curl -L -o /tmp/quarto.tar.gz \ + "https://github.com/quarto-dev/quarto-cli/releases/download/v1.9.38/quarto-1.9.38-linux-${quarto_arch}.tar.gz"; \ + tar -zxf /tmp/quarto.tar.gz \ + -C "/opt/quarto/1.9.38" \ + --strip-components=1; \ + test -x /opt/quarto/1.9.38/bin/quarto; \ + rm /tmp/quarto.tar.gz # Install VDJmatch and symlink RUN mkdir -p /opt/vdjmatch/1.3.1 \ @@ -68,7 +76,8 @@ RUN mkdir -p /opt/vdjmatch/1.3.1 \ && ln -s /opt/vdjmatch/1.3.1/vdjmatch-1.3.1/vdjmatch-1.3.1.jar /usr/local/bin/vdjmatch.jar # Add to PATH -ENV PATH="/opt/quarto/1.6.42/bin:${PATH}" +ENV PATH="/opt/quarto/1.9.38/bin:${PATH}" # Add LD_LIBRARY_PATH for pandas ENV LD_LIBRARY_PATH=/opt/conda/lib + diff --git a/conf/base.config b/conf/base.config index 0b8e8b3..6e5785a 100644 --- a/conf/base.config +++ b/conf/base.config @@ -9,7 +9,7 @@ */ process { - container = "ghcr.io/karchinlab/tcrtoolkit:main" + container = "${params.container}" publishDir = [ path: { "${params.outdir}/${task.process.tokenize(':')[-1].tokenize('_')[0].toLowerCase()}" }, diff --git a/env.yml b/env.yml index 331ef45..2e41aa9 100644 --- a/env.yml +++ b/env.yml @@ -11,9 +11,7 @@ dependencies: - seaborn=0.13.0 - dash=2.14.1 - matplotlib=3.8.1 - - pip=23.2.1 - - jupyterlab=4.0.8 - - notebook=7.0.6 + - pip=26.1.2 - fsspec=2024.3.1 - s3fs=2024.3.1 - python-igraph=0.11.8 @@ -26,6 +24,8 @@ dependencies: - rpy2=3.6.4 - unzip - openjdk=8 + - parasail-python + - jupyter # R and R packages - r-base=4.4.2 @@ -38,8 +38,8 @@ dependencies: # Pip packages - pip: + - papermill - csvkit==1.3.0 - - papermill==2.4.0 - plotly==5.18.0 - abydos==0.5.0 - fastcluster==1.2.6 diff --git a/nextflow.config b/nextflow.config index dd2aa28..09e2d3d 100644 --- a/nextflow.config +++ b/nextflow.config @@ -63,6 +63,9 @@ params { matrix_sparsity = "sparse" distance_metric = "tcrdist" db_path = "${projectDir}/assets/tcrdist3_files/alphabeta_gammadelta_db.tsv" + + //container to use for running the workflow + container = "ghcr.io/karchinlab/tcrtoolkit:main" } includeConfig 'conf/base.config' From f4b5b62399a21ebcaa7402ad37e6117484fc4fb2 Mon Sep 17 00:00:00 2001 From: dimalvovs Date: Mon, 22 Jun 2026 14:11:02 -0400 Subject: [PATCH 07/21] pick container for tests in ci --- .github/workflows/nextflow-test.yaml | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/.github/workflows/nextflow-test.yaml b/.github/workflows/nextflow-test.yaml index 46adbea..fcc538b 100644 --- a/.github/workflows/nextflow-test.yaml +++ b/.github/workflows/nextflow-test.yaml @@ -7,6 +7,9 @@ on: - 'main' workflow_dispatch: +env: + DEFAULT_CONTAINER: ${{ vars.CONTAINER || 'ghcr.io/karchinlab/tcrtoolkit:main' }} + jobs: unit-tests: name: Unit Tests @@ -17,6 +20,22 @@ jobs: with: submodules: recursive + - &resolve_test_container + name: Resolve test container + shell: bash + env: + IMAGE_BASE: ghcr.io/${{ github.repository }} + SOURCE_SHA: ${{ github.event.pull_request.head.sha || github.sha }} + run: | + container="$DEFAULT_CONTAINER" + if [[ "$container" == 'ghcr.io/karchinlab/tcrtoolkit:main' ]]; then + candidate="$IMAGE_BASE:${SOURCE_SHA::7}" + if docker manifest inspect "$candidate" >/dev/null 2>&1; then + container="$candidate" + fi + fi + echo "CONTAINER=$container" >> "$GITHUB_ENV" + - name: Setup Java uses: actions/setup-java@v2 with: @@ -62,6 +81,8 @@ jobs: with: submodules: recursive + - *resolve_test_container + - name: Setup Java uses: actions/setup-java@v2 with: From 41b5fb9bccd8155f08662c822ea6aaba64e89da4 Mon Sep 17 00:00:00 2001 From: dimalvovs Date: Mon, 22 Jun 2026 14:16:42 -0400 Subject: [PATCH 08/21] allow custom container in test --- tests/fixtures/test.qmd | 2 +- tests/main.nf.test | 3 ++- tests/modules/local/report/render_notebook.nf.test | 1 + 3 files changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/fixtures/test.qmd b/tests/fixtures/test.qmd index 2ab696a..1dee0ee 100644 --- a/tests/fixtures/test.qmd +++ b/tests/fixtures/test.qmd @@ -1,7 +1,7 @@ --- title: "test notebook" format: html -engine: knitr +jupyter: python3 --- ```{python} diff --git a/tests/main.nf.test b/tests/main.nf.test index 3b7fbe4..92bb2b9 100644 --- a/tests/main.nf.test +++ b/tests/main.nf.test @@ -5,7 +5,7 @@ nextflow_pipeline { test("Minimal example dataset") { - tag "basic-sd" + tag "minimal-example" when { params { @@ -14,6 +14,7 @@ nextflow_pipeline { input_format = "adaptive" max_memory = "8.GB" max_cpus = 4 + container = System.getenv("CONTAINER") ?: "ghcr.io/karchinlab/tcrtoolkit:main" } } diff --git a/tests/modules/local/report/render_notebook.nf.test b/tests/modules/local/report/render_notebook.nf.test index 6497f1d..4dd7980 100644 --- a/tests/modules/local/report/render_notebook.nf.test +++ b/tests/modules/local/report/render_notebook.nf.test @@ -14,6 +14,7 @@ nextflow_process { when { params { samplesheet = "${projectDir}/tests/fixtures/valid_samplesheet.csv" + container = System.getenv("CONTAINER") ?: "ghcr.io/karchinlab/tcrtoolkit:main" } process { From cc52095b75a1eb1348ff63dde4a764e6d78e71bd Mon Sep 17 00:00:00 2001 From: dimalvovs Date: Thu, 25 Jun 2026 09:03:11 -0400 Subject: [PATCH 09/21] update docker instructions --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index d6dbc44..5aafd0a 100644 --- a/README.md +++ b/README.md @@ -37,3 +37,5 @@ With the update to [Nextflow strict syntax](https://docs.seqera.io/nextflow/stri nextflow run KarchinLab/TCRtoolkit \ -params-file params.yml ``` +[!IMPORTANT] +If having an error similar to `qemu: uncaught target signal 11 (Segmentation fault) - core dumped` (observed on Mac OS Tahoe i.e. `26.5.1`), rebuild the Docker container locally, using the provided Dockerfile, and use that in the params, e.g. `--container `. It features arch selectors to pick the right `quarto` installation. \ No newline at end of file From 0819c7c2f30d2192ba598f2b8c2913fdfd385135 Mon Sep 17 00:00:00 2001 From: dimalvovs Date: Fri, 26 Jun 2026 09:17:53 -0400 Subject: [PATCH 10/21] init qc report --- conf/modules.config | 8 ++++++++ modules/local/report/render_notebook.nf | 14 +++++--------- nextflow.config | 3 +++ nextflow_schema.json | 4 ++++ params.yml | 5 +++-- subworkflows/local/sample.nf | 2 ++ workflows/tcrtoolkit.nf | 24 ++++++++++++++++++++++++ 7 files changed, 49 insertions(+), 11 deletions(-) diff --git a/conf/modules.config b/conf/modules.config index c697775..f664c9b 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -39,4 +39,12 @@ process { label = params.matrix_sparsity == 'sparse' ? 'process_medium' : ['process_high', 'process_high_memory'] } + withName: RENDER_NOTEBOOK { + publishDir = [ + path: { "${params.outdir}/report" }, + mode: params.publish_dir_mode, + saveAs: { filename -> filename.equals('versions.yml') ? null : filename } + ] + } + } \ No newline at end of file diff --git a/modules/local/report/render_notebook.nf b/modules/local/report/render_notebook.nf index 5b98d4d..5b15427 100644 --- a/modules/local/report/render_notebook.nf +++ b/modules/local/report/render_notebook.nf @@ -1,32 +1,28 @@ // Generic process to render a Quarto notebook to HTML process RENDER_NOTEBOOK { - tag "${report_name}" + tag "${notebook.getBaseName()}" label 'process_single' input: - tuple val(report_name), path(notebook), val(data_dir) + tuple path(notebook), path(files) // path(files) just stages files in root dir val project_name val workflow_cmd output: - path "${report_name}.html" + path "${notebook.getBaseName()}.html" script: """ - ## copy quarto notebook to working directory - cp $notebook ${report_name}.qmd - ## render qmd report to html - quarto render ${report_name}.qmd \\ + quarto render ${notebook} \\ -P project_name:${project_name} \\ -P workflow_cmd:'${workflow_cmd}' \\ - -P project_dir:${data_dir} \\ -P sample_table:${file(params.samplesheet)} \\ --to html """ stub: """ - touch ${report_name}.html + touch ${notebook.getBaseName()}.html """ } diff --git a/nextflow.config b/nextflow.config index 09e2d3d..226cdd2 100644 --- a/nextflow.config +++ b/nextflow.config @@ -66,6 +66,9 @@ params { //container to use for running the workflow container = "ghcr.io/karchinlab/tcrtoolkit:main" + + //reports + template_qc = "${projectDir}/notebooks/template_qc.qmd" } includeConfig 'conf/base.config' diff --git a/nextflow_schema.json b/nextflow_schema.json index ac90f76..48aae9e 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -111,6 +111,10 @@ "compare_stats_template": { "type": "string", "description": "Path to compare notebook template." + }, + "template_qc": { + "type": "string", + "description": "Path to QC notebook template." } } }, diff --git a/params.yml b/params.yml index af15fd8..6d7eec4 100644 --- a/params.yml +++ b/params.yml @@ -1,6 +1,7 @@ -samplesheet: "https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/samplesheet.csv" +samplesheet: "tests/test_data/minimal-example/samplesheet.csv" outdir: "out-minimal-dev" -workflow_level: "sample,patient,compare" +#workflow_level: "sample,compare,patient" - uncomment once patient is fixed +workflow_level: "sample,compare" input_format: "adaptive" max_memory: 10.GB max_cpus: 4 \ No newline at end of file diff --git a/subworkflows/local/sample.nf b/subworkflows/local/sample.nf index 52923bf..28f7b67 100644 --- a/subworkflows/local/sample.nf +++ b/subworkflows/local/sample.nf @@ -110,4 +110,6 @@ workflow SAMPLE { VDJDB_VDJMATCH (processed_samples, VDJDB_GET.out.ref_db) + emit: + sample_csv = SAMPLE_CALC.out.sample_csv } \ No newline at end of file diff --git a/workflows/tcrtoolkit.nf b/workflows/tcrtoolkit.nf index 7237e29..4406ccc 100644 --- a/workflows/tcrtoolkit.nf +++ b/workflows/tcrtoolkit.nf @@ -16,6 +16,7 @@ include { PATIENT } from '../subworkflows/local/patient' include { COMPARE } from '../subworkflows/local/compare' include { VALIDATE_PARAMS } from '../subworkflows/local/validate_params' include { ANNOTATE } from '../subworkflows/local/annotate' +include { REPORT } from '../subworkflows/local/report' include { PSEUDOBULK_PHENOTYPE }from '../subworkflows/local/pseudobulk_phenotype' @@ -105,6 +106,29 @@ workflow TCRTOOLKIT { ANNOTATE.out.cdr3_pgen ) } + + // Notebooks - works on channel tuples [report name, [report files]] + + // QC report requires sample-level aggregate outputs. + if (levels.contains('sample')) { + def sample_stats_agg = SAMPLE.out.sample_csv + .collectFile(name: "sample_stats.csv", keepHeader: true, skip: 1, sort: true) + + def ch_reports = sample_stats_agg + .combine(ANNOTATE.out.concat_cdr3_sorted) + .map { sample_stats_csv, concat_cdr3_sorted -> + tuple( + file(params.template_qc), + [sample_stats_csv, + concat_cdr3_sorted] + ) + } + + REPORT( ch_reports ) + } else { + println("\u001B[33m[WARN]\u001B[0m QC report skipped because workflow_level does not include 'sample'.") + } + } /* From 83be6078fbb54d25f22ab4a51bcb63a041ae913a Mon Sep 17 00:00:00 2001 From: dimalvovs Date: Fri, 26 Jun 2026 09:24:12 -0400 Subject: [PATCH 11/21] fix test case --- tests/modules/local/report/render_notebook.nf.test | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/modules/local/report/render_notebook.nf.test b/tests/modules/local/report/render_notebook.nf.test index 4dd7980..9004825 100644 --- a/tests/modules/local/report/render_notebook.nf.test +++ b/tests/modules/local/report/render_notebook.nf.test @@ -20,7 +20,6 @@ nextflow_process { process { """ input[0] = [ - "test_notebook", file("${projectDir}/tests/fixtures/test.qmd"), "${projectDir}/tests/fixtures" ] From c16e04ed630e06db39eb5992fa74fee06dea9abd Mon Sep 17 00:00:00 2001 From: dimalvovs Date: Fri, 26 Jun 2026 09:30:17 -0400 Subject: [PATCH 12/21] use inputs from current dir --- notebooks/template_qc.qmd | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/notebooks/template_qc.qmd b/notebooks/template_qc.qmd index 7928457..85d5241 100644 --- a/notebooks/template_qc.qmd +++ b/notebooks/template_qc.qmd @@ -48,7 +48,6 @@ Thank you for using TCRtoolkit! This report is generated from the data provided. # --------------------------------------------------------- workflow_cmd = '' project_name='' -project_dir='' sample_table='' timepoint_col = 'timepoint' @@ -66,9 +65,8 @@ subject_col = 'subject_id' # --------------------------------------------------------- # Define files -project_dir=f"{project_dir}/{project_name}" -sample_stats_csv = f"{project_dir}/sample/sample_stats.csv" -concat_file = f"{project_dir}/annotate/concatenated_cdr3_sorted.tsv" +sample_stats_csv = "sample_stats.csv" +concat_file = "concatenated_cdr3_sorted.tsv" ``` From 4fc6e9029e2a20ddacb53e0b441497e0e86450ca Mon Sep 17 00:00:00 2001 From: dimalvovs Date: Fri, 26 Jun 2026 09:31:32 -0400 Subject: [PATCH 13/21] add missing subworkflow --- subworkflows/local/report.nf | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 subworkflows/local/report.nf diff --git a/subworkflows/local/report.nf b/subworkflows/local/report.nf new file mode 100644 index 0000000..2bac961 --- /dev/null +++ b/subworkflows/local/report.nf @@ -0,0 +1,34 @@ +/* +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + IMPORT LOCAL MODULES +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +*/ + +include { RENDER_NOTEBOOK } from '../../modules/local/report/render_notebook' + +/* +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + RUN MAIN SUBWORKFLOW +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +*/ + +// Single-channel listener: receives tuples of +// (report_name, notebook_path, sample_stats_csv, concat_cdr3_sorted, data_dir) +// and renders each notebook to HTML. +workflow REPORT { + + take: + ch_reports // channel of tuples: tuple(val report_name, path notebook, path sample_stats_csv, path concat_cdr3_sorted, val data_dir) + + main: + RENDER_NOTEBOOK( + ch_reports, + params.project_name, + workflow.commandLine + ) + + emit: + + RENDER_NOTEBOOK.out + +} From d2ffbf984ed48f22c86b7f10366d83520252fa2b Mon Sep 17 00:00:00 2001 From: dimalvovs Date: Fri, 26 Jun 2026 12:15:09 -0400 Subject: [PATCH 14/21] extend to multiple reports --- workflows/tcrtoolkit.nf | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/workflows/tcrtoolkit.nf b/workflows/tcrtoolkit.nf index 4406ccc..f962d0b 100644 --- a/workflows/tcrtoolkit.nf +++ b/workflows/tcrtoolkit.nf @@ -107,27 +107,27 @@ workflow TCRTOOLKIT { ) } - // Notebooks - works on channel tuples [report name, [report files]] + // Report - works on channel of tuples [report name, [report files]] + ch_reports = channel.empty() // QC report requires sample-level aggregate outputs. - if (levels.contains('sample')) { - def sample_stats_agg = SAMPLE.out.sample_csv - .collectFile(name: "sample_stats.csv", keepHeader: true, skip: 1, sort: true) - - def ch_reports = sample_stats_agg - .combine(ANNOTATE.out.concat_cdr3_sorted) - .map { sample_stats_csv, concat_cdr3_sorted -> - tuple( - file(params.template_qc), - [sample_stats_csv, - concat_cdr3_sorted] - ) - } - - REPORT( ch_reports ) - } else { - println("\u001B[33m[WARN]\u001B[0m QC report skipped because workflow_level does not include 'sample'.") - } + ch_qc_report = SAMPLE.out.sample_csv + .collectFile(name: "sample_stats.csv", keepHeader: true, skip: 1, sort: true) + .combine(ANNOTATE.out.concat_cdr3_sorted) + .map { sample_stats_csv, concat_cdr3_sorted -> + tuple( + file(params.template_qc), + [sample_stats_csv, + concat_cdr3_sorted] + ) + } + ch_reports = ch_reports.mix(ch_qc_report) + + // Another report + // ch_reports = ch_reports.mix( ... ) + + REPORT( ch_reports ) + } From 4e34a1a159080dd4735b1d908fef3c8adf8656ee Mon Sep 17 00:00:00 2001 From: KevinMLanderos Date: Thu, 2 Jul 2026 17:55:58 -0400 Subject: [PATCH 15/21] Changed column names associated with file 'concatenated_cdr3_sorted' --- notebooks/template_qc.qmd | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/notebooks/template_qc.qmd b/notebooks/template_qc.qmd index 85d5241..5ac1f5e 100644 --- a/notebooks/template_qc.qmd +++ b/notebooks/template_qc.qmd @@ -187,8 +187,7 @@ For $q=0$, this simplifies to $R$ (Total unique clonotypes). For $q=1$, the limi # --- 1. CALCULATION FUNCTION --- def calculate_diversity_metrics(df): results = [] - # Check if 'counts' exists, if not try 'clone_count' or similar - count_col = 'counts' if 'counts' in df.columns else df.columns[0] # Fallback + count_col = 'duplicate_count' if 'duplicate_count' in df.columns else df.columns[0] # Fallback grouped = df.groupby('sample') for sample, data in grouped: @@ -1016,14 +1015,14 @@ Where $|A \cap B|$ is the count of unique clonotypes shared by both samples, and # --- 1. CALCULATION FUNCTION (Only Jaccard) --- def calculate_overlaps(df): samples = sorted(df['sample'].unique()) - clones = sorted(df['CDR3b'].unique()) + clones = sorted(df['junction_aa'].unique()) sample_map = {s: i for i, s in enumerate(samples)} clone_map = {c: i for i, c in enumerate(clones)} row_idx = df['sample'].map(sample_map).values - col_idx = df['CDR3b'].map(clone_map).values - values = df['counts'].values + col_idx = df['junction_aa'].map(clone_map).values + values = df['duplicate_count'].values n_samples = len(samples) n_clones = len(clones) @@ -1836,7 +1835,7 @@ In an ideal experiment, **we would sequence until we reach saturation** (the poi #| fig-height: 5 # --- 1. Merge the missing counts data --- -total_counts = concat_df.groupby('sample')['counts'].sum().to_frame(name='total_counts') +total_counts = concat_df.groupby('sample')['duplicate_count'].sum().to_frame(name='total_counts') bias_plot_df = plot_df.merge(total_counts, on='sample', how='inner') # --- 2. Setup Grouping Options --- @@ -1976,8 +1975,8 @@ You **cannot fairly compare the diversity of a sample with 1 million reads to on # --------------------------------------------------------- qc_stats = concat_df.groupby('sample').agg( - sequencing_depth=('counts', 'sum'), - clonotype_richness=('CDR3b', 'nunique') + sequencing_depth=('duplicate_count', 'sum'), + clonotype_richness=('junction_aa', 'nunique') ).reset_index() # We extract just the metadata columns we need from 'df' @@ -1998,7 +1997,7 @@ def get_rarefaction_curve(sample_df, sample, steps=50): """ # Flat array where each element is a clone ID, repeated by its count # Use integer IDs for memory efficiency instead of strings - counts = sample_df['counts'].values + counts = sample_df['duplicate_count'].values n_reads = counts.sum() if n_reads == 0: @@ -2132,9 +2131,9 @@ Metrics that measure **dominance** or **inequality** focus on the most abundant # --- 1. PCA Calculation Function (Unchanged) --- def analyze_v_usage(df): """ - Expects df with columns: ['sample', 'TRBV', 'counts'] + Expects df with columns: ['sample', 'v_call', 'duplicate_count'] """ - matrix = df.groupby(['sample', 'TRBV'])['counts'].sum().unstack(fill_value=0) + matrix = df.groupby(['sample', 'v_call'])['duplicate_count'].sum().unstack(fill_value=0) freq_matrix = matrix.div(matrix.sum(axis=1), axis=0) scaler = StandardScaler() @@ -2253,10 +2252,10 @@ This aggregates everything. We calculate Z-scores for every metric. If a sample # CDR3 Length Distribution & Distances def analyze_cdr3_lengths(df): """ - Expects df with columns: ['sample', 'CDR3b'] + Expects df with columns: ['sample', 'junction_aa'] """ # Calculate lengths - df['cdr3_len'] = df['CDR3b'].apply(len) + df['cdr3_len'] = df['junction_aa'].apply(len) # 1. Basic Stats per sample stats = df.groupby('sample')['cdr3_len'].agg( From 5bf055ea1906e085a7cdd4bb849b6b26b2865bfb Mon Sep 17 00:00:00 2001 From: KevinMLanderos Date: Fri, 3 Jul 2026 13:49:35 -0400 Subject: [PATCH 16/21] Update template_qc.qmd to use AIRR column names (junction_aa, v_call, duplicate_count). Update samplesheet to add 'alias' column, which is a sample's short name to be used for display purposes in notebooks. Added parameters to nextflow_schema.json. --- nextflow.config | 2 +- nextflow_schema.json | 27 ++++++- notebooks/template_qc.qmd | 2 +- .../test_data/minimal-example/samplesheet.csv | 72 +++++++++---------- 4 files changed, 64 insertions(+), 39 deletions(-) diff --git a/nextflow.config b/nextflow.config index dc21d9f..ccc4534 100644 --- a/nextflow.config +++ b/nextflow.config @@ -41,7 +41,7 @@ params { timepoint_col = 'timepoint' timepoint_order_col = 'timepoint_order' alias_col = 'alias' - subject_col = 'subject_id' + subject_col = 'patient' // OLGA parameters olga_chunk_length = 100000 // larger chunk size = less parallelization diff --git a/nextflow_schema.json b/nextflow_schema.json index 48aae9e..26d3c20 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -80,6 +80,11 @@ "default": "copy", "enum": ["copy", "move", "link", "symlink"], "description": "Method used by `publishDir` to save outputs." + }, + "container": { + "type": "string", + "description": "Docker/Singularity container image to use for pipeline processes.", + "fa_icon": "fas fa-box" } } }, @@ -128,7 +133,27 @@ "samplechart_x_col": { "type": "string", "default": "timepoint" }, "samplechart_color_col": { "type": "string", "default": "origin" }, "vgene_subject_col": { "type": "string", "default": "patient" }, - "vgene_x_cols": { "type": "string", "default": "origin,timepoint" } + "vgene_x_cols": { "type": "string", "default": "origin,timepoint" }, + "subject_col": { + "type": "string", + "default": "patient", + "description": "Samplesheet column identifying the subject/patient." + }, + "timepoint_col": { + "type": "string", + "default": "timepoint", + "description": "Samplesheet column identifying the timepoint." + }, + "timepoint_order_col": { + "type": "string", + "default": "timepoint_order", + "description": "Samplesheet column giving the sort order of timepoints." + }, + "alias_col": { + "type": "string", + "default": "alias", + "description": "Samplesheet column giving a display alias for each sample." + } } }, diff --git a/notebooks/template_qc.qmd b/notebooks/template_qc.qmd index 5ac1f5e..5c0d2b3 100644 --- a/notebooks/template_qc.qmd +++ b/notebooks/template_qc.qmd @@ -53,7 +53,7 @@ sample_table='' timepoint_col = 'timepoint' timepoint_order_col = 'timepoint_order' alias_col = 'alias' -subject_col = 'subject_id' +subject_col = 'patient' ``` diff --git a/tests/test_data/minimal-example/samplesheet.csv b/tests/test_data/minimal-example/samplesheet.csv index 79d1b0f..5f89b10 100644 --- a/tests/test_data/minimal-example/samplesheet.csv +++ b/tests/test_data/minimal-example/samplesheet.csv @@ -1,36 +1,36 @@ -sample,patient,timepoint,origin,file -Patient01_Base,Patient01,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient01_Base.tsv -Patient02_Base,Patient02,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient02_Base.tsv -Patient03_Base,Patient03,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient03_Base.tsv -Patient04_Base,Patient04,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient04_Base.tsv -Patient05_Base,Patient05,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient05_Base.tsv -Patient06_Base,Patient06,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient06_Base.tsv -Patient07_Base,Patient07,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient07_Base.tsv -Patient08_Base,Patient08,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient08_Base.tsv -Patient09_Base,Patient09,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient09_Base.tsv -Patient10_Base,Patient10,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient10_Base.tsv -Patient11_Base,Patient11,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient11_Base.tsv -Patient12_Base,Patient12,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient12_Base.tsv -Patient13_Base,Patient13,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient13_Base.tsv -Patient14_Base,Patient14,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient14_Base.tsv -Patient15_Base,Patient15,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient15_Base.tsv -Patient16_Base,Patient16,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient16_Base.tsv -Patient16_Post,Patient16,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient16_Post.tsv -Patient17_Base,Patient17,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient17_Base.tsv -Patient17_Post,Patient17,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient17_Post.tsv -Patient18_Base,Patient18,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient18_Base.tsv -Patient18_Post,Patient18,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient18_Post.tsv -Patient19_Base,Patient19,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient19_Base.tsv -Patient19_Post,Patient19,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient19_Post.tsv -Patient20_Base,Patient20,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient20_Base.tsv -Patient20_Post,Patient20,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient20_Post.tsv -Patient21_Base,Patient21,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient21_Base.tsv -Patient21_Post,Patient21,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient21_Post.tsv -Patient22_Base,Patient22,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient22_Base.tsv -Patient22_Post,Patient22,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient22_Post.tsv -Patient23_Base,Patient23,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient23_Base.tsv -Patient23_Post,Patient23,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient23_Post.tsv -Patient24_Base,Patient24,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient24_Base.tsv -Patient24_Post,Patient24,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient24_Post.tsv -Patient25_Base,Patient25,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient25_Base.tsv -Patient25_Post,Patient25,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient25_Post.tsv +sample,alias,patient,timepoint,origin,file +Patient01_Base,Patient01_Base,Patient01,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient01_Base.tsv +Patient02_Base,Patient02_Base,Patient02,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient02_Base.tsv +Patient03_Base,Patient03_Base,Patient03,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient03_Base.tsv +Patient04_Base,Patient04_Base,Patient04,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient04_Base.tsv +Patient05_Base,Patient05_Base,Patient05,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient05_Base.tsv +Patient06_Base,Patient06_Base,Patient06,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient06_Base.tsv +Patient07_Base,Patient07_Base,Patient07,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient07_Base.tsv +Patient08_Base,Patient08_Base,Patient08,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient08_Base.tsv +Patient09_Base,Patient09_Base,Patient09,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient09_Base.tsv +Patient10_Base,Patient10_Base,Patient10,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient10_Base.tsv +Patient11_Base,Patient11_Base,Patient11,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient11_Base.tsv +Patient12_Base,Patient12_Base,Patient12,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient12_Base.tsv +Patient13_Base,Patient13_Base,Patient13,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient13_Base.tsv +Patient14_Base,Patient14_Base,Patient14,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient14_Base.tsv +Patient15_Base,Patient15_Base,Patient15,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient15_Base.tsv +Patient16_Base,Patient16_Base,Patient16,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient16_Base.tsv +Patient16_Post,Patient16_Post,Patient16,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient16_Post.tsv +Patient17_Base,Patient17_Base,Patient17,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient17_Base.tsv +Patient17_Post,Patient17_Post,Patient17,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient17_Post.tsv +Patient18_Base,Patient18_Base,Patient18,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient18_Base.tsv +Patient18_Post,Patient18_Post,Patient18,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient18_Post.tsv +Patient19_Base,Patient19_Base,Patient19,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient19_Base.tsv +Patient19_Post,Patient19_Post,Patient19,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient19_Post.tsv +Patient20_Base,Patient20_Base,Patient20,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient20_Base.tsv +Patient20_Post,Patient20_Post,Patient20,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient20_Post.tsv +Patient21_Base,Patient21_Base,Patient21,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient21_Base.tsv +Patient21_Post,Patient21_Post,Patient21,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient21_Post.tsv +Patient22_Base,Patient22_Base,Patient22,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient22_Base.tsv +Patient22_Post,Patient22_Post,Patient22,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient22_Post.tsv +Patient23_Base,Patient23_Base,Patient23,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient23_Base.tsv +Patient23_Post,Patient23_Post,Patient23,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient23_Post.tsv +Patient24_Base,Patient24_Base,Patient24,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient24_Base.tsv +Patient24_Post,Patient24_Post,Patient24,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient24_Post.tsv +Patient25_Base,Patient25_Base,Patient25,Base,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient25_Base.tsv +Patient25_Post,Patient25_Post,Patient25,Post,tumor,https://raw.githubusercontent.com/KarchinLab/TCRtoolkit/refs/heads/main/tests/test_data/minimal-example/PD1_Patient25_Post.tsv From 9e56380d42bbaf55391c5cccbf76d3c8115b214e Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 4 Jul 2026 11:19:00 +0000 Subject: [PATCH 17/21] test: align minimal samplesheet header expectation with alias column --- tests/modules/local/samplesheet/samplesheet_check.nf.test | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/modules/local/samplesheet/samplesheet_check.nf.test b/tests/modules/local/samplesheet/samplesheet_check.nf.test index f405c4f..92328c9 100644 --- a/tests/modules/local/samplesheet/samplesheet_check.nf.test +++ b/tests/modules/local/samplesheet/samplesheet_check.nf.test @@ -55,7 +55,7 @@ nextflow_process { def lines = path(samplesheet_utf8.get(0)).readLines() assert lines.size() == 36 // Header + 35 samples - assert lines[0].replaceAll("^\uFEFF", "") == "sample,patient,timepoint,origin,file" + assert lines[0].replaceAll("^\uFEFF", "") == "sample,alias,patient,timepoint,origin,file" } } } From d628cdb99f9bc45aece71ce2e7d69f99fc5c871b Mon Sep 17 00:00:00 2001 From: dimalvovs Date: Sat, 4 Jul 2026 08:14:33 -0400 Subject: [PATCH 18/21] Apply suggestions from code review Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- README.md | 4 ++-- subworkflows/local/report.nf | 9 +++------ tests/modules/local/report/render_notebook.nf.test | 14 ++++++++++---- 3 files changed, 15 insertions(+), 12 deletions(-) diff --git a/README.md b/README.md index 40ea13b..ec9b44f 100644 --- a/README.md +++ b/README.md @@ -37,8 +37,8 @@ With the update to [Nextflow strict syntax](https://docs.seqera.io/nextflow/stri nextflow run KarchinLab/TCRtoolkit \ -params-file params.yml ``` -[!IMPORTANT] -If having an error similar to `qemu: uncaught target signal 11 (Segmentation fault) - core dumped` (observed on Mac OS Tahoe i.e. `26.5.1`), rebuild the Docker container locally, using the provided Dockerfile, and use that in the params, e.g. `--container `. It features arch selectors to pick the right `quarto` installation. +> [!IMPORTANT] +> If having an error similar to `qemu: uncaught target signal 11 (Segmentation fault) - core dumped` (observed on Mac OS Tahoe i.e. `26.5.1`), rebuild the Docker container locally, using the provided Dockerfile, and use that in the params, e.g. `--container `. It features arch selectors to pick the right `quarto` installation. ## Input Formats diff --git a/subworkflows/local/report.nf b/subworkflows/local/report.nf index 2bac961..b9ff66e 100644 --- a/subworkflows/local/report.nf +++ b/subworkflows/local/report.nf @@ -13,13 +13,12 @@ include { RENDER_NOTEBOOK } from '../../modules/local/report/render_notebook' */ // Single-channel listener: receives tuples of -// (report_name, notebook_path, sample_stats_csv, concat_cdr3_sorted, data_dir) +// (notebook_template, files_to_stage) // and renders each notebook to HTML. workflow REPORT { take: - ch_reports // channel of tuples: tuple(val report_name, path notebook, path sample_stats_csv, path concat_cdr3_sorted, val data_dir) - + ch_reports // channel of tuples: tuple(path notebook_template, path files_to_stage) main: RENDER_NOTEBOOK( ch_reports, @@ -28,7 +27,5 @@ workflow REPORT { ) emit: - - RENDER_NOTEBOOK.out - + rendered_html = RENDER_NOTEBOOK.out } diff --git a/tests/modules/local/report/render_notebook.nf.test b/tests/modules/local/report/render_notebook.nf.test index 9004825..aabc370 100644 --- a/tests/modules/local/report/render_notebook.nf.test +++ b/tests/modules/local/report/render_notebook.nf.test @@ -5,8 +5,8 @@ nextflow_process { process "RENDER_NOTEBOOK" test("Should render basic notebook") { - // note that on Mac M series, jupyter may fail in container - // so the test notebook uses knitr engine instead + // Note: this fixture uses the Jupyter `python3` engine (see tests/fixtures/test.qmd) + // If Jupyter kernels are not available in the container, this test will fail. tag "basic" tag "notebook" @@ -21,7 +21,7 @@ nextflow_process { """ input[0] = [ file("${projectDir}/tests/fixtures/test.qmd"), - "${projectDir}/tests/fixtures" + file("${projectDir}/tests/fixtures") ] input[1] = "TCRtoolkit" input[2] = "nextflow run main.nf" @@ -31,7 +31,13 @@ nextflow_process { then { assert process.success - assert process.out.size() == 1 + with(process.out) { + assert size() == 1 + def html_file = get(0) + def html_path = path(html_file) + assert html_path.exists() + assert html_path.getFileName().toString().endsWith('.html') + } } } From 0493dcdee19e98afb36e37bb6c7d5cc16a057de1 Mon Sep 17 00:00:00 2001 From: dimalvovs Date: Sat, 4 Jul 2026 09:01:18 -0400 Subject: [PATCH 19/21] adjust copilot suggestions --- modules/local/report/render_notebook.nf | 2 +- tests/modules/local/report/render_notebook.nf.test | 10 ++++------ 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/modules/local/report/render_notebook.nf b/modules/local/report/render_notebook.nf index 5b15427..a7f8136 100644 --- a/modules/local/report/render_notebook.nf +++ b/modules/local/report/render_notebook.nf @@ -9,7 +9,7 @@ process RENDER_NOTEBOOK { val workflow_cmd output: - path "${notebook.getBaseName()}.html" + path "${notebook.getBaseName()}.html", emit: report_html script: """ diff --git a/tests/modules/local/report/render_notebook.nf.test b/tests/modules/local/report/render_notebook.nf.test index aabc370..eba24f6 100644 --- a/tests/modules/local/report/render_notebook.nf.test +++ b/tests/modules/local/report/render_notebook.nf.test @@ -31,12 +31,10 @@ nextflow_process { then { assert process.success - with(process.out) { - assert size() == 1 - def html_file = get(0) - def html_path = path(html_file) - assert html_path.exists() - assert html_path.getFileName().toString().endsWith('.html') + with(process.out.report_html) { + def html_file = path(get(0)) + assert html_file.exists() + assert html_file.getFileName().toString().endsWith('.html') } } } From 353bbc4fffae28c852719d6b9972927e13394281 Mon Sep 17 00:00:00 2001 From: KevinMLanderos Date: Sat, 4 Jul 2026 12:29:06 -0400 Subject: [PATCH 20/21] added 'alias' column to samplesheet_check.nf.test --- tests/modules/local/samplesheet/samplesheet_check.nf.test | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/modules/local/samplesheet/samplesheet_check.nf.test b/tests/modules/local/samplesheet/samplesheet_check.nf.test index 92328c9..7a07c13 100644 --- a/tests/modules/local/samplesheet/samplesheet_check.nf.test +++ b/tests/modules/local/samplesheet/samplesheet_check.nf.test @@ -28,7 +28,7 @@ nextflow_process { def lines = path(samplesheet_utf8.get(0)).readLines() assert lines.size() == 5 // Header + 4 data rows // The output has UTF-8 BOM, so strip it for comparison - assert lines[0].replaceAll("^\uFEFF", "") == "sample,patient,timepoint,origin,file" + assert lines[0].replaceAll("^\uFEFF", "") == "sample,alias,patient,timepoint,origin,file" assert lines[1].startsWith("Patient01_Base,Patient01,Base,tumor,") } } @@ -79,7 +79,7 @@ nextflow_process { with(process.out) { def lines = path(samplesheet_utf8.get(0)).readLines() assert lines.size() == 1 // Only header - assert lines[0].replaceAll("^\uFEFF", "") == "sample,patient,timepoint,origin,file" + assert lines[0].replaceAll("^\uFEFF", "") == "sample,alias,patient,timepoint,origin,file" } } } From 08a89919ccc6484a78cf1137c0f3290193427259 Mon Sep 17 00:00:00 2001 From: KevinMLanderos Date: Sat, 4 Jul 2026 12:47:59 -0400 Subject: [PATCH 21/21] Revert recent commits to restore state at 0493dcd --- tests/modules/local/samplesheet/samplesheet_check.nf.test | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/modules/local/samplesheet/samplesheet_check.nf.test b/tests/modules/local/samplesheet/samplesheet_check.nf.test index 7a07c13..92328c9 100644 --- a/tests/modules/local/samplesheet/samplesheet_check.nf.test +++ b/tests/modules/local/samplesheet/samplesheet_check.nf.test @@ -28,7 +28,7 @@ nextflow_process { def lines = path(samplesheet_utf8.get(0)).readLines() assert lines.size() == 5 // Header + 4 data rows // The output has UTF-8 BOM, so strip it for comparison - assert lines[0].replaceAll("^\uFEFF", "") == "sample,alias,patient,timepoint,origin,file" + assert lines[0].replaceAll("^\uFEFF", "") == "sample,patient,timepoint,origin,file" assert lines[1].startsWith("Patient01_Base,Patient01,Base,tumor,") } } @@ -79,7 +79,7 @@ nextflow_process { with(process.out) { def lines = path(samplesheet_utf8.get(0)).readLines() assert lines.size() == 1 // Only header - assert lines[0].replaceAll("^\uFEFF", "") == "sample,alias,patient,timepoint,origin,file" + assert lines[0].replaceAll("^\uFEFF", "") == "sample,patient,timepoint,origin,file" } } }