-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathfamily_simulation.R
125 lines (91 loc) · 3.03 KB
/
family_simulation.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
# Generate a child from two parents (father=number of column in X (VCF without header) for the father, mother=number of column in VCF for the mother)
offspring <- function(X, father=NULL, mother=NULL) {
fatherGT <- X[,father]
motherGT <- X[,mother]
fatherAll <- str_split_fixed(fatherGT, "/", 2)
motherAll <- str_split_fixed(motherGT, "/", 2)
snps <- nrow(X)
f <- vector(length = snps)
m <- vector(length = snps)
f[1] <- sample(1:2, prob = c(0.5,0.5))[1]
m[1] <- sample(1:2, prob = c(0.5,0.5))[1]
for(n in 2:snps) {
d <- abs((X[n,2]-X[n-1,2])/1000)
if(d==0) {
prec <- 0
}else{
prec <- 0.5*(1-exp(-2*d))
}
if(f[n-1]==1) {
f[n] <- sample(1:2, prob = c(1-prec,prec))[1]
}else{
f[n] <- sample(1:2, prob = c(prec,1-prec))[1]
}
if(m[n-1]==1) {
m[n] <- sample(1:2, prob = c(1-prec,prec))[1]
} else {
m[n] <- sample(1:2, prob = c(prec,1-prec))[1]
}
}
child <- sapply(1:snps, function(i) {
paste(as.matrix(fatherAll[i, f[i]]), as.matrix(motherAll[i, m[i]]),
sep = "/")
})
X <- cbind(X, child)
colnames(X) <- 1:ncol(X)
return(X)
}
# Import the data from the founders
library(data.table)
library(stringr)
founders <- data.frame(fread("founders_sort.txt"))
colnames(founders) <- 1:ncol(founders)
# Run only one pedigree type
#type 1 pedigree (16 ind)
family <- offspring(founders, 10, 11)
family <- offspring(family, 12, 13)
family <- offspring(family, 14, 15)
family <- offspring(family, 14, 15)
family <- offspring(family, 14, 15)
family <- offspring(family, 14, 15)
family <- offspring(family, 14, 15)
family <- offspring(family, 14, 15)
family <- offspring(family, 14, 15)
family <- offspring(family, 14, 15)
family <- offspring(family, 14, 15)
family <- offspring(family, 14, 15)
#type 2 pedigree (16 ind)
family <- offspring(founders, 10, 11)
family <- offspring(family, 10, 11)
family <- offspring(family, 10, 11)
family <- offspring(family, 10, 11)
family <- offspring(family, 12, 14)
family <- offspring(family, 12, 14)
family <- offspring(family, 12, 14)
family <- offspring(family, 12, 14)
family <- offspring(family, 13, 14)
family <- offspring(family, 13, 14)
family <- offspring(family, 13, 14)
family <- offspring(family, 13, 14)
#type 3 pedigree (16 ind)
family <- offspring(founders, 10, 11)
family <- offspring(family, 10, 11)
family <- offspring(family, 14, 12)
family <- offspring(family, 14, 12)
family <- offspring(family, 14, 12)
family <- offspring(family, 14, 12)
family <- offspring(family, 14, 12)
family <- offspring(family, 13, 15)
family <- offspring(family, 13, 15)
family <- offspring(family, 13, 15)
family <- offspring(family, 13, 15)
family <- offspring(family, 13, 15)
#type 4 pedigree (12 ind)
family <- offspring(founders, 10, 11)
family <- offspring(family, 11, 12)
family <- offspring(family, 11, 12)
family <- offspring(family, 13, 17)
family <- offspring(family, 14, 18)
family <- offspring(family, 15, 19)
# Save simulated family
write.table(family, file = "simulated_family.txt", quote = F, sep = "\t", row.names = F, col.names = F)