Init channel-based reports

.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:

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
+

README.md

@@ -37,6 +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 <my image>`. It features arch selectors to pick the right `quarto` installation.
 
 ## Input Formats
 
@@ -92,4 +94,4 @@ Certain sub-reports are automatically appended based on input and workflow optio
 | `--sobject_gex` | — | Path to TSV file containing cell-barcode phenotypes for pseudo-bulk phenotyping |
 | `--max_memory` | `768.GB` | Maximum memory allocation |
 | `--max_cpus` | `192` | Maximum CPU allocation |
-
+

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()}" },

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 }
+        ]
+    }
+
 }

env.yml

@@ -11,9 +11,7 @@ dependencies:
   - seaborn=0.13.0
   - dash>=2.15.0
   - 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
   - upsetplot=0.9.0
 
   # R and R packages
@@ -39,8 +39,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

modules/local/report/render_notebook.nf

@@ -0,0 +1,28 @@
+// Generic process to render a Quarto notebook to HTML
+process RENDER_NOTEBOOK {
+    tag "${notebook.getBaseName()}"
+    label 'process_single'
+
+    input:
+    tuple path(notebook), path(files) // path(files) just stages files in root dir
+    val project_name
+    val workflow_cmd
+
+    output:
+    path "${notebook.getBaseName()}.html", emit: report_html
+
+    script:
+    """
+    ## render qmd report to html
+    quarto render ${notebook} \\
+        -P project_name:${project_name} \\
+        -P workflow_cmd:'${workflow_cmd}' \\
+        -P sample_table:${file(params.samplesheet)} \\
+        --to html
+    """
+
+    stub:
+    """
+    touch ${notebook.getBaseName()}.html
+    """
+}

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
@@ -69,6 +69,12 @@ 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"
+
+    //reports
+    template_qc = "${projectDir}/notebooks/template_qc.qmd"
 }
 
 includeConfig 'conf/base.config'

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"
         }
       }
     },
@@ -111,6 +116,10 @@
         "compare_stats_template": {
           "type": "string",
           "description": "Path to compare notebook template."
+        },
+        "template_qc": {
+          "type": "string",
+          "description": "Path to QC notebook template."
         }
       }
     },
@@ -124,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."
+        }
       }
     },
 

notebooks/template_qc.qmd

@@ -48,13 +48,12 @@ Thank you for using TCRtoolkit! This report is generated from the data provided.
 # ---------------------------------------------------------
 workflow_cmd = '<command used to run the pipeline>'
 project_name='<project_name>'
-project_dir='<path/to/project_dir>'
 sample_table='<path/to/sample_table.csv>'
 
 timepoint_col = 'timepoint'
 timepoint_order_col = 'timepoint_order'
 alias_col = 'alias'
-subject_col = 'subject_id'
+subject_col = 'patient'
 
 ```
 
@@ -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"
 
 
 ```
@@ -189,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:
@@ -1018,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)
@@ -1838,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 ---
@@ -1978,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'
@@ -2000,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:
@@ -2134,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()
@@ -2255,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(

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