Skip to content

Commit

Permalink
fix crosses
Browse files Browse the repository at this point in the history
  • Loading branch information
@cjGO
cjGO committed Jun 9, 2024
1 parent a4a95f4 commit a8980cd
Show file tree
Hide file tree
Showing 12 changed files with 322 additions and 275 deletions.
50 changes: 30 additions & 20 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,8 @@ pip install chewc
First, define the genome of your crop

``` python
import random

ploidy = 2
number_chromosomes = 10
loci_per_chromosome = 100
Expand All @@ -51,38 +53,46 @@ simparam.genome = crop_genome


#add a single additive trait
qtl_loci = 100
qtl_loci = 20
qtl_map = select_qtl_loci(qtl_loci,simparam.genome)

ta = TraitA(qtl_map,simparam,0, 1)
ta.sample_initial_effects()
ta.scale_genetic_effects()
ta.calculate_intercept()

# Create a 2x2 grid of plots
fig, axs = plt.subplots(2, 2, figsize=(10, 10)) # Adjust figsize as needed

# Plot 1: Scaled Effects Histogram
axs[0, 0].hist(ta.scaled_effects.flatten())
axs[0, 0].set_title('Marker Effects')

# Plot 2: True Genetic Values Histogram
genetic_values = ta.calculate_genetic_values(simparam.founder_pop)
axs[0, 1].hist(genetic_values)
axs[0, 1].set_title('True Genetic Values')

# Plot 3: Phenotype Histogram
phenotypes = ta.phenotype(simparam.founder_pop, h2=0.5) # Assuming 'h2' is a parameter for heritability in the phenotype method
axs[1, 0].hist(phenotypes)
axs[1, 0].set_title('Phenotype')

# Plot 4: Scatter Plot of Genetic Values vs Phenotypes
axs[1, 1].scatter(genetic_values, phenotypes)
axs[1, 1].set_title('Genetic Values vs Phenotype')
years = 50
current_pop = founder_pop
pmean = []
pvar = []
for i in range(years):
#phenotype current pop
TOPK = 2
new_pop=[]
pheno = ta.phenotype(current_pop,h2=.14)
topk = torch.topk(pheno,TOPK).indices

for i in range(200):
sampled_indices = torch.multinomial(torch.ones(topk.size(0)), 2, replacement=False)
sampled_parents = topk[sampled_indices]
m,f = current_pop[sampled_parents[0]], current_pop[sampled_parents[1]]
new_pop.append(make_cross(simparam.genome.genetic_map, m, f))

current_pop = torch.stack(new_pop)
pmean.append(ta.calculate_genetic_values(current_pop).mean())
pvar.append(ta.calculate_genetic_values(current_pop).var())


pmean_normalized = torch.tensor(pmean) / max(pmean)

pvar_normalized = torch.tensor(pvar) / max(pvar)

# Display the plots
plt.tight_layout() # Adjust layout to prevent overlap
plt.show()
plt.scatter(range(len(pmean_normalized)), pmean_normalized)
plt.scatter(range(len(pvar_normalized)), pvar_normalized)
```

![](index_files/figure-commonmark/cell-4-output-1.png)
1 change: 0 additions & 1 deletion chewc/_modidx.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,6 @@
'chewc.core.create_uniform_genetic_map': ('core.html#create_uniform_genetic_map', 'chewc/core.py')},
'chewc.crossing': { 'chewc.crossing.make_DH': ('crossing.html#make_dh', 'chewc/crossing.py'),
'chewc.crossing.make_cross': ('crossing.html#make_cross', 'chewc/crossing.py'),
'chewc.crossing.random_crosses': ('crossing.html#random_crosses', 'chewc/crossing.py'),
'chewc.crossing.truncate_select': ('crossing.html#truncate_select', 'chewc/crossing.py')},
'chewc.meiosis': { 'chewc.meiosis.gamma_interference_model': ('meiosis.html#gamma_interference_model', 'chewc/meiosis.py'),
'chewc.meiosis.simulate_gametes': ('meiosis.html#simulate_gametes', 'chewc/meiosis.py'),
Expand Down
59 changes: 2 additions & 57 deletions chewc/crossing.py
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,14 @@
# AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/04_crossing.ipynb.

# %% auto 0
__all__ = ['make_DH', 'make_cross', 'random_crosses', 'truncate_select']
__all__ = ['make_DH', 'make_cross', 'truncate_select']

# %% ../nbs/04_crossing.ipynb 4
from .core import *
from .meiosis import *
from .trait import *
import torch

import random
def make_DH(genetic_map, individual):
haploid_line = simulate_gametes(genetic_map, individual)
return gamete.repeat(2, 1, 1)
Expand All @@ -20,60 +20,5 @@ def make_cross(genetic_map, mother, father):
return torch.stack((egg, pollen), dim=0).squeeze(1)

# %% ../nbs/04_crossing.ipynb 6
def random_crosses(population: Population,
n_crosses: int,
) -> Population:
"""
Performs random crosses within a population.
Args:
----
population (Population): The population to perform crosses within.
n_crosses (int): The number of crosses to perform.
Returns:
-------
Population: A new population of offspring resulting from the crosses.
"""
# make sure there are enough individuals for the requested number of crosses.
if n_crosses > len(population.individuals) // 2:
raise ValueError(f"Not enough individuals in the population ({len(population.individuals)}) for {n_crosses} crosses. Need at least {2*n_crosses}.")

# sample parents randomly
selected_parents = random.sample(population.individuals, 2 * n_crosses)

# create offspring
offspring = []
for i in range(n_crosses):
mother = selected_parents[i * 2]
father = selected_parents[i * 2 + 1]
haplotypes = make_cross(population.individuals[0].genome.genetic_map, mother, father)
# genetic_value = trait.calculate_genetic_value(haplotypes.unsqueeze(0)).squeeze(0) # calculate GV
offspring.append(
Individual(
genome = mother.genome,
haplotypes=haplotypes,
mother_id=mother.id,
father_id=father.id,
# genetic_value=genetic_value,
)
)
return Population(individuals=offspring)

# %% ../nbs/04_crossing.ipynb 7
def random_crosses(parent_pop:Population, n_crosses,genetic_map)-> Population:
new_pop = []
total_cross=n_crosses
random_pairs = torch.randint(0,len(parent_pop), (total_cross*2,)).reshape(-1,2)
cross_pairs = parent_pop[random_pairs]
for pair_idx in range(len((random_pairs))):
egg = cross_pairs[pair_idx][0]
pollen = cross_pairs[pair_idx][1]
x = make_cross(genetic_map, egg , pollen)
new_pop.append(x)
new_pop = torch.stack(new_pop)
return new_pop

# %% ../nbs/04_crossing.ipynb 8
def truncate_select(tgv):
return torch.topk(tgv,20).indices
22 changes: 13 additions & 9 deletions chewc/trait.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -176,18 +176,19 @@ def calculate_intercept(self):
self.intercept = self.target_mean - current_mean


def phenotype(self, pop, h2=None, varE=None):
def phenotype(self, pop, h2=None, varE=None, reps=1):
"""
Generate phenotypes for the individuals in the population.
Args:
pop (Pop): The population object.
h2 (float, optional): Target narrow-sense heritability. Defaults to None.
varE (float, optional): Target environmental variance. Defaults to None.
reps (int, optional): Number of repetitions for averaging. Defaults to 1.
Returns:
torch.Tensor: A tensor of phenotype values for each individual.
Shape: (n_individuals,)
torch.Tensor: A tensor of average phenotype values for each individual.
Shape: (n_individuals,)
Raises:
ValueError: If both h2 and varE are provided or neither are provided.
Expand All @@ -204,13 +205,16 @@ def phenotype(self, pop, h2=None, varE=None):
# Calculate varE based on target heritability
varG = genetic_values.var()
varE = (varG / h2) - varG

# Generate environmental error
error = torch.randn(len(pop)) * torch.sqrt(varE)

# Calculate phenotypes
phenotypes = genetic_values + error

# Generate phenotypes for each repetition
phenotypes = torch.zeros(len(pop))
for _ in range(reps):
error = torch.randn(len(pop)) * torch.sqrt(varE)
phenotypes += genetic_values + error

# Average phenotypes across repetitions
phenotypes /= reps

return phenotypes


Expand Down
Binary file modified index_files/figure-commonmark/cell-4-output-1.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added index_files/figure-commonmark/cell-5-output-1.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added index_files/figure-commonmark/cell-6-output-1.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
189 changes: 179 additions & 10 deletions nbs/02_trait.ipynb
View file

Large diffs are not rendered by default.

77 changes: 0 additions & 77 deletions nbs/03_meiosis.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,83 +109,6 @@
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1ca41486",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"tensor([[[0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 1,\n",
" 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1,\n",
" 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1,\n",
" 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1,\n",
" 0, 0, 1, 1, 1, 0, 1, 1],\n",
" [1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1,\n",
" 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1,\n",
" 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1,\n",
" 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0,\n",
" 0, 0, 1, 0, 0, 1, 0, 0],\n",
" [0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0,\n",
" 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0,\n",
" 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0,\n",
" 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1,\n",
" 1, 0, 0, 0, 1, 0, 1, 1],\n",
" [1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1,\n",
" 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0,\n",
" 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0,\n",
" 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0,\n",
" 1, 1, 0, 1, 0, 0, 0, 0],\n",
" [1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0,\n",
" 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0,\n",
" 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1,\n",
" 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,\n",
" 1, 1, 0, 0, 1, 1, 1, 1],\n",
" [0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1,\n",
" 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0,\n",
" 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1,\n",
" 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0,\n",
" 1, 0, 1, 1, 0, 0, 1, 0],\n",
" [1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1,\n",
" 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1,\n",
" 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0,\n",
" 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0,\n",
" 0, 0, 0, 0, 0, 1, 0, 1],\n",
" [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0,\n",
" 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1,\n",
" 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0,\n",
" 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0,\n",
" 1, 1, 1, 0, 0, 0, 1, 1],\n",
" [1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1,\n",
" 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1,\n",
" 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0,\n",
" 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0,\n",
" 1, 1, 0, 1, 0, 0, 1, 1],\n",
" [1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0,\n",
" 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0,\n",
" 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1,\n",
" 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0,\n",
" 1, 1, 1, 0, 1, 0, 1, 1]]])"
]
},
"execution_count": null,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ploidy = 2\n",
"number_chromosomes = 10\n",
"loci_per_chromosome = 100\n",
"n_founders = 50\n",
"genetic_map = create_random_genetic_map(number_chromosomes,loci_per_chromosome)\n",
"crop_genome = Genome(ploidy, number_chromosomes, loci_per_chromosome, genetic_map)\n",
"founder_pop = create_random_founder_pop(crop_genome , n_founders)\n",
"simulate_gametes(genetic_map, founder_pop[0])"
]
},
{
"cell_type": "code",
"execution_count": null,
Expand Down
82 changes: 3 additions & 79 deletions nbs/04_crossing.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,12 +54,10 @@
"from chewc.meiosis import *\n",
"from chewc.trait import *\n",
"import torch\n",
"def double_haploid_operation(gamete):\n",
" return gamete.repeat(2, 1, 1)\n",
"\n",
"import random\n",
"def make_DH(genetic_map, individual):\n",
" haploid_line = simulate_gametes(genetic_map, individual)\n",
" return double_haploid_operation(haploid_line)"
" return gamete.repeat(2, 1, 1)"
]
},
{
Expand All @@ -71,84 +69,10 @@
"source": [
"#| export\n",
"\n",
"def cross_operation(egg, pollen):\n",
" return torch.stack((egg, pollen), dim=0).squeeze(1)\n",
"\n",
"def make_cross(genetic_map, mother, father):\n",
" egg = simulate_gametes(genetic_map, mother) # Simulating gametes for the mother\n",
" pollen = simulate_gametes(genetic_map, father) # Simulating gametes for the father\n",
" return cross_operation(egg, pollen)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d9670c3d",
"metadata": {},
"outputs": [],
"source": [
"#| export\n",
"def random_crosses(population: Population, \n",
" n_crosses: int, \n",
" ) -> Population:\n",
" \"\"\"\n",
" Performs random crosses within a population.\n",
"\n",
" Args:\n",
" ----\n",
" population (Population): The population to perform crosses within.\n",
" n_crosses (int): The number of crosses to perform.\n",
"\n",
" Returns:\n",
" -------\n",
" Population: A new population of offspring resulting from the crosses.\n",
" \"\"\"\n",
" # make sure there are enough individuals for the requested number of crosses.\n",
" if n_crosses > len(population.individuals) // 2:\n",
" raise ValueError(f\"Not enough individuals in the population ({len(population.individuals)}) for {n_crosses} crosses. Need at least {2*n_crosses}.\")\n",
"\n",
" # sample parents randomly\n",
" selected_parents = random.sample(population.individuals, 2 * n_crosses)\n",
"\n",
" # create offspring\n",
" offspring = []\n",
" for i in range(n_crosses):\n",
" mother = selected_parents[i * 2]\n",
" father = selected_parents[i * 2 + 1]\n",
" haplotypes = make_cross(population.individuals[0].genome.genetic_map, mother, father)\n",
"# genetic_value = trait.calculate_genetic_value(haplotypes.unsqueeze(0)).squeeze(0) # calculate GV\n",
" offspring.append(\n",
" Individual(\n",
" genome = mother.genome,\n",
" haplotypes=haplotypes,\n",
" mother_id=mother.id,\n",
" father_id=father.id,\n",
"# genetic_value=genetic_value,\n",
" )\n",
" )\n",
" return Population(individuals=offspring)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "72437879",
"metadata": {},
"outputs": [],
"source": [
"#| export\n",
"def random_crosses(parent_pop:Population, n_crosses,genetic_map)-> Population:\n",
" new_pop = []\n",
" total_cross=n_crosses\n",
" random_pairs = torch.randint(0,len(parent_pop), (total_cross*2,)).reshape(-1,2)\n",
" cross_pairs = parent_pop[random_pairs]\n",
" for pair_idx in range(len((random_pairs))):\n",
" egg = cross_pairs[pair_idx][0] \n",
" pollen = cross_pairs[pair_idx][1]\n",
" x = make_cross(genetic_map, egg , pollen)\n",
" new_pop.append(x)\n",
" new_pop = torch.stack(new_pop)\n",
" return new_pop"
" return torch.stack((egg, pollen), dim=0).squeeze(1)"
]
},
{
Expand Down
Loading

0 comments on commit a8980cd

Please sign in to comment.