diff --git a/.gitignore b/.gitignore
index 5f49094716f8f9f8824254962606fa0437ec99df..a704ac2f36b2df643b1fad857a209487fb37ef7b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,3 +4,4 @@
/dataset/
/test.py
/database/
+/dlomix/data/
diff --git a/dlomix/__init__.py b/dlomix/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..4f3165727c670bb8dbf089e6a78f00c7e5d25b05
--- /dev/null
+++ b/dlomix/__init__.py
@@ -0,0 +1,10 @@
+__version__ = "0.0.6"
+
+META_DATA = {
+ "author": "Omar Shouman",
+ "author_email": "o.shouman@tum.de",
+ "description": "Deep Learning for Proteomics",
+ "package_name": "DLOmix",
+ "copyright_text": "2023, Wilhelm Lab, TU Munich.",
+ "github_url": "https://github.com/wilhelm-lab/dlomix",
+}
diff --git a/dlomix/constants.py b/dlomix/constants.py
new file mode 100644
index 0000000000000000000000000000000000000000..e771b3c7f0ad71563c36e9c19b3e024b76807f0a
--- /dev/null
+++ b/dlomix/constants.py
@@ -0,0 +1,100 @@
+DEFAULT_PARQUET_ENGINE = "pyarrow"
+
+retention_time_pipeline_parameters = {
+ "model_params": {"seq_length": 30},
+ "data_params": {
+ "seq_length": 30,
+ },
+ "trained_model_path": "../pretrained_models/retention_time/example_rtmodel/",
+ "trained_model_zipfile_name": "rtmodel.zip",
+ "trained_model_stats": [0.0, 1.0],
+}
+
+retention_time_pipeline_parameters.update(
+ {
+ "trained_model_url": "https://raw.githubusercontent.com/wilhelm-lab/dlomix/develop"
+ + retention_time_pipeline_parameters["trained_model_path"].strip("..")
+ + retention_time_pipeline_parameters["trained_model_zipfile_name"]
+ }
+)
+
+ALPHABET_UNMOD = {
+ "A": 1,
+ "C": 2,
+ "D": 3,
+ "E": 4,
+ "F": 5,
+ "G": 6,
+ "H": 7,
+ "I": 8,
+ "K": 9,
+ "L": 10,
+ "M": 11,
+ "N": 12,
+ "P": 13,
+ "Q": 14,
+ "R": 15,
+ "S": 16,
+ "T": 17,
+ "V": 18,
+ "W": 19,
+ "Y": 20,
+}
+
+# relevant for feature extraction for PTMs, only for reference
+ALPHABET_PTMS = {
+ "A": 1,
+ "C": 2,
+ "D": 3,
+ "E": 4,
+ "F": 5,
+ "G": 6,
+ "H": 7,
+ "I": 8,
+ "K": 9,
+ "L": 10,
+ "M": 11, # amino acids
+ "N": 12,
+ "P": 13,
+ "Q": 14,
+ "R": 15,
+ "S": 16,
+ "T": 17,
+ "V": 18,
+ "W": 19,
+ "Y": 20,
+ "[]-": 21,
+ "-[]": 22, # termini
+ "M[UNIMOD:35]": 23,
+ "S[UNIMOD:21]": 24,
+ "T[UNIMOD:21]": 25,
+ "Y[UNIMOD:21]": 26,
+ "R[UNIMOD:7]": 27,
+ "K[UNIMOD:1]": 28,
+ "K[UNIMOD:121]": 29,
+ "Q(gl)": 30,
+ "R[UNIMOD:34]": 31,
+ "K[UNIMOD:34]": 32,
+ "T(ga)": 33,
+ "S(ga)": 34,
+ "T(gl)": 35,
+ "S(gl)": 36,
+ "C[UNIMOD:4]": 37,
+ "E(gl)": 39,
+ "[ac]-": 38,
+ "K[UNIMOD:36]": 40,
+ "K[UNIMOD:37]": 41,
+ "K[UNIMOD:122]": 42,
+ "K[UNIMOD:58]": 43,
+ "K[UNIMOD:1289]": 44,
+ "K[UNIMOD:747]": 45,
+ "K[UNIMOD:64]": 46,
+ "K[UNIMOD:1848]": 47,
+ "K[UNIMOD:1363]": 48,
+ "K[UNIMOD:1849]": 49,
+ "K[UNIMOD:3]": 50,
+ "R[UNIMOD:36]": 51,
+ "R[UNIMOD:36a]": 52,
+ "P[UNIMOD:35]": 53,
+ "Y[UNIMOD:354]": 54,
+}
diff --git a/dlomix/eval/__init__.py b/dlomix/eval/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..4df20838d7797f0f623af0f6904a50d85b91a39f
--- /dev/null
+++ b/dlomix/eval/__init__.py
@@ -0,0 +1,3 @@
+from .rt_eval import TimeDeltaMetric
+
+__all__ = [TimeDeltaMetric]
diff --git a/dlomix/eval/rt_eval.py b/dlomix/eval/rt_eval.py
new file mode 100644
index 0000000000000000000000000000000000000000..231cbf19846a8aea4c620bb168ce5b8fa49cf951
--- /dev/null
+++ b/dlomix/eval/rt_eval.py
@@ -0,0 +1,83 @@
+import tensorflow as tf
+import tensorflow.keras.backend as K
+
+
+class TimeDeltaMetric(tf.keras.metrics.Metric):
+ """Implementation of the time delta metric as a Keras Metric.
+
+ Parameters
+ ----------
+ mean (int, optional): Mean value of the targets in case normalization was performed. Defaults to 0.
+ std (int, optional): Standard deviation value of the targets in case normalization was performed. Defaults to 1.
+ percentage (float, optional): What percentage of the data points to consider, this is specific to the conmputation of the metric. Defaults to 0.95 which corresponds to 95% of the datapoints and is the mostly used value in papers.
+ name (str, optional): Name of the metric so that it can be reported and used later in Keras History objects. Defaults to 'timedelta'.
+ rescale_targets (bool, optional): Whether to rescale (denormalize) targets or not. Defaults to False.
+ rescale_predictions (bool, optional): Whether to rescale (denormalize) predictions or not. Defaults to False.
+ double_delta (bool, optional): Whether to multiple the computed delta by 2 in order to make it two-sided or not. Defaults to False.
+ """
+
+ def __init__(
+ self,
+ mean=0,
+ std=1,
+ percentage=0.95,
+ name="timedelta",
+ rescale_targets=False,
+ rescale_predictions=False,
+ double_delta=False,
+ **kwargs
+ ):
+
+ super(TimeDeltaMetric, self).__init__(name=name, **kwargs)
+ self.delta = self.add_weight(name="delta", initializer="zeros")
+ self.batch_count = self.add_weight(name="batch-count", initializer="zeros")
+ self.mean = mean
+ self.std = std
+ self.percentage = percentage
+ self.rescale_targets = rescale_targets
+ self.rescale_predictions = rescale_predictions
+ self.double_delta = double_delta
+
+ def update_state(self, y_true, y_pred, sample_weight=None):
+ # rescale
+ if self.rescale_targets:
+ y_true = y_true * self.std + self.mean
+
+ if self.rescale_predictions:
+ y_pred = y_pred * self.std + self.mean
+
+ # find position of the index
+ length = tf.shape(y_true)[0]
+ mark = tf.cast(length, dtype=tf.float32) * self.percentage
+ mark = tf.cast(mark, dtype=tf.int32)
+
+ # compute residuals and sort
+ abs_error = tf.abs(y_true - y_pred)
+ d = tf.sort(abs_error)[mark - 1]
+
+ # two-sided delta
+ if self.double_delta:
+ d = d * 2
+
+ # update count of batches
+ self.batch_count.assign_add(1.0)
+
+ # update delta
+ self.delta.assign_add(tf.math.reduce_sum(d))
+
+ def result(self):
+ # this is simple averaging over the batches, more complex reduction can be added based on domain expertises
+ # Examples are: take max or min of both deltas (translates to a strict or a relaxed metric)
+ return tf.math.divide(self.delta, self.batch_count)
+
+
+# code adopted and modified based on:
+# https://github.com/horsepurve/DeepRTplus/blob/cde829ef4bd8b38a216d668cf79757c07133b34b/RTdata_emb.py
+def delta95_metric(y_true, y_pred):
+ mark95 = tf.cast(
+ tf.cast(tf.shape(y_true)[0], dtype=tf.float32) * 0.95, dtype=tf.int32
+ )
+ abs_error = K.abs(y_true - y_pred)
+ delta = tf.sort(abs_error)[mark95 - 1]
+ norm_range = K.max(y_true) - K.min(y_true)
+ return (delta * 2) / (norm_range)
diff --git a/dlomix/layers/__init__.py b/dlomix/layers/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/dlomix/layers/attention.py b/dlomix/layers/attention.py
new file mode 100644
index 0000000000000000000000000000000000000000..74dc249fe255531df0d33ac2704e23c6a98051a8
--- /dev/null
+++ b/dlomix/layers/attention.py
@@ -0,0 +1,112 @@
+import tensorflow as tf
+import tensorflow.keras.backend as K
+from tensorflow.keras import constraints, initializers, regularizers
+
+
+class DecoderAttentionLayer(tf.keras.layers.Layer):
+ def __init__(self, time_steps):
+ super(DecoderAttentionLayer, self).__init__()
+ self.time_steps = time_steps
+
+ def build(self, input_shape):
+ self.permute = tf.keras.layers.Permute((2, 1))
+ self.dense = tf.keras.layers.Dense(self.time_steps, activation="softmax")
+ self.multiply = tf.keras.layers.Multiply()
+
+ def call(self, inputs):
+ x = self.permute(inputs)
+ x = self.dense(x)
+ x = self.permute(x)
+ x = self.multiply([inputs, x])
+ return x
+
+
+class AttentionLayer(tf.keras.layers.Layer):
+ def __init__(
+ self,
+ context=False,
+ W_regularizer=None,
+ b_regularizer=None,
+ u_regularizer=None,
+ W_constraint=None,
+ b_constraint=None,
+ u_constraint=None,
+ bias=True,
+ **kwargs
+ ):
+ self.supports_masking = True
+ self.init = initializers.get("glorot_uniform")
+ self.W_regularizer = regularizers.get(W_regularizer)
+ self.b_regularizer = regularizers.get(b_regularizer)
+ self.u_regularizer = regularizers.get(u_regularizer)
+ self.W_constraint = constraints.get(W_constraint)
+ self.b_constraint = constraints.get(b_constraint)
+ self.u_constraint = constraints.get(u_constraint)
+ self.bias = bias
+ self.context = context
+ super(AttentionLayer, self).__init__(**kwargs)
+
+ def build(self, input_shape):
+ assert len(input_shape) == 3
+ self.W = self.add_weight(
+ shape=(input_shape[-1],),
+ initializer=self.init,
+ name="{}_W".format(self.name),
+ regularizer=self.W_regularizer,
+ constraint=self.W_constraint,
+ )
+ if self.bias:
+ self.b = self.add_weight(
+ shape=(input_shape[1],),
+ initializer="zero",
+ name="{}_b".format(self.name),
+ regularizer=self.b_regularizer,
+ constraint=self.b_constraint,
+ )
+ else:
+ self.b = None
+ if self.context:
+ self.u = self.add_weight(
+ shape=(input_shape[-1],),
+ initializer=self.init,
+ name="{}_u".format(self.name),
+ regularizer=self.u_regularizer,
+ constraint=self.u_constraint,
+ )
+
+ self.built = True
+
+ def compute_mask(self, input, input_mask=None):
+ return None
+
+ def call(self, x, mask=None):
+ a = K.squeeze(K.dot(x, K.expand_dims(self.W)), axis=-1)
+ if self.bias:
+ a += self.b
+ a = K.tanh(a)
+ if self.context:
+ a = K.squeeze(K.dot(x, K.expand_dims(self.u)), axis=-1)
+ a = K.exp(a)
+ if mask is not None:
+ a *= K.cast(mask, K.floatx())
+ a /= K.cast(K.sum(a, axis=1, keepdims=True) + K.epsilon(), K.floatx())
+ a = K.expand_dims(a)
+ weighted_input = x * a
+ return K.sum(weighted_input, axis=1)
+
+ def compute_output_shape(self, input_shape):
+ return input_shape[0], input_shape[-1]
+
+ def get_config(self):
+ config = {
+ "bias": self.bias,
+ "context": self.context,
+ "W_regularizer": regularizers.serialize(self.W_regularizer),
+ "b_regularizer": regularizers.serialize(self.b_regularizer),
+ "u_regularizer": regularizers.serialize(self.u_regularizer),
+ "W_constraint": constraints.serialize(self.W_constraint),
+ "b_constraint": constraints.serialize(self.b_constraint),
+ "u_constraint": constraints.serialize(self.u_constraint),
+ }
+ base_config = super(AttentionLayer, self).get_config()
+ return dict(list(base_config.items()) + list(config.items()))
diff --git a/dlomix/losses/__init__.py b/dlomix/losses/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..c966525e400d2208418fd356b904461082b7736a
--- /dev/null
+++ b/dlomix/losses/__init__.py
@@ -0,0 +1,3 @@
+from .intensity import masked_pearson_correlation_distance, masked_spectral_distance
+
+__all__ = [masked_spectral_distance, masked_pearson_correlation_distance]
diff --git a/dlomix/losses/intensity.py b/dlomix/losses/intensity.py
new file mode 100644
index 0000000000000000000000000000000000000000..512233d655d2b59a8341c31dcd1abbb8863a5d37
--- /dev/null
+++ b/dlomix/losses/intensity.py
@@ -0,0 +1,69 @@
+import numpy as np
+import tensorflow as tf
+import tensorflow.keras.backend as K
+
+
+def masked_spectral_distance(y_true, y_pred):
+ """Masked, normalized spectral angles between true and pred vectors
+ > arccos(1*1 + 0*0) = 0 > SL = 0 > high correlation
+ > arccos(0*1 + 1*0) = pi/2 > SL = 1 > low correlation
+ """
+
+ # To avoid numerical instability during training on GPUs,
+ # we add a fuzzing constant epsilon of 1×10−7 to all vectors
+ epsilon = K.epsilon()
+
+ # Masking: we multiply values by (true + 1) because then the peaks that cannot
+ # be there (and have value of -1 as explained above) won't be considered
+ pred_masked = ((y_true + 1) * y_pred) / (y_true + 1 + epsilon)
+ true_masked = ((y_true + 1) * y_true) / (y_true + 1 + epsilon)
+
+ # L2 norm
+ pred_norm = K.l2_normalize(true_masked, axis=-1)
+ true_norm = K.l2_normalize(pred_masked, axis=-1)
+
+ # Spectral Angle (SA) calculation
+ # (from the definition below, it is clear that ions with higher intensities
+ # will always have a higher contribution)
+ product = K.sum(pred_norm * true_norm, axis=1)
+ arccos = tf.math.acos(product)
+ return 2 * arccos / np.pi
+
+
+def masked_pearson_correlation_distance(y_true, y_pred):
+ """
+ Calculates the masked Pearson correlation distance between true and predicted intensity vectors.
+
+ The masked Pearson correlation distance is a metric for comparing the similarity between two intensity vectors,
+ taking into account only the non-negative values in the true values tensor (which represent valid peaks).
+
+ Parameters:
+ -----------
+ y_true : tf.Tensor
+ A tensor containing the true values, with shape `(batch_size, num_values)`.
+ y_pred : tf.Tensor
+ A tensor containing the predicted values, with the same shape as `y_true`.
+
+ Returns:
+ --------
+ tf.Tensor
+ A tensor containing the masked Pearson correlation distance between `y_true` and `y_pred`.
+
+ Raises:
+ -------
+ ValueError
+ If `y_true` and `y_pred` have different shapes.
+ """
+ epsilon = K.epsilon()
+
+ # Masking: we multiply values by (true + 1) because then the peaks that cannot
+ # be there (and have value of -1 as explained above) won't be considered
+ pred_masked = ((y_true + 1) * y_pred) / (y_true + 1 + epsilon)
+ true_masked = ((y_true + 1) * y_true) / (y_true + 1 + epsilon)
+
+ mx = tf.math.reduce_mean(true_masked)
+ my = tf.math.reduce_mean(pred_masked)
+ xm, ym = true_masked - mx, pred_masked - my
+ r_num = tf.math.reduce_mean(tf.multiply(xm, ym))
+ r_den = tf.math.reduce_std(xm) * tf.math.reduce_std(ym)
+ return 1 - (r_num / r_den)
diff --git a/dlomix/models/__init__.py b/dlomix/models/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..5f6c4df63bcefeede6e3a74c83e6a4bd8700eb51
--- /dev/null
+++ b/dlomix/models/__init__.py
@@ -0,0 +1,10 @@
+from .base import *
+from .deepLC import *
+from .prosit import *
+
+__all__ = [
+ "RetentionTimePredictor",
+ "PrositRetentionTimePredictor",
+ "DeepLCRetentionTimePredictor",
+ "PrositIntensityPredictor",
+]
diff --git a/dlomix/models/base.py b/dlomix/models/base.py
new file mode 100644
index 0000000000000000000000000000000000000000..632da9bb62b5f938a5ad6c78158722fb78eaf20b
--- /dev/null
+++ b/dlomix/models/base.py
@@ -0,0 +1,81 @@
+import tensorflow as tf
+from tensorflow.keras.layers.experimental import preprocessing
+
+from ..constants import ALPHABET_UNMOD
+
+
+class RetentionTimePredictor(tf.keras.Model):
+ """A simple class for Retention Time prediction models.
+
+ Parameters
+ ----------
+ embedding_dim (int, optional): Dimensionality of the embeddings to be used for representing the Amino Acids. Defaults to 16.
+ seq_length (int, optional): Sequence length of the peptide sequences. Defaults to 30.
+ encoder (str, optional): String for specifying the decoder to use, either based on 1D conv-layers or LSTMs. Defaults to "conv1d".
+ vocab_dict (dict, optional): Dictionary mapping for the vocabulary (the amino acids in this case). Defaults to ALPHABET_UNMOD.
+ """
+
+ def __init__(
+ self,
+ embedding_dim=16,
+ seq_length=30,
+ encoder="conv1d",
+ vocab_dict=ALPHABET_UNMOD,
+ ):
+ super(RetentionTimePredictor, self).__init__()
+
+ # tie the count of embeddings to the size of the vocabulary (count of amino acids)
+ self.embeddings_count = len(vocab_dict) + 2
+
+ self.string_lookup = preprocessing.StringLookup(
+ vocabulary=list(vocab_dict.keys())
+ )
+
+ self.embedding = tf.keras.layers.Embedding(
+ input_dim=self.embeddings_count,
+ output_dim=embedding_dim,
+ input_length=seq_length,
+ )
+
+ self._build_encoder(encoder)
+
+ self.flatten = tf.keras.layers.Flatten()
+ self.regressor = tf.keras.Sequential(
+ [
+ tf.keras.layers.Dense(128, activation="relu"),
+ tf.keras.layers.Dense(64, activation="relu"),
+ ]
+ )
+
+ self.output_layer = tf.keras.layers.Dense(1)
+
+ def _build_encoder(self, encoder_type):
+ if encoder_type.lower() == "conv1d":
+ self.encoder = tf.keras.Sequential(
+ [
+ tf.keras.layers.Conv1D(
+ filters=256, kernel_size=3, padding="same", activation="relu"
+ ),
+ tf.keras.layers.Conv1D(
+ filters=512, kernel_size=3, padding="valid", activation="relu"
+ ),
+ tf.keras.layers.MaxPooling1D(pool_size=2),
+ ]
+ )
+ else:
+ self.encoder = tf.keras.Sequential(
+ [
+ tf.keras.layers.LSTM(256, return_sequences=True),
+ tf.keras.layers.LSTM(256),
+ ]
+ )
+
+ def call(self, inputs, **kwargs):
+ x = self.string_lookup(inputs)
+ x = self.embedding(x)
+ x = self.encoder(x)
+ x = self.flatten(x)
+ x = self.regressor(x)
+ x = self.output_layer(x)
+
+ return x
diff --git a/dlomix/models/deepLC.py b/dlomix/models/deepLC.py
new file mode 100644
index 0000000000000000000000000000000000000000..3c6b6a2b6727d86d71c3b4cd226fd184538a2b6b
--- /dev/null
+++ b/dlomix/models/deepLC.py
@@ -0,0 +1,131 @@
+import tensorflow as tf
+from tensorflow.keras.layers.experimental import preprocessing
+
+from ..constants import ALPHABET_UNMOD
+
+
+class DeepLCRetentionTimePredictor(tf.keras.Model):
+ def __init__(
+ self, seq_length=60, vocab_dict=ALPHABET_UNMOD, use_global_features=False
+ ):
+ super(DeepLCRetentionTimePredictor, self).__init__()
+ self.seq_length = seq_length
+ self._use_global_features = use_global_features
+
+ self.leaky_relu = tf.keras.layers.ReLU(max_value=20, negative_slope=0.1)
+ self.string_lookup = preprocessing.StringLookup(
+ vocabulary=list(vocab_dict.keys())
+ )
+
+ self._build_aminoacid_branch()
+ self._build_diaminoacid_branch()
+ self._build_onehot_encoding_branch()
+ self._build_regressor()
+ self.output_layer = tf.keras.layers.Dense(1)
+
+ if self._use_global_features:
+ self._build_global_features_branch()
+
+ def _build_aminoacid_branch(self):
+ self.aminoacid_branch = tf.keras.Sequential(
+ [
+ self._build_conv_pool_block(n_filters=256, kernel=8, padding="same"),
+ self._build_conv_pool_block(n_filters=128, kernel=8, padding="same"),
+ self._build_conv_pool_block(
+ n_filters=64, kernel=8, padding="same", pool=False
+ ),
+ tf.keras.layers.Flatten(),
+ ]
+ )
+
+ def _build_diaminoacid_branch(self):
+ self.diaminoacid_branch = tf.keras.Sequential(
+ [
+ self._build_conv_pool_block(n_filters=128, kernel=2, padding="same"),
+ self._build_conv_pool_block(n_filters=64, kernel=2, padding="same"),
+ tf.keras.layers.Flatten(),
+ ]
+ )
+
+ def _build_global_features_branch(self):
+ self.global_features_branch = tf.keras.Sequential(
+ [
+ tf.keras.layers.Dense(16, activation=self.leaky_relu),
+ tf.keras.layers.Dense(16, activation=self.leaky_relu),
+ tf.keras.layers.Dense(16, activation=self.leaky_relu),
+ ]
+ )
+
+ def _build_onehot_encoding_branch(self):
+ self.onehot_encoding_branch = tf.keras.Sequential(
+ [
+ self._build_conv_pool_block(
+ n_filters=2,
+ kernel=2,
+ padding="same",
+ activation="tanh",
+ pool_strides=10,
+ pool_size=10,
+ ),
+ tf.keras.layers.Flatten(),
+ ]
+ )
+
+ def _build_regressor(self):
+ self.regressor = tf.keras.Sequential(
+ [tf.keras.layers.Dense(128, activation=self.leaky_relu) for _ in range(5)]
+ )
+
+ def _build_conv_pool_block(
+ self,
+ n_conv_layers=2,
+ n_filters=256,
+ kernel=8,
+ padding="same",
+ activation="leaky_relu",
+ pool=True,
+ pool_strides=2,
+ pool_size=2,
+ ):
+ # leaky relu by default
+ activation_fn = self.leaky_relu
+
+ if activation in ["tanh", "relu"]:
+ activation_fn = activation
+
+ block = tf.keras.Sequential(
+ [
+ tf.keras.layers.Conv1D(
+ filters=n_filters,
+ kernel_size=kernel,
+ padding=padding,
+ activation=activation_fn,
+ )
+ for _ in range(n_conv_layers)
+ ]
+ )
+ if pool:
+ block.add(tf.keras.layers.MaxPooling1D(pool_size=2, strides=2))
+
+ return block
+
+ def call(self, inputs, **kwargs):
+ outputs = {}
+
+ integer_encoded = self.string_lookup(inputs["seq"])
+ onehot_encoded = tf.one_hot(integer_encoded, depth=self.seq_length)
+
+ if self._use_global_features:
+ outputs["global_features_output"] = self.global_features_branch(
+ inputs["global_features"]
+ )
+
+ outputs["onehot_branch_output"] = self.onehot_encoding_branch(onehot_encoded)
+ outputs["aminoacids_branch_output"] = self.aminoacid_branch(inputs["counts"])
+ outputs["diaminoacids_branch_output"] = self.diaminoacid_branch(
+ inputs["di_counts"]
+ )
+
+ concatenated_output = tf.concat(outputs.values(), axis=1)
+ concatenated_output = self.regressor(concatenated_output)
+ return self.output_layer(concatenated_output)
diff --git a/dlomix/models/prosit.py b/dlomix/models/prosit.py
new file mode 100644
index 0000000000000000000000000000000000000000..5fa1b2cde246954a37f15f2e833176e5535b1381
--- /dev/null
+++ b/dlomix/models/prosit.py
@@ -0,0 +1,333 @@
+import tensorflow as tf
+from tensorflow.keras.layers.experimental import preprocessing
+
+from ..constants import ALPHABET_UNMOD
+from ..data.feature_extractors import (
+ ModificationGainFeature,
+ ModificationLocationFeature,
+ ModificationLossFeature,
+)
+from ..layers.attention import AttentionLayer, DecoderAttentionLayer
+
+
+class PrositRetentionTimePredictor(tf.keras.Model):
+ """Implementation of the Prosit model for retention time prediction.
+
+ Parameters
+ -----------
+ embedding_output_dim (int, optional): Size of the embeddings to use. Defaults to 16.
+ seq_length (int, optional): Sequence length of the peptide sequences. Defaults to 30.
+ vocab_dict (dict, optional): Dictionary mapping for the vocabulary (the amino acids in this case). Defaults to None, which is mapped to `ALPHABET_UNMOD`.
+ dropout_rate (float, optional): Probability to use for dropout layers in the encoder. Defaults to 0.5.
+ latent_dropout_rate (float, optional): Probability to use for dropout layers in the regressor layers after encoding. Defaults to 0.1.
+ recurrent_layers_sizes (tuple, optional): A tuple of 2 values for the sizes of the two GRU layers in the encoder. Defaults to (256, 512).
+ regressor_layer_size (int, optional): Size of the dense layer in the regressor after the encoder. Defaults to 512.
+ """
+
+ DEFAULT_INPUT_KEYS = {
+ "SEQUENCE_KEY": "sequence",
+ }
+
+ def __init__(
+ self,
+ embedding_output_dim=16,
+ seq_length=30,
+ vocab_dict=None,
+ dropout_rate=0.5,
+ latent_dropout_rate=0.1,
+ recurrent_layers_sizes=(256, 512),
+ regressor_layer_size=512,
+ ):
+ super(PrositRetentionTimePredictor, self).__init__()
+
+ self.dropout_rate = dropout_rate
+ self.latent_dropout_rate = latent_dropout_rate
+ self.regressor_layer_size = regressor_layer_size
+ self.recurrent_layers_sizes = recurrent_layers_sizes
+
+ if vocab_dict:
+ self.vocab_dict = vocab_dict
+ else:
+ self.vocab_dict = ALPHABET_UNMOD
+
+ # tie the count of embeddings to the size of the vocabulary (count of amino acids)
+ self.embeddings_count = len(self.vocab_dict) + 2
+
+ self.string_lookup = preprocessing.StringLookup(
+ vocabulary=list(self.vocab_dict.keys())
+ )
+
+ self.embedding = tf.keras.layers.Embedding(
+ input_dim=self.embeddings_count,
+ output_dim=embedding_output_dim,
+ input_length=seq_length,
+ )
+ self._build_encoder()
+
+ self.attention = AttentionLayer()
+
+ self.regressor = tf.keras.Sequential(
+ [
+ tf.keras.layers.Dense(self.regressor_layer_size, activation="relu"),
+ tf.keras.layers.Dropout(rate=self.latent_dropout_rate),
+ ]
+ )
+
+ self.output_layer = tf.keras.layers.Dense(1)
+
+ def _build_encoder(self):
+ self.encoder = tf.keras.Sequential(
+ [
+ tf.keras.layers.Bidirectional(
+ tf.keras.layers.GRU(
+ units=self.recurrent_layers_sizes[0], return_sequences=True
+ )
+ ),
+ tf.keras.layers.Dropout(rate=self.dropout_rate),
+ tf.keras.layers.GRU(
+ units=self.recurrent_layers_sizes[1], return_sequences=True
+ ),
+ tf.keras.layers.Dropout(rate=self.dropout_rate),
+ ]
+ )
+
+ def call(self, inputs, **kwargs):
+ if isinstance(inputs, dict):
+ x = inputs.get(
+ PrositRetentionTimePredictor.DEFAULT_INPUT_KEYS["SEQUENCE_KEY"]
+ )
+ else:
+ x = inputs
+ x = self.string_lookup(x)
+ x = self.embedding(x)
+ x = self.encoder(x)
+ x = self.attention(x)
+ x = self.regressor(x)
+ x = self.output_layer(x)
+ return x
+
+
+class PrositIntensityPredictor(tf.keras.Model):
+ """Implementation of the Prosit model for intensity prediction.
+
+ Parameters
+ -----------
+ embedding_output_dim (int, optional): Size of the embeddings to use. Defaults to 16.
+ seq_length (int, optional): Sequence length of the peptide sequences. Defaults to 30.
+ vocab_dict (dict, optional): Dictionary mapping for the vocabulary (the amino acids in this case). Defaults to None, which is mapped to `ALPHABET_UNMOD`.
+ dropout_rate (float, optional): Probability to use for dropout layers in the encoder. Defaults to 0.5.
+ latent_dropout_rate (float, optional): Probability to use for dropout layers in the regressor layers after encoding. Defaults to 0.1.
+ recurrent_layers_sizes (tuple, optional): A tuple of 2 values for the sizes of the two GRU layers in the encoder. Defaults to (256, 512).
+ regressor_layer_size (int, optional): Size of the dense layer in the regressor after the encoder. Defaults to 512.
+ use_ptm_counts (boolean, optional): Whether to use PTM counts and create corresponding layers, has to be aligned with input_keys. Defaults to False.
+ input_keys (dict, optional): dict of string keys and values mapping a fixed key to a value key in the inputs dict from the dataset class. Defaults to None, which corresponds then to the required default input keys `DEFAULT_INPUT_KEYS`.
+ meta_data_keys (list, optional): list of string values corresponding to fixed keys in the inputs dict that are considered meta data. Defaults to None, which corresponds then to the default meta data keys `META_DATA_KEYS`.
+ """
+
+ # consider using kwargs in the call function instead !
+
+ DEFAULT_INPUT_KEYS = {
+ "SEQUENCE_KEY": "sequence",
+ "COLLISION_ENERGY_KEY": "collision_energy",
+ "PRECURSOR_CHARGE_KEY": "precursor_charge",
+ "FRAGMENTATION_TYPE_KEY": "fragmentation_type",
+ }
+
+ # can be extended to include all possible meta data
+ META_DATA_KEYS = [
+ "COLLISION_ENERGY_KEY",
+ "PRECURSOR_CHARGE_KEY",
+ "FRAGMENTATION_TYPE_KEY",
+ ]
+ PTM_INPUT_KEYS = [
+ ModificationLossFeature.__name__.lower(),
+ ModificationGainFeature.__name__.lower(),
+ ModificationLocationFeature.__name__.lower(),
+ ]
+
+ def __init__(
+ self,
+ embedding_output_dim=16,
+ seq_length=30,
+ len_fion=6,
+ vocab_dict=None,
+ dropout_rate=0.2,
+ latent_dropout_rate=0.1,
+ recurrent_layers_sizes=(256, 512),
+ regressor_layer_size=512,
+ use_ptm_counts=False,
+ input_keys=None,
+ meta_data_keys=None,
+ ):
+ super(PrositIntensityPredictor, self).__init__()
+
+ self.dropout_rate = dropout_rate
+ self.latent_dropout_rate = latent_dropout_rate
+ self.regressor_layer_size = regressor_layer_size
+ self.recurrent_layers_sizes = recurrent_layers_sizes
+ self.seq_length = seq_length
+ self.len_fion = len_fion
+ self.use_ptm_counts = use_ptm_counts
+ self.input_keys = input_keys
+ self.meta_data_keys = meta_data_keys
+
+ # maximum number of fragment ions
+ self.max_ion = self.seq_length - 1
+
+ if vocab_dict:
+ self.vocab_dict = vocab_dict
+ else:
+ self.vocab_dict = ALPHABET_UNMOD
+
+ # tie the count of embeddings to the size of the vocabulary (count of amino acids)
+ self.embeddings_count = len(self.vocab_dict) + 2
+
+ self.string_lookup = preprocessing.StringLookup(
+ vocabulary=list(self.vocab_dict.keys())
+ )
+
+ self.embedding = tf.keras.layers.Embedding(
+ input_dim=self.embeddings_count,
+ output_dim=embedding_output_dim,
+ input_length=seq_length,
+ )
+
+ if self.input_keys is None:
+ self.input_keys = PrositIntensityPredictor.DEFAULT_INPUT_KEYS
+
+ if self.meta_data_keys is None:
+ self.meta_data_keys = PrositIntensityPredictor.META_DATA_KEYS
+
+ self._build_encoders()
+ self._build_decoder()
+
+ self.attention = AttentionLayer(name="encoder_att")
+
+ self.meta_data_fusion_layer = tf.keras.Sequential(
+ [
+ tf.keras.layers.Multiply(name="add_meta"),
+ tf.keras.layers.RepeatVector(self.max_ion, name="repeat"),
+ ]
+ )
+
+ self.regressor = tf.keras.Sequential(
+ [
+ tf.keras.layers.TimeDistributed(
+ tf.keras.layers.Dense(self.len_fion), name="time_dense"
+ ),
+ tf.keras.layers.LeakyReLU(name="activation"),
+ tf.keras.layers.Flatten(name="out"),
+ ]
+ )
+
+ def _build_encoders(self):
+ self.meta_encoder = tf.keras.Sequential(
+ [
+ tf.keras.layers.Concatenate(name="meta_in"),
+ tf.keras.layers.Dense(
+ self.recurrent_layers_sizes[1], name="meta_dense"
+ ),
+ tf.keras.layers.Dropout(self.dropout_rate, name="meta_dense_do"),
+ ]
+ )
+
+ self.sequence_encoder = tf.keras.Sequential(
+ [
+ tf.keras.layers.Bidirectional(
+ tf.keras.layers.GRU(
+ units=self.recurrent_layers_sizes[0], return_sequences=True
+ )
+ ),
+ tf.keras.layers.Dropout(rate=self.dropout_rate),
+ tf.keras.layers.GRU(
+ units=self.recurrent_layers_sizes[1], return_sequences=True
+ ),
+ tf.keras.layers.Dropout(rate=self.dropout_rate),
+ ]
+ )
+ if not self.use_ptm_counts:
+ self.ptm_encoder, self.ptm_aa_fusion = None, None
+ else:
+ self.ptm_encoder = tf.keras.Sequential(
+ [
+ tf.keras.layers.Concatenate(name="ptm_ac_loss_gain"),
+ tf.keras.layers.Bidirectional(
+ tf.keras.layers.GRU(
+ units=self.recurrent_layers_sizes[0], return_sequences=True
+ )
+ ),
+ tf.keras.layers.Dropout(rate=self.dropout_rate),
+ tf.keras.layers.GRU(
+ units=self.recurrent_layers_sizes[1], return_sequences=True
+ ),
+ tf.keras.layers.Dropout(rate=self.dropout_rate),
+ ]
+ )
+
+ self.ptm_aa_fusion = tf.keras.layers.Multiply(name="aa_ptm_in")
+
+ def _build_decoder(self):
+ self.decoder = tf.keras.Sequential(
+ [
+ tf.keras.layers.GRU(
+ units=self.regressor_layer_size,
+ return_sequences=True,
+ name="decoder",
+ ),
+ tf.keras.layers.Dropout(rate=self.dropout_rate),
+ DecoderAttentionLayer(self.max_ion),
+ ]
+ )
+
+ def call(self, inputs, **kwargs):
+ peptides_in = inputs.get(self.input_keys["SEQUENCE_KEY"])
+
+ # read meta data from the input dict
+ meta_data = []
+ # note that the value here is the key to use in the inputs dict passed from the dataset
+ for meta_key, key_in_inputs in self.input_keys.items():
+ if meta_key in PrositIntensityPredictor.META_DATA_KEYS:
+ # get the input under the specified key if exists
+ meta_in = inputs.get(key_in_inputs, None)
+ if meta_in is not None:
+ # add the input to the list of meta data inputs
+ meta_data.append(meta_in)
+
+ if self.meta_encoder and len(meta_data) > 0:
+ encoded_meta = self.meta_encoder(meta_data)
+ else:
+ raise ValueError(
+ f"Following metadata keys are expected in the model for Prosit Intesity: {PrositIntensityPredictor.META_DATA_KEYS}. The actual input passed to the model contains the following keys: {list(inputs.keys())}"
+ )
+
+ # read PTM atom count features from the input dict
+ ptm_ac_features = []
+ for ptm_key in PrositIntensityPredictor.PTM_INPUT_KEYS:
+ ptm_ac_f = inputs.get(ptm_key, None)
+ if ptm_ac_f is not None:
+ ptm_ac_features.append(ptm_ac_f)
+
+ if self.ptm_encoder and len(ptm_ac_features) > 0:
+ encoded_ptm = self.ptm_encoder(ptm_ac_features)
+ elif self.use_ptm_counts:
+ raise ValueError(
+ f"PTM features enabled and following PTM features are expected in the model for Prosit Intesity: {PrositIntensityPredictor.PTM_INPUT_KEYS}. The actual input passed to the model contains the following keys: {list(inputs.keys())}"
+ )
+ else:
+ encoded_ptm = None
+
+ x = self.string_lookup(peptides_in)
+ x = self.embedding(x)
+ x = self.sequence_encoder(x)
+
+ if self.use_ptm_counts and self.ptm_aa_fusion and encoded_ptm is not None:
+ x = self.ptm_aa_fusion([x, encoded_ptm])
+
+ x = self.attention(x)
+
+ x = self.meta_data_fusion_layer([x, encoded_meta])
+
+ x = self.decoder(x)
+ x = self.regressor(x)
+
+ return x
diff --git a/dlomix/pipelines/__init__.py b/dlomix/pipelines/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..8d4a3a851f12410fc2245ff8b42d405b64111bac
--- /dev/null
+++ b/dlomix/pipelines/__init__.py
@@ -0,0 +1,3 @@
+from .pipeline import RetentionTimePipeline
+
+__all__ = [RetentionTimePipeline]
diff --git a/dlomix/pipelines/pipeline.py b/dlomix/pipelines/pipeline.py
new file mode 100644
index 0000000000000000000000000000000000000000..926ad54388ef0a4c3208daa8842c9ce0b760d001
--- /dev/null
+++ b/dlomix/pipelines/pipeline.py
@@ -0,0 +1,101 @@
+import zipfile
+from os import makedirs
+from os.path import dirname, splitext
+
+import numpy as np
+import requests
+
+from ..constants import retention_time_pipeline_parameters
+from ..data.RetentionTimeDataset import RetentionTimeDataset
+from ..models.base import RetentionTimePredictor
+from ..reports import RetentionTimeReport
+
+# pipelines can be used to train the model further or from scratch given a dataset
+# add string arguments (e.g. prosit to create the model, data source to create the dataset)
+
+# if neither train nor test are provided --> use toy datasets to (train if necessary or load pre-trained weights), predict on test, and generate report
+# if test only --> load pre-trained weights, predict and generate report
+# if train and test --> do what you have to do
+
+
+class RetentionTimePipeline:
+ def __init__(self, pre_trained=True):
+ super(RetentionTimePipeline, self).__init__()
+ self.model = None
+ self.test_dataset = None
+ self.pre_trained = pre_trained
+
+ # pass the config in the constructor
+ # refactor to have a base class Pipeline
+
+ self._build_model()
+
+ def _build_model(self):
+ self.model = RetentionTimePredictor(
+ **retention_time_pipeline_parameters["model_params"]
+ )
+
+ if self.pre_trained:
+ self._download_unzip_pretrained_model(
+ retention_time_pipeline_parameters["trained_model_url"],
+ retention_time_pipeline_parameters["trained_model_path"]
+ + retention_time_pipeline_parameters["trained_model_zipfile_name"],
+ )
+
+ self.model.load_weights(
+ retention_time_pipeline_parameters["trained_model_path"]
+ + splitext(
+ retention_time_pipeline_parameters["trained_model_zipfile_name"]
+ )[0]
+ )
+
+ def _download_unzip_pretrained_model(self, model_url, save_path):
+ makedirs(model_url)
+ r = requests.get(model_url)
+
+ with open(save_path, "wb") as f:
+ f.write(r.content)
+
+ self._unzip_model(save_path)
+
+ def _unzip_model(self, model_zipfile_path):
+ zip_ref = zipfile.ZipFile(model_zipfile_path)
+ model_folder = dirname(model_zipfile_path)
+ zip_ref.extractall(model_folder)
+ zip_ref.close()
+
+ """
+
+ Predict retention times given data either as numpy array of sequences or a filepath to a csv file
+
+ """
+
+ def predict(self, data=None):
+ if not (isinstance(data, str) or isinstance(data, np.ndarray)):
+ raise ValueError(
+ "Dataset should be provided either as a numpy array or a string pointing to a file."
+ )
+
+ self.test_dataset = RetentionTimeDataset(
+ data,
+ **retention_time_pipeline_parameters["data_params"],
+ val_ratio=0,
+ test=True
+ )
+ (
+ self.test_dataset.data_mean,
+ self.test_dataset.data_std,
+ ) = retention_time_pipeline_parameters["trained_model_stats"]
+
+ predictions = self.model.predict(self.test_dataset.test_data)
+ predictions = self.test_dataset.denormalize_targets(predictions)
+ predictions = predictions.ravel()
+
+ return predictions
+
+ def predict_report(self, data, output_path="./") -> None:
+ predictions = self.predict(data)
+ report = RetentionTimeReport(output_path=output_path, history=None)
+
+ test_targets = self.test_dataset.get_split_targets(split="test")
+ report.generate_report(test_targets, predictions)
diff --git a/dlomix/reports/IntensityReport.py b/dlomix/reports/IntensityReport.py
new file mode 100644
index 0000000000000000000000000000000000000000..21d0196b2ee26998415dd0fa8a9b458b88f06500
--- /dev/null
+++ b/dlomix/reports/IntensityReport.py
@@ -0,0 +1,74 @@
+from os.path import join
+
+import pandas as pd
+import seaborn as sns
+
+from .postprocessing import normalize_intensity_predictions
+from .Report import PDFFile, Report
+
+
+class IntensityReport(Report):
+ """Report generation for Fragment Ion Intensity Prediction tasks."""
+
+ TARGETS_LABEL = "x"
+ PREDICTIONS_LABEL = "y"
+ DEFAULT_BATCH_SIZE = 600
+
+ def __init__(self, output_path, history, figures_ext="png", batch_size=0):
+ super(IntensityReport, self).__init__(output_path, history, figures_ext)
+
+ self.pdf_file = PDFFile("DLOmix - Fragment Ion Intensity Report")
+
+ if batch_size:
+ self.batch_size = batch_size
+ else:
+ self.batch_size = IntensityReport.DEFAULT_BATCH_SIZE
+
+ def generate_report(self, dataset, predictions):
+ self._init_report_resources()
+
+ predictions_df = self.generate_intensity_results_df(dataset, predictions)
+ self.plot_all_metrics()
+
+ # make custom plots
+ self.plot_spectral_angle(predictions_df)
+
+ self._compile_report_resources_add_pdf_pages()
+ self.pdf_file.output(join(self._output_path, "intensity_Report.pdf"), "F")
+
+ def generate_intensity_results_df(self, dataset, predictions):
+ predictions_df = pd.DataFrame()
+
+ predictions_df["sequences"] = dataset.sequences
+ predictions_df["intensities_pred"] = predictions.tolist()
+ predictions_df["precursor_charge_onehot"] = dataset.precursor_charge.tolist()
+ predictions_df["intensities_raw"] = dataset.intensities.tolist()
+
+ return predictions_df
+
+ def plot_spectral_angle(self, predictions_df):
+ """Create spectral plot
+
+ Arguments
+ ---------
+ predictions_df: dataframe with raw intensities, predictions, sequences, precursor_charges
+ """
+
+ predictions_acc = normalize_intensity_predictions(
+ predictions_df, self.batch_size
+ )
+ violin_plot = sns.violinplot(predictions_acc["spectral_angle"])
+
+ save_path = join(
+ self._output_path, "violin_spectral_angle_plot" + self._figures_ext
+ )
+
+ fig = violin_plot.get_figure()
+ fig.savefig(save_path)
+
+ self._add_report_resource(
+ "spectral_angle_plot",
+ "Spectral angle plot",
+ "The following figure shows the spectral angle plot for the test data.",
+ save_path,
+ )
diff --git a/dlomix/reports/Report.py b/dlomix/reports/Report.py
new file mode 100644
index 0000000000000000000000000000000000000000..8b40cb3391e37d05a4ed93117b7294f1298c01a4
--- /dev/null
+++ b/dlomix/reports/Report.py
@@ -0,0 +1,238 @@
+import abc
+import glob
+import warnings
+from os import makedirs
+from os.path import join
+
+import tensorflow as tf
+from fpdf import FPDF
+from matplotlib import pyplot as plt
+
+
+class Report(abc.ABC):
+ """Base class for reports, child classes should implement the abstract method generate_report.
+
+ Parameters
+ ----------
+ output_path: path to save output files and figures.
+ history : reference to a Keras History object or its history dict attribute (History.history).
+ figures_ext: File extension and format for saving figures.
+ """
+
+ VALID_FIGURE_FORMATS = ["pdf", "jpeg", "jpg", "png"]
+
+ def __init__(self, output_path, history, figures_ext):
+ self._output_path = output_path
+ makedirs(self._output_path, exist_ok=True)
+
+ if history is None:
+ warnings.warn(
+ "The passed History object is None, no training/validation data can be reported."
+ )
+ self._history_dict = {}
+ else:
+ self._set_history_dict(history)
+ self._set_figures_format(figures_ext)
+
+ # an empty dict to use to list the report resources
+ self._init_report_resources()
+
+ # empty pdf file
+ self.pdf_file = None
+
+ def _set_history_dict(self, history):
+ if isinstance(history, dict):
+ self._history_dict = history
+ elif not isinstance(history, tf.keras.callbacks.History):
+ raise ValueError(
+ "Reporting requires a History object (tf.keras.callbacks.History) or its history dict attribute (History.history), which is returned from a call to "
+ f"model.fit(). Passed history argument is of type {type(history)} ",
+ )
+ elif not hasattr(history, "history"):
+ raise ValueError(
+ "The passed History object does not have a history attribute, which is a dict with results."
+ )
+ else:
+ self._history_dict = history.history
+
+ if len(self._history_dict.keys()) == 0:
+ warnings.warn(
+ "The passed History object contains an empty history dict, no training was done."
+ )
+
+ def _set_figures_format(self, figures_ext):
+ figures_ext = figures_ext.lower()
+ if figures_ext.startswith("."):
+ figures_ext = figures_ext[1:]
+ if figures_ext not in Report.VALID_FIGURE_FORMATS:
+ raise ValueError(
+ f"Allowed figure formats are: {Report.VALID_FIGURE_FORMATS}"
+ )
+ self._figures_ext = "." + figures_ext
+
+ def _get_all_saved_plots(self):
+ all_plots = glob.glob(join(self._output_path, "*" + self._figures_ext))
+ return all_plots
+
+ def _add_report_resource(self, key, title, paragraph_text, value):
+ self._report_resources[key] = (title, paragraph_text, value)
+
+ def _init_report_resources(self):
+ self._report_resources = {}
+
+ def _compile_report_resources_add_pdf_pages(self):
+ for key, resource in self._report_resources.items():
+ value_is_fig_path = self._figures_ext in str(resource[2])
+ plot_word_is_in_key = "plot" in key
+ if value_is_fig_path or plot_word_is_in_key:
+ self.pdf_file.add_content_plot_page(
+ section_title=resource[0],
+ section_body=resource[1],
+ plot_filepath=resource[2],
+ )
+ else:
+ self.pdf_file.add_content_text_page(
+ section_title=resource[0], section_body=resource[1]
+ )
+
+ def plot_keras_metric(self, metric_name, save_plot=True):
+ """Plot a keras metric given its name and the history object returned by model.fit()
+
+ Arguments
+ ---------
+ metric_name: String with the name of the metric.
+ save_plot (bool, optional): whether to save plot to disk or not. Defaults to True.
+ """
+
+ if metric_name.lower() not in self._history_dict.keys():
+ raise ValueError(
+ f"Metric name to plot is not available in the history dict. Available metrics to plot are {self._history_dict.keys()}",
+ )
+
+ if (
+ "val_" + metric_name.lower() not in self._history_dict.keys()
+ and metric_name.lower() not in ["lr"]
+ ):
+ raise ValueError(
+ f"""No validation epochs were run during training, the metric name to plot is not available in the history dict.
+ Available metrics to plot are {self._history_dict.keys()}
+ """,
+ )
+ plt.plot(self._history_dict[metric_name])
+ plt.plot(self._history_dict["val_" + metric_name])
+ plt.title(metric_name)
+ plt.ylabel(metric_name)
+ plt.xlabel("epoch")
+ plt.legend(["train", "val"], loc="upper left")
+ if save_plot:
+ save_path = join(self._output_path, metric_name + self._figures_ext)
+ plt.savefig(save_path)
+ plt.show()
+ plt.close()
+ metric_name_spaced = metric_name.replace("_", " ")
+ self._add_report_resource(
+ metric_name + "_plot",
+ metric_name_spaced.title(),
+ f"The following figure shows the {metric_name_spaced} for training and validation.",
+ save_path,
+ )
+
+ def plot_all_metrics(self):
+ """Plot all available Keras metrics in the History object."""
+ metrics = self._history_dict.keys()
+ metrics = filter(lambda x: not x.startswith(tuple(["val_", "_"])), metrics)
+ print("Plotting all metrics: ", list(metrics))
+ for metric in metrics:
+ self.plot_keras_metric(metric)
+
+ @abc.abstractmethod
+ def generate_report(self, targets, predictions, **kwargs):
+ """Abstract method to generate a complete report. Child classes need to implement this method.
+
+ Arguments
+ ---------
+ targets: Array with target values.
+ predictions: Array with prediction values.
+ """
+
+
+class PDFFile(FPDF):
+ """PDF file template class.
+
+ Parameters
+ ----------
+ title: Title for the pdf file
+ """
+
+ PAGE_WIDTH = 210
+ PAGE_HEIGHT = 297
+
+ SECTION_PARAGRAPH_FONT = ["Arial", "", 11]
+ SECTION_TITLE_FONT = ["Arial", "B", 13]
+ LINE_HEIGHT = 5
+
+ def __init__(self, title):
+ super().__init__()
+ self.title = title
+ self.width = PDFFile.PAGE_WIDTH
+ self.height = PDFFile.PAGE_HEIGHT
+
+ self.set_auto_page_break(True)
+ self.document_empty = True
+
+ def header(self):
+ self.set_font("Arial", "B", 11)
+ self.cell(self.width - 80)
+ self.cell(60, 1, self.title, 0, 0, "R")
+ self.ln(20)
+
+ def footer(self):
+ # Page numbers in the footer
+ self.set_y(-15)
+ self.set_font("Arial", "I", 8)
+ self.set_text_color(128)
+ self.cell(0, 10, "Page " + str(self.page_no()), 0, 0, "C")
+
+ def _add_plot(self, plot_filepath):
+ self.image(plot_filepath)
+ self.ln(3 * PDFFile.LINE_HEIGHT)
+
+ def _add_section_content(self, section_title, section_body):
+ if section_title != "":
+ self.set_font(*PDFFile.SECTION_TITLE_FONT)
+ self.cell(w=0, txt=section_title)
+ self.ln(PDFFile.LINE_HEIGHT)
+ if section_body != "":
+ self.set_font(*PDFFile.SECTION_PARAGRAPH_FONT)
+ self.multi_cell(w=0, h=PDFFile.LINE_HEIGHT, txt=section_body)
+ self.ln(PDFFile.LINE_HEIGHT)
+
+ def _create_first_page_if_document_empty(self):
+ if self.document_empty:
+ self.add_page()
+ self.document_empty = False
+
+ def add_content_text_page(self, section_title, section_body):
+ """Add a section title and a paragraph.
+
+ Arguments
+ ---------
+ section_title: title for the section.
+ section_body: paragraph text to add.
+ """
+ self._create_first_page_if_document_empty()
+ self._add_section_content(section_title, section_body)
+
+ def add_content_plot_page(self, plot_filepath, section_title="", section_body=""):
+ """Add a new page with a section title, a paragraph, and a plot. At least a plot has to be provided.
+
+ Arguments
+ ---------
+ plot_filepath (str): filepath of the plot to be inserted in the new page.
+ section_title (str, optional): title for the section. Defaults to "".
+ section_body (str, optional): paragraph text to add. Defaults to "".
+ """
+
+ self._create_first_page_if_document_empty()
+ self._add_section_content(section_title, section_body)
+ self._add_plot(plot_filepath)
diff --git a/dlomix/reports/RetentionTimeReport.py b/dlomix/reports/RetentionTimeReport.py
new file mode 100644
index 0000000000000000000000000000000000000000..5842635928d9af73726a9530526238c0b305c75d
--- /dev/null
+++ b/dlomix/reports/RetentionTimeReport.py
@@ -0,0 +1,169 @@
+from os.path import join
+from warnings import warn
+
+import numpy as np
+from matplotlib import pyplot as plt
+from matplotlib.colors import LogNorm
+from matplotlib.ticker import LogLocator
+
+from .Report import PDFFile, Report
+
+
+class RetentionTimeReport(Report):
+ """Report generation for Retention Time Prediction tasks."""
+
+ TARGETS_LABEL = "iRT (measured)"
+ PREDICTIONS_LABEL = "iRT (predicted)"
+
+ def __init__(self, output_path, history, figures_ext="png"):
+ super(RetentionTimeReport, self).__init__(output_path, history, figures_ext)
+
+ warn(f"{self.__class__.__name__} This class is deprecated and will not further developed. Use RetentionTimeReportWandb instead for creating a report with the Weights & Biases Report API.",
+ DeprecationWarning,
+ stacklevel=2
+ )
+
+
+ self.pdf_file = PDFFile("DLOmix - Retention Time Report")
+
+ def generate_report(self, targets, predictions, **kwargs):
+ self._init_report_resources()
+
+ _ = self.calculate_r2(targets, predictions)
+ self.plot_all_metrics()
+ self.plot_residuals(targets, predictions)
+ self.plot_density(targets, predictions)
+
+ self._compile_report_resources_add_pdf_pages()
+
+ self.pdf_file.output(join(self._output_path, "iRT_Report.pdf"), "F")
+
+ def calculate_r2(self, targets, predictions):
+ """Calculate R-squared using sklearn given true targets and predictions
+
+ Arguments
+ ---------
+ targets: Array with target values
+ predictions: Array with prediction values
+
+ Returns:
+ r_squared (float): float value of R squared
+ """
+ from sklearn.metrics import r2_score
+
+ r2 = r2_score(np.ravel(targets), np.ravel(predictions))
+
+ self._add_report_resource(
+ "r2",
+ "R-Squared",
+ f"The R-squared value for the predictions is {round(r2, 4)}",
+ r2,
+ )
+
+ return r2
+
+ def plot_residuals(self, targets, predictions, xrange=(0, 0)):
+ """Plot histogram of residuals
+
+ Argsuments
+ ----------
+ targets: Array with target values
+ predictions: Array with prediction values
+ xrange (tuple, optional): X-axis range for plotting the histogram. Defaults to (-10, 10).
+ """
+ error = np.ravel(predictions) - np.ravel(targets)
+
+ x_min, x_max = xrange
+ if xrange == (0, 0):
+ mean, std_dev = np.mean(error), np.std(error)
+ x_min, x_max = mean - (3 * std_dev), mean + (3 * std_dev)
+
+ bins = np.linspace(x_min, x_max, 200)
+
+ plt.hist(error, bins, alpha=0.5, color="orange")
+ plt.title("Historgram of Residuals")
+ plt.xlabel("Residual value")
+ plt.ylabel("Count")
+ save_path = join(self._output_path, "histogram_residuals" + self._figures_ext)
+ plt.savefig(save_path)
+ plt.show()
+ plt.close()
+
+ self._add_report_resource(
+ "residuals_plot",
+ "Error Residuals",
+ "The following plot shows a historgram of residuals for the test data.",
+ save_path,
+ )
+
+ def plot_density(
+ self,
+ targets,
+ predictions,
+ irt_delta95=5,
+ palette="Reds_r",
+ delta95_line_color="#36479E",
+ nbins=1000,
+ ):
+ """Create density plot
+
+ Arguments
+ ---------
+ targets: Array with target values
+ predictions: Array with prediction values
+ irt_delta95 (int, optional): iRT Value of the delta 95% . Defaults to 5.
+ palette (str, optional): Color palette from matplotlib. Defaults to 'Reds_r'.
+ delta95_line_color (str, optional): Color for the delta 95% line. Defaults to '#36479E'.
+ nbins (int, optional): Number of bins to use for creating the 2D histogram. Defaults to 1000.
+ """
+
+ H, xedges, yedges = np.histogram2d(targets, predictions, bins=nbins)
+
+ x_min = np.min(targets)
+ x_max = np.max(targets)
+
+ # H needs to be rotated and flipped
+ H = np.rot90(H)
+ H = np.flipud(H)
+
+ # Mask zeros
+ Hmasked = np.ma.masked_where(H == 0, H) # Mask pixels with a value of zero
+
+ # Plot 2D histogram using pcolor
+ cm = plt.cm.get_cmap(palette)
+ plt.pcolormesh(
+ xedges, yedges, Hmasked, cmap=cm, norm=LogNorm(vmin=1e0, vmax=1e2)
+ )
+
+ plt.xlabel(RetentionTimeReport.TARGETS_LABEL, fontsize=18)
+ plt.ylabel(RetentionTimeReport.PREDICTIONS_LABEL, fontsize=18)
+
+ cbar = plt.colorbar(ticks=LogLocator(subs=range(5)))
+ cbar.ax.set_ylabel("Counts", fontsize=14)
+
+ plt.plot([x_min, x_max], [x_min, x_max], c="black")
+ plt.plot(
+ [x_min, x_max],
+ [x_min - irt_delta95, x_max - irt_delta95],
+ color=delta95_line_color,
+ )
+ plt.plot(
+ [x_min, x_max],
+ [x_min + irt_delta95, x_max + irt_delta95],
+ color=delta95_line_color,
+ )
+
+ font_size = 14 # Adjust as appropriate.
+ cbar.ax.tick_params(labelsize=font_size)
+ cbar.ax.minorticks_on()
+ save_path = join(self._output_path, "density_plot" + self._figures_ext)
+ plt.savefig(save_path)
+ plt.show()
+ plt.close()
+
+ self._add_report_resource(
+ "density_plot",
+ "Density Plot",
+ "The following figure shows the density plot with the delta-95 highlighted for the test data.",
+ save_path,
+ )
diff --git a/dlomix/reports/RetentionTimeReportModelComparisonWandb.py b/dlomix/reports/RetentionTimeReportModelComparisonWandb.py
new file mode 100644
index 0000000000000000000000000000000000000000..0aa736342da85c45f61cae37fd4dd12102919b80
--- /dev/null
+++ b/dlomix/reports/RetentionTimeReportModelComparisonWandb.py
@@ -0,0 +1,321 @@
+import os
+import re
+
+import numpy as np
+import pandas as pd
+import wandb
+import wandb.apis.reports as wr
+
+from ..data.RetentionTimeDataset import RetentionTimeDataset
+
+
+class RetentionTimeReportModelComparisonWandb:
+
+ # Wilhelmlab WandB account that has all VEGA presets required for the reports
+ VEGA_LITE_PRESETS_ID = "prosit-compms"
+
+ def __init__(
+ self,
+ models: dict,
+ project: str,
+ title: str,
+ description: str,
+ test_dataset: RetentionTimeDataset,
+ ):
+ """Creates WandB report for comparing models.
+
+ Args:
+ models (dict): keys are model names, values are model objects
+ project (str): Name of the project.
+ title (str): Title of the report.
+ description (str): Description of the report.
+ test_dataset (RetentionTimeDataset): Test dataset object to compare predictions of models on.
+ """
+ self.project = project
+ self.title = title
+ self.description = description
+ self.models = models
+ self.test_dataset = test_dataset
+ self.entity = wandb.apis.PublicApi().default_entity
+ self.api = wandb.Api()
+
+ def create_report(
+ self,
+ add_data_section=True,
+ add_residuals_section=True,
+ add_r2_section=True,
+ add_density_section=True,
+ ):
+ """Creates the report in wandb.
+
+ Args:
+ add_data_section (bool, optional): Add a section for input data to the report. Defaults to True.
+ add_residuals_section (bool, optional): Add a section for residual plots. Defaults to True.
+ add_r2_section (bool, optional): Add a section for the R2 metric. Defaults to True.
+ add_density_section (bool, optional): Add a section for the density plot. Defaults to True.
+ """
+ report = wr.Report(
+ project=self.project, title=self.title, description=self.description
+ )
+
+ report.blocks = [wr.TableOfContents()]
+
+ if add_data_section:
+ report.blocks += self._build_data_section()
+ if add_residuals_section:
+ report.blocks += self._build_residuals_section()
+ if add_r2_section:
+ report.blocks += self._build_r2_section()
+ if add_density_section:
+ report.blocks += self._build_density_section()
+
+ report.save()
+
+ def calculate_r2(self, targets, predictions):
+ from sklearn.metrics import r2_score
+
+ r2 = r2_score(targets, predictions)
+ return r2
+
+ def calculate_residuals(self, targets, predictions):
+ residuals = predictions - targets
+ return residuals
+
+ def _build_data_section(self):
+ data_block = [
+ wr.H1(text="Data"),
+ wr.P(
+ "The following section is showing a simple explorative data analysis of the used dataset. The first histogram shows the distribution of peptide lengths in the data set, while the second histogram shows the distribution of indexed retention times."
+ ),
+ wr.PanelGrid(
+ runsets=[
+ wr.Runset(self.entity, self.project),
+ ],
+ panels=[
+ wr.CustomChart(
+ query={"summaryTable": {"tableKey": self.table_key_len}},
+ chart_name=f"{RetentionTimeReportModelComparisonWandb.VEGA_LITE_PRESETS_ID}/histogram_peptide_length",
+ chart_fields={"value": self.test_dataset.sequence_col},
+ ),
+ wr.CustomChart(
+ query={"summaryTable": {"tableKey": self.table_key_rt}},
+ chart_name=f"{RetentionTimeReportModelComparisonWandb.VEGA_LITE_PRESETS_ID}/histogram_irt",
+ chart_fields={"value": self.test_dataset.target_col},
+ ),
+ ],
+ ),
+ wr.HorizontalRule(),
+ ]
+ return data_block
+
+ def _build_residuals_section(self):
+ panel_list_models = []
+ for model in self.models:
+ panel_list_models.append(
+ wr.CustomChart(
+ query={"summaryTable": {"tableKey": f"results_table_{model}"}},
+ chart_name=f"{RetentionTimeReportModelComparisonWandb.VEGA_LITE_PRESETS_ID}/histogram_residuals_irt",
+ chart_fields={"value": "residuals", "name": model},
+ )
+ )
+
+ residuals_block = [
+ wr.H1(text="Residuals"),
+ wr.P(
+ "This section shows the residuals histograms. Each plot shows the residuals of each of the compared models"
+ ),
+ wr.PanelGrid(
+ runsets=[
+ wr.Runset(self.entity, self.project),
+ ],
+ panels=panel_list_models,
+ ),
+ wr.HorizontalRule(),
+ ]
+
+ return residuals_block
+
+ def _build_r2_section(self):
+ r2_block = [
+ wr.H1(text="R2"),
+ wr.P("The following plot displays the R2 score for all the compared models."),
+ wr.PanelGrid(
+ runsets=[
+ wr.Runset(self.entity, self.project),
+ ],
+ panels=[
+ wr.BarPlot(
+ title="R2",
+ metrics=["r2"],
+ orientation="h",
+ title_x="R2",
+ max_runs_to_show=20,
+ max_bars_to_show=20,
+ font_size="auto",
+ ),
+ ],
+ ),
+ wr.HorizontalRule(),
+ ]
+ return r2_block
+
+ def _build_density_section(self, irt_delta95=5):
+ panel_list_models = []
+ targets = self.test_dataset.get_split_targets(
+ split=self.test_dataset.main_split
+ )
+ for model in self.models:
+ panel_list_models.append(
+ wr.CustomChart(
+ query={"summaryTable": {"tableKey": f"results_table_{model}"}},
+ chart_name=f"{RetentionTimeReportModelComparisonWandb.VEGA_LITE_PRESETS_ID}/density_plot_irt",
+ chart_fields={
+ "measured": "irt",
+ "predicted": "predicted_irt",
+ "name": model,
+ "irt_delta95": irt_delta95,
+ },
+ )
+ )
+
+ density_block = [
+ wr.H1(text="Density"),
+ wr.P(
+ "This section displays the density plots for all compared models."
+ ),
+ wr.PanelGrid(
+ runsets=[
+ wr.Runset(self.entity, self.project),
+ ],
+ panels=panel_list_models,
+ ),
+ wr.HorizontalRule(),
+ ]
+
+ return density_block
+
+ def compare_models(self):
+ for model in self.models:
+ # initialize WANDB
+ current_model = model
+ wandb.init(project=self.project, name=current_model)
+
+ # predict on test_dataset
+ predictions = self.models[model].predict(self.test_dataset.test_data)
+ predictions = predictions.ravel()
+ targets = self.test_dataset.get_split_targets(
+ split=self.test_dataset.main_split
+ )
+ # create result df
+ results_df = pd.DataFrame(
+ {
+ "sequence": self.test_dataset.sequences,
+ "irt": targets,
+ "predicted_irt": predictions,
+ "residuals": self.calculate_residuals(targets, predictions),
+ }
+ )
+ # log df as table to wandb
+ table = wandb.Table(dataframe=results_df)
+ wandb.log({f"results_table_{current_model}": table})
+
+ # log r2 to wandb
+ r2 = self.calculate_r2(targets, predictions)
+ wandb.log({"r2": r2})
+
+ # finish run
+ wandb.finish()
+
+ # function to log sequence length table to wandb
+ def log_sequence_length_table(
+ self, data: pd.DataFrame, seq_col: str = "modified_sequence"
+ ):
+ name_hist = "counts_hist"
+ counts = self.count_seq_length(data, seq_col)
+ # convert to df for easier handling
+ counts_df = counts.to_frame()
+ table = wandb.Table(dataframe=counts_df)
+ # log to wandb
+ hist = wandb.plot_table(
+ vega_spec_name=f"{RetentionTimeReportModelComparisonWandb.VEGA_LITE_PRESETS_ID}/histogram_peptide_length",
+ data_table=table,
+ fields={"value": seq_col},
+ )
+ wandb.log({name_hist: hist})
+ name_hist_table = name_hist + "_table"
+ return name_hist_table
+
+ # function to count sequence length
+ def count_seq_length(self, data: pd.DataFrame, seq_col: str) -> pd.Series:
+ pattern = re.compile(r"\[UNIMOD:.*\]", re.IGNORECASE)
+ data[seq_col].replace(pattern, "", inplace=True)
+ return data[seq_col].str.len()
+
+ # function to log retention time table to wandb
+ def log_rt_table(self, data: pd.DataFrame, rt_col: str = "indexed_retention_time"):
+ name_hist = "rt_hist"
+ rt = data.loc[:, rt_col]
+ # convert to df for easier handling
+ rt_df = rt.to_frame()
+ table = wandb.Table(dataframe=rt_df)
+ # log to wandb
+ hist = wandb.plot_table(
+ vega_spec_name=f"{RetentionTimeReportModelComparisonWandb.VEGA_LITE_PRESETS_ID}/histogram_irt",
+ data_table=table,
+ fields={"value": rt_col},
+ )
+ wandb.log({name_hist: hist})
+ name_hist_table = name_hist + "_table"
+ return name_hist_table
+
+ def log_data(self):
+ wandb.init(project=self.project, name="data_run")
+ # check if datasource is a string
+ if isinstance(self.test_dataset.data_source, str):
+ # read corresponding file
+ file_extension = os.path.splitext(self.test_dataset.data_source)[-1].lower()
+
+ if file_extension == ".csv":
+ data = pd.read_csv(self.test_dataset.data_source)
+ if file_extension == ".json":
+ data = pd.read_json(self.test_dataset.data_source)
+ if file_extension == ".parquet":
+ data = pd.read_parquet(
+ self.test_dataset.data_source, engine="fastparquet"
+ )
+
+ self.table_key_len = self.log_sequence_length_table(
+ data, self.test_dataset.sequence_col
+ )
+ self.table_key_rt = self.log_rt_table(data, self.test_dataset.target_col)
+
+ # check if datasource is a tuple of two ndarrays or two lists
+ if (
+ isinstance(self.test_dataset.data_source, tuple)
+ and all(
+ isinstance(item, (np.ndarray, list))
+ for item in self.test_dataset.data_source
+ )
+ and len(self.test_dataset.data_source) == 2
+ ):
+ data = pd.DataFrame(
+ {
+ self.test_dataset.sequence_col: self.test_dataset.data_source[0],
+ self.test_dataset.target_col: self.test_dataset.data_source[1],
+ }
+ )
+ self.table_key_len = self.log_sequence_length_table(
+ data, self.test_dataset.sequence_col
+ )
+ self.table_key_rt = self.log_rt_table(data, self.test_dataset.target_col)
+
+ # check if datasource is a single ndarray or list
+ # does not work? maybe error in RetentionTimeDataset
+ if isinstance(self.test_dataset.data_source, (np.ndarray, list)):
+ data = pd.DataFrame(
+ {self.test_dataset.sequence_col: self.test_dataset.data_source}
+ )
+ self.table_key_len = self.log_sequence_length_table(
+ data, self.test_dataset.sequence_col
+ )
+ wandb.finish()
diff --git a/dlomix/reports/RetentionTimeReportRunComparisonWandb.py b/dlomix/reports/RetentionTimeReportRunComparisonWandb.py
new file mode 100644
index 0000000000000000000000000000000000000000..3921f553cde3cc6fdcce76d12a76d5f90bf688f7
--- /dev/null
+++ b/dlomix/reports/RetentionTimeReportRunComparisonWandb.py
@@ -0,0 +1,430 @@
+import os
+import re
+
+import numpy as np
+import pandas as pd
+import wandb
+import wandb.apis.reports as wr
+from wandb.keras import WandbCallback, WandbMetricsLogger
+
+from ..data import RetentionTimeDataset
+
+# ToDo: add R2 plot, TimeDelta plot, residuals
+
+
+class RetentionTimeReportRunComparisonWandb:
+
+ METRICS_TO_EXCLUDE = [
+ "epoch/learning_rate",
+ "epoch/epoch",
+ "batch/learning_rate",
+ "batch/batch_step",
+ ]
+
+ # Wilhelmlab WandB account that has all VEGA presets required for the reports
+ VEGA_LITE_PRESETS_ID = "prosit-compms"
+
+ def __init__(
+ self,
+ project: str,
+ title: str,
+ description: str,
+ dataset: RetentionTimeDataset = None,
+ ):
+ """Create WandB report for comparing runs.
+
+ Args:
+ project (str): Name of the project to be used in wandb.
+ title (str): Title of the report in wandb.
+ description (str): Description of the report in wandb.
+ dataset (RetentionTimeDataset, optional): The retention time dataset if logging the data is desired. Defaults to None, no logging of input data.
+ """
+ self.project = project
+ self.title = title
+ self.description = description
+ self.dataset = dataset
+ self.entity = wandb.apis.PublicApi().default_entity
+ self.wandb_api = wandb.Api()
+ self.table_key_len = ""
+ self.table_key_rt = ""
+ self.model_info = []
+
+ def create_report(
+ self,
+ add_config_section=True,
+ add_data_section=True,
+ add_train_section=True,
+ add_val_section=True,
+ add_train_val_section=True,
+ add_model_section=True,
+ ):
+ """Create a report in wandb.
+
+ Args:
+ add_config_section (bool, optional): Add a section for config parameters and the run to the report. Defaults to True.
+ add_data_section (bool, optional): Add a section for input data to the report. Defaults to True.
+ add_train_section (bool, optional): Add a section for training metrics to the report. Defaults to True.
+ add_val_section (bool, optional): Add a section for validation metrics to the report. Defaults to True.
+ add_train_val_section (bool, optional): Add a section for train-val metrics to the report. Defaults to True.
+ add_model_section (bool, optional): Add a section for model summary and number of parameters to the report. Defaults to True.
+ """
+ report = wr.Report(
+ project=self.project, title=self.title, description=self.description
+ )
+
+ report.blocks = [wr.TableOfContents()]
+
+ if add_model_section:
+ report.blocks += self._build_model_section()
+ if add_config_section:
+ report.blocks += self._build_config_section()
+ if add_data_section and self.dataset is not None:
+ report.blocks += self._build_data_section()
+ if add_train_section:
+ report.blocks += self._build_train_section()
+ if add_val_section:
+ report.blocks += self._build_val_section()
+ if add_train_val_section:
+ report.blocks += self._build_train_val_section()
+
+ report.save()
+
+ # get metrics of last run in project or from specified run_id
+ def _get_metrics(self, run_id=None):
+ if run_id:
+ # run is specified by <entity>/<project>/<run_id>
+ run = self.wandb_api.run(path=f"{self.entity}/{self.project}/{run_id}")
+ metrics_dataframe = run.history()
+ return metrics_dataframe
+ else:
+ # get metrics of latest run
+ runs = self.wandb_api.runs(path=f"{self.entity}/{self.project}")
+ run = runs[0]
+ metrics_dataframe = run.history()
+ return metrics_dataframe
+
+ # get metric names split into train/val, train is further split into batch/epoch
+ def _get_metrics_names(self):
+ metrics = self._get_metrics()
+ # filter strings from list that are not starting with "_" and do not contain "val"
+ pre_filter = [string for string in metrics if not string.startswith("_")]
+ batch_train_metrics_names = [
+ string
+ for string in pre_filter
+ if ("val" not in string.lower())
+ & ("epoch" not in string.lower())
+ & ("table" not in string.lower())
+ ]
+ epoch_train_metrics_names = [
+ string
+ for string in pre_filter
+ if ("val" not in string.lower())
+ & ("batch" not in string.lower())
+ & ("table" not in string.lower())
+ ]
+ # filter strings from list that contain "val"
+ epoch_val_metrics_names = list(filter(lambda x: "val" in x.lower(), metrics))
+
+ # filter strings from train metrics that are 'epoch/learning_rate' and 'epoch/epoch'
+ strings_to_filter = RetentionTimeReportRunComparisonWandb.METRICS_TO_EXCLUDE
+ batch_train_metrics_names = [
+ string
+ for string in batch_train_metrics_names
+ if string not in strings_to_filter
+ ]
+ epoch_train_metrics_names = [
+ string
+ for string in epoch_train_metrics_names
+ if string not in strings_to_filter
+ ]
+ batch_train_metrics_names.sort()
+ epoch_train_metrics_names.sort()
+ return (
+ batch_train_metrics_names,
+ epoch_train_metrics_names,
+ epoch_val_metrics_names,
+ )
+
+ def get_train_val_metrics_names(self):
+ (
+ _,
+ epoch_train_metrics_names,
+ epoch_val_metrics_names,
+ ) = self._get_metrics_names()
+ epoch_train_metrics_names.sort()
+ epoch_val_metrics_names.sort()
+ return list(zip(epoch_train_metrics_names, epoch_val_metrics_names))
+
+ def _build_config_section(self):
+ config_block = [
+ wr.H1(text="Config"),
+ wr.PanelGrid(
+ runsets=[
+ wr.Runset(self.entity, self.project),
+ ],
+ panels=[wr.RunComparer(layout={"w": 24})],
+ ),
+ wr.HorizontalRule(),
+ ]
+ return config_block
+
+ def _build_data_section(self):
+ data_block = [
+ wr.H1(text="Data"),
+ wr.P(
+ "The following section is showing a simple explorative data analysis of the used dataset. The first histogram shows the distribution of peptide lengths in the data set, while the second histogram shows the distribution of indexed retention times."
+ ),
+ wr.PanelGrid(
+ runsets=[
+ wr.Runset(self.entity, self.project),
+ ],
+ panels=[
+ wr.CustomChart(
+ query={"summaryTable": {"tableKey": self.table_key_len}},
+ chart_name=f"{RetentionTimeReportRunComparisonWandb.VEGA_LITE_PRESETS_ID}/histogram_peptide_length",
+ chart_fields={"value": self.dataset.sequence_col},
+ ),
+ wr.CustomChart(
+ query={"summaryTable": {"tableKey": self.table_key_rt}},
+ chart_name=f"{RetentionTimeReportRunComparisonWandb.VEGA_LITE_PRESETS_ID}/histogram_irt",
+ chart_fields={"value": self.dataset.target_col},
+ ),
+ ],
+ ),
+ wr.HorizontalRule(),
+ ]
+ return data_block
+
+ def _build_train_section(self):
+ (
+ batch_train_metrics_names,
+ epoch_train_metrics_names,
+ _,
+ ) = self._get_metrics_names()
+ panel_list_batch = []
+ panel_list_epoch = []
+ for name in batch_train_metrics_names:
+ panel_list_batch.append(wr.LinePlot(x="Step", y=[name]))
+ for name in epoch_train_metrics_names:
+ panel_list_epoch.append(wr.LinePlot(x="Step", y=[name]))
+ train_block = [
+ wr.H1(text="Training metrics"),
+ wr.P(
+ "The following section shows the different metrics that were used to track the training. All used metrics are added by default. The first subsection shows the metrics per epoch, whereas the second subsection show the metrics per batch."
+ ),
+ wr.H2(text="per batch"),
+ wr.PanelGrid(
+ runsets=[
+ wr.Runset(self.entity, self.project),
+ ],
+ panels=panel_list_batch,
+ ),
+ wr.H2(text="per epoch"),
+ wr.PanelGrid(
+ runsets=[
+ wr.Runset(self.entity, self.project),
+ ],
+ panels=panel_list_epoch,
+ ),
+ wr.HorizontalRule(),
+ ]
+ return train_block
+
+ def _build_val_section(self):
+ _, _, epoch_val_metrics_names = self._get_metrics_names()
+ panel_list_epoch = []
+ for name in epoch_val_metrics_names:
+ panel_list_epoch.append(wr.LinePlot(x="Step", y=[name]))
+ val_block = [
+ wr.H1(text="Validation metrics"),
+ wr.P(
+ "The following section shows the different metrics that were used to track the validation. All used metrics are added by default. The metrics are shown per epoch."
+ ),
+ wr.H2(text="per epoch"),
+ wr.PanelGrid(
+ runsets=[
+ wr.Runset(self.entity, self.project),
+ ],
+ panels=panel_list_epoch,
+ ),
+ wr.HorizontalRule(),
+ ]
+ return val_block
+
+ def _build_model_section(self):
+ model_block = [
+ wr.H1(text="Model information"),
+ wr.P(
+ "The following section shows information about the model. The table below contains information about the models' layers."
+ ),
+ wr.UnorderedList(items=self.model_info),
+ wr.PanelGrid(
+ runsets=[
+ wr.Runset(self.entity, self.project),
+ ],
+ panels=[wr.WeavePanelSummaryTable("layer_table")],
+ ),
+ wr.HorizontalRule(),
+ ]
+ return model_block
+
+ def _build_train_val_section(self):
+ train_val_metrics_names = self.get_train_val_metrics_names()
+ panel_list_epoch = []
+ for name in train_val_metrics_names:
+ panel_list_epoch.append(wr.LinePlot(x="Step", y=list(name)))
+ train_val_block = [
+ wr.H1(text="Train - Validation metrics"),
+ wr.P(
+ "The following section shows the different metrics for both training and validation in comaprison. All used metrics are added by default. The metrics are shown per epoch."
+ ),
+ wr.H2(text="per epoch"),
+ wr.PanelGrid(
+ runsets=[
+ wr.Runset(self.entity, self.project),
+ ],
+ panels=panel_list_epoch,
+ ),
+ wr.HorizontalRule(),
+ ]
+ return train_val_block
+
+ def log_sequence_length_table(
+ self, data: pd.DataFrame, seq_col: str = "modified_sequence"
+ ):
+ """Log sequence length table to wandb
+
+ Args:
+ data (pd.DataFrame): input data
+ seq_col (str, optional): Name of the column containing the sequences in the data frame. Defaults to "modified_sequence".
+
+ Returns:
+ str: Name of the histogram created by wandb after logging the data.
+ """
+ name_hist = "counts_hist"
+ counts = self.count_seq_length(data, seq_col)
+ # convert to df for easier handling
+ counts_df = counts.to_frame()
+ table = wandb.Table(dataframe=counts_df)
+ # log to wandb
+ hist = wandb.plot_table(
+ vega_spec_name=f"{RetentionTimeReportRunComparisonWandb.VEGA_LITE_PRESETS_ID}/histogram_peptide_length",
+ data_table=table,
+ fields={"value": seq_col},
+ )
+ wandb.log({name_hist: hist})
+ name_hist_table = name_hist + "_table"
+ return name_hist_table
+
+ # function to count sequence length
+ def count_seq_length(self, data: pd.DataFrame, seq_col: str):
+ pattern = re.compile(r"\[UNIMOD:.*\]", re.IGNORECASE)
+ data[seq_col].replace(pattern, "", inplace=True)
+ return data[seq_col].str.len()
+
+ # function to log retention time table to wandb
+ def log_rt_table(self, data: pd.DataFrame, rt_col: str = "indexed_retention_time"):
+ name_hist = "rt_hist"
+ rt = data.loc[:, rt_col]
+ # convert to df for easier handling
+ rt_df = rt.to_frame()
+ table = wandb.Table(dataframe=rt_df)
+ # log to wandb
+ hist = wandb.plot_table(
+ vega_spec_name=f"{RetentionTimeReportRunComparisonWandb.VEGA_LITE_PRESETS_ID}/histogram_irt",
+ data_table=table,
+ fields={"value": rt_col},
+ )
+ wandb.log({name_hist: hist})
+ name_hist_table = name_hist + "_table"
+ return name_hist_table
+
+ def log_data(self):
+ # check if datasource is a string
+ if isinstance(self.dataset.data_source, str):
+ # read corresponding file
+ file_extension = os.path.splitext(self.dataset.data_source)[-1].lower()
+
+ if file_extension == ".csv":
+ data = pd.read_csv(self.dataset.data_source)
+ if file_extension == ".json":
+ data = pd.read_json(self.dataset.data_source)
+ if file_extension == ".parquet":
+ data = pd.read_parquet(self.dataset.data_source, engine="fastparquet")
+ self.table_key_len = self.log_sequence_length_table(
+ data, self.dataset.sequence_col
+ )
+ self.table_key_rt = self.log_rt_table(data, self.dataset.target_col)
+
+ # check if datasource is a tuple of two ndarrays or two lists
+ if (
+ isinstance(self.dataset.data_source, tuple)
+ and all(
+ isinstance(item, (np.ndarray, list))
+ for item in self.dataset.data_source
+ )
+ and len(self.dataset.data_source) == 2
+ ):
+ data = pd.DataFrame(
+ {
+ self.dataset.sequence_col: self.dataset.data_source[0],
+ self.dataset.target_col: self.dataset.data_source[1],
+ }
+ )
+ self.table_key_len = self.log_sequence_length_table(
+ data, self.dataset.sequence_col
+ )
+ self.table_key_rt = self.log_rt_table(data, self.dataset.target_col)
+
+ # check if datasource is a single ndarray or list
+ # does not work? maybe error in RetentionTimeDataset
+ if isinstance(self.dataset.data_source, (np.ndarray, list)):
+ data = pd.DataFrame({self.dataset.sequence_col: self.dataset.data_source})
+ self.table_key_len = self.log_sequence_length_table(
+ data, self.dataset.sequence_col
+ )
+
+ def log_model_data(self, model):
+ import io
+ model_summary_buffer = io.StringIO()
+ model.summary(print_fn=lambda x: model_summary_buffer.write(x + "<br>"))
+ model_summary_lines = model_summary_buffer.getvalue().split("<br>")
+
+ lines = [line.rstrip() for line in model_summary_lines]
+
+ # remove formatting lines
+ strings_to_remove = ["____", "===="]
+ cleaned_list = [
+ item
+ for item in lines
+ if not any(string in item for string in strings_to_remove)
+ ]
+
+ # split into words by splitting if there are more than two whitespaces
+ words = []
+ for line in cleaned_list:
+ words.append(re.split(r"\s{2,}", line))
+
+ # remove lines that contain less than 3 characters
+ filtered_list_of_lists = [
+ sublist for sublist in words if all(len(item) > 3 for item in sublist)
+ ]
+
+ # extract layer info and model info
+ layer_info = [sublist for sublist in filtered_list_of_lists if len(sublist) > 2]
+ model_info = [sublist for sublist in filtered_list_of_lists if len(sublist) < 2]
+
+ # flatten model_info and filter entries with length smaller than 5
+ model_info_flat = [item for sublist in model_info for item in sublist]
+ model_info_flat_filtered = [item for item in model_info_flat if len(item) >= 5]
+
+ # create layer_info_df
+ column_names = layer_info[0]
+ layer_info_df = pd.DataFrame(layer_info[1:], columns=column_names)
+
+ # log layer_table to wandb
+ layer_table = wandb.Table(dataframe=layer_info_df)
+ wandb.log({"layer_table": layer_table})
+
+ # attach model_info to object
+ self.model_info = model_info_flat_filtered
diff --git a/dlomix/reports/__init__.py b/dlomix/reports/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..8a71bc47eaf4ede089e5759380f7e81ea1989961
--- /dev/null
+++ b/dlomix/reports/__init__.py
@@ -0,0 +1,12 @@
+from .IntensityReport import IntensityReport
+from .RetentionTimeReport import RetentionTimeReport
+from .RetentionTimeReportModelComparisonWandb import (
+ RetentionTimeReportModelComparisonWandb,
+)
+from .RetentionTimeReportRunComparisonWandb import RetentionTimeReportRunComparisonWandb
+
+__all__ = ["RetentionTimeReport",
+ "IntensityReport",
+ "RetentionTimeReportRunComparisonWandb",
+ "RetentionTimeReportModelComparisonWandb",
+ ]
diff --git a/dlomix/reports/postprocessing.py b/dlomix/reports/postprocessing.py
new file mode 100644
index 0000000000000000000000000000000000000000..a75316fe28b09bae36dcf48b1ea142989351a4d6
--- /dev/null
+++ b/dlomix/reports/postprocessing.py
@@ -0,0 +1,101 @@
+import functools
+
+import numpy as np
+import tensorflow as tf
+
+from ..losses import masked_spectral_distance
+
+
+def reshape_dims(array):
+ n, dims = array.shape
+ assert dims == 174
+ nlosses = 1
+ return array.reshape([n, 30 - 1, 2, nlosses, 3])
+
+
+def reshape_flat(array):
+ s = array.shape
+ flat_dim = [s[0], functools.reduce(lambda x, y: x * y, s[1:], 1)]
+ return array.reshape(flat_dim)
+
+
+def normalize_base_peak(array):
+ # flat
+ maxima = array.max(axis=1)
+ array = array / maxima[:, np.newaxis]
+ return array
+
+
+def mask_outofrange(array, lengths, mask=-1.0):
+ # dim
+ for i in range(array.shape[0]):
+ array[i, lengths[i] - 1 :, :, :, :] = mask
+ return array
+
+
+def mask_outofcharge(array, charges, mask=-1.0):
+ # dim
+ for i in range(array.shape[0]):
+ if charges[i] < 3:
+ array[i, :, :, :, charges[i] :] = mask
+ return array
+
+
+def get_spectral_angle(true, pred, batch_size=600):
+ n = true.shape[0]
+ sa = np.zeros([n])
+
+ def iterate():
+ if n > batch_size:
+ for i in range(n // batch_size):
+ true_sample = true[i * batch_size : (i + 1) * batch_size]
+ pred_sample = pred[i * batch_size : (i + 1) * batch_size]
+ yield i, true_sample, pred_sample
+ i = n // batch_size
+ yield i, true[(i) * batch_size :], pred[(i) * batch_size :]
+ else:
+ yield 0, true, pred
+
+ for i, t_b, p_b in iterate():
+ tf.compat.v1.reset_default_graph()
+ with tf.compat.v1.Session() as s:
+ sa_graph = masked_spectral_distance(t_b, p_b)
+ sa_b = 1 - s.run(sa_graph)
+ sa[i * batch_size : i * batch_size + sa_b.shape[0]] = sa_b
+ sa = np.nan_to_num(sa)
+ return sa
+
+
+def normalize_intensity_predictions(data, batch_size=600):
+ assert (
+ "sequences" in data
+ ), "Key sequences is missing in the data provided for post-processing"
+ assert (
+ "intensities_pred" in data
+ ), "Key intensities_pred is missing in the data provided for post-processing"
+ assert (
+ "precursor_charge_onehot" in data
+ ), "Key precursor_charge_onehot is missing in the data provided for post-processing"
+
+ sequence_lengths = data["sequences"].apply(lambda x: len(x))
+ intensities = np.stack(data["intensities_pred"].to_numpy()).astype(np.float32)
+ precursor_charge_onehot = np.stack(data["precursor_charge_onehot"].to_numpy())
+ charges = list(precursor_charge_onehot.argmax(axis=1) + 1)
+
+ intensities[intensities < 0] = 0
+ intensities = reshape_dims(intensities)
+ intensities = mask_outofrange(intensities, sequence_lengths)
+ intensities = mask_outofcharge(intensities, charges)
+ intensities = reshape_flat(intensities)
+ m_idx = intensities == -1
+ intensities = normalize_base_peak(intensities)
+ intensities[m_idx] = -1
+ data["intensities_pred"] = intensities
+
+ if "intensities_raw" in data:
+ data["spectral_angle"] = get_spectral_angle(
+ np.stack(data["intensities_raw"].to_numpy()).astype(np.float32),
+ intensities,
+ batch_size=batch_size,
+ )
+ return data
diff --git a/dlomix/utils.py b/dlomix/utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..e2dbb7756ebcb8eda78537ea805cd3915b3cbe49
--- /dev/null
+++ b/dlomix/utils.py
@@ -0,0 +1,48 @@
+import pickle
+
+import numpy as np
+
+
+def save_obj(obj, name):
+ with open(name + ".pkl", "wb") as f:
+ pickle.dump(obj, f, pickle.HIGHEST_PROTOCOL)
+
+
+def load_obj(name):
+ with open(name + ".pkl", "rb") as f:
+ return pickle.load(f)
+
+
+def convert_nested_list_to_numpy_array(nested_list, dtype=np.float32):
+ return np.array([np.array(x, dtype=dtype) for x in nested_list])
+
+
+def lower_and_trim_strings(strings):
+ return [s.lower().trim() for s in strings]
+
+
+def get_constructor_call_object_creation(object_instance):
+ members = [
+ attr
+ for attr in vars(object_instance)
+ if not callable(getattr(object_instance, attr))
+ and not attr.startswith(("_", "__"))
+ ]
+ values = [object_instance.__getattribute__(m) for m in members]
+
+ repr_str = ", ".join([f"{m}={v}" for m, v in zip(members, values)])
+
+ return f"{object_instance.__class__.__name__}({repr_str})"
+
+
+def flatten_dict_for_values(d):
+ if not isinstance(d, dict):
+ return d
+ else:
+ items = []
+ for v in d.values():
+ if isinstance(v, dict):
+ return flatten_dict_for_values(v)
+ else:
+ items.append(v)
+ return items