Skip to content

Conversation

@andrewkern
Copy link
Member

PR for the SPIDNA network and associated simulator / processor for genotypes. Also includes testing.

Here are representative results from a 3 epoch model where I trained the embedder and normalizing flow separately.

Posterior at prior high:
image

Posterior at prior low:
image

Posterior at prior mean:
image

Calibration (looks great!):
image

Concentration:
image

Posterior Expectation:
image

Overall-- I think this is ready to go! We need to decide on parameterization for a 'production' run still, but this run was done with the parameters specified in workflow/config/variable_popnSize_spidna.yaml

@andrewkern andrewkern requested a review from nspope February 13, 2025 19:07
@andrewkern
Copy link
Member Author

one decision I punted on here: what to do about the summary stats embeddings that @ningyuxin1999 did

@fbaumdicker
Copy link
Collaborator

I just discussed this with @ningyuxin1999 today. We can write the summary stats processor for the summary statistics we used, a combination of the SFS and mean LD in different distance bins.

@andrewkern
Copy link
Member Author

the big question to me here is do we do mean aggregation of those summaries as you were doing originally, and if so how given our workflow.

a few alternatives to consider:

  1. we could not do mean aggregation and instead just get one summary stat vector (LD+SFS stats)
  2. we could simulate a very long stretch of sequence, chop it in to windows, calculate stats from each window, and then collapse to means
  3. we could write a specialized simulator and associate processor that for a given parameter set simulates multiple reps (as you are doing) and then does summary stats / mean aggregation

Copy link
Contributor

@nspope nspope left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM!

# MAF filtering
if self.maf != 0:
num_sample = snp.shape[1] # Now using shape[1] since matrix isn't transposed
if (snp==2).any():
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is detecting ploidy, right? It's a little messy to do this on the fly (e.g. would fail if there were no homozygous derived alleles, which probably won't happen in practice, but still ...)

A better way is to do this directly from the individuals in the ts, e.g.

ploidy = [ind.nodes.size for ind in ts.individuals()]

or if we only need the number of haploids

samples = ts.num_samples

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this MAF filtering step was part of Yuxin's pipeline that I just copied over. I agree it's a bit messy to do it at this step. let me push a change to use samples from ts.num_samples.

@andrewkern
Copy link
Member Author

okay pushed a change @nspope and it still passes tests

@nspope nspope merged commit d69877c into main Feb 14, 2025
1 check passed
@andrewkern andrewkern deleted the spidna branch February 15, 2025 01:06
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants