aegis_sim.submodels.genetics.modifying.gpm

  1import logging
  2import numpy as np
  3from aegis_sim import parameterization
  4from numba import njit
  5
  6
  7class GPM:
  8    """Genotype-phenotype map
  9
 10    Order of elements in the vector does not matter. # TODO Explain better
 11
 12    ### GENOTYPE-PHENOTYPE MAP (GPM) ###
 13    In AEGIS, every individual carries a genome which encodes an intrinsic phenotype.
 14    A genome can be converted into an intrinsic phenotype using the genotype-phenotype map (GPM).
 15    Conceptually, the GPM contains the information on how each site affects the intrinsic phenotype
 16    of the individual (e.g. the first site decreases fertility by 0.15% at age class 28).
 17
 18    The GPM can be saved in two formats: a list or a matrix.
 19
 20    If it is a list, it will be a list of quadruple (4-tuple) with the following structure: `index`, `trait`, `age`, `magnitude`.
 21    Thus a single quadruple encodes an effect of a single site at the index `index` (e.g. 1)
 22    on the trait `trait` (e.g. fertility) expressed at the age `age` (e.g. 28). The change to the trait is of magnitude `magnitude` (0.85).
 23    When a site is pleiotropic, there will be multiple quadruples with the same `index`.
 24    We distringuish between age-pleiotropy (a single site affecting at least one trait at multiple ages) and trait-pleiotropy (a single site affecting multiple traits).
 25
 26    If the GPM is encoded as a matrix, it is a 3D matrix where dimensions encode `index`, `trait` and `age`,
 27    while the matrix values encode the `magnitude`s.
 28
 29    When most sites are age-pleiotropic and trait-pleiotropic, the optimal encoding format is a matrix.
 30    When most sites have age-specific and trait-specific effects, the optimal encoding format is a list
 31    rather than a matrix because the matrix will be very sparse (it will carry a lot of 0's).
 32    """
 33
 34    def __init__(self, phenomatrix, phenolist):
 35        self.phenomatrix = phenomatrix
 36        self.phenolist = phenolist
 37
 38        self.dummy = len(self.phenolist) == 0 and self.phenomatrix is None
 39        if self.dummy:
 40            logging.info("Phenomap inactive.")
 41
 42    def phenodiff(self, vectors, zeropheno):
 43        """
 44        vectors .. haploidized genomes of all individuals; shape is (n_individuals, ?)
 45        phenomatrix ..
 46        phenolist .. list of (bit_index, trait, age, magnitude)
 47        zeropheno .. phenotypes of all-zero genomes (i.e. how phenotypes would be if all bits were 0)
 48        """
 49
 50        if self.phenomatrix is not None:
 51            # TODO BUG resolve phenomatrix
 52            return vectors.dot(self.phenomatrix)
 53
 54        elif self.phenolist is not None:
 55            phenodiff = zeropheno.copy()
 56            for vec_index, trait, age, magnitude in self.phenolist:
 57                vec_state = vectors[:, vec_index]
 58                phenotype_change = vec_state * magnitude
 59                phenotype_index = parameterization.traits[trait].start + age
 60                phenodiff[:, phenotype_index] += phenotype_change
 61            return phenodiff
 62
 63        else:
 64            raise Exception("Neither phenomatrix nor phenolist has been provided.")
 65
 66    def phenodiff_accelerated(self, vectors, zeropheno):
 67        if self.phenomatrix is not None:
 68            return vectors.dot(self.phenomatrix)
 69
 70        elif self.phenolist is not None:
 71            phenodiff = zeropheno.copy()
 72            vec_indices, traits, ages, magnitudes = zip(*self.phenolist)
 73            vec_indices = np.array(vec_indices)
 74            magnitudes = np.array(magnitudes)
 75            vec_states = vectors[:, vec_indices]
 76            phenotype_indices = np.array(
 77                [parameterization.traits[trait].start + age for trait, age in zip(traits, ages)]
 78            )
 79            phenodiff = apply_phenolist_numba(phenodiff, vec_states, phenotype_indices, magnitudes)
 80            return phenodiff
 81
 82        else:
 83            raise Exception("Neither phenomatrix nor phenolist has been provided.")
 84
 85    def __call__(self, interpretome, zeropheno):
 86        if self.dummy:
 87            return zeropheno
 88        else:
 89            # phenodiff_old = self.phenodiff(vectors=interpretome, zeropheno=zeropheno)
 90            phenodiff = self.phenodiff_accelerated(vectors=interpretome, zeropheno=zeropheno)
 91            return phenodiff
 92
 93
 94@njit
 95def apply_phenolist_numba(phenodiff, vec_states, phenotype_indices, magnitudes):
 96    """
 97    phenodiff .. modification to phenotype vector for a specific trait; shape is (individual, phenotypic site)
 98    vec_states .. state of each genomic site for each individual; shape is (individual, genomic site)
 99    phenotype_indices .. position in the phenotype vector that is modified by the genomic site; shape is (genomic site)
100    magnitudes .. effect size for each site; shape is (genomic site)
101    """
102    n_individuals, n_phenos = vec_states.shape
103    for i in range(n_phenos):
104        for j in range(n_individuals):
105            phenodiff[j, phenotype_indices[i]] += vec_states[j, i] * magnitudes[i]
106    return phenodiff