OSCR

Neural posterior estimation for population genetics.

Code ↔ Paper

21 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 21 matches · 1 of them tie a paragraph to a whole file, not to given lines: a weak match, whose lines are not tinted
  1. [1] § Methods › Population genetic tasks › Comparisons to MSMC2 ↔ experiments/variable-popn-size/baselines/msmc2_utils.py, lines 279–351 · score 0.91 · EM iterations, segment pattern, diploid individuals, fixed recombination, mutation rate, tree sequences
  2. [2] § Methods › Population genetic tasks › Embedding networks and summary statistics ↔ workflow/scripts/ts_processors.py, lines 325–410 · score 0.89 · Rogers Huff, distance bin, SNP pairs, gap, logarithmically, r2
  3. [3] § Methods › Population genetic tasks › Embedding networks and summary statistics ↔ workflow/scripts/embedding_networks.py, lines 328–397 · score 0.75 · bi directional GRU, dropout, PyTorch, vector, concatenated, layer
  4. [4] § Methods › Population genetic tasks › Embedding networks and summary statistics ↔ workflow/scripts/embedding_networks.py, lines 71–179 · score 0.71 · permutation invariance, exchangeable CNN, haplotype, LD, SNPs, genotypes
  5. [5] § Results › Inference of population bottleneck parameters ↔ workflow/scripts/ts_simulators.py, lines 163–212 · score 0.70 · South Middle Atlas, Arabidopsis thaliana, demographic model, stdpopsim, event, population
  6. [6] § Methods › Population genetic tasks › Comparisons to ABC ↔ experiments/aratha-2epoch-extra/abc_utils.py, lines 274–334 · score 0.66 · closest simulations, ABC rejection, training simulations, posterior sample, quantile, NPE
  7. [7] § Results › Inference of population bottleneck parameters ↔ experiments/aratha-2epoch-extra/plot-posterior-surfaces.py, lines 20–73 · score 0.63 · npe cnn, npe rnn, npe sfs, curvature, surfaces, match
  8. [8] § Results › Inference of population bottleneck parameters ↔ experiments/aratha-2epoch-extra/plot-posterior-surfaces.py, lines 20–73 · score 0.63 · npe cnn, npe rnn, npe sfs, surface, rows, likelihood
  9. [9] § Methods › Population genetic tasks › Application to Drosophila melanogaster ↔ workflow/scripts/ts_simulators.py, lines 95–160 · score 0.62 · ancestral population, isolation, haploid, CO, FR, migration
  10. [10] § Results › D. melanogaster out-of-Africa model ↔ experiments/variable-popn-size/compare-posteriors.py, lines 325–370 · score 0.59 · marginal distributions, generations ago, credible intervals, posterior samples, model, population
  11. [11] § Results › D. melanogaster out-of-Africa model ↔ experiments/dependent-variable-popn-size/compare-posteriors.py, lines 333–377 · score 0.59 · marginal distributions, generations ago, credible intervals, posterior samples, model, population
  12. [12] § Results › Inference of historical population sizes in a one-population-model ↔ experiments/dependent-variable-popn-size/compare-posteriors.py, lines 488–543 · score 0.57 · credible interval width, CI width, geometric, recombination rate, log, models
  13. [13] § Results › Inference of historical population sizes in a one-population-model ↔ experiments/variable-popn-size/compare-posteriors.py, lines 467–522 · score 0.57 · credible interval width, CI width, geometric, recombination rate, log, models
  14. [14] § Results › D. melanogaster out-of-Africa model ↔ depr/amortized_dadi_workflow/scripts/dadi_simulators.py, lines 15–84 · score 0.56 · migration rates, ancestral population, demographic model, split, growth
  15. [15] § Results › Inference of historical population sizes in a one-population-model ↔ workflow/scripts/ts_simulators.py, lines 380–519 · score 0.56 · 100–100000, uniform distribution, tree sequences, log10, msprime, spaced
  16. [16] § Results › Inference of population bottleneck parameters ↔ depr/amortized_dadi_workflow/scripts/dadi_simulators.py, lines 86–125 · score 0.56 · Arabidopsis thaliana, catalog, demographic model, stdpopsim, event, population
  17. [17] § Methods › Population genetic tasks › Application to Drosophila melanogaster ↔ depr/amortized_dadi_workflow/scripts/dadi_simulators.py, lines 15–84 · score 0.54 · migration rates, ancestral population, split, growth, uniform, Demographic
  18. [18] § Methods › NPE workflow ↔ experiments/dependent-variable-popn-size/simulate-for-posterior.py, lines 98–220 · score 0.52 · simulating tree sequences, model parameters, Python, recombination rate, raw, flow
  19. [19] § Methods › NPE workflow ↔ experiments/variable-popn-size/simulate-for-posterior.py, lines 224–366 · score 0.52 · simulating tree sequences, model parameters, Python, recombination rate, raw, flow
  20. [20] § Results › Inference of historical population sizes in a one-population-model ↔ experiments/variable-popn-size/run-ne-comparison-three-scenarios.sh, the whole file · a weak match · score 0.51 · linear decline, linear growth, CNN, RNN, population, model
  21. [21] § Methods › Population genetic tasks › Application to Drosophila melanogaster ↔ workflow/scripts/ts_simulators.py, lines 95–160 · score 0.50 · mutation rate, recombination rate, Drosophila, melanogaster, uniform, demographic

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 519 lines · 20 KB · MIT · 4 matches

  1. ## ts_simulators outputs tree sequence. This cannot be used as a simulator for simulate_for_sbi!
  3. import tskit
  4. import msprime
  5. import demes
  6. import torch
  7. import numpy as np
  8. import stdpopsim
  9. from sbi.utils import BoxUniform
  12. class BaseSimulator:
  13. def __init__(self, config: dict, default: dict):
  14. for key in config:
  15. if key == "class_name": continue
  16. assert key in default, f"Option {key} not available for simulator"
  17. for key, default in default.items():
  18. setattr(self, key, config.get(key, default))
  21. class YRI_CEU(BaseSimulator):
  22. """
  23. Simulate a model defined in dadi's manual
  24. (https://dadi.readthedocs.io/en/latest/examples/YRI_CEU/YRI_CEU/).
  25. Ancestral population changes in size; splits into YRI/CEU; CEU population
  26. undergoes a bottleneck upon splitting and subsequently grows exponentially.
  27. There is continuous symmetric migration between YRI and CEU.
  28. """
  30. default_config = {
  31. # FIXED PARAMETERS
  32. "samples": {"YRI": 10, "CEU": 10},
  33. "sequence_length": 10e6,
  34. "recombination_rate": 1.5e-8,
  35. "mutation_rate": 1.5e-8,
  36. # RANDOM PARAMETERS (UNIFORM)
  37. "N_A": [1e2, 1e5],
  38. "N_YRI": [1e2, 1e5],
  39. "N_CEU_initial": [1e2, 1e5],
  40. "N_CEU_final": [1e2, 1e5],
  41. "M": [0, 5e-4],
  42. "Tp": [0, 6e4],
  43. "T": [0, 6e4],
  44. }
  46. def __init__(self, config: dict):
  47. super().__init__(config, self.default_config)
  48. self.parameters = ["N_A", "N_YRI", "N_CEU_initial", "N_CEU_final", "M", "Tp", "T"]
  49. self.prior = BoxUniform(
  50. low=torch.tensor([getattr(self, p)[0] for p in self.parameters]),
  51. high=torch.tensor([getattr(self, p)[1] for p in self.parameters]),
  52. )
  54. def __call__(self, seed: int = None) -> (tskit.TreeSequence, np.ndarray):
  56. torch.manual_seed(seed)
  57. theta = self.prior.sample().numpy()
  58. N_A, N_YRI, N_CEU_initial, N_CEU_final, M, Tp, T = theta
  60. graph = demes.Builder()
  61. graph.add_deme(
  62. "ancestral",
  63. epochs=[dict(start_size=N_A, end_time=Tp + T)]
  64. )
  65. graph.add_deme(
  66. "AMH",
  67. ancestors=["ancestral"],
  68. epochs=[dict(start_size=N_YRI, end_time=T)],
  69. )
  70. graph.add_deme(
  71. "CEU",
  72. ancestors=["AMH"],
  73. epochs=[dict(start_size=N_CEU_initial, end_size=N_CEU_final)],
  74. )
  75. graph.add_deme(
  76. "YRI",
  77. ancestors=["AMH"],
  78. epochs=[dict(start_size=N_YRI)],
  79. )
  80. graph.add_migration(demes=["CEU", "YRI"], rate=M)
  82. demog = msprime.Demography.from_demes(graph.resolve())
  83. ts = msprime.sim_ancestry(
  84. self.samples,
  85. demography=demog,
  86. sequence_length=self.sequence_length,
  87. recombination_rate=self.recombination_rate,
  88. random_seed=seed,
  89. )
  90. ts = msprime.sim_mutations(ts, rate=self.mutation_rate, random_seed=seed)
  92. return ts, theta
  95. class DroMel_CO_FR(BaseSimulator):
  96. """
  97. Simulate a two-population isolation-with-migration model for
  98. Drosophila melanogaster, Congolese and French populations,
  99. with growth in the French population
  100. """
  102. default_config = {
  103. # FIXED PARAMETERS
  104. "samples": {"CO": 10, "FR": 10},
  105. "sequence_length": 1e6,
  106. "recombination_rate": [0.5e-8, 3.0e-8],
  107. "mutation_rate": 5.49e-9,
  108. "haploid": True,
  109. # RANDOM PARAMETERS (UNIFORM IN LOGSPACE)
  110. "N_ANC": [3.5, 6.5], # log10 ancestral population size
  111. "N_CO": [3.5, 6.5], # log10 Congolese population size
  112. "N_FR0": [3.5, 6.5], # log10 French population size after split
  113. "N_FR1": [3.5, 6.5], # log10 French population size after growth
  114. "T": [3, 6], # log10 split time
  115. "m_CO_FR": [-8, -3], # log10 migration Congolese to French
  116. "m_FR_CO": [-8, -3], # log10 migration French to Congolese
  117. }
  119. def __init__(self, config: dict):
  120. super().__init__(config, self.default_config)
  121. self.parameters = ["N_ANC", "N_CO", "N_FR0", "N_FR1", "T", "m_CO_FR", "m_FR_CO"]
  122. self.prior = BoxUniform(
  123. low=torch.tensor([getattr(self, p)[0] for p in self.parameters]),
  124. high=torch.tensor([getattr(self, p)[1] for p in self.parameters]),
  125. )
  127. def __call__(self, seed: int = None) -> (tskit.TreeSequence, np.ndarray):
  128. torch.manual_seed(seed)
  129. theta = self.prior.sample().numpy()
  130. seeds = torch.randint(2 ** 32 - 1, (2, )).numpy()
  131. r = self.recombination_rate[0] + \
  132. (self.recombination_rate[1] - self.recombination_rate[0]) * \
  133. torch.rand(1).item()
  134. N_ANC, N_CO, N_FR0, N_FR1, T, m_CO_FR, m_FR_CO = (10 ** theta)
  135. G_FR = np.log(N_FR1 / N_FR0) / T
  137. demogr = msprime.Demography()
  138. demogr.add_population(name="CO", initial_size=N_CO)
  139. demogr.add_population(name="FR", initial_size=N_FR1, growth_rate=G_FR)
  140. demogr.add_population(name="ANC", initial_size=N_ANC)
  141. demogr.migration_matrix = np.array([[0, m_CO_FR, 0], [m_FR_CO, 0, 0], [0, 0, 0]])
  142. demogr.add_population_split(time=T, derived=["CO", "FR"], ancestral="ANC")
  143. samples = [
  144. msprime.SampleSet(n, population=p, ploidy=1 if self.haploid else 2)
  145. for p, n in self.samples.items()
  146. ]
  147. ts = msprime.sim_ancestry(
  148. samples,
  149. demography=demogr,
  150. sequence_length=self.sequence_length,
  151. recombination_rate=r,
  152. random_seed=seeds[0],
  153. )
  154. ts = msprime.sim_mutations(
  155. ts,
  156. rate=self.mutation_rate,
  157. random_seed=seeds[1],
  158. )
  160. return ts, theta
  163. class AraTha_2epoch(BaseSimulator):
  164. """
  165. Simulate the African2Epoch_1H18 model from stdpopsim for Arabidopsis thaliana.
  166. The model consists of a single population that undergoes a size change.
  167. """
  169. species = stdpopsim.get_species("AraTha")
  170. model = species.get_demographic_model("African2Epoch_1H18")
  171. default_config = {
  172. # FIXED PARAMETERS
  173. "samples": {"SouthMiddleAtlas": 10},
  174. "sequence_length": 10e6,
  175. # RANDOM PARAMETERS (UNIFORM)
  176. "nu": [0.01, 1], # Ratio of current to ancestral population size
  177. "T": [0.01, 1.5], # Time of size change (scaled)
  178. }
  180. def __init__(self, config: dict):
  181. super().__init__(config, self.default_config)
  182. self.parameters = ["nu", "T"]
  183. self.prior = BoxUniform(
  184. low=torch.tensor([getattr(self, p)[0] for p in self.parameters]),
  185. high=torch.tensor([getattr(self, p)[1] for p in self.parameters]),
  186. )
  188. def __call__(self, seed: int = None) -> (tskit.TreeSequence, np.ndarray):
  189. torch.manual_seed(seed)
  190. theta = self.prior.sample().numpy()
  191. nu, T = theta
  193. species = self.species
  194. contig = species.get_contig(
  195. length=self.sequence_length,
  196. )
  197. model = self.model
  199. # Scale the population size and time parameters
  200. N_A = model.model.events[0].initial_size # ancestral population size
  201. model.populations[0].initial_size = nu * N_A # current population size
  202. model.model.events[0].time = T * 2 * N_A # time of size change
  204. engine = stdpopsim.get_engine("msprime")
  205. ts = engine.simulate(
  206. model,
  207. contig,
  208. samples=self.samples,
  209. random_seed=seed
  210. )
  212. return ts, theta
  215. class VariablePopulationSize(BaseSimulator):
  216. """
  217. Simulate a model with varying population size across multiple time windows.
  218. The model consists of a single population that undergoes multiple size changes.
  219. """
  221. default_config = {
  222. # FIXED PARAMETERS
  223. "samples": {"pop": 10},
  224. "sequence_length": 10e6,
  225. "mutation_rate": 1.5e-8,
  226. "num_time_windows": 3,
  227. # RANDOM PARAMETERS (UNIFORM)
  228. "pop_sizes": [1e2, 1e5], # Range for population sizes (log10 space)
  229. "recomb_rate": [1e-9, 2e-8], # Range for recombination rate
  230. # TIME PARAMETERS
  231. "max_time": 100000, # Maximum time for population events
  232. "time_rate": 0.1, # Rate at which time changes across windows
  233. }
  235. def __init__(self, config: dict):
  236. super().__init__(config, self.default_config)
  237. # Set up parameters list
  238. self.parameters = [f"N_{i}" for i in range(self.num_time_windows)] + ["recomb_rate"]
  240. # Create parameter ranges in the same format as AraTha_2epoch
  241. # Population sizes (in log10 space)
  242. pop_size_ranges = [[np.log10(self.pop_sizes[0]), np.log10(self.pop_sizes[1])]
  243. for _ in range(self.num_time_windows)]
  244. # Add recombination rate range
  245. param_ranges = pop_size_ranges + [self.recomb_rate]
  247. # Set up prior using BoxUniform
  248. self.prior = BoxUniform(
  249. low=torch.tensor([r[0] for r in param_ranges]),
  250. high=torch.tensor([r[1] for r in param_ranges])
  251. )
  253. # Calculate fixed time points for population size changes
  254. self.change_times = self._calculate_change_times()
  256. def _calculate_change_times(self) -> np.ndarray:
  257. """Calculate the times at which population size changes occur using an exponential spacing."""
  258. #times = [(np.exp(np.log(1 + self.time_rate * self.max_time) * i /
  259. # (self.num_time_windows - 1)) - 1) / self.time_rate
  260. # for i in range(self.num_time_windows)]
  261. win = np.logspace(2, np.log10(self.max_time), self.num_time_windows)
  262. win[0] = 0
  263. times = win
  264. return np.around(times).astype(int)
  266. def __call__(self, seed: int = None) -> (tskit.TreeSequence, np.ndarray):
  267. if seed is not None:
  268. torch.manual_seed(seed)
  270. min_snps = 400
  271. max_attempts = 100
  272. attempt = 0
  274. while attempt < max_attempts:
  275. # Sample parameters directly from prior (like AraTha_2epoch)
  276. theta = self.prior.sample().numpy()
  278. # Convert population sizes from log10 space
  279. pop_sizes = 10 ** theta[:-1] # All but last element are population sizes
  280. recomb_rate = theta[-1] # Last element is recombination rate
  282. # Create demography
  283. demography = msprime.Demography()
  284. demography.add_population(name="pop0", initial_size=float(pop_sizes[0]))
  286. # Add population size changes at calculated time intervals
  287. for i in range(1, len(pop_sizes)):
  288. demography.add_population_parameters_change(
  289. time=self.change_times[i],
  290. initial_size=float(pop_sizes[i]),
  291. growth_rate=0,
  292. population="pop0"
  293. )
  295. # Simulate ancestry
  296. ts = msprime.sim_ancestry(
  297. samples={"pop0": self.samples["pop"]},
  298. demography=demography,
  299. sequence_length=self.sequence_length,
  300. recombination_rate=recomb_rate,
  301. random_seed=seed
  302. )
  304. # Add mutations
  305. ts = msprime.sim_mutations(ts, rate=self.mutation_rate, random_seed=seed)
  307. # Check if we have enough SNPs after MAF filtering
  308. geno = ts.genotype_matrix().T
  309. num_sample = geno.shape[0]
  310. if (geno==2).any():
  311. num_sample *= 2
  313. row_sum = np.sum(geno, axis=0)
  314. keep = np.logical_and.reduce([
  315. row_sum != 0,
  316. row_sum != num_sample,
  317. row_sum > num_sample * 0.05,
  318. num_sample - row_sum > num_sample * 0.05
  319. ])
  321. if np.sum(keep) >= min_snps:
  322. return ts, theta
  324. attempt += 1
  325. if seed is not None:
  326. seed += 1
  328. raise RuntimeError(f"Failed to generate tree sequence with at least {min_snps} SNPs after {max_attempts} attempts")
  331. class recombination_rate(BaseSimulator):
  332. """
  333. Simulate a one population model where recombination rate varies
  334. among replicates. The prior is a beta distribution shifted/scaled
  335. to a given interval (by default, a noninformative beta).
  336. """
  338. default_config = {
  339. # FIXED PARAMETERS
  340. "samples": {0: 10},
  341. "sequence_length": 1e6,
  342. "mutation_rate": 1.5e-8,
  343. "pop_size": 1e4,
  344. "mean_and_dispersion": [0.5, 0.5], # beta mean and dispersion
  345. # RANDOM PARAMETERS (UNIFORM)
  346. "recombination_rate": [0, 1e-8], # bounds on recombination rate
  347. }
  349. def __init__(self, config: dict):
  350. super().__init__(config, self.default_config)
  351. self.parameters = ["recombination_rate"]
  352. self.prior = BoxUniform(
  353. low=torch.tensor([getattr(self, p)[0] for p in self.parameters]),
  354. high=torch.tensor([getattr(self, p)[1] for p in self.parameters]),
  355. )
  358. def __call__(self, seed: int = None) -> (tskit.TreeSequence, np.ndarray):
  360. torch.manual_seed(seed)
  361. mean, dispersion = self.mean_and_dispersion
  362. alpha, beta = mean / dispersion, (1 - mean) / dispersion
  363. prior = torch.distributions.Beta(torch.FloatTensor([alpha]), torch.FloatTensor([beta]))
  364. low, high = self.recombination_rate
  365. theta = low + (high - low) * prior.sample().numpy()
  366. recombination_rate = theta.item()
  368. ts = msprime.sim_ancestry(
  369. self.samples,
  370. population_size=self.pop_size,
  371. sequence_length=self.sequence_length,
  372. recombination_rate=recombination_rate,
  373. random_seed=seed,
  374. discrete_genome=False, # don't want overlapping mutations
  375. )
  376. ts = msprime.sim_mutations(ts, rate=self.mutation_rate, random_seed=seed)
  378. return ts, theta
  380. class DependentVariablePopulationSize(BaseSimulator):
  381. """
  382. Simulate a population with variable population size across multiple time windows, with each
  383. population size dependent on the previous one.
  384. The model consists of a single population that undergoes multiple size changes.
  385. """
  386. default_config = {
  387. # FIXED PARAMETERS
  388. "samples": {"pop0": 25},
  389. "sequence_length": 2e6,
  390. "mutation_rate": 1e-8,
  391. "num_time_windows": 21,
  392. "maf": 0.05,
  393. # RANDOM PARAMETERS (UNIFORM)
  394. "pop_sizes": [1e2, 1e5], # Range for population sizes (log10 space)
  395. "pop_changes": [-1, 1], # Range for population size changes (* 10 ** beta)
  396. "recomb_rate": [1e-9, 1e-8], # Range for recombination rate
  397. # TIME PARAMETERS
  398. "max_time": 130000, # Maximum time for population events
  399. "time_rate": 0.1, # Rate at which time changes across windows
  400. }
  402. def __init__(self, config: dict):
  403. super().__init__(config, self.default_config)
  404. # Set up parameters list
  405. self.parameters = [f"N_{i}" for i in range(self.num_time_windows)] + ["recomb_rate"]
  407. # Create parameter ranges in the same format as AraTha_2epoch
  408. # Population sizes (in log10 space)
  409. pop_size_range = [[np.log10(self.pop_sizes[0]), np.log10(self.pop_sizes[1])]
  410. for _ in range(self.num_time_windows)]
  411. pop_change_ranges = [[self.pop_changes[0], self.pop_changes[1]]
  412. for _ in range(1,self.num_time_windows)]
  413. # Add recombination rate range
  414. param_ranges = pop_size_range + [self.recomb_rate]
  416. # Set up prior using BoxUniform
  417. self.prior = BoxUniform(
  418. low=torch.tensor([r[0] for r in param_ranges]),
  419. high=torch.tensor([r[1] for r in param_ranges])
  420. )
  421. self.beta_prior = BoxUniform(
  422. low=torch.tensor([b[0] for b in pop_change_ranges]),
  423. high=torch.tensor([b[1] for b in pop_change_ranges])
  424. )
  426. # Calculate fixed time points for population size changes
  427. self.change_times = self._calculate_change_times()
  429. def _calculate_change_times(self) -> np.ndarray:
  430. """Calculate the times at which population size changes occur using an exponential spacing."""
  431. times = [(np.exp(np.log(1 + self.time_rate * self.max_time) * i /
  432. (self.num_time_windows - 1)) - 1) / self.time_rate
  433. for i in range(self.num_time_windows)]
  434. return np.around(times).astype(int)
  436. def _generate_dependent_pop_sizes(self) -> np.ndarray:
  437. """
  438. Generate a sequence of population sizes where the first population size N_1
  439. is sampled from a uniform distribution corresponding to the most recent
  440. time window. The following population sizes are generated following
  441. N_i = N_{i-1} * 10 ^ β for i in [2,...,num_time_windows], unless N_i
  442. is outside of pop_ranges. If so N_i is set to the max/min population size
  443. """
  444. prior_sample = self.prior.sample().numpy()
  445. # only the first sample from prior will be used
  446. beta_sample = self.beta_prior.sample().numpy()
  447. modified_prior = prior_sample.copy()
  448. # The first value is uniformly sampled within the log10 bounds
  449. for i in range(self.num_time_windows-1):
  450. # For subsequent time windows, calculate the new value based on the previous one and beta
  451. new_value = modified_prior[i] + beta_sample[i]
  452. # If the new value is outside the bounds, set it to the max/min
  453. if new_value > np.log10(self.pop_sizes[1]):
  454. new_value = np.log10(self.pop_sizes[1])
  455. if new_value < np.log10(self.pop_sizes[0]):
  456. new_value = np.log10(self.pop_sizes[0])
  457. modified_prior[i+1] = new_value
  459. # Return the sampled prior and recombination rate
  460. return modified_prior
  462. def __call__(self, seed: int = None) -> (tskit.TreeSequence, np.ndarray):
  463. if seed is not None:
  464. torch.manual_seed(seed)
  466. min_snps = 400
  467. max_attempts = 100
  468. attempt = 0
  470. while attempt < max_attempts:
  471. # Sample parameters directly from prior (like AraTha_2epoch)
  472. theta = self._generate_dependent_pop_sizes()
  474. pop_sizes = 10**theta[:-1] # All but last element are for the population sizes
  475. recomb_rate = theta[-1] # Last element is recombination rate
  477. # Create demography
  478. demography = msprime.Demography()
  479. demography.add_population(name="pop0", initial_size=float(pop_sizes[0]))
  481. # Add population size changes at calculated time intervals
  482. for i in range(1, len(pop_sizes)):
  483. demography.add_population_parameters_change(
  484. time=self.change_times[i],
  485. initial_size=float(pop_sizes[i]),
  486. growth_rate=0,
  487. population="pop0"
  488. )
  490. # Simulate ancestry
  491. ts = msprime.sim_ancestry(
  492. samples=self.samples,
  493. demography=demography,
  494. sequence_length=self.sequence_length,
  495. recombination_rate=recomb_rate,
  496. random_seed=seed
  497. )
  499. # Add mutations
  500. ts = msprime.sim_mutations(ts, rate=self.mutation_rate, random_seed=seed)
  502. # Check if we have enough SNPs after MAF filtering
  503. geno = ts.genotype_matrix().T
  504. num_sample = ts.num_samples
  506. row_sum = np.sum(geno, axis=0)
  507. keep = np.logical_and.reduce([
  508. row_sum > num_sample * self.maf,
  509. num_sample - row_sum > num_sample * self.maf
  510. ])
  512. if np.sum(keep) >= min_snps:
  513. return ts, theta
  515. attempt += 1
  516. if seed is not None:
  517. seed += 1
  519. raise RuntimeError(f"Failed to generate tree sequence with at least {min_snps} SNPs after {max_attempts} attempts")

ts_simulators.py at commit e81dfe6, under MIT · at the source

Overview

Authors: Jiseon Min1, Yuxin Ning2, Nathaniel S Pope1, Franz Baumdicker2,3, Andrew D Kern1
ORCID iDs: Nathaniel S Pope, Franz Baumdicker, Andrew D Kern
  1. Institute of Ecology and Evolution, University of Oregon, Eugene, OR 97403, United States
  2. Institute for Bioinformatics and Medical Informatics (IBMI), University of Tübingen, Tübingen 72076, Germany
  3. Big Data Analytics in Bioinformatics, Justus-Liebig University, Gießen 35390, Germany
Institutions: University of Oregon (United States); University of Tübingen (Germany); Justus-Liebig-Universität Gießen (Germany)
Journal: Genetics, volume 233, issue 3, article iyag107
Dates: published online 6 July 2026; in print July 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1093/genetics/iyag107 · PMID 42032815 · PMCID PMC13334119 · OpenAlex W4417116550
Open access: hybrid, a free copy (OpenAlex)
Preprint: doi:10.64898/2025.12.01.691638
Status: code verified
Categories: drosophila (organism), computational (subfield)
Methods: Statistics, Preprocessing, Machine learning
Keywords: population genetics, deep learning, machine learning, Bayesian inference, simulation
MeSH: Genetics, Population*, Models, Genetic*, Neural Networks, Computer*, Animals, Bayes Theorem, Computer Simulation, Drosophila melanogaster, Machine Learning (* major topic)
Topic: Markov Chains and Monte Carlo Methods (Statistics and Probability, Mathematics), according to OpenAlex
Funding: NIH (R01HG012473, R01HG010774, R35GM148253); NIGMS NIH HHS (R35 GM148253); NIH HHS (R01HG012473, R01HG010774, R35GM148253); DFG; German Research Foundation; NHGRI NIH HHS (R01 HG010774, R01 HG012473); Deutsche Forschungsgemeinschaft
Citations: cited by 2 papers (Europe PMC); 74 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

Its files are read in the Code ↔ Paper reader above, with 21 matches between paragraphs and lines of code.

kr-colab/popgen-npe

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: e81dfe66d32ff194e42c3c68941e647f824e0a32, 22 June 2026
Languages: Python (52), Shell (9), Jupyter (1)
Size: 183 files, 62 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README, license file, environment (environment.yaml, pyproject.toml, docs/requirements.txt), tests, documentation, 1 notebook
Not found: CITATION.cff, continuous integration
Tools: NumPy (47 files), PyTorch (25 files), Matplotlib (22 files), seaborn (13 files), PyTorch Lightning (6 files), pandas (5 files), Snakemake (5 files), SciPy (3 files), PyMC (2 files), pysam (2 files), BCFtools (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
64 files
At the source: github.com/kr-colab/popgen-npe

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 62 scripts, each with its path and the digest of its content;
  • 21 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Code and data availability statement

The paper has a code and data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1093/genetics/iyag107.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

  • Publisher: n/a → Oxford University Press

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 5 keywords, 8 MeSH terms, 7 funders, 61 references.

Cite

This paper

Min, J., Ning, Y., Pope, N. S., Baumdicker, F., & Kern, A. D. (2026). Neural posterior estimation for population genetics. Genetics, 233(3), iyag107. https://doi.org/10.1093/genetics/iyag107

BibTeX

@article{min2026neural,
author = {Min, Jiseon and Ning, Yuxin and Pope, Nathaniel S and Baumdicker, Franz and Kern, Andrew D},
title = {{Neural posterior estimation for population genetics}},
journal = {Genetics},
year = {2026},
month = jul,
volume = {233},
number = {3},
pages = {iyag107},
publisher = {Oxford University Press},
issn = {0016-6731},
doi = {10.1093/genetics/iyag107},
url = {https://doi.org/10.1093/genetics/iyag107},
pmid = {42032815},
pmcid = {PMC13334119}
}

RIS

TY - JOUR
AU - Min, Jiseon
AU - Ning, Yuxin
AU - Pope, Nathaniel S
AU - Baumdicker, Franz
AU - Kern, Andrew D
TI - Neural posterior estimation for population genetics
T2 - Genetics
J2 - Genetics
PY - 2026
DA - 2026/07/01
VL - 233
IS - 3
SP - iyag107
SN - 0016-6731
PB - Oxford University Press
DO - 10.1093/genetics/iyag107
UR - https://doi.org/10.1093/genetics/iyag107
LA - en
ER -

CSL-JSON

{
"id": "10.1093/genetics/iyag107",
"type": "article-journal",
"title": "Neural posterior estimation for population genetics",
"container-title": "Genetics",
"author": [
{
"family": "Min",
"given": "Jiseon"
},
{
"family": "Ning",
"given": "Yuxin"
},
{
"family": "Pope",
"given": "Nathaniel S"
},
{
"family": "Baumdicker",
"given": "Franz"
},
{
"family": "Kern",
"given": "Andrew D"
}
],
"container-title-short": "Genetics",
"volume": "233",
"issue": "3",
"page": "iyag107",
"DOI": "10.1093/genetics/iyag107",
"PMID": "42032815",
"PMCID": "PMC13334119",
"ISSN": "0016-6731",
"publisher": "Oxford University Press",
"URL": "https://doi.org/10.1093/genetics/iyag107",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
1
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s44400-026-00094-8 [code]
Haplotype-resolved DNA methylation at the &lt;i&gt;APOE&lt;/i&gt; locus identifies allele-specific epigenetic signatures relevant to Alzheimer's disease risk.
Journal: NPJ dementia
In common: Snakemake, BCFtools, pysam, 6 other tools
[2] doi:10.1016/j.celrep.2026.117110 [code]
Single-nucleus multiome analysis in the human prefrontal cortex identifies gene expression and cis-regulatory elements associated with aging.
Journal: Cell reports
In common: Snakemake, BCFtools, pysam, 6 other tools
[3] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: PyMC, BCFtools, PyTorch Lightning, 6 other tools
[4] doi:10.1093/evlett/qrag004 [code]
AI solutions for evolutionary genomics of nonmodel species.
Journal: Evolution letters
In common: NumPy, 6 references
[5] doi:10.1038/s41467-026-71790-5 [code]
Recurrent DNA break clusters drive replication-stress-induced copy number variants and genome diversification.
Journal: Nature communications
In common: Snakemake, BCFtools, pysam, 5 other tools
[6] doi:10.1016/j.xcrm.2026.102904 [code]
High-dose furmonertinib as first-line treatment for untreated EGFR-mutated advanced NSCLC with central nervous system metastases: A phase 2 trial.
Journal: Cell reports. Medicine
In common: PyMC, pysam, PyTorch Lightning, 4 other tools
[7] doi:10.1371/journal.pcbi.1014364 [code]
A comparative study of simulation-based inference methods for epidemic models with identifiability considerations.
Journal: PLoS computational biology
In common: PyTorch, pandas, SciPy, 2 other tools, computational, 3 references
[8] doi:10.1093/bioinformatics/btag652 [code]
mmVelo: a deep generative model for estimating cell state-dependent dynamics across multiple modalities.
Journal: Bioinformatics (Oxford, England)
In common: pysam, PyTorch Lightning, PyTorch, 5 other tools
[9] doi:10.1186/s13059-026-04177-w [code]
Genomic sequence evolution underlying human neocortical interareal diversification.
Journal: Genome biology
In common: Snakemake, pysam, seaborn, 4 other tools
[10] doi:10.1016/j.celrep.2026.117073 [code]
Single-cell epigenomics uncovers heterochromatin instability and transcription factor dysfunction during mouse brain aging.
Journal: Cell reports
In common: Snakemake, pysam, seaborn, 4 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.