-
Notifications
You must be signed in to change notification settings - Fork 0
/
pedtest.jl
121 lines (95 loc) · 3.6 KB
/
pedtest.jl
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
112
113
114
115
116
117
118
119
120
121
#------------------------------------------------------------------------------#
# makeA4.jl
#------------------------------------------------------------------------------#
# Austin: Austin Putz
# Created: 2017-09-27
# Updated: 2023-02-18
# License: MIT
# Version: 0.0.3
#------------------------------------------------------------------------------#
# Description
#------------------------------------------------------------------------------#
# Create A with the tabular method
# Make Sure:
# - Animals are ordered 1,2,...,n
# - Missing parents are coded as 0
# Example:
#ped = DataFrame(
# animal = [1, 2, 3, 4, 5, 6],
# sire = [0, 0, 1, 1, 4, 5],
# dam = [0, 0, 2, 0, 3, 2]
#)
#First press ']' and enter Pkg menu. Then write 'generate project_name', 'instantiate' respectively.
#use git clone or copy the files from Pedigree.jl package and paste and replace in the newly created project folder.
using Pkg
using Pedigree
using DataFrames
Pkg.add("CSV")
using CSV
working_dir = "C://Users/USER/gg/"
data_file = "swine_ped.csv"
ped_MSU = CSV.read(working_dir * data_file, # this will just combine the 2 strings
DataFrame,
header=true,
delim=',',
missingstring="NA")
#------------------------------------------------------------------------------#
# Function
#------------------------------------------------------------------------------#
# Begin Function
function makeA(ped::DataFrame)
# extract 1st 3 columns (if more from renumered pedigree)
ped = ped[:, 1:3]
# check eltypes (should be Int64, not a string)
if eltype(ped[:,1]) == Int64
@info "Gave the correct type (Int64)"
else
@info "Gave wrong column type for this algorithm"
error("Need to provide Int64 numbered 1 to n")
end
# pull out animal, sire, dam
a = ped[:,1]
s = ped[:,2]
d = ped[:,3]
# add 1 to the sire and dam if ordered 1,2,3,...,n
sp1 = s .+ 1
dp1 = d .+ 1
# I do this because I will add an extra row and column so that
# when parents are referenced it refers to row 1 or column 1 and
# I won't need to write a bunch of if/then statements
# set number of animals
n = size(ped, 1)
# add 1 for use in loop (will add a padding of 0's in row 1 and column 1)
N = n + 1
# check to make sure the pedigree is ordered 1 to n
if a[1] != 1
# throw an error if the first animal isn't listed as animal #1
error("Pedigree must be sequential 1 to n")
elseif a[n] != n
error("Pedigree must be sequential 1 to n")
end
# allocate the whole A Matrix
# This matrix will be padded with a row and column of 0's
A = zeros(N, N)
# Begin FOR loop
for i in 2:N
# calculate the diagonal element of A
# diag = 1 + half the relationship of the sire and dam
A[i, i] = 1 + (A[dp1[i-1], sp1[i-1]] / 2)
# loop to calculate off-diagonal elements of A
# Loop down a column for memory efficiency
for j in (i+1):N
# bottom left triangle of A
# average the relationship between animal and parents of other individual
A[j, i] = 0.5 * A[sp1[j-1], i] +
0.5 * A[dp1[j-1], i]
# Symmetric
A[i, j] = A[j, i]
end
end
# return the A matrix (remove first row and column)
return A[2:end, 2:end]
end # end function
ped_MSU_sort = sort_ped(ped_MSU)
ped_MSU_renum = renum_ped(ped_MSU_sort)
A = makeA(ped_MSU_renum)