import itertools
import warnings
from pathlib import Path
from typing import List
import numpy as np
import pandas as pd
from immuneML.data_model.dataset.RepertoireDataset import RepertoireDataset
from immuneML.data_model.receptor.receptor_sequence.Chain import Chain
from immuneML.reports.ReportOutput import ReportOutput
from immuneML.reports.ReportResult import ReportResult
from immuneML.reports.encoding_reports.EncodingReport import EncodingReport
from immuneML.util.PathBuilder import PathBuilder
[docs]class Matches(EncodingReport):
"""
Reports the number of matches that were found when using one of the following encoders:
* :ref:`MatchedSequences` encoder
* :ref:`MatchedReceptors` encoder
* :ref:`MatchedRegex` encoder
Report results are:
* A table containing all matches, where the rows correspond to the Repertoires, and the
columns correspond to the objects to match (regular expressions or receptor sequences).
* The repertoire sizes (read frequencies and the number of unique sequences per repertoire), for each of the chains.
This can be used to calculate the percentage of matched sequences in a repertoire.
* When using :ref:`MatchedSequences` encoder or
:ref:`MatchedReceptors` encoder, tables describing
the chains and receptors (ids, chains, V and J genes and sequences).
* When using :ref:`MatchedReceptors` encoder or using
:ref:`MatchedRegex` encoder with chain pairs, tables describing
the paired matches (where a match was found in both chains) per repertoire.
YAML Specification:
.. indent with spaces
.. code-block:: yaml
my_match_report: Matches
"""
[docs] @classmethod
def build_object(cls, **kwargs):
return Matches(**kwargs)
def __init__(self, dataset: RepertoireDataset = None, result_path: Path = None, name: str = None, number_of_processes: int = 1):
super().__init__(dataset=dataset, result_path=result_path, name=name, number_of_processes=number_of_processes)
def _generate(self) -> ReportResult:
PathBuilder.build(self.result_path)
return self._write_reports()
def _write_reports(self) -> ReportResult:
all_matches_table = self._write_match_table()
repertoire_sizes = self._write_repertoire_sizes()
output_tables = [all_matches_table, repertoire_sizes]
if self.dataset.encoded_data.encoding == "MatchedSequencesEncoder":
output_tables += self._write_sequence_info(self.result_path / "sequence_info")
else:
if len(self.dataset.encoded_data.feature_annotations["chain"].unique()) == 2:
output_tables += self._write_paired_matches(self.result_path / "paired_matches")
if self.dataset.encoded_data.encoding == "MatchedReceptorsEncoder":
output_tables += self._write_receptor_info(self.result_path / "receptor_info")
return ReportResult(self.name,
info=f"Reports the number of matches that were found when using {self.dataset.encoded_data.encoding}",
output_tables=output_tables)
def _write_match_table(self):
id_df = pd.DataFrame({"repertoire_id": self.dataset.encoded_data.example_ids})
label_df = pd.DataFrame(self.dataset.encoded_data.labels)
matches_df = pd.DataFrame(self.dataset.encoded_data.examples, columns=self.dataset.encoded_data.feature_names)
result_path = self.result_path / "complete_match_count_table.csv"
id_df.join(label_df).join(matches_df).to_csv(result_path, index=False)
return ReportOutput(result_path, "All matches")
def _write_paired_matches(self, paired_matches_path: Path) -> List[ReportOutput]:
PathBuilder.build(paired_matches_path)
report_outputs = []
for i in range(0, len(self.dataset.encoded_data.example_ids)): # todo don't mention subject in the name twice
file_name = "example_{}_".format(self.dataset.encoded_data.example_ids[i])
file_name += "_".join(["{label}_{value}".format(label=label, value=values[i]) for
label, values in self.dataset.encoded_data.labels.items()])
file_name += ".csv"
file_path = paired_matches_path / file_name
if self.dataset.encoded_data.encoding == "MatchedReceptorsEncoder":
self._write_paired_receptor_matches_for_repertoire(self.dataset.encoded_data.examples[i], file_path)
elif self.dataset.encoded_data.encoding == "MatchedRegexEncoder":
self._write_paired_regex_matches_for_repertoire(self.dataset.encoded_data.examples[i], file_path)
report_outputs.append(ReportOutput(file_path, f"Example {self.dataset.encoded_data.example_ids[i]} paired matches"))
return report_outputs
def _write_paired_receptor_matches_for_repertoire(self, matches, filename):
match_identifiers = []
match_values = []
for i in range(0, int(len(matches) / 2)):
first_match_idx = i * 2
second_match_idx = i * 2 + 1
if matches[first_match_idx] > 0 and matches[second_match_idx] > 0:
match_identifiers.append(self.dataset.encoded_data.feature_names[first_match_idx])
match_identifiers.append(self.dataset.encoded_data.feature_names[second_match_idx])
match_values.append(matches[first_match_idx])
match_values.append(matches[second_match_idx])
results_df = pd.DataFrame([match_values], columns=match_identifiers)
results_df.to_csv(filename, index=False)
def _write_paired_regex_matches_for_repertoire(self, matches, filename):
match_identifiers = []
match_values = []
annotation_df = self.dataset.encoded_data.feature_annotations
for receptor_id in sorted(set(annotation_df["receptor_id"])):
chain_ids = list(annotation_df.loc[annotation_df["receptor_id"] == receptor_id]["chain_id"])
if len(chain_ids) == 2:
first_match_idx = self.dataset.encoded_data.feature_names.index(chain_ids[0])
second_match_idx = self.dataset.encoded_data.feature_names.index(chain_ids[1])
if matches[first_match_idx] > 0 and matches[second_match_idx] > 0:
match_identifiers.append(chain_ids[0])
match_identifiers.append(chain_ids[1])
match_values.append(matches[first_match_idx])
match_values.append(matches[second_match_idx])
results_df = pd.DataFrame([match_values], columns=match_identifiers)
results_df.to_csv(filename, index=False)
def _write_repertoire_sizes(self):
"""
Writes the repertoire sizes (# clones & # reads) per subject, per chain.
"""
all_subjects = self.dataset.encoded_data.example_ids
all_chains = sorted(set(self.dataset.encoded_data.feature_annotations["chain"]))
results_df = pd.DataFrame(list(itertools.product(all_subjects, all_chains)),
columns=["subject_id", "chain"])
results_df["n_reads"] = 0
results_df["n_clones"] = 0
for repertoire in self.dataset.repertoires:
rep_counts = repertoire.get_counts()
rep_chains = repertoire.get_chains()
for chain in all_chains:
indices = rep_chains == Chain.get_chain(chain.upper())
results_df.loc[(results_df.subject_id == repertoire.metadata["subject_id"]) & (results_df.chain == chain),
'n_reads'] += np.sum(rep_counts[indices])
results_df.loc[(results_df.subject_id == repertoire.metadata["subject_id"]) & (results_df.chain == chain),
'n_clones'] += len(rep_counts[indices])
results_path = self.result_path / "repertoire_sizes.csv"
results_df.to_csv(results_path, index=False)
return ReportOutput(results_path, "Repertoire sizes")
def _write_receptor_info(self, receptor_info_path) -> List[ReportOutput]:
PathBuilder.build(receptor_info_path)
receptor_chains = self.dataset.encoded_data.feature_annotations
chain_types = receptor_chains["chain"].unique()
first_chains = receptor_chains.loc[receptor_chains.chain == chain_types[0]]
second_chains = receptor_chains.loc[receptor_chains.chain == chain_types[1]]
first_chains.drop(columns=["chain"], inplace=True)
second_chains.drop(columns=["chain"], inplace=True)
on_cols = ["receptor_id"]
if "clonotype_id" in second_chains.columns and first_chains.columns:
on_cols += ["clonotype_id"]
receptors = pd.merge(first_chains, second_chains,
on=on_cols,
suffixes=(f"_{chain_types[0]}", f"_{chain_types[1]}"))
unique_alpha_chains = first_chains.drop_duplicates(subset=["sequence", "v_gene", "j_gene"])
unique_beta_chains = second_chains.drop_duplicates(subset=["sequence", "v_gene", "j_gene"])
unique_receptors = receptors.drop_duplicates(subset=[f"sequence_{chain_types[0]}", f"v_gene_{chain_types[0]}", f"j_gene_{chain_types[0]}",
f"sequence_{chain_types[1]}", f"v_gene_{chain_types[1]}", f"j_gene_{chain_types[1]}"])
receptor_chains_path = receptor_info_path / "all_chains.csv"
receptor_chains.to_csv(receptor_chains_path, index=False)
receptors_path = receptor_info_path / "all_receptors.csv"
receptors.to_csv(receptors_path, index=False)
unique_chain1_path = receptor_info_path / f"unique_{chain_types[0]}_chains.csv"
unique_alpha_chains.to_csv(unique_chain1_path, index=False)
unique_chain2_path = receptor_info_path / f"unique_{chain_types[1]}_chains.csv"
unique_beta_chains.to_csv(unique_chain2_path, index=False)
unique_receptors_path = receptor_info_path / "unique_receptors.csv"
unique_receptors.to_csv(unique_receptors_path, index=False)
return [ReportOutput(path=path, name=name) for path, name in [(receptors_path, "All receptors info"),
(receptor_chains_path, "All receptor chains info"),
(unique_receptors_path, "Unique receptors info"),
(unique_chain1_path, "Unique chain 1 info"),
(unique_chain2_path, "Unique chain 2 info")]]
def _write_sequence_info(self, sequence_info_path) -> List[ReportOutput]:
PathBuilder.build(sequence_info_path)
chains = self.dataset.encoded_data.feature_annotations
unique_chains = chains.drop_duplicates(subset=["sequence", "v_gene", "j_gene"])
chains_path = sequence_info_path / "all_chains.csv"
chains.to_csv(chains_path, index=False)
unique_chains_path = sequence_info_path / "unique_chains.csv"
unique_chains.to_csv(unique_chains_path, index=False)
return [ReportOutput(path=path, name=name) for path, name in [(chains_path, "All chains info"), (unique_chains_path, "Unique chains info")]]
[docs] def check_prerequisites(self):
if self.dataset.encoded_data is None or self.dataset.encoded_data.examples is None:
warnings.warn(f"No encoding was specified for dataset {self.dataset.identifier}. Please use one of the following encodings: MatchedReceptorsEncoder, MatchedSequencesEncoder, MatchedRegexEncoder. Matches report will not be created.")
return False
if self.dataset.encoded_data.encoding not in ("MatchedReceptorsEncoder", "MatchedSequencesEncoder", "MatchedRegexEncoder"):
warnings.warn(f"Encoding {self.dataset.encoded_data.encoding} is not compatible with this report type. Matches report will not be created.")
return False
else:
return True