L!M>($Rmu_e9EXANWpck;LX>6vv{E}~FnW!ei}#47%e>o9Q?nRi0_0cRawrDo
z&VoLs&
z_u%|dUofh9m+2x#X6|#2r1l!^kaN$-iO5sf=4m+O#LM+O*GWlE*KnlzRN;@)5SM7k
zXMFY*lTb}dx8BX5pD=sysenk7F5|meQEHdh7c-nLi-;oh||Hba!bI)d;HclQdFD7RDP5dlZ-~R_GOvkzFY$_UdQF
z3-MVmLi|0;=0`U~ENUiSxDB>K;!_64g4EFb+r}>f=xE_x0q3ti@C`bI0`&0{zc4Tx
zD=#4DiG`7p8t{A=yooDkmb%ZYr1(aIE-L_7Un#2w7mL}F%
z$?cGtpX80#SFQ4Ef2(#?44q7q?t@>2d?gaA?dV`l!!suWKl5W`C@uZiMTW*PLORhQ
zPcdj>Y16B9yRQ>I=&DkKLHuV8i|2Mcy*_&s)xiNt`!2ze(uqlL@OjUgwe#)Du(!0(
zH1^n;Ti|Q+J6LVK{(%8(rWc*KwJ<`UcQ&(8y$c8fJ^gDv)7y{^_K^}Wz@g--K~?vT
zXU*4!q37HSyU8VYfJJ2U&eB(#PGLXOlsiwA=%L^ABmOfw@%P&&L<#jBJAzRr|{
zYU8Ma)#fh6*%o*d4wlz=)gM!r26^RhoZ3EqAjovBKMH1`+S^rTy#lQ?Ua?qx^Tn3w
z>xUz92zqknswQZfgEz!27l^_dzu-sOVCwpkIlnfKf=WA0X9_xq-Oiv;|Fh2>ObjbZWMul%YzRY-Iatr@E@N8g>-(4d
zA0XzIa*l@k_%crKFd}Z@{w=f0^xwtQgV%#b93UXD!84Jut&yUGt(_xyS?CSW(8j>R
zn$ZPlI-3cy#W!>RY(^c%tC3`(nS@oQRzlFg@UwpUl#)xlMy5axB^;u;)mDqVnURyB
zvR|%3IyQa#H239MRD?ph5eAI#E54{q;#bweF)>Ly-0-U}h?&C2)LYGQ?6Md5GDv8Y
zjy_-b6BkDGdWw1N*|w#UTnLhzec?_E&n;Rx6C3W>6e9Hw8+px-FSq3Y@Uyae8o0Bm
zutL*_gnzJe{y9>rz09in4m&C5=LS+q^wU6)_x|AB@uq0Bn2#sNM4b2NRWtq_y7>4H
zz%66@-DPt24U&ui?O^qKb^hl-d@h@*Lqz)yMjr(uTd@wQX#yMg4#l-Uh|$
z*Y&d1W%2m~rWwN=>GMhvpli*?A=Z@Xwr>AC2tI}}hm5ryH2JIb|oNM;^YP`ZaYPf{Pcp0SdQK*8z$!laZ
z#M-zgdl=bn>!}2WQZVLhlr0LLyNw$TA#9{C(>cI!M}tDyq8V6_GZ|n$fa_8q!Ey$|
z(E2^zayZjE>V~*){ovcPP^Q*nWn}wogduEM@TJ2(rVUgaG2i+H@TDDRzTl4l81ry2
zIKRuT4zj%xj5$1*&Gq*Bf?ydE?Oh@m&UL9XkibwAeB~Z{#;Rfs6(tM;_!SIL(SGZW
zY!6m1_#XQ+zU?Y1N+<+y8hpvN?ZOn<9xN4nFLesC-9QD$9|X|Lr-9@$na$o9H0<_qj#^Ng3MJBBP^Le2z%%~PStuhHGSF+M2hn
zFSC@=mVSo-KpLo5&BGF-4}lwRK69I{ryWT`r|!HGV*;r~s$o`9KYwi-rB@$VnGX}T
zFF4!>W=<5N4^2{Z+vQ20Lrgi^tU(5|iq9xRiM#rMIm?~{)w(r@-?@YCnPck{UZTxY
zxlh+y*1Dl>hzI=ae@uo&!x8%zvN#&fzj0TGn1HbzkE?;I?-?i*!DfHBrBpY?zmYY&
z|CL+38|x4E8?(PaRJ&FGkpE!Ib#we1S-1PmzmY||k^YUW(k=IgY@+>kAw7lOVkYT)
z+!z#esHE})a|wTU#=Qc>8KH(+MGQ9j71$Rf%J&z%a+87h5_Dl2m`KH3l1v~Gb6P*r
z76q^gy38riWO6CSDy>1NMMYo--4e7~gP3PZLuzt`PwMl)a-Get
zI_ST&;}o;}qh22=V>sM+`eQ|~(A0ms{67>a4k`IJMgG|T?~2qjE;pGgw(@{O>+fa)
zq1ubG7tY_?;g?D7V-_h781|9okJ{zp=dGatuX-TOU_ID8ZENLa$v|xhx)2RiBun#p
zVCb0I@17+Mhmm572Ru&$lGOfXNR!&Q5Q<@ROO3x6@53I5_bcp!OhHMnfPYyL52)Ae
zHvH}}QMeoSpTq}yKjJK-H`HV=%$csIP}cq(3EC3c==M_oBx2kD3n!Yr_&+(Z>{b6s
z#ITq9ClSY9_aDSx27aUx2-9BWFV05nr87$Jlj-_Aevm4Z_<2-wyp
zcuaiA9*>C+x#JgcO$(ekz$0#zzT5{o@{6wRw3+~`+gFq1
zePM9Pd1mWKgRsW3?xxGbu4j#ICA`VZbJ%=?J>hwLh66n^x3}U^!h&vx?~}5}2F*OU
z1lo7iZg$jw0U+6Pr|oV1?Fe1H=#9knI5tX!Y{5*K2yvQ_ko~Wkam~@0gy@pdekts|
z{?sbDB%TgQW<&=+3bF?>R}!ejLFI}$3T<;J8|T|KgX`$wxw4>(q-};xBKzZzKq)H$
z25Bj=w+d17U$wV}SDsRW)86-~)J(-34@ZLjqWBQD-0w4u`58=zb7b{eWwSfnbl);U
z5VGfy&)Ni{h>(~qq8xL-hX?I-56T+mzv*s?DJ&7rifZEk?Ib}r3c>xUUxZPCfT83C5}2+O!S;0+4bUZ
zO=p=zDG3ZJ%4pAyB&W4yqJ#`qv5)bf3q|OA=8A!m8kn_<@n9udXn>+T#6gVxExesrYA3E
zE^@|e)wd`;%{PM8{{NW9(*GMwpKDTnWB<3r{<{k3w%c_2c=7ruWVRK|&vGq^@(I;u
zxjQbrMcKI?a{t-yg&|mbB)V0^btZQ{G=Ng}{~Xg6Sn&#ljf&0wC}{%8UfG9UpyP;%
zJv4F$WB5%Ue2LZCOoQd*3@!@+uKh&En
z*b9sDmS&zzWKUjgit}zRaZ1iwqq@l9-eic-U;X~B%i$}P;zILY+onPOW`OncX8ND0
z-L|P4E?zb@nX1O`?f1
z9i5DABX*0VvQaaSSq0rZ5}uV63ibkjl*eT$P6OAXigKfP(vpz3G(iSOGMu@JdPXm3
zVe^H4ED2p_4M`+3%Vy5PQQ424q_(sL#!i1IC7-jCAt+lqD&3B=d
zs$s-u^dr-r+U4C3?Brbm#tZtir833|I-xx6WRI)&gN7N`A5*rM*w-i<;#!6JGOzCT
zbz7!kaRejQZMsx@tnM0hyWZ`bb}Wpp^!?NX3S91FXw93=YF^cp
zo|?JK%I)zqC-V`quwHmH&GuH$^yb=bjA4^WL`@22jSCXZC93dYu{gAxZfNcd=J7dN
zK1CzzwepD_za?!VE8vX6O&z0ra@BOYx7sS$_Li!ef`4gS3M`!U#ZrOWBQAEOj;dxR
zejv+u)pe(+^faM{3Xez6gSSLr=NH>o%fVl4?cXeY*1$CHSOGUq?M_WI1y%=7=Z)$a
z$6vy9nmcsAS{J3KEe|kPmTFhUzW_-}S{EMx^}bs&@q+IQuvVBiyI59SvOZSlkz5=Q
zGpiwU8kKBxY@&2C7R0$Rkq6Qkt_}#KPPC!lh|BLd9e*pIIhuADl=R6`(T>1vbX54z
z=f%+^RNzeO{ziQ#(7F*PSE#RuK`dgnk!o~JZ>*t<+y``1M#k^+F=c#yHg79`yWZhz
zC)zaw#7ra;dT8g8ohb%eOvNbY}h6G?s)ixts0T%G09KU
z3sVNo?tyCsHBT>9<5krBD5clJOu1PlExu30GD6|%K?`leB)yNdm4FUD?fEHdRO!;R
zttw-?^T!|yhP{4@x)#yinZ5a$uMVT_!8^mojU0p)VAWJ?o
zZHhX+22}=-8VFp0R5aIEX6fF;n|U-3{?VGM{bgYCDdKRZJ%_fHj};04-ojM7E)`U^
z8-d5urV7#of{IL5U4dLBJG^xfvX*D=H#8SpezCnF@2w6BE3e)Z9@yfp@}U&HxD)Kj
zb8ojg4`B=Cyqr2Qnaf)uZ%Cq+fr}rx_B!?7IX@+sKpr%k*7f+>By;&H&RSvO+Bp%N
zSr0bZJSrfu@nU53&U(qU>Qqg>#(TTAlZ4%ZPbvI!&k3udn79^7Sh@H%+Z}MwEgbI{qX$}`NgOss}zu?)zPg;Cr?;Xd6p>o2O+f-!x`6TMGVlB!f7`?2_xw{oWn&%gR!hSSKaiSkdQi?Ic7
zUuW-toSB{W*EwTlgy<6@&%gAuq8mv{1tjbdwg)S(e|h&)+Do+*0PfpAd*7%1tg#(@
zB?#Pm@Fvpz{=1QZwUMZ~frZWe7Q9Xj2;9|YL3^+ZIjeengQQz6s6U@>LMJfV1Q1DR
zi!IF4f?{T24%f)YecZe?yRvW%H!^+4t098qubLbqeXrygqdplbm#cl|!`%Ojuwf}k
zSI2reQan)A(jg8Bl)ec)D%7r}W%P<>R3H-6NOZt|PITp49K~Q*9aIr+I0)2(P%Y9#
zC7%4$%$R`-uy>+gJ<@C6
zI0n^=QU&-ycZ4l#zuKKJOiv(h+F3}i*DjOylpSaU3X2)*H4U#$%V#cA7o@gpI0nr=eHd9t*5{pjD4TeBBf*i6lLLT@l8wp
z#Pc45OM3EqQ{Mb?(1A-mECCv84~GV2ipfXV%l?;QlJ6Ixt1pIB;ZX@hVL4t|+2OuD
z2U!VEt`GXj(5mhH|iVxm2V*&SV0@f={6E
zS1I??%&%!tEQj0U+{X@NjD)p#>0YAhFZrQ<&!ZGn=7aWj
zV`Rj;A^s4ZY4?cY8!7sLr^8CTlsu`s_~z?asM2A22xK@L&5Qg1}=a>{d+tb*xCJWID+Hw_a`Gx*cQb6EN~C{
z^2M_`e4aVQXMUxIuhfq#4Ity{&E!59`}Quf
z&q~J|5?X44&XX<+YIqE-3(vhP7ZWqW!+zB31FG}#7!xxZEUhAgIpDWGyTVW`!2Klp
zURfJyk5~88Vmu%}<*Ad8D^NZ69Eb{+A$?5AhY?Z*2O5j0Q+TB*)N>ud0Fe`Cof2fk
zb2eagCe9=YjkEEYL?6mviNJp23-qleC*>$#cOTEe7v=+bZjhfXuj2rhUaG<+l?P$7
zac69x%1KwCndhr!S1Sa)NjJgQpb6u`3t0IrL8LJOTh=QQf5J(g)bx&1GqBd-_O1O<*vWJ|K$|CqRf8+{O8fzzXb2U+rVSu
zZwGN73qC%4_(${vyr1{`z{+1o5g!Zx=NYL#q7V?_h!4X5A19_B<2*hm^9O10*?%eV
zk5e;`Q65*I|3PU+dq8`UBX+
z{Q&qMPfCxa|MOwvk1Pa40ngvM^w*=uWAXo-8~!fN#s4?)-!p}x3@kY8As~>!f2rV+
Kb1HcM>;C{s!H#VJ
literal 0
HcmV?d00001
diff --git a/Julia/beefblup.jl b/Julia/beefblup.jl
index 547cfb4..16c0c6e 100644
--- a/Julia/beefblup.jl
+++ b/Julia/beefblup.jl
@@ -7,6 +7,9 @@
# Import the required packages
using XLSX
+using LinearAlgebra
+using Dates
+using Gtk
# Display stuff
println("beefblup v 0.0.0.1")
@@ -16,19 +19,23 @@ print("\n")
### Prompt User
# Ask for an input spreadsheet
-# print("Enter the full filename of a beefblup worksheet> ")
-# path = readline(stdin)
-path = "C:\\Users\\cclea\\source\\repos\\beefblup\\Excel\\Master BLUP Worksheet.xlsx"
+path = open_dialog_native(
+ "Select a beefblup worksheet",
+ GtkNullContainer(),
+ ("*.xlsx", GtkFileFilter("*.xlsx", name="beefblup worksheet"))
+)
# Ask for an output text filename
-# print("Enter the full filename of your desired results> ")
-# savepath = readline(stdin)
-savepath = "C:\\Users\\cclea\\source\\repos\\beefblup\\results.txt"
+savepath = save_dialog_native(
+ "Save your beefblup results",
+ GtkNullContainer(),
+ (GtkFileFilter("*.txt", name="Results file"),
+ "*.txt")
+)
# Ask for heritability
-# print("What is the heritability for this trait?> ")
-# h2 = parse(Float64, readline(stdin))
-h2 = 0.4
+print("What is the heritability for this trait?> ")
+h2 = parse(Float64, readline(stdin))
### Import input filename
print("[🐮]: Importing Excel file...")
@@ -39,7 +46,7 @@ data = XLSX.readxlsx(path)[1][:]
print("Done!\n")
### Process input file
-print("[🐮]: Processing and formatting data ...")
+print("[🐮]: Processing and formatting data...")
# Extract the headers into a separate array
headers = data[1,:]
@@ -49,9 +56,9 @@ data = data[2:end,:]
data = sortslices(data, dims=1, lt=(x,y)->isless(x[2],y[2]))
# Define fields to hold id values for animals and their parents
-ids = data[:,1]
-damids = data[:,3]
-sireids = data[:,4]
+ids = string.(data[:,1])
+damids = string.(data[:,3])
+sireids = string.(data[:,4])
numanimals = length(ids)
# Find the index values for animals and their parents
@@ -98,23 +105,183 @@ normal = Array{String}(undef,1,length(headers)-5)
for i in 6:length(headers)
for j in numanimals:-1:1
if !ismissing(data[j,i])
- normal[i-5] = data[j,i]
+ normal[i-5] = string(data[j,i])
break
end
end
end
-# Print the results of the "normal" definition
-println(" ")
-println("For the purposes of this analysis, a 'normal' animal will be defined by")
-println("the following traits:")
-for i in 6:length(headers)
- print(headers[i])
- print(": ")
- print(normal[i-5])
- print("\n")
+print("Done!\n")
+
+### Create the fixed-effect matrix
+print("[🐮]: Creating the fixed-effect matrix...")
+
+# Form the fixed-effect matrix
+X = zeros(Int8, numanimals, floor(Int,sum(numgroups))-length(numgroups)+1)
+X[:,1] = ones(Int8, 1, numanimals)
+
+# Create an external counter that will increment through both loops
+counter = 2
+
+# Store the traits in a string array
+adjustedtraits =
+Array{String}(undef,floor(Int,sum(numgroups))-length(numgroups))
+# Iterate through each group
+for i in 1:length(normal)
+ # Find the traits that are present in this trait
+ localdata = string.(data[:,i+5])
+ traits = unique(localdata)
+ # Remove the normal version from the analysis
+ effecttraits = traits[findall(x -> x != normal[i], traits)]
+ # Iterate inside of the group
+ for j in 1:length(effecttraits)
+ matchedindex = findall(x -> x != effecttraits[j], localdata)
+ X[matchedindex, counter] .= 1
+ # Add this trait to the string
+ adjustedtraits[counter - 1] = traits[j]
+ # Increment the big counter
+ global counter = counter + 1
+ end
end
-println("If no animal matching this description exists, the results may appear")
-println("outlandish, but are still as correct as the accuracy suggests")
+
+print("Done!\n")
+
+### Additive relationship matrix
+print("[🐮]: Creating additive relationship matrix...")
+
+# Create an empty matrix for the additive relationship matrix
+A = zeros(numanimals, numanimals)
+
+# Create the additive relationship matrix by the FORTRAN method presented by
+# Henderson
+for i in 1:numanimals
+ if !isnothing(dam[i]) && !isnothing(sire[i])
+ for j in 1:(i-1)
+ A[j,i] = 0.5*(A[j,sire[i]] + A[j,dam[i]])
+ A[i,j] = A[j,i]
+ end
+ A[i,i] = 1 + 0.5*A[sire[i], dam[i]]
+ elseif !isnothing(dam[i]) && isnothing(sire[i])
+ for j in 1:(i-1)
+ A[j,i] = 0.5*A[j,dam[i]]
+ A[i,j] = A[j,i]
+ end
+ A[i,i] = 1
+ elseif isnothing(dam[i]) && !isnothing(sire[i])
+ for j in 1:(i-1)
+ A[j,i] = 0.5*A[j,sire[i]]
+ A[i,j] = A[j,i]
+ end
+ A[i,i] = 1
+ else
+ for j in 1:(i-1)
+ A[j,i] = 0
+ A[i,j] = 0
+ end
+ A[i,i] = 1
+ end
+end
+
+print("Done!\n")
+
+### Perform BLUP
+print("[🐮]: Solving the mixed-model equations...")
+
+# Extract the observed data
+Y = convert(Array{Float64}, data[:,5])
+
+# The random effects matrix
+Z = Matrix{Int}(I, numanimals, numanimals)
+
+# Remove items where there is no data
+nullobs = findall(isnothing, Y)
+Z[nullobs, nullobs] .= 0
+
+# Calculate heritability
+λ = (1-h2)/h2
+
+# Use the mixed-model equations
+MME = [X'*X X'*Z; Z'*X (Z'*Z)+(inv(A).*λ)]
+MMY = [X'*Y; Z'*Y]
+solutions = MME\MMY
+
+# Find the accuracies
+diaginv = diag(inv(MME))
+reliability = ones(Float64, length(diaginv)) - diaginv.*λ
+
+print("Done!\n")
+
+### Output the results
+print("[🐮]: Saving results...")
+
+# Find how many traits we found BLUE for
+numgroups = numgroups .- 1
+
+# Start printing results to output
+fileID = open(savepath, "w")
+write(fileID, "beefblup Results Report\n")
+write(fileID, "Produced using beefblup for Julia (")
+write(fileID, "https://github.com/millironx/beefblup")
+write(fileID, ")\n\n")
+write(fileID, "Input:\t")
+write(fileID, path)
+write(fileID, "\nAnalysis performed:\t")
+write(fileID, string(Dates.today()))
+write(fileID, "\nTrait examined:\t")
+write(fileID, headers[5])
+write(fileID, "\n\n")
+
+# Print base population stats
+write(fileID, "Base Population:\n")
+for i in 1:length(numgroups)
+ write(fileID, "\t")
+ write(fileID, headers[i+5])
+ write(fileID, ":\t")
+ write(fileID, normal[i])
+ write(fileID, "\n")
+end
+write(fileID, "\tMean ")
+write(fileID, headers[5])
+write(fileID, ":\t")
+write(fileID, string(solutions[1]))
+write(fileID, "\n\n")
+
+# Contemporary group adjustments
+counter = 2
+write(fileID, "Contemporary Group Effects:\n")
+for i in 1:length(numgroups)
+ write(fileID, "\t")
+ write(fileID, headers[i+5])
+ write(fileID, "\tEffect\tReliability\n")
+ for j in 1:numgroups[i]
+ write(fileID, "\t")
+ write(fileID, adjustedtraits[counter - 1])
+ write(fileID, "\t")
+ write(fileID, string(solutions[counter]))
+ write(fileID, "\t")
+ write(fileID, string(reliability[counter]))
+ write(fileID, "\n")
+
+ global counter = counter + 1
+ end
+ write(fileID, "\n")
+end
+write(fileID, "\n")
+
+# Expected breeding values
+write(fileID, "Expected Breeding Values:\n")
+write(fileID, "\tID\tEBV\tReliability\n")
+for i in 1:numanimals
+ write(fileID, "\t")
+ write(fileID, ids[i])
+ write(fileID, "\t")
+ write(fileID, string(solutions[i+counter-1]))
+ write(fileID, "\t")
+ write(fileID, string(reliability[i+counter-1]))
+ write(fileID, "\n")
+end
+
+write(fileID, "\n - END REPORT -")
+close(fileID)
print("Done!\n")
diff --git a/Julia/install.jl b/Julia/install.jl
index 3276384..1a3698d 100644
--- a/Julia/install.jl
+++ b/Julia/install.jl
@@ -9,4 +9,5 @@
using Pkg
# Install requisite packages
-Pkg.add("XLSX")
\ No newline at end of file
+Pkg.add("XLSX")
+Pkg.add("Gtk")
\ No newline at end of file