nf-core/epitopeprediction

Edit

Introduction

This document describes the output produced by the pipeline. The version of all tools used in the pipeline are summarized in a MultiQC report which is generated at the end of the pipeline.

The directories listed below will be created in the results directory after the pipeline has finished. All paths are relative to the top-level results directory.

Variant prediction

Variant (VCF) input is processed with an offline chain of bcftools, Ensembl VEP and pVACtools (see usage). Only peptides that overlap the mutation are kept, within the length bounds set by --min_peptide_length_class[I|II] and --max_peptide_length_class[I|II]. That means the mutated residue for missense, the junction for in-frame indels, and the novel C-terminal tail for frameshifts. Each peptide carries provenance (gene, transcript, consequence, HGVSp, genomic anchor, UniProt).

Example: for the missense mutation p.Cys138Tyr with min_peptide_length_classI = max_peptide_length_classI = 9, the length-9 table looks like this (WT counterpart shown when --wild_type is set):

Tables are written per peptide length as a tsv, then passed to the MHC binding prediction subworkflow where they are scored against the sample’s individual MHC alleles.

Output directories:

• variant_peptides/[sample]_length_[k].tsv — mutation-overlapping peptides with provenance. As in pvacseq run, windows are cut with k - 1 residues on each side of the mutation for each peptide length k, so every k-mer covers the mutation; nearby somatic missense variants are folded in (see usage), and a combined window keeps the identity of the variant it was built for, so its k-mers that also occur in the single-variant window are counted twice in counts
• variant_fasta/[sample].annotated.fasta — WT/MT protein windows with --mutation_flanking_aas residues on each side of the mutation (frameshifts to the new stop) and provenance-annotated headers (schema below), e.g. as a search database for nf-core/mhcquant

Each pvacseq defline is rewritten into a fixed, pipe-delimited schema (NA for any missing value). The values come from pVACtools’ own variant table, joined to the FASTA records on its index:

>{kind}|{numbering}|{genomic_anchor}|{gene}|{transcript}|{uniprot}|{consequence}|{aa_change}|{hgvs}

Example: >MT|170|3:126730598:G:C|CHCHD6|ENST00000290913.8|Q9BRQ6|missense|78Q/H|p.Gln78His

One record is written per variant and transcript, so identical windows can recur across isoforms; deduplicate by protein grouping downstream if you use this FASTA as a search database.

Epitopeprediction

Depending on the specified predictor(s) in --tools, the tools individual binding prediction files are written in the respective directories. The number of input peptides for the MHC binding subworkflow is splitted into chunks to enable scalability. The chunksize is controlled by --peptides_split_minchunksize and --peptides_split_maxchunks.

Tools output directory:

• mhcflurry/[sample]_[split]_c[0-9]_predicted_mhcflurry.csv
• mhcnuggets/[sample]_[split]_c[0-9]_predicted_mhcnuggets.csv
• mhcnuggetsii/[sample]_[split]_c[0-9]_predicted_mhcnuggetsii.csv
• netmhcpan/[sample]_[split]_c[0-9]_predicted_netmhcpan.xls
• netmhciipan/[sample]_[split]_c[0-9]_predicted_netmhciipan.xls
• mixmhcpred/[sample]_[split]_c[0-9]_predicted_mixmhcpred.txt
• mixmhciipred/[sample]_[split]_c[0-9]_predicted_mixmhciipred.txt

The name is built from the sample and its split coordinates, so it stays short no matter how many stages ran:

• [split] is the peptide length for variant and protein input (length_9). It is omitted for peptide input.
• _c[0-9] is the peptide chunk, controlled by --peptides_split_minchunksize and --peptides_split_maxchunks.
• _a[0-9] is appended when a sample has more alleles than a predictor accepts per call and they have to be chunked too, e.g. netmhcpan/[sample]_length_9_c0_a3_predicted_netmhcpan.xls.

These predictor-specific output files are harmonized and chunks are merged on the sample information of your samplesheet.

Output directory: predictions/[sample].tsv.

Output files always contain the columns --peptide_col_name (default:‘sequence’), allele, BA, rank, binder, predictor. All further metadata columns are parsed into the output files.

An example prediction result looks like this in TSV format:

The prediction results are given as allele-specific Binding Affinity (BA) and percentile ranks (rank) per peptide. The computation of these values depends on the applied prediction method. Binding Affinity represents the predicted strength of the interaction between a peptide and an MHC molecule. It is derived from the predicted IC50 value (in nanomolar, nM) and normalized to a scale between 0 and 1 using the formula:

$BA = 1 - \frac{\log_{10}(\text{aff})}{\log_{10}(50000)}$

where aff is the predicted IC50 binding affinity. Lower IC50 values indicate stronger binding, with peptides having IC50 values below 500 nM typically considered strong binders.

Percentile rank (rank) indicates the relative binding strength of a peptide compared to a large set of random natural peptides. This measure is not affected by inherent biases of certain MHC molecules towards higher or lower mean predicted affinities. Strong binders are defined as having rank < 0.5, and weak binders with rank < 2. For example, a peptide with a rank of 0.1 is among the top 0.1% of best binders. This approach ensures a more consistent selection across different MHC alleles, as it accounts for variability in binding thresholds. It is advised to select candidate binders based on rank rather than binding affinities. Consequently, the binder column is defined based on the rank. An exception to this is the percentile rank computation of MHCnuggets, which is considered experimental and therefore it is implemented and advised to use the BA column for the binder definition.

For netMHCpan and netMHCIIpan predictions specifically, the pipeline uses EL_Rank (Eluted Ligand Rank) by default, which is the rank metric recommended by the developers. EL_Rank is computed against a reference set of eluted ligands rather than binding affinity measurements. If you prefer to use BA_Rank (Binding Affinity Rank) instead, which correlates more directly with the BA column, you can enable the --use_ba_rank parameter. Note that EL_Rank and BA_Rank can differ significantly.

For MixMHCpred (Class I) and MixMHCIIpred (Class II), the BA column is set to na since these tools output likelihood-based scores rather than IC50 binding affinities. The rank column contains the %Rank percentile rank, with binders defined as %Rank <= 2. Both tools are licensed for academic non-commercial research only and run only with --accept_mixmhcpred_license; see Usage.

Optionally you can provide --wide_format_output to obtain your results in wide format.

An example of the wide format looks like this:

MultiQC

Binding prediction results are summarized into tables, such as the number of binders/non-binders. Binding prediction score distributions are also highlighted to give the user an appropriate overview of the binding prediction results.

Output directory: multiqc/

• multiqc_data/
  • Underlying data to generate MultiQC plots
• multiqc_plots/
  • Plots in pdf, png, and svg format that are part of the MultiQC report
• multiqc_report.html
  • The main multiQC report comprising statistics and distributions of the binding prediction results.

For more information about how to use MultiQC reports, see http://multiqc.info.

Pipeline information

Nextflow provides excellent functionality for generating various reports relevant to the running and execution of the pipeline. This will allow you to troubleshoot errors with the running of the pipeline, and also provide you with other information such as launch commands, run times and resource usage.