Réunion 18/02/26

Lionel Viguier

Ordre du jour

  • Previously
  • Popfreq
  • ReadVCF
    • Parallélisation
  • Soutenance de stage
  • Discussion

Previously:

  • Changer le nom des individus générer (quand un des deux parents est absent du pédigree) nouvel id: indiv.indiv + “_fa”. ✔️
  • Changer l’issue 11: BlackList -> Pruning
    • tester le pruning avec le pédigree complet d’Elise
  • Vérifier le calcul des frequences allelique.
    • on peut regarder l’évolution de l’estimation de popFreq si on a la vrais popFreq
  • Mettre a jours les Anterior des Founders a chaque Cycle ie recalculer les frequences allelique ✔️

Popfreq

La fonction ReEstimateFoundersAnterior a pour but de recalculer la PopFreq a chaque fin de cycle et de l’associer aux anterior des fondateurs.

def ReEstimateFoundersAnterior(self):
  fndrIdxList = [self.d.ped.idx_dict[i.indiv] for i in self.d.ped.founders]
  fndrProbalog = self.d.peeling[..., fndrIdxList] - self.d.L[..., fndrIdxList]
  # Estimation of MAF
  estimateMAF = fndrProbalog * np.array([0, 0.5, 1])[None, :, None]
  estimateMAF = np.mean(estimateMAF, axis=(1, 2))

  log_1_minus_MAF = np.log(1 - np.exp(estimateMAF.clip(max=np.log(1 - 1e-10))))
  # Reestimate founders anterior probabilities
  reestimated_anterior_freq = np.vstack(
      (
          2 * log_1_minus_MAF,  # log((1-MAF)^2)
          np.log(2) + estimateMAF + log_1_minus_MAF,  # log(2 * MAF * (1-MAF))
          2 * estimateMAF,  # log(MAF^2)
      )
  )
  self.d.anterior[..., self.d.founder_idx] = reestimated_anterior_freq.reshape(
      reestimated_anterior_freq.shape[::-1] + (1,)
  )
  1. On calcule la MAF(en réalité c’est la fréquence allèlique de l’allele alternatif)
  2. A partir de la “MAF” on calcule les anteriors des fondateurs: \(p^2,2pq,q^2\)
  3. On associe cette mesure a tous les anteriors des fondateurs

ReadVCF

La procédure de lecture de fichier .vcf est séparer en 3 fonctions:

  • Une de lecture du vcf au global, et un partitionnement du traitement (chunksize)
  • Une lecture des variant de chaque sous région et écriture dans une archive zarr
  • Une de lecture de l’archive zarr en vue de créer la np.Array de Penetrance

Le partitionnement se fait par région, ce qui signifie que pour une région donné, on manipule des np.Array de forme (nb_variant_de_région,nbGénotype = 3, nb_individus_de_ped).

L’architecture de l’archive zarr (qui peut être discuter) est : “variants/{region}/LL”. Avec LL pour logarithm likelihood.

ReadVCF

ReadVCF

Fonction de lecture du vcf au global, et un partitionnement du traitement: découpage des régions en chunks.

def ReadVCF(
  fname: str,
  samples: list,
  idx_dict: dict = None,
  mode: str = "likelihood",
  baseValueLikelihood: float = 0.0,
  MissingGenotypeValue: int = -1,
  chunksize=-1,
):
  
  assert mode in [
      "likelihood",
      "genotype",
  ], "mode can only be likelihood or genotype."

  if idx_dict is None:
      idx_dict = {}
      for idx, id in enumerate(samples):
          idx_dict[id] = idx

  output_prefix = fname.split(".")[0]
  fzname = output_prefix + ".zarr"
  fz = zarr.open(fzname, "w")
  fz.attrs["archive"] = fzname

  v = vcf.VCF(fname, gts012=True, strict_gt=True, samples=samples, lazy=True)
  regions = vcf.vcf_chunk_regions(fname, chunksize)
  missigSamples = len(samples) - len(v.samples)
  if missigSamples > 0:
      logger.warning(missigSamples, "samples are missing in the .vcf file.")

  if mode == "likelihood":
      for r in regions:
          VCF_reg_to_Likelihood(
              v, fz, r, samples, idx_dict, baseValueLikelihood, MissingGenotypeValue
          )
  if mode == "genotype":
      logger.error("'genotype' mode not implemented yet")
      # Work in progress

  return fz

ReadVCF

VCF_reg_to_Likelihood

Fonction lecture des variant de chaque sous région et écriture du likelihood, genotype et ID dans une archive zarr.

def VCF_reg_to_Likelihood(
   v,
   zarrFile,
   region: str,
   samples: list,
   idx_dict,
   baseValueLikelihood: float = 0.0,
   MissingGenotypeValue: int = -1,
):
   samples_idx = [idx_dict[key] for key in v.samples]

   # Get Sample ID
   try:
       ID = sorted(samples, key=lambda x: idx_dict[x])
   except KeyError:
       raise Exception(
           "at least one sample's individual is missing in idx_dict"
       )  # ici orthographe!

   likelihood_all_var = []
   genotypes_all_var = []

   # Read from .vcf
   for s in v(region):

       # Get Likelihood
       likeli = np.full((3, len(samples)), baseValueLikelihood)
       likeli[0, samples_idx] = s.gt_phred_ll_homref
       likeli[1, samples_idx] = s.gt_phred_ll_het
       likeli[2, samples_idx] = s.gt_phred_ll_homalt
       likelihood_all_var.append(vcf.phred2ln(likeli))
       # Get Genotype
       genotypes = np.full((len(samples)), MissingGenotypeValue)
       genotypes_np = np.sum(
           np.array(s.genotypes)[..., 0:2], axis=1
       )  # transorm "[0, 1, False]" into "1"
       genotypes_np[genotypes_np < 0] = (
           MissingGenotypeValue  # missing value is "[-1, -1, False]" so we want to encode MissingGenotypeValue insted of -2
       )

       genotypes[samples_idx] = genotypes_np
       genotypes_all_var.append(genotypes)

   # Write in .zarr
   try:
       # try if the region got variants
       _ = genotypes_all_var[0]

       zarrFile[f"variants/{region}/ID"] = ID
       zarrFile[f"variants/{region}/LL"] = np.stack(likelihood_all_var, axis=0)
       zarrFile[f"variants/{region}/GT"] = np.stack(genotypes_all_var, axis=0)

   except IndexError:
       pass

ReadVCF

VCF_reg_to_Likelihood

Une fonction de lecture de l’archive zarr en vue de créer la np.Array de Penetrance directement utilisable dans l’algo de peeling.

def zarr_to_genoArray(
    fz, region: str
):  # Evolution possible: region est une liste de str
    return np.array(fz[f"variants/{region}/LL"])

Temp de calcul et parallélisation L’idée d’utiliser des archives zarr émergeai à la base de l’incapacité de l’ordi du stagiere à lire un fichier .vcf de test (~768 000 loci x 1047 individus). En effet la mémoire RAM de 16G étais saturée. L’archive zarr permet d’écrire directement sur le disque et donc de pallié ce problème.

L’idée ici est de lire les variants régions par région avec la fonction dédiée VCF_reg_to_Likelihood puis de stocker les données par région dans une archive zarr. Permettant ainsi la parallélisation à la lecture du vcf et au Peeling.

Soutenance de stage

Dates: 📆

  • rendu du rapport le 29 mai
  • soutenance prévue le 9 ou 10 juin en salle de conférences de l’IBCG.

⚠️ choisir un rapporteur ⚠️

  • Doit être sur Toulouse.
  • Peut-être
    • membre de l’équipe enseignante le plus proche du sujet
    • un chercheur ou enseignant-chercheur, ingénieur, hors équipe enseignante.
  • Doit être proche du sujet pour pouvoir juger au mieux le travail.
  • Ne doit pas collaborer directement avec vous sur le sujet et ne doit pas assister aux réunions de travail ni répétitions du stagiaire.

Discussion

  1. Comment mesurer l’écart entre ma popFreq estimée et la vrais popFreq
  2. Architecture et noms de l’archive zarr
  3. Vincent me manque déjà… que faire? 🤔

Retour:

  • [X]Modification de la structure du zarr afin de matcher la structure de yapp:
    • Stocker genotype dans /genotype/{chr}/{chunk}/GT
    • Stocker ID dans /samples/ID
    • Stocker la log vraisemblance dans /penetrance/{chr}/{chunk}/LL
  • [X]Modification sur la fonction ReEstimateFoundersAnterior:
    • [X]Refaire le calcul des MAF (OSKOUR!)
    • [X]Modifier la fonction pour qu’elle ne dépende plus de la classe “IterativePeeling”
      • ça permet de la tester plus facilement.
  • [X]Pour le Pruning on peut avoir la liste des individus génotypé dans le VCF sans avoir à lire les génotypes.
  • Pour les rapporteurs, les noms qui ont été mentionnés sont : Jean-Michel, Tom, Isabelle. Aucun ne semblait faire consensus.
  • [X]Ne pas utiliser de idx_dict dans la fonction readVCF -> plutot utiliser une liste de sample trié

Point hébdomadaire semaine 7.