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 peptides

For variant (VCF) input, the pipeline makes peptides with bcftools, Ensembl VEP and pVACtools. See Variant input for the variants that the pipeline uses.

The pipeline keeps only peptides that overlap the mutation:

The peptide lengths are set by --min_peptide_length_classI, --max_peptide_length_classI and the class II equivalents.

Output directories:

• variant_peptides/[sample]_length_[k].tsv: the peptides of length k with their variant information.
• variant_fasta/[sample].annotated.fasta: the wild-type and mutant protein sequence around each variant. See Variant FASTA.
• references/: the VEP cache and the genome FASTA. The pipeline writes this directory only with --vep_download_cache.

Peptide tables

Each peptide table has these columns:

A column contains values separated by ; if a peptide comes from more than one variant or transcript.

For each peptide length k, the pipeline takes k - 1 residues on each side of the mutation, as pvacseq run does. Thus each peptide of length k contains the mutation.

The pipeline also makes sequences that combine nearby somatic missense variants (see Nearby somatic variants). A combined sequence belongs to one of its variants. A peptide that occurs in the combined sequence and in the single-variant sequence therefore gets a counts value of 2.

Example: the missense mutation p.Cys138Tyr with --min_peptide_length_classI 9 --max_peptide_length_classI 9 and --wild_type gives this table (some columns are left out):

Without --wild_type, the table has no WT rows.

Variant FASTA

You can use this file as a search database in proteogenomics approaches, for example in nf-core/mhcquant. The search can then identify mutant peptides in immunopeptidomics data.

Each sequence in variant_fasta/ contains --mutation_flanking_aas residues on each side of the mutation (default: 25). For frameshifts, the sequence continues to the new stop codon. This parameter changes only the FASTA file. It does not change the predicted peptides.

The header of each sequence has this format. NA replaces missing values.

>{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

The file contains one record for each variant and transcript. Thus the same sequence can occur several times for different isoforms. If you use the file as a search database, group the proteins in your downstream analysis.

Wild-type peptides

The wildtype column contains the wild-type peptide of each missense mutant peptide. With --wild_type, the pipeline also predicts these wild-type peptides. It adds each one as a separate row with the same alleles.

• peptide_origin is MT;WT if a sequence is mutant for one variant and wild type for another.
• Wild-type rows have the variant information of their mutant peptide. Their protein_ids have the format WT.<index>.
• Indels, frameshifts and peptides without variant information have wildtype NA. They get no wild-type row. The same applies to wild-type peptides with non-standard residues.
• The pipeline applies --proteome_reference before it adds the wild-type rows. If a mutant peptide occurs in the reference proteome, the pipeline removes it and its wild-type row.
• The MultiQC binder statistics do not include wild-type rows. predictions/[sample].tsv includes them.
• With --binder_only, the pipeline keeps the wild-type row of each mutant binder, also if the wild-type peptide does not bind. In long format, the pipeline matches the rows for each tool and allele.

Binding predictions

Each prediction tool in --tools writes its results to its own directory. To run in parallel, the pipeline splits the peptides into chunks. --peptides_split_minchunksize and --peptides_split_maxchunks control the chunk size.

Output directories:

• 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 parts of the file name are:

• [split]: the peptide length for variant and protein input, for example length_9. Peptide input has no [split].
• _c[0-9]: the peptide chunk.
• _a[0-9]: the allele chunk. The pipeline adds it if a sample has more alleles than a tool accepts in one call, for example netmhcpan/[sample]_length_9_c0_a3_predicted_netmhcpan.xls.

The pipeline converts the results of all tools to one format. It then merges the chunks into one file for each sample.

Output directory: predictions/[sample].tsv

Each file contains these columns:

The file also contains all other columns of the input file.

An example prediction result looks like this:

Binding affinity

The BA column is calculated from the predicted IC50 value in nM (aff):

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

A low IC50 value means strong binding. Peptides with an IC50 below 500 nM are usually considered binders, and peptides below 50 nM strong binders.

MixMHCpred and MixMHC2pred do not predict an IC50 value. For these tools, BA is empty.

Percentile rank

The percentile rank compares the score of a peptide with the scores of a large set of random natural peptides. A rank of 0.1 means that the peptide is in the best 0.1%. The rank is not affected by alleles that have higher or lower mean affinities. Thus you can compare ranks between alleles.

We recommend that you select binders by rank, not by BA.

For NetMHCpan and NetMHCIIpan, rank is the eluted ligand rank (EL_Rank), which the tool developers recommend. The hidden parameter --use_ba_rank selects the binding affinity rank (BA_Rank) instead. The two ranks can be very different.

For MixMHCpred and MixMHC2pred, rank is the %Rank output of the tool.

Binder definition

The binder column uses these thresholds:

MHCnuggets uses BA, because its percentile rank is experimental.

Missing values

A row without prediction values (peptide3 in the example) is a peptide that no tool could predict. The tools do not support the allele or the peptide length. assets/supported_alleles.json lists the supported alleles. Peptide lengths lists the supported lengths. The MultiQC report shows the number of peptides that the pipeline could not predict.

Wide format

With --wide_format_output, the pipeline writes one row for each peptide (wide format). The file has these columns in addition to the input columns:

An example with MHCflurry, NetMHCpan and one allele looks like this:

MultiQC

The MultiQC report shows the number of binders and non-binders for each sample, allele and tool. It also shows the distributions of the prediction scores.

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 MultiQC report with binding statistics and score distributions.

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.