# Simulate genotype data of family

```###############
### 1
###############

# When we have one marker
# No. alleles
Na<-4
# Allele frequency vector
Pa<-c(0.3,0.2,0.1,0.4)
# Sum of allele frequency should be 1.
Pa<-Pa/sum(Pa)

###############
### 2
###############

# Generate allele frequency values with uniform distribution
Pa<-runif(Na)
# When you use uniform distribution, allele frequency tends to be 1/Na,
# you can use dirichlet distribution to avoid this phenomenon.
install.packages("MCMCpack")
library(MCMCpack)
# Generate 1 set of allele frequency accoring to dirichlet distribution.
Pa<-rdirichlet(1,rep(1,Na))
Pa
# Double-check sum of allele freq being 1.
sum(Pa)

###############
### 3
###############

# Now we have multiple markers on one chromosome.
Ns<-4
# NaS is a vector of No. alleles of markers
NaS<-c(3,4,5,2)
# Allele frequency
PaS<-list()
for(i in 1:Ns){
PaS[[i]]<-c(rdirichlet(1,rep(1,NaS[i])))
}
NaS
PaS

hap1<-hap2<-rep(0,Ns)
for(i in 1:Ns){
hap1[i]<-sample(1:NaS[i],1,prob=PaS[[i]])
hap2[i]<-sample(1:NaS[i],1,prob=PaS[[i]])
}
hap1
hap2

###############
### 4
###############

# Function to generate haplotypes according to No. allele vector and allele freq. vector.
MakeHaplotype<-function(NaS,PaS,N){
# haplotypeはN x length(NaS)の行列に納める
H<-matrix(0,N,length(NaS))
for(i in 1:length(NaS)){
H[,i]<-sample(1:NaS[i],N,prob=PaS[[i]])
}
H
}
# Use this function
N<-2 # Generate 2 haplotypes for 1 indivudual.
Haps<-MakeHaplotype(NaS,PaS,N)
Haps

###############
### 5
###############

# When you simulate recombination, you use two haplotypes.
# Define genetic distance between neighboring markers.
Dg<-runif(Ns-1)
# Let them cross-over at random.
# Decide which inbetweens cross-over.
rands.crossing<-runif(length(Dg))
# Take out inbetween with cross-over.
Cross<-which((Dg-rands.crossing)>0)
Cross
# Make recombinants.
hap1<-Haps[1,]
hap2<-Haps[2,]
rec.hap1<-hap1
rec.hap2<-hap2
if(length(Cross)>0){
for(i in 1:length(Cross)){
# Exchange segments downstream from cross-over points
tmp.tail1<-rec.hap1[(Cross[i]+1):Ns]
tmp.tail2<-rec.hap2[(Cross[i]+1):Ns]
rec.hap1[(Cross[i]+1):Ns]<-tmp.tail2
rec.hap2[(Cross[i]+1):Ns]<-tmp.tail1
}
}
Cross
hap1
hap2
rec.hap1
rec.hap2

###############
### 6
###############

# Make function to generate recombinants.
Rec.Haplotype<-function(hap1,hap2,d){
Ns<-length(hap1)
rands.crossing<-runif(length(d))
# Take out inbetweens with crossover.
Cross<-which((d-rands.crossing)>0)
rec.hap1<-hap1
rec.hap2<-hap2
if(length(Cross)>0){
for(i in 1:length(Cross)){
tmp.tail1<-rec.hap1[(Cross[i]+1):Ns]
tmp.tail2<-rec.hap2[(Cross[i]+1):Ns]
rec.hap1[(Cross[i]+1):Ns]<-tmp.tail2
rec.hap2[(Cross[i]+1):Ns]<-tmp.tail1
}
}
matrix(c(rec.hap1,rec.hap2),byrow=TRUE,nrow=2)
}
# Use the function

rec.haps<-Rec.Haplotype(Haps[1,],Haps[2,],Dg)
rec.haps

###############
### 7
###############

# Generate parents' haplotype pairs from population.
father<-MakeHaplotype(NaS,PaS,2)
mother<-MakeHaplotype(NaS,PaS,2)
father
mother
# Make gametes' haplotype
father.gamates<-Rec.Haplotype(father[1,],father[2,],Dg)
mother.gamates<-Rec.Haplotype(mother[1,],mother[2,],Dg)
# Select gamates and decide offspring's haplotype pair.
from.father<-father.gamates[sample(1:2,1),]
from.mother<-mother.gamates[sample(1:2,1),]
offspring<-matrix(c(from.father,from.mother),byrow=TRUE,nrow=2)
offspring

father
mother
offspring

###############
### 8
###############
# Make the job into a function
MakeOffspring<-function(father,mother,d){
father.gamates<-Rec.Haplotype(father[1,],father[2,],d)
mother.gamates<-Rec.Haplotype(mother[1,],mother[2,],d)
from.father<-father.gamates[sample(1:2,1),]
from.mother<-mother.gamates[sample(1:2,1),]
offspring<-matrix(c(from.father,from.mother),byrow=TRUE,nrow=2)
offspring
}
offspring<-MakeOffspring(father,mother,Dg)
father
mother
offspring

###############
### 9
###############
# We make matrix to handle familial relationship.
# Trio with parents and an offspring
M.Oyako.Trio<-matrix(c(0,0,0,0,1,2),byrow=TRUE,ncol=2)
M.Oyako.Trio
# Somebody with its parents and 4 grandparents
M.3gen<-matrix(c(0,0,0,0,0,0,0,0,1,2,3,4,5,6),byrow=TRUE,ncol=2)
M.3gen
# A pair of cousins.
M.cousins<-matrix(c(0,0,0,0,1,2,0,0,1,2,0,0,3,4,5,6),byrow=TRUE,ncol=2)
M.cousins

###############
### 10
###############

# Generate haplotype pairs for all the individuals in the pedigree.
sim.G<-list()

# From the 1-st to the last member in the pedigree.
for(i in 1:n.rows){

tmp.father<-M.cousins[i,1]
tmp.mother<-M.cousins[i,2]
if(tmp.father==0){
father.genotype<-MakeHaplotype(NaS,PaS,2)
}else{
father.genotype<-sim.G[[tmp.father]]
}
if(tmp.mother==0){
mother.genotype<-MakeHaplotype(NaS,PaS,2)
}else{
mother.genotype<-sim.G[[tmp.mother]]
}
sim.G[[i]]<-MakeOffspring(father.genotype,mother.genotype,Dg)
}
sim.G
M.cousins

###############
### 11
###############

# Make the job into a function/
Simulate.Genotype<-function(M,NaS,PaS,Dg){
sim.G<-list()

for(i in 1:n.rows){
tmp.father<-M.cousins[i,1]
tmp.mother<-M.cousins[i,2]
if(tmp.father==0){
father.genotype<-MakeHaplotype(NaS,PaS,2)
}else{
father.genotype<-sim.G[[tmp.father]]
}
if(tmp.mother==0){
mother.genotype<-MakeHaplotype(NaS,PaS,2)
}else{
mother.genotype<-sim.G[[tmp.mother]]
}
sim.G[[i]]<-MakeOffspring(father.genotype,mother.genotype,Dg)
}
sim.G
}
sim.G<-Simulate.Genotype(M.cousins,NaS,PaS,Dg)
sim.G

###############
### 12
###############

# We now handle multiple chromosomes
num.Chromosomes<-22
NaS.list<-list()
PaS.list<-list()
Dg.list<-list()
for(i in 1:num.Chromosomes){
NaS.list[[i]]<-NaS
PaS.list[[i]]<-PaS
Dg.list[[i]]<-Dg
}
sim.Gs<-list()
for(i in 1:num.Chromosomes){
sim.Gs[[i]]<-Simulate.Genotype(M.cousins,NaS.list[[i]],PaS.list[[i]],Dg.list[[i]])
}

###############
### 13
###############

# Make a function for the job.
Simulate.Genotype.Multi<-function(M,NaS.list,PaS.list,Dg.list){
num.Chromosomes<-length(NaS.list)
sim.Gs<-list()
for(i in 1:num.Chromosomes){
sim.Gs[[i]]<-Simulate.Genotype(M.cousins,NaS.list[[i]],PaS.list[[i]],Dg.list[[i]])
}
sim.Gs
}
sim.Gs<-Simulate.Genotype.Multi(M.cousins,NaS.list,PaS.list,Dg.list)

###############
### 14
###############

# We make multiple families.
Niter<-10
sim.Gs.list<-list()
for(i in 1:Niter){
sim.Gs.list[[i]]<-Simulate.Genotype.Multi(M.cousins,NaS.list,PaS.list,Dg.list)
}

```