diff --git a/modules/msk/generatemutfasta/main.nf b/modules/msk/generatemutfasta/main.nf index b6e3c155..20f3c203 100644 --- a/modules/msk/generatemutfasta/main.nf +++ b/modules/msk/generatemutfasta/main.nf @@ -1,7 +1,7 @@ process GENERATEMUTFASTA { tag "$meta.id" label 'process_single' - container "ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.4.0" + container "ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.6.0" input: tuple val(meta), path(inputMaf) @@ -32,13 +32,13 @@ process GENERATEMUTFASTA { mkdir ${prefix}_out - MUTALYZER_SETTINGS="\$(pwd)/config.txt" generateMutFasta.py --sample_id ${prefix} \ + MUTALYZER_SETTINGS="\$(pwd)/config.txt" generate_mut_fasta.py --sample_id ${prefix} \ --output_dir ${prefix}_out \ --maf_file ${inputMaf} cat <<-END_VERSIONS > versions.yml "${task.process}": - generateMutFasta: \$(echo \$(generateMutFasta.py -v)) + generateMutFasta: \$(echo \$(generate_mut_fasta.py -v)) mutalyzer: \$(echo \$(mutalyzer_normalizer -v | tr '\n' ' ' | awk '{print \$3}')) END_VERSIONS """ @@ -54,7 +54,7 @@ process GENERATEMUTFASTA { touch ${prefix}_out/${prefix}_generate_mut_fasta.log cat <<-END_VERSIONS > versions.yml "${task.process}": - generateMutFasta: \$(echo \$(generateMutFasta.py -v)) + generateMutFasta: \$(echo \$(generate_mut_fasta.py -v)) mutalyzer: \$(echo \$(mutalyzer_normalizer -v | tr '\n' ' ' | awk '{print \$3}')) END_VERSIONS """ diff --git a/modules/msk/mutalyzer/normalizer/main.nf b/modules/msk/mutalyzer/normalizer/main.nf index c70f60c6..21abe610 100644 --- a/modules/msk/mutalyzer/normalizer/main.nf +++ b/modules/msk/mutalyzer/normalizer/main.nf @@ -1,7 +1,7 @@ process MUTALYZER_NORMALIZER { tag "$meta.id" label 'process_single' - container "ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.4.0" + container "ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.6.0" input: tuple val(meta), val(hgvs_description) diff --git a/modules/msk/mutalyzer/retriever/main.nf b/modules/msk/mutalyzer/retriever/main.nf index 33cb7199..bdc2bfa3 100644 --- a/modules/msk/mutalyzer/retriever/main.nf +++ b/modules/msk/mutalyzer/retriever/main.nf @@ -1,7 +1,7 @@ process MUTALYZER_RETRIEVER { tag "$meta.id" label 'process_medium' - container "ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.4.0" + container "ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.6.0" input: tuple val(meta), path(fasta), path(gff3) diff --git a/modules/msk/neoantigenediting/aligntoiedb/resources/usr/bin/align_neoantigens_to_IEDB.py b/modules/msk/neoantigenediting/aligntoiedb/resources/usr/bin/align_neoantigens_to_IEDB.py deleted file mode 100755 index 69d8d224..00000000 --- a/modules/msk/neoantigenediting/aligntoiedb/resources/usr/bin/align_neoantigens_to_IEDB.py +++ /dev/null @@ -1,293 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -"""script for computing alignments of neoantigens to IEDB peptides - Copyright (C) 2022 Marta Luksza -""" - - -import json -import os -import subprocess -import tempfile -from collections import defaultdict -import numpy as np -import argparse - -import pandas as pd -from Bio import SeqIO - -try: - from Bio.Align import PairwiseAligner, substitution_matrices - _USE_PAIRWISE_ALIGNER = True -except ImportError: - from Bio.pairwise2 import align as _pairwise2_align - _USE_PAIRWISE_ALIGNER = False - - -def load_blosum62_mat(): - raw_blosum62_mat_str = """ - A R N D C Q E G H I L K M F P S T W Y V B Z X * -A 4 -1 -2 -2 0 -1 -1 0 -2 -1 -1 -1 -1 -2 -1 1 0 -3 -2 0 -2 -1 0 -4 -R -1 5 0 -2 -3 1 0 -2 0 -3 -2 2 -1 -3 -2 -1 -1 -3 -2 -3 -1 0 -1 -4 -N -2 0 6 1 -3 0 0 0 1 -3 -3 0 -2 -3 -2 1 0 -4 -2 -3 3 0 -1 -4 -D -2 -2 1 6 -3 0 2 -1 -1 -3 -4 -1 -3 -3 -1 0 -1 -4 -3 -3 4 1 -1 -4 -C 0 -3 -3 -3 9 -3 -4 -3 -3 -1 -1 -3 -1 -2 -3 -1 -1 -2 -2 -1 -3 -3 -2 -4 -Q -1 1 0 0 -3 5 2 -2 0 -3 -2 1 0 -3 -1 0 -1 -2 -1 -2 0 3 -1 -4 -E -1 0 0 2 -4 2 5 -2 0 -3 -3 1 -2 -3 -1 0 -1 -3 -2 -2 1 4 -1 -4 -G 0 -2 0 -1 -3 -2 -2 6 -2 -4 -4 -2 -3 -3 -2 0 -2 -2 -3 -3 -1 -2 -1 -4 -H -2 0 1 -1 -3 0 0 -2 8 -3 -3 -1 -2 -1 -2 -1 -2 -2 2 -3 0 0 -1 -4 -I -1 -3 -3 -3 -1 -3 -3 -4 -3 4 2 -3 1 0 -3 -2 -1 -3 -1 3 -3 -3 -1 -4 -L -1 -2 -3 -4 -1 -2 -3 -4 -3 2 4 -2 2 0 -3 -2 -1 -2 -1 1 -4 -3 -1 -4 -K -1 2 0 -1 -3 1 1 -2 -1 -3 -2 5 -1 -3 -1 0 -1 -3 -2 -2 0 1 -1 -4 -M -1 -1 -2 -3 -1 0 -2 -3 -2 1 2 -1 5 0 -2 -1 -1 -1 -1 1 -3 -1 -1 -4 -F -2 -3 -3 -3 -2 -3 -3 -3 -1 0 0 -3 0 6 -4 -2 -2 1 3 -1 -3 -3 -1 -4 -P -1 -2 -2 -1 -3 -1 -1 -2 -2 -3 -3 -1 -2 -4 7 -1 -1 -4 -3 -2 -2 -1 -2 -4 -S 1 -1 1 0 -1 0 0 0 -1 -2 -2 0 -1 -2 -1 4 1 -3 -2 -2 0 0 0 -4 -T 0 -1 0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1 1 5 -2 -2 0 -1 -1 0 -4 -W -3 -3 -4 -4 -2 -2 -3 -2 -2 -3 -2 -3 -1 1 -4 -3 -2 11 2 -3 -4 -3 -2 -4 -Y -2 -2 -2 -3 -2 -1 -2 -3 2 -1 -1 -2 -1 3 -3 -2 -2 2 7 -1 -3 -2 -1 -4 -V 0 -3 -3 -3 -1 -2 -2 -3 -3 3 1 -2 1 -1 -2 -2 0 -3 -1 4 -3 -2 -1 -4 -B -2 -1 3 4 -3 0 1 -1 0 -3 -4 0 -3 -3 -2 0 -1 -4 -3 -3 4 1 -1 -4 -Z -1 0 0 1 -3 3 4 -2 0 -3 -3 1 -1 -3 -1 0 -1 -3 -2 -2 1 4 -1 -4 -X 0 -1 -1 -1 -2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -2 0 0 -2 -1 -1 -1 -1 -1 -4 -* -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 1 -""" - amino_acids = "ACDEFGHIKLMNPQRSTVWY" - blosum62_mat_str_list = [ - l.split() for l in raw_blosum62_mat_str.strip().split("\n") - ] - blosum_aa_order = [blosum62_mat_str_list[0].index(aa) for aa in amino_acids] - - blosum62_mat = np.zeros((len(amino_acids), len(amino_acids))) - for i, bl_ind in enumerate(blosum_aa_order): - blosum62_mat[i] = np.array( - [int(x) for x in blosum62_mat_str_list[bl_ind + 1][1:]] - )[blosum_aa_order] - blosum62 = { - (aaA, aaB): blosum62_mat[i, j] - for i, aaA in enumerate(amino_acids) - for j, aaB in enumerate(amino_acids) - } - return blosum62 - - -class _AlignmentResult: - """Minimal alignment result with a score attribute.""" - def __init__(self, score=0.0): - self.score = score - - -def align_peptides(seq1, seq2, matrix): - gap_open = -11 - gap_extend = -1 - s1 = seq1.upper() - s2 = seq2.upper() - if _USE_PAIRWISE_ALIGNER: - blosum62 = substitution_matrices.load("BLOSUM62") - aligner = PairwiseAligner() - aligner.mode = "local" - aligner.substitution_matrix = blosum62 - aligner.open_gap_score = gap_open - aligner.extend_gap_score = gap_extend - alignments = aligner.align(s1, s2) - if not alignments: - return _AlignmentResult(0.0) - return alignments[0] - else: - aln = _pairwise2_align.localds(s1, s2, matrix, gap_open, gap_extend) - if not aln: - return _AlignmentResult(0.0) - return aln[0] - - -def run_blastp_n(pep_list, blastdb): - """ - Run BLASTP on the given n neoantigens - - :param pep_list: list - list of peptides (neoantigen sequnces) - - :param blastdb: str - fasta file with IEDB peptides - - :return: dict - str (peptide) -> list of IEDB identifiers - """ - - if blastdb is None: - raise ValueError("No BLAST database specified") - os_fid, fa_file = tempfile.mkstemp(suffix=".fa", dir=os.getcwd()) - os.close(os_fid) - os_fid, txt_file = tempfile.mkstemp(suffix=".txt", dir=os.getcwd()) - os.close(os_fid) - id2seq = {} - seq2id = {} - with open(fa_file, "w") as fh: - for seqid, neoseq in enumerate(pep_list): - seqid2 = "seq_" + str(seqid) - # seqid2id[seqid2] = seqid - id2seq[seqid2] = neoseq - seq2id[neoseq] = seqid2 - fh.write(">seq_{}\n{}\n".format(seqid, neoseq)) - # run BLASTP - blastpexe = "blastp" - blast_args = [ - blastpexe, - "-db", - blastdb, - "-query", - fa_file, - "-outfmt", - "6 qseqid sacc score", - "-gapopen", - "32767", - "-gapextend", - "32767", - "-evalue", - "1e6", - "-max_hsps", - "1", - "-matrix", - "BLOSUM62", - "-max_target_seqs", - "10000000", - "-out", - txt_file, - ] - - subprocess.check_call(blast_args) - os.unlink(fa_file) - alignments = defaultdict(list) - with open(txt_file) as fh: - for line in fh: - S = line.split() - epi_id = int(S[1].split("|")[0]) - seq_id = id2seq[S[0]] - alignments[seq_id].append(epi_id) - os.unlink(txt_file) - return alignments - - -def run_blastp(peplist, blastdb, n=1000): - """ - Blast peptides in neolist against peptides in blastdb. - - :param peplist: list - list of peptides (neoantigens) - - :param blastdb: - iedb fasta file - - :param n: int - run blastp in batches of size n - - :return: dict - dictionary mapping neoantigen peptide sequences to alignment candidates - """ - - alignments = defaultdict(set) - for i in range(0, len(peplist) + n, n): # run blastp in batches of size n - peplist0 = peplist[i : (i + n)] - if len(peplist0) == 0: - continue - alignments0 = run_blastp_n(peplist0, blastdb) - for pepseq in alignments0: - for epi in alignments0[pepseq]: - alignments[pepseq].add(epi) - return alignments - - -def prepare_blastdb(peptidesfasta): - """ - Builds BLAST database - - :param peptidesfasta: str - path to the IEDB.fasta file - """ - instr = ["makeblastdb", "-in", peptidesfasta, "-dbtype", "prot", ">", "/dev/null"] - instr = "\t".join(instr) - os.system(instr) - - -def load_epitopes(iedbfasta): - """ - Load IEDB epitopes from fasta file - - :param iedbfasta: str - - :return: dict - IEDB epitope identifiers mapped to epitope sequence - """ - epitopes = {} - with open(iedbfasta) as f: - seqs = SeqIO.parse(f, "fasta") - for seq in seqs: - seqid = int((seq.id).split("|")[0]) - epitopes[seqid] = str(seq.seq) - return epitopes - - -if __name__ == "__main__": - - """ - - Aligns neoantigens peptides of all patients to IEDB - Requirement: blastp installed and in the PATH - - run as: - python align_neoantigens_to_IEDB.py - - """ - - parser = argparse.ArgumentParser(prog="align_neoantigens_to_IEDB") - parser.add_argument("--fasta", help="IEDB fasta file", required=True) - parser.add_argument("--input", help="patient_data file", required=True) - - args = parser.parse_args() - - iedb_file = args.fasta - patient_file = args.input - - # blosum62 - - blosum62 = load_blosum62_mat() - - # prepare blast database - prepare_blastdb(iedb_file) - epitopes = load_epitopes(iedb_file) - - with open(patient_file) as f: - pjson = json.load(f) - patient = pjson["patient"] - neoantigens = pjson["neoantigens"] - peptides = set( - [("_".join(neo["id"].split("_")[:-1]), neo["sequence"]) for neo in neoantigens] - ) - pepseq2pepid = defaultdict(set) - for pep_id, pep_seq in peptides: - pepseq2pepid[pep_seq].add(pep_id) - - seqlist = list(set([pep_seq for pep_id, pep_seq in peptides])) - alignments = run_blastp(seqlist, iedb_file, n=100) - scores = [] - aln_data = [] - for pep_seq in alignments: - for epitope_id in alignments[pep_seq]: - episeq = epitopes[epitope_id] - score = align_peptides(pep_seq, episeq, blosum62).score - pep_ids = pepseq2pepid[pep_seq] - for pep_id in pep_ids: - aln_data.append([pep_id, pep_seq, epitope_id, score]) - if len(aln_data): - aln_data = pd.DataFrame(aln_data) - aln_data.columns = [ - "Peptide_ID", - "Peptide", - "Epitope_ID", - "Alignment_score", - ] - else: - aln_data = pd.DataFrame( - columns=["Peptide_ID", "Peptide", "Epitope_ID", "Alignment_score"] - ) - aln_data.to_csv("iedb_alignments_" + patient + ".txt", sep="\t", index=False) diff --git a/modules/msk/neoantigenediting/computefitness/main.nf b/modules/msk/neoantigenediting/computefitness/main.nf index 62bf35a1..a39e029c 100644 --- a/modules/msk/neoantigenediting/computefitness/main.nf +++ b/modules/msk/neoantigenediting/computefitness/main.nf @@ -3,8 +3,8 @@ process NEOANTIGENEDITING_COMPUTEFITNESS { label 'process_medium' container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? - 'docker://ghcr.io/mskcc-omics-workflows/neoantigen-editing:1.3': - 'ghcr.io/mskcc-omics-workflows/neoantigen-editing:1.3' }" + 'docker://ghcr.io/mskcc-omics-workflows/neoantigen-editing:1.1': + 'ghcr.io/mskcc-omics-workflows/neoantigen-editing:1.1' }" input: tuple val(meta), path(patient_data), path(alignment_file) diff --git a/modules/msk/neoantigenediting/computefitness/resources/usr/bin/EpitopeDistance.py b/modules/msk/neoantigenediting/computefitness/resources/usr/bin/EpitopeDistance.py deleted file mode 100755 index ea090734..00000000 --- a/modules/msk/neoantigenediting/computefitness/resources/usr/bin/EpitopeDistance.py +++ /dev/null @@ -1,86 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -"""Class for computing the crossreactivity distance between two epitopes - Copyright (C) 2022 Zachary Sethna - - Use is subject to the included term of use found at - https://github.com/LukszaLab/NeoantigenEditing -""" -import numpy as np -import json -import os - -#% -class EpitopeDistance(object): - """Base class for epitope crossreactivity. - - Model: - dist({a_i}, {b_i}) = \sum_i d_i M_ab(a_i, b_i) - - Attributes - ---------- - amino_acids : str - Allowed amino acids in specified order. - - amino_acid_dict : dict - Dictionary of amino acids and corresponding indicies - - d_i : ndarray - Position scaling array d_i. - d_i.shape == (9,) - - M_ab : ndarray - Amino acid substitution matrix. Indexed by the order of amino_acids. - M_ab.shape == (20, 20) - - - """ - - def __init__( - self, - model_file=os.path.join( - os.path.dirname(__file__), "distance_data", "epitope_distance_model_parameters.json" - ), - amino_acids="ACDEFGHIKLMNPQRSTVWY", - ): - """Initialize class and compute M_ab.""" - - self.amino_acids = amino_acids - # self.amino_acid_dict = {aa: i for i, aa in enumerate(self.amino_acids)} - self.amino_acid_dict = {} - for i, aa in enumerate(self.amino_acids): - self.amino_acid_dict[aa.upper()] = i - self.amino_acid_dict[aa.lower()] = i - - self.set_model(model_file) - - def set_model(self, model_file): - """Load model and format substitution matrix M_ab.""" - with open(model_file, "r") as modelf: - c_model = json.load(modelf) - self.d_i = c_model["d_i"] - self.M_ab_dict = c_model["M_ab"] - M_ab = np.zeros((len(self.amino_acids), len(self.amino_acids))) - for i, aaA in enumerate(self.amino_acids): - for j, aaB in enumerate(self.amino_acids): - M_ab[i, j] = self.M_ab_dict[aaA + "->" + aaB] - self.M_ab = M_ab - - def epitope_dist(self, epiA, epiB): - """Compute the model difference between the 9-mers epiA and epiB. - - Ignores capitalization. - - Model: - dist({a_i}, {b_i}) = \sum_i d_i M_ab(a_i, b_i) - """ - - return sum( - [ - self.d_i[i] - * self.M_ab[ - self.amino_acid_dict[epiA[i]], self.amino_acid_dict[epiB[i]] - ] - for i in range(9) - ] - ) diff --git a/modules/msk/neoantigenediting/computefitness/resources/usr/bin/compute_fitness.py b/modules/msk/neoantigenediting/computefitness/resources/usr/bin/compute_fitness.py deleted file mode 100755 index 77ce239f..00000000 --- a/modules/msk/neoantigenediting/computefitness/resources/usr/bin/compute_fitness.py +++ /dev/null @@ -1,355 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -"""script for computing neoantigen qualities and fitness of tumor clones - Copyright (C) 2022 Marta Luksza -""" - -import argparse -import json - -from collections import defaultdict - -import numpy as np -import pandas as pd - -from EpitopeDistance import EpitopeDistance - - -def fill_up_clone_mutations(tree, mut2missense): - """ - Fills up the field with all mutations for each clone on the tree - - :param tree: dict - imported from json file - - :param mut2missense: dict - maps mutation identifiers to 1 if missense else to 0 - - :return: dict - json dictionary with filled up mutations - """ - - nodes = [(tree["topology"], [])] - while len(nodes) > 0: - (node, anc_mutations) = nodes[0] - nodes = nodes[1:] - cmutations = node["clone_mutations"] - node["all_mutations"] = list(set(cmutations + anc_mutations)) - node["TMB"] = sum([mut2missense[mid] for mid in node["all_mutations"]]) - if "children" in node: - for child in node["children"]: - nodes.append((child, node["all_mutations"])) - - -def fill_up_clone_neoantigens(tree, mut2neo): - """ - Adds neoantigen field for each clone on the tree - - :param tree: dict - tree imported from json file - - :param mut2neo: dict - str -> list - mapts mutation identifiers to neoantigen entries - - :return: dict - annotated json dictionary - """ - - nodes = [tree["topology"]] - while len(nodes) > 0: - node = nodes[0] - nodes = nodes[1:] - node["neoantigens"] = [ - neo["id"] for mid in node["all_mutations"] for neo in mut2neo[mid] - ] - node["neoantigen_load"] = len(node["neoantigens"]) - node["NA_Mut"] = sum([len(mut2neo[mid]) > 0 for mid in node["all_mutations"]]) - if "children" in node: - for child in node["children"]: - nodes.append(child) - - -def mark_driver_gene_mutations(pjson): - """ - Create a dictionary mapping mutation identifiers to their driver gene status - (1 - in a driver gene, 0 - not in a driver gene) - - :param pjson: dict - json representation of a sample - - :return dict - (str -> int) - """ - - dg_genes = set(["TP53", "KRAS", "CDKN2A", "SMAD4"]) - mutid2dg = {} - for mut in pjson["mutations"]: - mutid2dg[mut["id"]] = mut["gene"] in dg_genes - return mutid2dg - - -def mark_missense_mutations(pjson): - """ - Create a dictionary mapping mutation identifiers to their missense status - (1 - is a missense mutation, 0 - not a missense mutation) - - :param pjson: dict - json representation of a sample - - :return dict - (str -> int) - """ - - mut2missense = {} - for mut in pjson["mutations"]: - mut2missense[mut["id"]] = mut["missense"] - return mut2missense - - -def map_neoantigen_qualities(pjson): - neoid2quality = {} - for neo in pjson["neoantigens"]: - neoid2quality[neo["id"]] = neo["quality"] - return neoid2quality - - -def get_property(tree, property): - """ - Auxiliary function to extract a clone attribute values into a list - :param tree: dict - json representation of a tree - - :param property: str - name of the attribute - - :return list - """ - - nodes = [tree["topology"]] - vals = [] - while nodes: - node = nodes[0] - nodes = nodes[1:] - if "children" in node: - for child in node["children"]: - nodes.append(child) - vals.append((node["clone_id"], node[property])) - vals.sort(key=lambda x: x[0]) - vals = [x for (_, x) in vals] - return vals - - -def compute_effective_sample_size(sample_json): - """ - Computes the effective cancer cell population size for each sample (see Methods, p.10) - - :param sample_json: dict - json representation of a sample - - :return float - """ - - mut_freqs = {} - for mut in sample_json["mutations"]: - mut_freqs[mut["id"]] = [] - for tree in sample_json["sample_trees"]: - clone_muts_list = get_property(tree, "clone_mutations") - freqs = get_property(tree, "X") - for clone_muts, X in zip(clone_muts_list, freqs): - for mid in clone_muts: - mut_freqs[mid].append(X) - variances = [np.var(mut_freqs[mid]) if mut_freqs[mid] else 0 for mid in mut_freqs] - if not variances: - return 0 - avev = np.mean(variances) - if avev == 0 or np.isnan(avev): - return 0 - n = 1 / avev - return n - - -INF = float("inf") - - -def log_sum2(v1, v2): - ma = max(v1, v2) - if ma == -INF: - return -INF - return ma + np.log(np.exp(v1 - ma) + np.exp(v2 - ma)) - - -def log_sum(v): - if len(v): - ma = max(v) - if ma == -INF: - return -INF - return np.log(sum(list(map(lambda x: np.exp(x - ma), v)))) + ma - return -INF - - -def compute_R(scores, a, k): - """ - Computes the value of the R component given the alignment scores and parameters a and k - - :param scores: list - list of alignment scores for a given neoantigen and IEDB epitopes - - :param a: float - shift parameter of the sigmoid function - - :param k: float - slope parameter of the sigmoid function - - :return float - - """ - v = [-k * (a - score) for score in scores] - lgb = log_sum(v) - lZ = log_sum2(0, lgb) - bindProb = np.exp(lgb - lZ) - return bindProb - - -def set_immune_fitness(tree, neo2qualities): - """ - Sets the value of the immune fitness component for each clone in the tree - as the negative of the max quality neoantigen. - - :param tree: dict - json representation of a tree - - :param neo2qualities: dict - mapping neoantigen identifiers to their qualities (str->float) - """ - - nodes = [tree["topology"]] - while len(nodes) > 0: - node = nodes[0] - nodes = nodes[1:] - if node["neoantigens"]: - node["F_I"] = -max([neo2qualities[neoid] for neoid in node["neoantigens"]]) - else: - node["F_I"] = 0 - if "children" in node: - for child in node["children"]: - nodes.append(child) - - -def set_driver_gene_fitness(tree, mut2dg): - """ - Sets the driver-gene fitness component for each clone in the tree - as the number of driver gene mutations in that clone. - - :param tree: dict - json representation of a tree - - :param mut2dg: dict - mapping mutation identifiers to 1 if this is a mutation in a driver gene or to 0. - """ - nodes = [tree["topology"]] - while len(nodes) > 0: - node = nodes[0] - nodes = nodes[1:] - if node["all_mutations"]: - node["F_P"] = sum([mut2dg[mut_id] for mut_id in node["all_mutations"]]) - else: - node["F_P"] = 0 - if "children" in node: - for child in node["children"]: - nodes.append(child) - - -def clean_data(tree): - """ - Removes no longer needed clone attributes. - """ - nodes = [tree["topology"]] - while len(nodes) > 0: - node = nodes[0] - nodes = nodes[1:] - del node["all_mutations"] - del node["clone_mutations"] - del node["neoantigens"] - if "children" in node: - for child in node["children"]: - nodes.append(child) - - -if __name__ == "__main__": - - """ - Computes components contributing to neoantigen quality score, fitness of clones and annotates clones with neoantigens in *_annotated.json - files. - - Run as: - - python compute_fitness.py - - """ - - parser = argparse.ArgumentParser(prog="align_neoantigens_to_IEDB") - parser.add_argument("--alignment", help="neoantigen alignment file", required=True) - parser.add_argument("--input", help="patient_data file", required=True) - parser.add_argument("--a_param", help="weight corresponding to a", default = 22.897590714815188) - parser.add_argument("--k_param", help="weight corresponding to k", default = 1) - parser.add_argument("--w_param", help="weight corresponding to w", default = 0.22402192838740312) - - args = parser.parse_args() - - alignment_file = args.alignment - patient_file = args.input - - a = float(args.a_param) - k = float(args.k_param) - w = float(args.w_param) - - - epidist = EpitopeDistance() - - sample_file = patient_file - output_file = patient_file.replace(".json", "_annotated.json") - - norm = 1 - - with open(sample_file) as f: - sjson = json.load(f) - patient = sjson["patient"] - neoantigens = sjson["neoantigens"] - nalist = [neo["sequence"] for neo in neoantigens] - - alignments = pd.read_csv(alignment_file, sep="\t") - naseq2scores = defaultdict(list) - for r in alignments.itertuples(): - naseq2scores[r.Peptide].append(r.Alignment_score) - - mut2neo = defaultdict(list) - for neo in neoantigens: - score_list = naseq2scores[neo["sequence"]] - neo["R"] = compute_R(score_list, a, k) if score_list else 0.0 - if neo["Kd"] == 0 or neo["KdWT"] == 0: - neo["logC"] = neo["logA"] = neo["quality"] = 0.0 - else: - neo["logC"] = epidist.epitope_dist(neo["sequence"], neo["WT_sequence"]) - neo["logA"] = np.log(neo["KdWT"] / neo["Kd"]) - neo["quality"] = (w * neo["logC"] + (1 - w) * neo["logA"]) * neo["R"] - mut2neo[neo["mutation_id"]].append(neo) - - mut2dg = mark_driver_gene_mutations(sjson) - mut2missense = mark_missense_mutations(sjson) - neo2qualities = map_neoantigen_qualities(sjson) - - for tree in sjson["sample_trees"]: - fill_up_clone_mutations(tree, mut2missense) - fill_up_clone_neoantigens(tree, mut2neo) - set_immune_fitness(tree, neo2qualities) - set_driver_gene_fitness(tree, mut2dg) - - neff = compute_effective_sample_size(sjson) - sjson["Effective_N"] = neff / norm - - for tree in sjson["sample_trees"]: - clean_data(tree) - - with open(output_file, "w") as of: - json.dump(sjson, of, indent=True) diff --git a/modules/msk/neoantigenediting/computefitness/resources/usr/bin/distance_data/epitope_distance_model_parameters.json b/modules/msk/neoantigenediting/computefitness/resources/usr/bin/distance_data/epitope_distance_model_parameters.json deleted file mode 100644 index c2c3261d..00000000 --- a/modules/msk/neoantigenediting/computefitness/resources/usr/bin/distance_data/epitope_distance_model_parameters.json +++ /dev/null @@ -1,432 +0,0 @@ -{ - "terms_of_use": "Use of these data is subject to the terms of use found at https://github.com/LukszaLab/NeoantigenEditing", - "d_i": [ - 0.4265316417189761, 0.7769577054740209, 1.1858350996115459, 1.2011211326454208, 1.6962418186727721, - 1.4795549981576157, 1.38899685274111, 0.7822985074612253, 0.9749012286522828 - ], - "euclid_coords": { - "A": [3.5312316688405128, -4.038507673277327], - "C": [7.272490502087666, -3.5962613422282153], - "D": [-3.3267651105489104, -4.134381581418677], - "E": [-2.791220079193741, -1.9929601506592383], - "F": [3.7349203453934803, 4.27588200724904], - "G": [1.0828701007327262, -7.270285095315686], - "H": [-1.6594024299600756, 1.7960492980168123], - "I": [6.94579315511744, -0.0012785445892339589], - "K": [-2.4118355486039613, -0.9046522878973625], - "L": [6.454412382783402, 0.772812962193054], - "M": [4.769499796942178, 0.5052384809555668], - "N": [-2.1016216436806427, -3.1474943437753105], - "P": [-0.9992387821861454, -6.803979888540144], - "Q": [-1.5702795695104803, -0.9551312866936055], - "R": [-2.6578083660129788, -0.06342218626699866], - "S": [0.8142701713381114, -3.464839915690124], - "T": [2.606749707538103, -2.6220589959640557], - "V": [6.508726391948886, -0.5670130202032231], - "W": [2.7313743259108967, 5.737800116112936], - "Y": [2.2822689847157145, 4.075146976868077] - }, - "M_ab": { - "A->A": 0.0, - "A->C": 3.7673066605568826, - "A->D": 6.858666899068506, - "A->E": 6.645123081893246, - "A->F": 8.316884322665505, - "A->G": 4.054486363743763, - "A->H": 7.80927252686884, - "A->I": 5.287574981132023, - "A->K": 6.718712490316743, - "A->L": 5.629723949153348, - "A->M": 4.709451886491528, - "A->N": 5.702888846341898, - "A->P": 5.307824326502934, - "A->Q": 5.96092500006195, - "A->R": 7.355645531294955, - "A->S": 2.7768641441391444, - "A->T": 1.6914472360189268, - "A->V": 4.573483349970818, - "A->W": 9.808973736427841, - "A->Y": 8.209220399536266, - "C->A": 3.7673066605568826, - "C->C": 0.0, - "C->D": 10.612906903098432, - "C->E": 10.19062536720303, - "C->F": 8.630471802181635, - "C->G": 7.19790603239027, - "C->H": 10.433394720342184, - "C->I": 3.609796735529243, - "C->K": 10.05141434612365, - "C->L": 4.445004171745916, - "C->M": 4.804920734998789, - "C->N": 9.384847912477928, - "C->P": 8.871919951434778, - "C->Q": 9.22876755094634, - "C->R": 10.540008923700857, - "C->S": 6.459557371202782, - "C->T": 4.766362069061748, - "C->V": 3.124048817228584, - "C->W": 10.380098238149005, - "C->Y": 9.151656483431939, - "D->A": 6.858666899068506, - "D->C": 10.612906903098432, - "D->D": 0.0, - "D->E": 2.2073726972863037, - "D->F": 10.981800176179599, - "D->G": 5.410986374502014, - "D->H": 6.160365957017682, - "D->I": 11.07284940900755, - "D->K": 3.3568210278104593, - "D->L": 10.943125305182395, - "D->M": 9.331429674784955, - "D->N": 1.5731887789560899, - "D->P": 3.541769886738501, - "D->Q": 3.6321995942385272, - "D->R": 4.125556147119299, - "D->S": 4.194813374630591, - "D->T": 6.123211387738817, - "D->V": 10.462457242314292, - "D->W": 11.582790030932431, - "D->Y": 9.94271706482523, - "E->A": 6.645123081893246, - "E->C": 10.19062536720303, - "E->D": 2.2073726972863037, - "E->E": 0.0, - "E->F": 9.049248081592172, - "E->G": 6.546658177551038, - "E->H": 3.9544410215949877, - "E->I": 9.93862275896961, - "E->K": 1.1525392080966916, - "E->L": 9.650451830548775, - "E->M": 7.9627558953165725, - "E->N": 1.344803035889008, - "E->P": 5.133917401663634, - "E->Q": 1.6024308032067456, - "E->R": 1.9341446277992693, - "E->S": 3.8943536035377697, - "E->T": 5.434504869409326, - "E->V": 9.408630590313022, - "E->W": 9.500721197130058, - "E->Y": 7.909628018608502, - "F->A": 8.316884322665505, - "F->C": 8.630471802181635, - "F->D": 10.981800176179599, - "F->E": 9.049248081592172, - "F->F": 0.0, - "F->G": 11.846828489539009, - "F->H": 5.937026904981663, - "F->I": 5.3482526666590955, - "F->K": 8.038690415955843, - "F->L": 4.434763767762229, - "F->M": 3.910000926335655, - "F->N": 9.44307888557397, - "F->P": 12.04888386013437, - "F->Q": 7.4504124864668855, - "F->R": 7.726353684727768, - "F->S": 8.273389470312676, - "F->T": 6.989589334998694, - "F->V": 5.581006381543492, - "F->W": 1.7732200005198253, - "F->Y": 1.4664550890160293, - "G->A": 4.054486363743763, - "G->C": 7.19790603239027, - "G->D": 5.410986374502014, - "G->E": 6.546658177551038, - "G->F": 11.846828489539009, - "G->G": 0.0, - "G->H": 9.471983845230948, - "G->I": 9.338753823510922, - "G->K": 7.26183507210032, - "G->L": 9.671860878370772, - "G->M": 8.605231281152383, - "G->N": 5.2094521066253945, - "G->P": 2.1336864685786554, - "G->Q": 6.849844582124883, - "G->R": 8.11982440582287, - "G->S": 3.8149127037464723, - "G->T": 4.891647465499837, - "G->V": 8.624023017473489, - "G->W": 13.112125954558945, - "G->Y": 11.408654021725596, - "H->A": 7.80927252686884, - "H->C": 10.433394720342184, - "H->D": 6.160365957017682, - "H->E": 3.9544410215949877, - "H->F": 5.937026904981663, - "H->G": 9.471983845230948, - "H->H": 0.0, - "H->I": 8.790891788166004, - "H->K": 2.803559283159802, - "H->L": 8.178080668130237, - "H->M": 6.557207973559407, - "H->N": 4.963283345859569, - "H->P": 8.625330025655463, - "H->Q": 2.7526237472532045, - "H->R": 2.110556470225403, - "H->S": 5.8134337062799055, - "H->T": 6.141639435657454, - "H->V": 8.503081322149713, - "H->W": 5.900535571611977, - "H->Y": 4.553137354727778, - "I->A": 5.287574981132023, - "I->C": 3.609796735529243, - "I->D": 11.07284940900755, - "I->E": 9.93862275896961, - "I->F": 5.3482526666590955, - "I->G": 9.338753823510922, - "I->H": 8.790891788166004, - "I->I": 0.0, - "I->K": 9.401132861352892, - "I->L": 0.9168820667305412, - "I->M": 2.234460175971931, - "I->N": 9.578851100030663, - "I->P": 10.4594587842663, - "I->Q": 8.569324927015641, - "I->R": 9.60380258068956, - "I->S": 7.042146794266393, - "I->T": 5.069101322156581, - "I->V": 0.7149005891484997, - "I->W": 7.120277378143191, - "I->Y": 6.194005385847643, - "K->A": 6.718712490316743, - "K->C": 10.05141434612365, - "K->D": 3.3568210278104593, - "K->E": 1.1525392080966916, - "K->F": 8.038690415955843, - "K->G": 7.26183507210032, - "K->H": 2.803559283159802, - "K->I": 9.401132861352892, - "K->K": 0.0, - "K->L": 9.02353822212172, - "K->M": 7.318426697404776, - "K->N": 2.264193709562544, - "K->P": 6.066093946206223, - "K->Q": 0.843068559055228, - "K->R": 0.8764534846716978, - "K->S": 4.118533574677532, - "K->T": 5.304308039139411, - "K->V": 8.926949401129649, - "K->W": 8.400879820164631, - "K->Y": 6.843465356607772, - "L->A": 5.629723949153348, - "L->C": 4.445004171745916, - "L->D": 10.943125305182395, - "L->E": 9.650451830548775, - "L->F": 4.434763767762229, - "L->G": 9.671860878370772, - "L->H": 8.178080668130237, - "L->I": 0.9168820667305412, - "L->K": 9.02353822212172, - "L->L": 0.0, - "L->M": 1.706026531134751, - "L->N": 9.4114041266561, - "L->P": 10.628485573776425, - "L->Q": 8.208621824427695, - "L->R": 9.15051125338444, - "L->S": 7.0547080789185435, - "L->T": 5.131243872079878, - "L->V": 1.3409264240426382, - "L->W": 6.205812582692975, - "L->Y": 5.3209200781931925, - "M->A": 4.709451886491528, - "M->C": 4.804920734998789, - "M->D": 9.331429674784955, - "M->E": 7.9627558953165725, - "M->F": 3.910000926335655, - "M->G": 8.605231281152383, - "M->H": 6.557207973559407, - "M->I": 2.234460175971931, - "M->K": 7.318426697404776, - "M->L": 1.706026531134751, - "M->M": 0.0, - "M->N": 7.781694348961095, - "M->P": 9.31144553586053, - "M->Q": 6.505803737707016, - "M->R": 7.449045677262714, - "M->S": 5.604048881547064, - "M->T": 3.8022989701449266, - "M->V": 2.0431917263232804, - "M->W": 5.615483656942405, - "M->Y": 4.350926772824042, - "N->A": 5.702888846341898, - "N->C": 9.384847912477928, - "N->D": 1.5731887789560899, - "N->E": 1.344803035889008, - "N->F": 9.44307888557397, - "N->G": 5.2094521066253945, - "N->H": 4.963283345859569, - "N->I": 9.578851100030663, - "N->K": 2.264193709562544, - "N->L": 9.4114041266561, - "N->M": 7.781694348961095, - "N->N": 0.0, - "N->P": 3.8190489015448565, - "N->Q": 2.2558324791172084, - "N->R": 3.133822704432521, - "N->S": 2.9331098323975704, - "N->T": 4.737598873449183, - "N->V": 8.988713876632943, - "N->W": 10.114658060508463, - "N->Y": 8.448967078194036, - "P->A": 5.307824326502934, - "P->C": 8.871919951434778, - "P->D": 3.541769886738501, - "P->E": 5.133917401663634, - "P->F": 12.04888386013437, - "P->G": 2.1336864685786554, - "P->H": 8.625330025655463, - "P->I": 10.4594587842663, - "P->K": 6.066093946206223, - "P->L": 10.628485573776425, - "P->M": 9.31144553586053, - "P->N": 3.8190489015448565, - "P->P": 0.0, - "P->Q": 5.876658706110999, - "P->R": 6.941611570958817, - "P->S": 3.799825059499118, - "P->T": 5.521921344947621, - "P->V": 9.760599201522258, - "P->W": 13.084866061501023, - "P->Y": 11.363260737036889, - "Q->A": 5.96092500006195, - "Q->C": 9.22876755094634, - "Q->D": 3.6321995942385272, - "Q->E": 1.6024308032067456, - "Q->F": 7.4504124864668855, - "Q->G": 6.849844582124883, - "Q->H": 2.7526237472532045, - "Q->I": 8.569324927015641, - "Q->K": 0.843068559055228, - "Q->L": 8.208621824427695, - "Q->M": 6.505803737707016, - "Q->N": 2.2558324791172084, - "Q->P": 5.876658706110999, - "Q->Q": 0.0, - "Q->R": 1.4063655296564266, - "Q->S": 3.4618946935227064, - "Q->T": 4.497357176081819, - "Q->V": 8.088323257269053, - "Q->W": 7.956101872064475, - "Q->Y": 6.3360736873501935, - "R->A": 7.355645531294955, - "R->C": 10.540008923700857, - "R->D": 4.125556147119299, - "R->E": 1.9341446277992693, - "R->F": 7.726353684727768, - "R->G": 8.11982440582287, - "R->H": 2.110556470225403, - "R->I": 9.60380258068956, - "R->K": 0.8764534846716978, - "R->L": 9.15051125338444, - "R->M": 7.449045677262714, - "R->N": 3.133822704432521, - "R->P": 6.941611570958817, - "R->Q": 1.4063655296564266, - "R->R": 0.0, - "R->S": 4.860552637259271, - "R->T": 5.853391669257117, - "R->V": 9.180357465640842, - "R->W": 7.918173418697156, - "R->Y": 6.444541795134483, - "S->A": 2.7768641441391444, - "S->C": 6.459557371202782, - "S->D": 4.194813374630591, - "S->E": 3.8943536035377697, - "S->F": 8.273389470312676, - "S->G": 3.8149127037464723, - "S->H": 5.8134337062799055, - "S->I": 7.042146794266393, - "S->K": 4.118533574677532, - "S->L": 7.0547080789185435, - "S->M": 5.604048881547064, - "S->N": 2.9331098323975704, - "S->P": 3.799825059499118, - "S->Q": 3.4618946935227064, - "S->R": 4.860552637259271, - "S->S": 0.0, - "S->T": 1.9807227383836576, - "S->V": 6.3893843494236595, - "S->W": 9.400205949574959, - "S->Y": 7.681563828806445, - "T->A": 1.6914472360189268, - "T->C": 4.766362069061748, - "T->D": 6.123211387738817, - "T->E": 5.434504869409326, - "T->F": 6.989589334998694, - "T->G": 4.891647465499837, - "T->H": 6.141639435657454, - "T->I": 5.069101322156581, - "T->K": 5.304308039139411, - "T->L": 5.131243872079878, - "T->M": 3.8022989701449266, - "T->N": 4.737598873449183, - "T->P": 5.521921344947621, - "T->Q": 4.497357176081819, - "T->R": 5.853391669257117, - "T->S": 1.9807227383836576, - "T->T": 0.0, - "T->V": 4.410060771483332, - "T->W": 8.360787981361641, - "T->Y": 6.705061937224871, - "V->A": 4.573483349970818, - "V->C": 3.124048817228584, - "V->D": 10.462457242314292, - "V->E": 9.408630590313022, - "V->F": 5.581006381543492, - "V->G": 8.624023017473489, - "V->H": 8.503081322149713, - "V->I": 0.7149005891484997, - "V->K": 8.926949401129649, - "V->L": 1.3409264240426382, - "V->M": 2.0431917263232804, - "V->N": 8.988713876632943, - "V->P": 9.760599201522258, - "V->Q": 8.088323257269053, - "V->R": 9.180357465640842, - "V->S": 6.3893843494236595, - "V->T": 4.410060771483332, - "V->V": 0.0, - "V->W": 7.349765799987526, - "V->Y": 6.2779448590733224, - "W->A": 9.808973736427841, - "W->C": 10.380098238149005, - "W->D": 11.582790030932431, - "W->E": 9.500721197130058, - "W->F": 1.7732200005198253, - "W->G": 13.112125954558945, - "W->H": 5.900535571611977, - "W->I": 7.120277378143191, - "W->K": 8.400879820164631, - "W->L": 6.205812582692975, - "W->M": 5.615483656942405, - "W->N": 10.114658060508463, - "W->P": 13.084866061501023, - "W->Q": 7.956101872064475, - "W->R": 7.918173418697156, - "W->S": 9.400205949574959, - "W->T": 8.360787981361641, - "W->V": 7.349765799987526, - "W->W": 0.0, - "W->Y": 1.7222401310301723, - "Y->A": 8.209220399536266, - "Y->C": 9.151656483431939, - "Y->D": 9.94271706482523, - "Y->E": 7.909628018608502, - "Y->F": 1.4664550890160293, - "Y->G": 11.408654021725596, - "Y->H": 4.553137354727778, - "Y->I": 6.194005385847643, - "Y->K": 6.843465356607772, - "Y->L": 5.3209200781931925, - "Y->M": 4.350926772824042, - "Y->N": 8.448967078194036, - "Y->P": 11.363260737036889, - "Y->Q": 6.3360736873501935, - "Y->R": 6.444541795134483, - "Y->S": 7.681563828806445, - "Y->T": 6.705061937224871, - "Y->V": 6.2779448590733224, - "Y->W": 1.7222401310301723, - "Y->Y": 0.0 - }, - "blosum62_reg": 0.1 -} diff --git a/modules/msk/neoantigenutils/convertannotjson/main.nf b/modules/msk/neoantigenutils/convertannotjson/main.nf index 981e6a85..a21122ae 100644 --- a/modules/msk/neoantigenutils/convertannotjson/main.nf +++ b/modules/msk/neoantigenutils/convertannotjson/main.nf @@ -2,8 +2,8 @@ process NEOANTIGENUTILS_CONVERTANNOTJSON { tag "$meta.id" label 'process_single' container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? - 'docker://ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.4.0': - 'ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.4.0' }" + 'docker://ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.6.0': + 'ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.6.0' }" input: tuple val(meta), path(annotatedJSON) diff --git a/modules/msk/neoantigenutils/convertannotjson/resources/usr/bin/convertannotjson.py b/modules/msk/neoantigenutils/convertannotjson/resources/usr/bin/convertannotjson.py deleted file mode 100755 index c927b52d..00000000 --- a/modules/msk/neoantigenutils/convertannotjson/resources/usr/bin/convertannotjson.py +++ /dev/null @@ -1,63 +0,0 @@ -#!/usr/bin/env python3 - -import json -import argparse - -VERSION = 1.1 - - -def process_json_file(json_file_path, output_file_path): - with open(json_file_path, "r") as json_file: - data = json.load(json_file) - - # Define the TSV header - tsv_header = [ - "id", - "mutation_id", - "Gene", - "HLA_gene_id", - "sequence", - "WT_sequence", - "mutated_position", - "Kd", - "KdWT", - "R", - "logC", - "logA", - "quality", - "NMD", - "git_branch", - ] - - # Convert JSON to TSV - tsv_lines = [] - tsv_lines.append("\t".join(tsv_header)) - - for neoantigen in data["neoantigens"]: - tsv_lines.append( - "\t".join(str(neoantigen.get(field, "")) for field in tsv_header[:-1]) - ) - - tsv_output = "\n".join(tsv_lines) - - # Write the TSV output to a file - with open(output_file_path, "w") as tsv_file: - tsv_file.write(tsv_output) - - -def main(): - parser = argparse.ArgumentParser( - description="Process an annotated JSON file and output TSV format." - ) - parser.add_argument("--json_file", help="Path to the annotated JSON file") - parser.add_argument("--output_file", help="Path to the output TSV file") - parser.add_argument( - "-v", "--version", action="version", version="v{}".format(VERSION) - ) - args = parser.parse_args() - - process_json_file(args.json_file, args.output_file) - - -if __name__ == "__main__": - main() diff --git a/modules/msk/neoantigenutils/formatnetmhcpan/main.nf b/modules/msk/neoantigenutils/formatnetmhcpan/main.nf index 58efeeaa..2f2b94c9 100644 --- a/modules/msk/neoantigenutils/formatnetmhcpan/main.nf +++ b/modules/msk/neoantigenutils/formatnetmhcpan/main.nf @@ -2,8 +2,8 @@ process NEOANTIGENUTILS_FORMATNETMHCPAN { tag "$meta.id" label 'process_single' container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? - 'docker://ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.4.0': - 'ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.4.0' }" + 'docker://ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.6.0': + 'ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.6.0' }" input: tuple val(meta), path(netmhcPanOutput) diff --git a/modules/msk/neoantigenutils/formatnetmhcpan/resources/usr/bin/format_netmhcpan_output.py b/modules/msk/neoantigenutils/formatnetmhcpan/resources/usr/bin/format_netmhcpan_output.py deleted file mode 100755 index 1c773de6..00000000 --- a/modules/msk/neoantigenutils/formatnetmhcpan/resources/usr/bin/format_netmhcpan_output.py +++ /dev/null @@ -1,116 +0,0 @@ -#!/usr/bin/env python3 - -import argparse - -VERSION = 1.0 - -PAN_HEADER = [ - "pos", - "MHC", - "peptide", - "core", - "OF", - "Gp", - "Gl", - "Ip", - "Il", - "icore", - "Identity", - "score_el", - "rank_el", - "score_ba", - "rank_ba", - "affinity", -] -STAB_PAN_HEADER = [ - "pos", - "HLA", - "peptide", - "Identity", - "Pred", - "Thalf(h)", - "%Rank_Stab", -] -NETMHC3_HEADER = ["pos", "peptide", "score_el", "affinity", "Identity", "MHC"] - - -def netMHCpan_out_reformat(netMHCoutput, mut, stab, netmhc3, prefix): - file_li = [] - stab_prefix = "" - type_prefix = "WT" - pan_prefix = "pan" - if stab: - stab_prefix = "stab" - if mut: - type_prefix = "MUT" - if netmhc3: - pan_prefix = "" - outfilename = "{}_netmhc{}{}.output.{}.tsv".format( - prefix, stab_prefix, pan_prefix, type_prefix - ) - with open(netMHCoutput, "r") as file: - # data = file.read() - for line in file: - # Remove leading whitespace - line = line.lstrip() - # Check if the line starts with a digit - if line == "": - pass - elif line[0].isdigit(): - # Print or process the line as needed - match = ( - line.strip() - .replace(" <= WB", "") - .replace(" <= SB", "") - .replace(" WB ", " ") - .replace(" SB ", " ") - ) # strip to remove leading/trailing whitespace - splititem = match.split() - tab_separated_line = "\t".join(splititem) - file_li.append(tab_separated_line) - if stab: - header = "\t".join(STAB_PAN_HEADER) + "\n" - elif netmhc3: - header = "\t".join(NETMHC3_HEADER) + "\n" - else: - header = "\t".join(PAN_HEADER) + "\n" - with open(outfilename, "w") as file: - file.writelines(header) - for item in file_li: - file.writelines(item) - file.writelines("\n") - - -def parse_args(): - parser = argparse.ArgumentParser(description="Process input files and parameters") - parser.add_argument( - "--netMHCpan_output", required=True, help="Path to netMHC output" - ) - parser.add_argument("--type_MUT", action="store_true", help="Output is a MUT type") - parser.add_argument( - "--from_STAB", - action="store_true", - help="Output is from netmhcstab", - ) - parser.add_argument( - "--from_NETMHC3", - action="store_true", - help="Output is from the older netmhc version 3.4", - ) - parser.add_argument("--id", required=True, help="Prefix to label the output") - parser.add_argument( - "-v", "--version", action="version", version="%(prog)s {}".format(VERSION) - ) - - return parser.parse_args() - - -def main(args): - netMHCpan_out_reformat( - args.netMHCpan_output, args.type_MUT, args.from_STAB, args.from_NETMHC3, args.id - ) - - -if __name__ == "__main__": - args = parse_args() - main(args) diff --git a/modules/msk/neoantigenutils/generatehlastring/main.nf b/modules/msk/neoantigenutils/generatehlastring/main.nf index e1c7ad2a..19914a71 100644 --- a/modules/msk/neoantigenutils/generatehlastring/main.nf +++ b/modules/msk/neoantigenutils/generatehlastring/main.nf @@ -2,8 +2,8 @@ process NEOANTIGENUTILS_GENERATEHLASTRING { tag "$meta.id" label 'process_single' container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? - 'docker://ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.4.0': - 'ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.4.0' }" + 'docker://ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.6.0': + 'ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.6.0' }" input: tuple val(meta), path(inputHLA) diff --git a/modules/msk/neoantigenutils/generatehlastring/resources/usr/bin/generateHLAString.sh b/modules/msk/neoantigenutils/generatehlastring/resources/usr/bin/generateHLAString.sh deleted file mode 100755 index 262494fc..00000000 --- a/modules/msk/neoantigenutils/generatehlastring/resources/usr/bin/generateHLAString.sh +++ /dev/null @@ -1,39 +0,0 @@ -#!/bin/bash - -### Versioning - -VERSION=1.0.0 - -get_help() { echo "USAGE: generateHLASTRING.sh -f [HLA_FILE]"; exit 0; } -get_version() { echo $VERSION; exit 0; } - -while (( "$#" )); do - case $1 in - -h|--help) get_help ;; - -f) file=$2; shift ;; - -v) get_version ;; - *) get_help ;; - esac - shift -done - -cat $file | tr "\t" "\n" | grep -v "HLA" | tr "\n" "," > massaged.winners.hla.txt - -input_string=`head -n 1 massaged.winners.hla.txt` - -IFS=',' read -ra items <<< "$input_string" - -for item in "${items[@]}"; do - - # Append the transformed item to the output string - truncated_value=$(echo "$item" | cut -c 1-11) - - # Replace the first '_', the next '_', and remaining '_' with '-', '*', and ':', respectively - modified_value=$(echo "$truncated_value" | tr '[:lower:]' '[:upper:]' | sed 's/_/-/; s/_//; s/_/:/g') - output_hla+=",$modified_value" - -done - -# Remove leading comma -output_hla="${output_hla:1}" -echo $output_hla diff --git a/modules/msk/neoantigenutils/neoantigeninput/main.nf b/modules/msk/neoantigenutils/neoantigeninput/main.nf index 08454756..1e03d58c 100644 --- a/modules/msk/neoantigenutils/neoantigeninput/main.nf +++ b/modules/msk/neoantigenutils/neoantigeninput/main.nf @@ -1,7 +1,7 @@ process NEOANTIGENUTILS_NEOANTIGENINPUT { tag "$meta.id" label 'process_medium' - container "ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.4.0" + container "ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.6.0" input: tuple val(meta), path(inputMaf), path(inputBedpe, arity: '0..*'), path(hlaFile) diff --git a/modules/msk/neoantigenutils/neoantigeninput/resources/usr/bin/generate_input.py b/modules/msk/neoantigenutils/neoantigeninput/resources/usr/bin/generate_input.py deleted file mode 100755 index 35323413..00000000 --- a/modules/msk/neoantigenutils/neoantigeninput/resources/usr/bin/generate_input.py +++ /dev/null @@ -1,1160 +0,0 @@ -#!/usr/bin/env python3 - -import json -import pandas as pd -import argparse -import os -import Bio -from Bio import pairwise2 -from Bio.pairwise2 import format_alignment -import numpy as np -from pyensembl.genome import Genome -from pyensembl import EnsemblRelease - -VERSION = 1.9 - - -def main(args): - - def makeChild(subTree, start): - if start: - subTree = 0 - - newsubtree = { - "clone_id": int(subTree), - "clone_mutations": [], - "children": [], - "X": 0, - "x": 0, - "new_x": 0, - } - - if str(subTree) in trees[tree]["structure"]: - for item in trees[tree]["structure"][str(subTree)]: - - child_dict = makeChild(item, False) - - newsubtree["children"].append(child_dict) - - try: - ssmli = [] - if start: - pass - else: - for ssm in treefile["mut_assignments"][str(subTree)]["ssms"]: - ssmli.append( - chrom_pos_dict[mut_data["ssms"][ssm]["name"]]["id"] - ) - newsubtree["clone_mutations"] = ssmli - newsubtree["X"] = trees[tree]["populations"][str(subTree)][ - "cellular_prevalence" - ][0] - newsubtree["x"] = trees[tree]["populations"][str(subTree)][ - "cellular_prevalence" - ][0] - newsubtree["new_x"] = 0.0 - except Exception as e: - print("Error in adding new subtree. Error not in base case**") - print(subTree) - print(e) - pass - - return newsubtree - - else: - # Base Case - # make childrendict and return it - ssmli = [] - - for ssm in treefile["mut_assignments"][str(subTree)]["ssms"]: - try: - ssmli.append(chrom_pos_dict[mut_data["ssms"][ssm]["name"]]["id"]) - except Exception as e: - print( - "Error in appending to mutation list. Error in base case appending ssm to ssmli" - ) - print(e) - # print(str(subTree)) - pass - - try: - newsubtree["clone_mutations"] = ssmli - newsubtree["X"] = trees[tree]["populations"][str(subTree)][ - "cellular_prevalence" - ][0] - newsubtree["x"] = trees[tree]["populations"][str(subTree)][ - "cellular_prevalence" - ][0] - newsubtree["new_x"] = 0.0 - except Exception as e: - print("Error in adding new subtree. Error in base case") - print(e) - pass - return newsubtree - - with open(args.summary_file, "r") as f: - # Load the JSON data into a dictionary - summ_data = json.load(f) - - with open(args.mutation_file, "r") as f: - # Load the JSON data into a dictionary - mut_data = json.load(f) - - chrom_pos_dict = {} # Just used for mapping right now - mutation_list = [] # Used as the output for mutations - mutation_dict = ( - {} - ) # Used for matching mutation without the subsititution information from netMHCpan to phyloWGS output - gene_dict = {} - - mafdf = pd.read_csv(args.maf_file, delimiter="\t") - - for index, row in mafdf.iterrows(): - if ( - # We - row["Variant_Type"] == "SNP" - or row["Variant_Type"] == "DEL" - or row["Variant_Type"] == "INS" - or row["Variant_Type"] == "DNP" - or row["Variant_Type"] == "TNP" - ): - if row["Variant_Classification"] == "Missense_Mutation": - missense = 1 - - else: - missense = 0 - - if ( - row["Variant_Type"] == "SNP" - or row["Variant_Type"] == "DNP" - or row["Variant_Type"] == "TNP" - ): - chrom_pos_dict[ - str(row["Chromosome"]) - + "_" - + str(row["Start_Position"]) - + "_" - + row["Reference_Allele"] - + "_" - + row["Tumor_Seq_Allele2"] - ] = { - "id": str(row["Chromosome"]) - + "_" - + str(row["Start_Position"]) - + "_" - + row["Reference_Allele"] - + "_" - + row["Tumor_Seq_Allele2"], - "gene": row["Hugo_Symbol"], - "missense": missense, - "transcript": row["Feature"] - } - - mutation_list.append( - { - "id": str(row["Chromosome"]) - + "_" - + str(row["Start_Position"]) - + "_" - + row["Reference_Allele"] - + "_" - + row["Tumor_Seq_Allele2"], - "gene": row["Hugo_Symbol"], - "missense": missense, - } - ) - - mutation_dict[makeID(row)] = ( - str(row["Chromosome"]) - + "_" - + str(row["Start_Position"]) - + "_" - + row["Reference_Allele"] - + "_" - + row["Tumor_Seq_Allele2"] - ) - - elif row["Variant_Type"] == "DEL": - chrom_pos_dict[ - str(row["Chromosome"]) - + "_" - + str(row["Start_Position"] - 1) - + "_" - + row["Reference_Allele"] - + "_" - + "D" - ] = { - "id": str(row["Chromosome"]) - + "_" - + str(row["Start_Position"]) - + "_" - + row["Reference_Allele"] - + "_" - + "D", - "gene": row["Hugo_Symbol"], - "missense": missense, - "transcript": row["Feature"] - } - - mutation_list.append( - { - "id": str(row["Chromosome"]) - + "_" - + str(row["Start_Position"]) - + "_" - + row["Reference_Allele"] - + "_" - + "D", - "gene": row["Hugo_Symbol"], - "missense": missense, - } - ) - mutation_dict[makeID(row)] = ( - str(row["Chromosome"]) - + "_" - + str(row["Start_Position"] - 1) - + "_" - + row["Reference_Allele"] - + "_" - + "D" - ) - - elif row["Variant_Type"] == "INS": - # print( - # str(row["Chromosome"]) - # + "_" - # + str(row["Start_Position"]) - # + "_" - # + "I" - # + "_" - # + row["Tumor_Seq_Allele2"] - # ) - chrom_pos_dict[ - str(row["Chromosome"]) - + "_" - + str(row["Start_Position"]) - + "_" - + "I" - + "_" - + row["Tumor_Seq_Allele2"] - ] = { - "id": str(row["Chromosome"]) - + "_" - + str(row["Start_Position"]) - + "_" - + "I" - + "_" - + row["Tumor_Seq_Allele2"], - "gene": row["Hugo_Symbol"], - "missense": missense, - "transcript": row["Feature"] - } - - mutation_list.append( - { - "id": str(row["Chromosome"]) - + "_" - + str(row["Start_Position"]) - + "_" - + "I" - + "_" - + row["Tumor_Seq_Allele2"], - "gene": row["Hugo_Symbol"], - "missense": missense, - } - ) - mutation_dict[makeID(row)] = ( - str(row["Chromosome"]) - + "_" - + str(row["Start_Position"]) - + "_" - + "I" - + "_" - + row["Tumor_Seq_Allele2"] - ) - - outer_dict = {"id": args.id, "sample_trees": []} - - trees = summ_data["trees"] - - for tree in trees: - - inner_sample_tree_dict = {"topology": [], "score": trees[tree]["llh"]} - with open("./" + args.tree_directory + "/" + str(tree) + ".json", "r") as f: - # Load the JSON data into a dictionary - treefile = json.load(f) - - bigtree = makeChild(tree, True) - - inner_sample_tree_dict["topology"] = bigtree - - outer_dict["sample_trees"].append(inner_sample_tree_dict) - - outer_dict["mutations"] = mutation_list - - # TODO format HLA_gene input data, depending on format inputted. They should look like this A*02:01 - # this will be setup for polysolver winners output - def convert_polysolver_hla(polyHLA): - allele = polyHLA[4] - shortHLA = polyHLA.split("_")[2:4] - return allele.upper() + "*" + shortHLA[0] + ":" + shortHLA[1] - - HLA_gene_li = [] - with open(args.HLA_genes, "r") as f: - for line in f: - line = line.split("\t") - - HLA_gene_li.append(convert_polysolver_hla(line[1])) - HLA_gene_li.append(convert_polysolver_hla(line[2])) - - outer_dict["HLA_genes"] = HLA_gene_li - - if args.patient_data_file: - # TODO format optional input data, depending on format of inputted file. I was imagining a tsv, but can be anything - status = 0 - OS_tmp = 0 - PFS = 0 - - outer_dict["status"] = status - outer_dict["OS"] = OS_tmp - outer_dict["PFS"] = PFS - else: - outer_dict["status"] = 0 - outer_dict["OS"] = 0 - outer_dict["PFS"] = 0 - - outer_dict["id"] = args.id - outer_dict["patient"] = args.patient_id - outer_dict["cohort"] = args.cohort - - outer_dict["neoantigens"] = [] - - if args.bedpe_file: - bedpe_list, bedpe_dict = bedpe_load(args.bedpe_file) - - bedpe_match_dict = {} - - ensembl = ensembl_load(args.release, args.gtf_file, args.cdna_file) - neoantigen_mut_in = pd.read_csv(args.netMHCpan_MUT_input, sep="\t") - neoantigen_WT_in = pd.read_csv(args.netMHCpan_WT_input, sep="\t") - - def find_first_difference_index(str1, str2): - min_length = min(len(str1), len(str2)) - for i in range(min_length): - if str1[i] != str2[i]: - return i - # If no difference found in the common length, return the length of the shorter string - return min_length - - WTdict = {} - SVWTdict = {} - for index_WT, row_WT in neoantigen_WT_in.iterrows(): - no_positon_ID = "" - id = "" - wtsvid = "" - row_WT_identity = trim_id(row_WT["Identity"]) - IDsplit = row_WT_identity.split("_") - print(row_WT_identity) - if len(IDsplit[0]) < 3: - # it is from neoSV - IDsplit = row_WT_identity.split("_") - wtsvid = ( - IDsplit[0] - + IDsplit[1][0:7] - + "_" - + str(len(row_WT["peptide"])) - + "_" - + row_WT["MHC"].split("-")[1].replace(":", "").replace("*", "") - + "_" - + str(row_WT["pos"]) - ) - no_positon_ID = ( - IDsplit[0] - + "_" - + IDsplit[1][0:7] - + "_" - + str(len(row_WT["peptide"])) - + "_" - + row_WT["MHC"].split("-")[1].replace(":", "").replace("*", "") - ) - WTdict[wtsvid] = { - "affinity": row_WT["affinity"], - "peptide": row_WT["peptide"], - } - id = wtsvid - if no_positon_ID not in WTdict: - WTdict[no_positon_ID] = { - "peptides": { - row_WT["peptide"]: id - }, # This is a dict so we can match the peptide with the actual ID later - "affinity": row_WT["affinity"], - } - - else: - WTdict[no_positon_ID]["peptides"][row_WT["peptide"]] = id - - else: - id = ( - row_WT_identity[:-2] - + "_" - + str(len(row_WT["peptide"])) - + "_" - + row_WT["MHC"].split("-")[1].replace(":", "").replace("*", "") - + "_" - + str(row_WT["pos"]) - ) - - no_positon_ID = ( - row_WT_identity[:-2] - + "_" - + str(len(row_WT["peptide"])) - + "_" - + row_WT["MHC"].split("-")[1].replace(":", "").replace("*", "") - ) - WTdict[id] = {"affinity": row_WT["affinity"], "peptide": row_WT["peptide"]} - - # This is used as last resort for the matching. We will preferentially find the peptide matching in length as well as POS. Worst case we will default to the WT pos 0 - if no_positon_ID not in WTdict: - WTdict[no_positon_ID] = { - "peptides": { - row_WT["peptide"]: id - }, # This is a dict so we can match the peptide with the ID later - "affinity": row_WT["affinity"], - } - - else: - WTdict[no_positon_ID]["peptides"][row_WT["peptide"]] = id - - def find_most_similar_string(target, strings): - max_score = -1 - max_score2 = -2 - most_similar_string = None - most_similar_string2 = None - first_AA_same = None - first_AA_same_score = -1 - len_target = len(target) - for s in strings: - if len(s) == len_target: - alignments = pairwise2.align.globalxx(target, s) - score = alignments[0][2] # The third element is the score - - if score > max_score2: - - if score > max_score: - max_score2 = max_score - most_similar_string2 = most_similar_string - max_score = score - most_similar_string = s - - else: - max_score2 = score - most_similar_string2 = s - - if target[0] == s[0]: - if score > first_AA_same_score: - first_AA_same = s - first_AA_same_score = score - - return ( - most_similar_string, - most_similar_string2, - first_AA_same, - first_AA_same_score, - max_score, - ) - NMD_dict = {} - - for index_mut, row_mut in neoantigen_mut_in.iterrows(): - row_MUT_identity = trim_id(row_mut["Identity"]) - IDsplit = row_MUT_identity.split("_") - SV = False - if row_mut["affinity"] < float(args.kD_cutoff): - peplen = len(row_mut["peptide"]) - matchfound = False - frameshift= False - print(row_MUT_identity) - if IDsplit[1][0] == "S" and IDsplit[1][1] != "p": - # If it is a silent mutation. Silent mutations can either be S or SY. These include intron mutations. Splices can be Sp - continue - if row_MUT_identity.count("_") == 1: - # its an SV - SV = True - WTid = ( - IDsplit[0] - + IDsplit[1][0:8] - + "_" - + str(len(row_mut["peptide"])) - + "_" - + row_mut["MHC"].split("-")[1].replace(":", "").replace("*", "") - + "_" - + str(row_mut["pos"]) - ) - no_positon_ID = ( - IDsplit[0] - + "_" - + IDsplit[1][0:8] - + "_" - + str(len(row_mut["peptide"])) - + "_" - + row_mut["MHC"].split("-")[1].replace(":", "").replace("*", "") - ) - # this part makes the dict that matches this to the bedpe - bedpe_match_dict[row_MUT_identity] = ( - IDsplit[0] + "_" + IDsplit[1][0:4] - ) - else: - # first find match in WT - WTid = ( - row_MUT_identity[:-2] - + "_" - + str(peplen) - + "_" - + row_mut["MHC"].split("-")[1].replace(":", "").replace("*", "") - + "_" - + str(row_mut["pos"]) - ) - no_positon_ID = ( - row_MUT_identity[:-2] - + "_" - + str(peplen) - + "_" - + row_mut["MHC"].split("-")[1].replace(":", "").replace("*", "") - ) - if ( - ("M" == IDsplit[1][0] and "Sp" not in row_MUT_identity) - or SV == False - ): - # match - if ( - (WTid in WTdict) - and IDsplit[1][0] not in ("I", "F") - ): - #This block takes care of Missense mutations caused by polymorphisims - matchfound = True - best_pepmatch = WTdict[WTid]["peptide"] - - else: - # Here we take care of INDELS and everything else - - if IDsplit[1].startswith("Fi") or IDsplit[1].startswith("Fd"): - frameshift = True - ( - best_pepmatch, - best_pepmatch2, - first_AA_same, - first_AA_same_score, - match_score, - ) = find_most_similar_string( - row_mut["peptide"], list(WTdict[no_positon_ID]["peptides"].keys()) - ) - if ( - best_pepmatch == row_mut["peptide"] - or best_pepmatch2 == row_mut["peptide"] - ): - # it seems this can happen where the row_mut is actually the canonical sequence. - # In this case we don't want to report the peptide as a neoantigen, its not neo - continue - - elif best_pepmatch2 is not None and ( - ( - best_pepmatch[0] != row_mut["peptide"][0] - and best_pepmatch2[0] == row_mut["peptide"][0] - ) or ( - best_pepmatch[-1] != row_mut["peptide"][-1] - and best_pepmatch2[-1] == row_mut["peptide"][-1] - ) - ): - # We should preferentially match the first AA if we can. Sometimes the pairwise alignment isnt the best at this so we do a little check here. - # It will also do this when the last AA of the best match doesnt match but the last A of the second best match does - best_pepmatch = best_pepmatch2 - - WTid = WTdict[no_positon_ID]["peptides"][best_pepmatch] - matchfound = True - - if matchfound == True and best_pepmatch != row_mut["peptide"]: - mut_pos = ( - find_first_difference_index( - row_mut["peptide"], best_pepmatch # WTdict[WTid]["peptide"] - ) - + 1 - ) - - - chrom, pos = mutation_dict[row_MUT_identity].split("_")[0:2] - - if frameshift: - mut_pos = "Frameshifted peptide" - num_windows = len(list(WTdict[no_positon_ID]["peptides"].keys())) - num_windows_li = len(list(WTdict[no_positon_ID]["peptides"].values())) - - print("FRAMESHIFT") - - if no_positon_ID in NMD_dict: - #NMD must only be calculated once per mutation - pass - else: - split_mutation_dict_ID = mutation_dict[row_MUT_identity].split("_") - if split_mutation_dict_ID[2] == "I": - len_indel = len(split_mutation_dict_ID[3]) - elif split_mutation_dict_ID[3] == "D": - len_indel = 0 - len(split_mutation_dict_ID[2]) - else: - len_indel = 0 - transcriptID = chrom_pos_dict[mutation_dict[row_MUT_identity]]["transcript"] - NMD_dict[no_positon_ID] = determine_NMD(chrom, pos,num_windows,len_indel,ensembl,transcriptID) - - else: - NMD_dict[no_positon_ID] = "False" - - if SV: - neo_dict = { - "id": row_MUT_identity - + "_" - + str(peplen) - + "_" - + str(row_mut["pos"]) - + "_" - + row_mut["MHC"].split("-")[1].replace(":", "").replace("*", "") - , - "mutation_id": bedpe_dict[ - bedpe_match_dict[row_MUT_identity] - ].id, - "Gene": chrom_pos_dict[mutation_dict[row_MUT_identity]]["gene"], - "HLA_gene_id": row_mut["MHC"], - "sequence": row_mut["peptide"], - "WT_sequence": best_pepmatch, # WTdict[WTid]["peptide"], - "mutated_position": mut_pos, - "Kd": float(row_mut["affinity"]), - "KdWT": float(WTdict[WTid]["affinity"]), - "NMD" : NMD_dict[no_positon_ID] - } - else: - neo_dict = { - "id": row_MUT_identity - + "_" - + str(peplen) - + "_" - + str(row_mut["pos"]) - + "_" - + row_mut["MHC"].split("-")[1].replace(":", "").replace("*", ""), - "mutation_id": mutation_dict[row_MUT_identity], - "Gene": chrom_pos_dict[mutation_dict[row_MUT_identity]]["gene"], - "HLA_gene_id": row_mut["MHC"], - "sequence": row_mut["peptide"], - "WT_sequence": best_pepmatch, # WTdict[WTid]["peptide"], - "mutated_position": mut_pos, - "Kd": float(row_mut["affinity"]), - "KdWT": float(WTdict[WTid]["affinity"]), - "NMD" : NMD_dict[no_positon_ID] - } - outer_dict["neoantigens"].append(neo_dict) - - outjson = args.patient_id + "_" + args.id + "_" + "input.json" - with open(outjson, "w") as tstout: - json.dump(outer_dict, tstout, indent=1) - -# Sometimes the id is set as .*_M_1 and we want to make sure its _M, otherwise it will not match -def trim_id(id_string): - if "_M_" in id_string: - return id_string.partition("_M_")[0]+"_M" - elif "_V_" in id_string: - return id_string.partition("_V_")[0]+"_V" - else: - return id_string - - -def makeID(maf_row): - ##ENCODING FASTA ID FOR USE IN MATCHING LATER - ALPHABET = [ - "A", - "B", - "C", - "D", - "E", - "F", - "G", - "H", - "I", - "J", - "K", - "L", - "M", - "N", - "O", - "P", - "Q", - "R", - "S", - "T", - "U", - "V", - "W", - "X", - "Y", - "Z", - ] - - variant_type_map = { - "missense_mutation": "M", - "nonsense_mutation": "X", - "silent_mutation": "S", - "silent": "S", - "frame_shift_ins": "Fi", - "frame_shift_del": "Fd", - "in_frame_ins": "If", - "in_frame_del": "Id", - "splice_site": "Sp", - } - - position = int(str(maf_row["Start_Position"])[0:2]) - - if position < 26: - encoded_start = ALPHABET[position] - elif position < 100: - encoded_start = ALPHABET[position // 4] - - position = int(str(maf_row["Start_Position"])[-2:]) - - if position < 26: - encoded_end = ALPHABET[position] - elif position < 100: - encoded_end = ALPHABET[position // 4] - - sum_remaining = sum(int(d) for d in str(maf_row["Start_Position"])[2:-2]) - - encoded_position = encoded_start + ALPHABET[sum_remaining % 26] + encoded_end - - if maf_row["Tumor_Seq_Allele2"] == "-": - # handles deletion - if len(maf_row["Reference_Allele"]) > 3: - Allele2code = maf_row["Reference_Allele"][0:3] - else: - Allele2code = maf_row["Reference_Allele"] - - elif len(maf_row["Tumor_Seq_Allele2"]) > 1: - # handles INS and DNP - if len(maf_row["Tumor_Seq_Allele2"]) > 3: - Allele2code = maf_row["Tumor_Seq_Allele2"][0:3] - else: - Allele2code = maf_row["Tumor_Seq_Allele2"] - - else: - # SNPs - Allele2code = maf_row["Tumor_Seq_Allele2"] - - if maf_row["Variant_Classification"].lower() in variant_type_map: - identifier_key = ( - str(maf_row["Chromosome"]) - + encoded_position - + "_" - + variant_type_map[maf_row["Variant_Classification"].lower()] - + Allele2code - + "_M" # This indicates mutated. It is added in the generateMutFasta script as well but not in this function. - ) - else: - identifier_key = ( - str(maf_row["Chromosome"]) - + encoded_position - + "_" - + "SY" - + Allele2code - + "_M" - ) - return identifier_key - - -class VariantCallingFormat(object): - """ - Class for storing SV information in VCF format, - all components are in string format - """ - - def __init__(self, chrom, pos, ref, alt): - self.chrom = chrom - self.pos = pos - self.ref = ref - self.alt = alt - - def __str__(self): - return "%s(chrom = %s, pos = %s, ref = %s, alt = %s)" % ( - self.__class__.__name__, - self.chrom, - self.pos, - self.ref, - self.alt, - ) - - def __repr__(self): - return "%s(%s, %s, %s, %s)" % ( - self.__class__.__name__, - self.chrom, - self.pos, - self.ref, - self.alt, - ) - - -class BedpeFormat(object): - """ - Class for storing SV information in BEDPE format, - all components are in string format - """ - - def __init__(self, chrom1, pos1, strand1, chrom2, pos2, strand2, id): - self.chrom1 = chrom1 - self.pos1 = pos1 - self.strand1 = strand1 - self.chrom2 = chrom2 - self.pos2 = pos2 - self.strand2 = strand2 - self.id = id - - def __str__(self): - return ( - "%s(chrom1 = %s, pos1 = %s, strand1 = %s, chrom2 = %s, pos2 = %s, strand2 = %s, id = %s)" - % ( - self.__class__.__name__, - self.chrom1, - self.pos1, - self.strand1, - self.chrom2, - self.pos2, - self.strand2, - self.id, - ) - ) - - def __repr__(self): - return "%s(%s, %s, %s, %s, %s, %s, %s)" % ( - self.__class__.__name__, - self.chrom1, - self.pos1, - self.strand1, - self.chrom2, - self.pos2, - self.strand2, - self.id, - ) - - -def bedpe_load(filepath): - """ - :param filepath: the absolute path of a BEDPE file - :return: a list of BEDPE objects - """ - bedpe_list = [] - bedpedict = {} - filename = os.path.basename(filepath) - line_num = 0 - print("Loading SVs from {0}.".format(filename)) - with open(filepath, "r") as f: - header = next(f) - header = header.rstrip().split("\t") - for line in f: - line_num += 1 - tmpline = line.rstrip().split("\t") - chrom1 = tmpline[header.index("chrom1")].replace("chr", "") - pos1 = tmpline[header.index("start1")] - chrom2 = tmpline[header.index("chrom2")].replace("chr", "") - pos2 = tmpline[header.index("start2")] - strand1 = tmpline[header.index("strand1")] - strand2 = tmpline[header.index("strand2")] - svclass = tmpline[header.index("svclass")] - sv_bedpe_id = tmpline[header.index("sv_id")] - custom_id = makeID_bedpe(chrom1, pos1, svclass) - bedpe = BedpeFormat( - chrom1, pos1, strand1, chrom2, pos2, strand2, sv_bedpe_id - ) - bedpe_list.append(bedpe) - bedpedict[custom_id] = bedpe - - return bedpe_list, bedpedict - - -def makeID_bedpe(chrom1, pos1, svclass): - ##ENCODING FASTA ID FOR USE IN MATCHING LATER - ALPHABET = [ - "A", - "B", - "C", - "D", - "E", - "F", - "G", - "H", - "I", - "J", - "K", - "L", - "M", - "N", - "O", - "P", - "Q", - "R", - "S", - "T", - "U", - "V", - "W", - "X", - "Y", - "Z", - ] - - position = int(str(pos1)[0:2]) - - if position < 26: - encoded_start = ALPHABET[position] - elif position < 100: - encoded_start = ALPHABET[position // 4] - - position = int(str(pos1)[-2:]) - - if position < 26: - encoded_end = ALPHABET[position] - elif position < 100: - encoded_end = ALPHABET[position // 4] - sum_remaining = sum(int(d) for d in str(pos1)[2:-2]) - - encoded_position = encoded_start + ALPHABET[sum_remaining % 26] + encoded_end - - identifier_key = ( - str(chrom1) - + "_" - + encoded_position - + "V" # This indicates structural variant. It is added in the generateMutFasta script as well but not in this function. - ) - - return identifier_key - - -def get_exon_range(transcript): - """ - :param transcript: transcript instance in pyensembl - :return: exon intervals of this transcript - from 5' to 3', exon 1, exon 2, ... - [start, end], start < end - """ - exon_ranges = [] - for exon in transcript.exons: - exon_ranges.append((exon.start, exon.end)) - return exon_ranges - -def get_longest_transcript(transcripts): - """ - :param transcripts: a list of Transcript(pyensembl) instances - :return: the longest transcript - """ - transcripts = sorted(transcripts, key=lambda t: t.end-t.start) - transcript = transcripts[-1] - return transcript - - -def get_transcript(chrom, pos, ensembl, complete=True): - """ - :param chrom: chromosome with no chr - :param pos: position of mutation - :param ensembl: Genome instance in pyensembl - :param complete: only consider complete transcripts - :return: firstly return the longest complete transcript, - if there is no complete transcript and - complete = False, return the longest transcript - """ - transcripts = ensembl.transcripts_at_locus(contig=str(chrom), position=int(pos)) - transcripts_comp = [transcript for transcript in transcripts if transcript.complete] - if transcripts_comp: - return get_longest_transcript(transcripts_comp) - else: - if complete: - return None - else: - if transcripts: - return get_longest_transcript(transcripts) - else: - return None - - -def ensembl_load(release, gtf_file, cdna_file): - """ - :param release: the release number in EMSEMBL, could be custom - :param gtf_file: the path of gtf file if release == custom - :param cdna_file: the path of cdna file if release == custom - :param cache_dir: directory for pyensembl downloading - :return: a Genome class in pyensembl - """ - # if release != 'custom': - # print("doing a new one") - # ensembl = EnsemblRelease(int(release)) - # ensembl.download() - # ensembl.index() - - # else: - ensembl = Genome(gtf_path_or_url=gtf_file, - transcript_fasta_paths_or_urls=cdna_file, - reference_name='User-defined', - annotation_name='User-defined') - ensembl.index() - return ensembl - - - -def get_exons_from_transcriptID(transcriptid, ensembl, complete=True): - """ - :param transcriptid: transcriptid from MAF col - :param ensembl: Genome instance in pyensembl - :param complete: only consider complete transcripts - :return: a list of tuples of exon ranges - """ - - transcripts= ensembl.exon_ids_of_transcript_id(transcriptid) - - exon_ranges = [] - for exonid in transcripts: - exon = ensembl.exon_by_id(exonid) - exon_ranges.append((exon.start, exon.end)) - - return exon_ranges - - -def determine_NMD(chrom, pos,num_windows,len_indel, ensembl, transcriptID=None): - """ - :param chrom: chromosome where alteration takes place - :param pos: position where alteration takes place - :pos transcriptID: transcriptID from the MAF. If it isnt annotated, then use longest transcript - :num_windows: number of windows created by the peptide - :len_indel: Length of the indel. Negative if del, positive if ins - :param gtf_file: the path of gtf file if release == custom - :param cdna_file: the path of cdna file if release == custom - :return: NMD value - """ - - if transcriptID == None: - #do these if maf was not annotated - transcript = get_transcript(chrom, pos, ensembl) - exon_ranges = get_exon_range(transcript) - else: - exon_ranges = get_exons_from_transcriptID(transcriptID,ensembl) - - NMD = "False" - - PTC_exon = None - pos = int(pos) + 1 - for i in range(0,len(exon_ranges)): - if pos>=exon_ranges[i][0] and pos<=exon_ranges[i][1]: - exon_ranges_dist = [exon_ranges[p][1]-exon_ranges[p][0] for p in range(0,len(exon_ranges))] - mut_to_stop_dist = (num_windows*3)+len_indel+1 - - for d in range(i,len(exon_ranges_dist)): - if exon_ranges_dist[d] == exon_ranges_dist[i]: - dist = exon_ranges_dist[d] - (exon_ranges[i][1]-pos) - else: - dist = exon_ranges_dist[d] - - if dist - mut_to_stop_dist >= 0: - PTC_exon = exon_ranges[d] - PTC_pos = exon_ranges[d][0] + mut_to_stop_dist - break - elif len(exon_ranges_dist)-1 == d and dist - mut_to_stop_dist < 0: - PTC_exon = exon_ranges[d] - PTC_pos = exon_ranges[d][0] + mut_to_stop_dist - else: - mut_to_stop_dist = mut_to_stop_dist - dist - - if PTC_exon is not None and PTC_exon == exon_ranges[-1]: - # "on the last exon" - NMD = "Last Exon" - else: - if exon_ranges[0][0] - PTC_pos < 150: - # less than 150 nt away from the start exon - NMD = "Start-proximal" - else: - if (PTC_exon[1] - PTC_exon[0]) > 407: - # in a long exon with more than 407 nt - NMD = "Long Exon" - else: - # it is in the last 50 nt of the penultimate exon - if PTC_exon == exon_ranges[-2] and (exon_ranges[-2][0] - PTC_pos) < 50 : - NMD = "50nt Rule" - else: - NMD = "Trigger NMD" - - return NMD - - - -def parse_args(): - parser = argparse.ArgumentParser(description="Process input files and parameters") - parser.add_argument("--maf_file", required=True, help="Path to the MAF file") - parser.add_argument("--bedpe_file", required=False, help="Path to the bedpe file") - parser.add_argument( - "--summary_file", required=True, help="Path to the summary file" - ) - parser.add_argument( - "--mutation_file", required=True, help="Path to the mutation file" - ) - parser.add_argument( - "--tree_directory", - required=True, - help="Path to the tree directory containing json files", - ) - parser.add_argument('-r', '--release', dest='release', metavar='RELEASE', default='75', - help='Which reference (ENSEMBL release) you want to use. Ensembl releases that' - 'correspond to hg18/NCBI36, hg19/GRCh37, hg38/GRCh38 are 54, 75, 95.' - 'If your data are from other species(custom), please download the gtf ' - 'file and the cdna file from ENSEMBL website ftp://ftp.ensembl.org/pub' - ' and specify them using --gtf-file and --cdna-file.') - parser.add_argument('-gf', '--gtf-file', dest='gtf_file', metavar='GTF_FILE', default=None, - help='GTF file for the reference.') - parser.add_argument('-cf', '--cdna-file', dest='cdna_file', metavar='CDNA_FILE', default=None, - help='cDNA file for the reference.') - - parser.add_argument("--id", required=True, help="ID") - parser.add_argument("--patient_id", required=True, help="Patient ID") - parser.add_argument("--cohort", required=True, help="Cohort") - parser.add_argument( - "--HLA_genes", required=True, help="Path to the file containing HLA genes" - ) - parser.add_argument( - "--netMHCpan_MUT_input", - required=True, - help="Path to the file containing MUT netmhcpan results", - ) - parser.add_argument( - "--netMHCpan_WT_input", - required=True, - help="Path to the file containing WT netmhcpan results", - ) - - parser.add_argument( - "--patient_data_file", - help="Path to the optional file containing status, overall survival, and PFS", - ) - parser.add_argument( - "-v", "--version", action="version", version="%(prog)s {}".format(VERSION) - ) - - parser.add_argument( - "--kD_cutoff", default=500, help="Cutoff value for the kD, default is 500", - ) - - return parser.parse_args() - - -if __name__ == "__main__": - args = parse_args() - print("MAF File:", args.maf_file) - print("Summary File:", args.summary_file) - print("Mutation File:", args.mutation_file) - print("Tree directory:", args.tree_directory) - print("ID:", args.id) - print("Patient ID:", args.patient_id) - print("Cohort:", args.cohort) - print("HLA Genes File:", args.HLA_genes) - print("netMHCpan Files:", args.netMHCpan_MUT_input, args.netMHCpan_WT_input) - print("kD Cutoff Value:", args.kD_cutoff) - print("Ensembl files:", args.gtf_file, args.cdna_file) - if args.patient_data_file: - print("patient_data_file File:", args.patient_data_file) - - main(args) diff --git a/modules/msk/neosv/main.nf b/modules/msk/neosv/main.nf index 89746364..851cf7c3 100644 --- a/modules/msk/neosv/main.nf +++ b/modules/msk/neosv/main.nf @@ -4,8 +4,8 @@ process NEOSV { conda "${moduleDir}/environment.yml" container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? - 'docker://ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.4.0': - 'ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.4.0' }" + 'docker://ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.6.0': + 'ghcr.io/mskcc-omics-workflows/neoantigen-utils-base:1.6.0' }" input: tuple val(meta), path(inputBedpe), val(hlaString)