diff --git a/.idea/dictionaries/cclea.xml b/.idea/dictionaries/cclea.xml deleted file mode 100644 index 8c49cbc..0000000 --- a/.idea/dictionaries/cclea.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - beefblup - blup - matlab - - - \ No newline at end of file diff --git a/.idea/workspace.xml b/.idea/workspace.xml deleted file mode 100644 index bbe788f..0000000 --- a/.idea/workspace.xml +++ /dev/null @@ -1,167 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1567039339175 - - - 1567130353481 - - - 1567131158833 - - - 1567131230009 - - - 1567131261939 - - - 1567133774405 - - - - - - - - - - - - - - - - - file://$PROJECT_DIR$/Python/beefblup.py - 44 - - - file://$PROJECT_DIR$/Python/beefblup.py - 59 - - - file://$PROJECT_DIR$/Python/beefblup.py - 71 - - - - - \ No newline at end of file diff --git a/Excel/van der Werf2.xlsx b/Excel/van der Werf2.xlsx new file mode 100644 index 0000000..e788425 Binary files /dev/null and b/Excel/van der Werf2.xlsx differ 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