1
0
Fork 0
mirror of https://github.com/MillironX/beefblup.git synced 2024-09-20 21:12:03 +00:00

Complete port of beefblup to Julia

This commit is contained in:
Thomas A. Christensen II 2020-02-03 18:11:55 -08:00
parent 44ec5ae1da
commit fc771d338a
Signed by: millironx
GPG key ID: 139C07724802BC5D
5 changed files with 194 additions and 202 deletions

View file

@ -1,9 +0,0 @@
<component name="ProjectDictionaryState">
<dictionary name="cclea">
<words>
<w>beefblup</w>
<w>blup</w>
<w>matlab</w>
</words>
</dictionary>
</component>

View file

@ -1,167 +0,0 @@
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ChangeListManager">
<list default="true" id="816df031-a834-466c-aafb-fce0c72415a4" name="Default Changelist" comment="">
<change beforePath="$PROJECT_DIR$/.gitignore" beforeDir="false" afterPath="$PROJECT_DIR$/.gitignore" afterDir="false" />
<change beforePath="$PROJECT_DIR$/.idea/.gitignore" beforeDir="false" afterPath="$PROJECT_DIR$/.idea/.gitignore" afterDir="false" />
<change beforePath="$PROJECT_DIR$/.idea/beefblup.iml" beforeDir="false" afterPath="$PROJECT_DIR$/.idea/beefblup.iml" afterDir="false" />
<change beforePath="$PROJECT_DIR$/.idea/inspectionProfiles/profiles_settings.xml" beforeDir="false" afterPath="$PROJECT_DIR$/.idea/inspectionProfiles/profiles_settings.xml" afterDir="false" />
<change beforePath="$PROJECT_DIR$/.idea/misc.xml" beforeDir="false" afterPath="$PROJECT_DIR$/.idea/misc.xml" afterDir="false" />
<change beforePath="$PROJECT_DIR$/.idea/modules.xml" beforeDir="false" afterPath="$PROJECT_DIR$/.idea/modules.xml" afterDir="false" />
<change beforePath="$PROJECT_DIR$/.idea/vcs.xml" beforeDir="false" afterPath="$PROJECT_DIR$/.idea/vcs.xml" afterDir="false" />
<change beforePath="$PROJECT_DIR$/.vscode/launch.json" beforeDir="false" afterPath="$PROJECT_DIR$/.vscode/launch.json" afterDir="false" />
<change beforePath="$PROJECT_DIR$/Python/beefblup.py" beforeDir="false" afterPath="$PROJECT_DIR$/Python/beefblup.py" afterDir="false" />
</list>
<option name="EXCLUDED_CONVERTED_TO_IGNORED" value="true" />
<option name="SHOW_DIALOG" value="false" />
<option name="HIGHLIGHT_CONFLICTS" value="true" />
<option name="HIGHLIGHT_NON_ACTIVE_CHANGELIST" value="false" />
<option name="LAST_RESOLUTION" value="IGNORE" />
</component>
<component name="Git.Settings">
<option name="RECENT_BRANCH_BY_REPOSITORY">
<map>
<entry key="$PROJECT_DIR$" value="develop" />
</map>
</option>
<option name="RECENT_GIT_ROOT_PATH" value="$PROJECT_DIR$" />
</component>
<component name="IgnoredFileRootStore">
<option name="generatedRoots">
<set>
<option value="C:\Users\cclea\source\repos\beefblup\.idea" />
</set>
</option>
</component>
<component name="ProjectId" id="1Q4ogIZMNQ0ZN4Ces5FEsodDyuh" />
<component name="PropertiesComponent">
<property name="last_opened_file_path" value="$PROJECT_DIR$/Python/beefblup.py" />
<property name="settings.editor.selected.configurable" value="configurable.group.language" />
</component>
<component name="RunDashboard">
<option name="ruleStates">
<list>
<RuleState>
<option name="name" value="ConfigurationTypeDashboardGroupingRule" />
</RuleState>
<RuleState>
<option name="name" value="StatusDashboardGroupingRule" />
</RuleState>
</list>
</option>
</component>
<component name="RunManager">
<configuration name="beefblup" type="PythonConfigurationType" factoryName="Python" nameIsGenerated="true">
<module name="beefblup" />
<option name="INTERPRETER_OPTIONS" value="" />
<option name="PARENT_ENVS" value="true" />
<envs>
<env name="PYTHONUNBUFFERED" value="1" />
</envs>
<option name="SDK_HOME" value="C:\tools\Anaconda3\python.exe" />
<option name="WORKING_DIRECTORY" value="$PROJECT_DIR$/Python" />
<option name="IS_MODULE_SDK" value="false" />
<option name="ADD_CONTENT_ROOTS" value="true" />
<option name="ADD_SOURCE_ROOTS" value="true" />
<option name="SCRIPT_NAME" value="$PROJECT_DIR$/Python/beefblup.py" />
<option name="PARAMETERS" value="" />
<option name="SHOW_COMMAND_LINE" value="false" />
<option name="EMULATE_TERMINAL" value="false" />
<option name="MODULE_MODE" value="false" />
<option name="REDIRECT_INPUT" value="false" />
<option name="INPUT_FILE" value="" />
<method v="2" />
</configuration>
</component>
<component name="SvnConfiguration">
<configuration />
</component>
<component name="TaskManager">
<task active="true" id="Default" summary="Default task">
<changelist id="816df031-a834-466c-aafb-fce0c72415a4" name="Default Changelist" comment="" />
<created>1567039339175</created>
<option name="number" value="Default" />
<option name="presentableId" value="Default" />
<updated>1567039339175</updated>
</task>
<task id="LOCAL-00001" summary="Added all results files to gitignore">
<created>1567130353481</created>
<option name="number" value="00001" />
<option name="presentableId" value="LOCAL-00001" />
<option name="project" value="LOCAL" />
<updated>1567130353481</updated>
</task>
<task id="LOCAL-00002" summary="Add new gitignore to avoid uploading junk from PyCharm">
<created>1567131158833</created>
<option name="number" value="00002" />
<option name="presentableId" value="LOCAL-00002" />
<option name="project" value="LOCAL" />
<updated>1567131158833</updated>
</task>
<task id="LOCAL-00003" summary="Added PyCharm's config files">
<created>1567131230009</created>
<option name="number" value="00003" />
<option name="presentableId" value="LOCAL-00003" />
<option name="project" value="LOCAL" />
<updated>1567131230009</updated>
</task>
<task id="LOCAL-00004" summary="Create initial python script">
<created>1567131261939</created>
<option name="number" value="00004" />
<option name="presentableId" value="LOCAL-00004" />
<option name="project" value="LOCAL" />
<updated>1567131261939</updated>
</task>
<task id="LOCAL-00005" summary="Added Julius van der Werf's example to test agains">
<created>1567133774405</created>
<option name="number" value="00005" />
<option name="presentableId" value="LOCAL-00005" />
<option name="project" value="LOCAL" />
<updated>1567133774405</updated>
</task>
<option name="localTasksCounter" value="6" />
<servers />
</component>
<component name="Vcs.Log.Tabs.Properties">
<option name="TAB_STATES">
<map>
<entry key="MAIN">
<value>
<State>
<option name="COLUMN_ORDER" />
</State>
</value>
</entry>
</map>
</option>
</component>
<component name="VcsManagerConfiguration">
<MESSAGE value="Added all results files to gitignore" />
<MESSAGE value="Add new gitignore to avoid uploading junk from PyCharm" />
<MESSAGE value="Added PyCharm's config files" />
<MESSAGE value="Create initial python script" />
<MESSAGE value="Added Julius van der Werf's example to test agains" />
<option name="LAST_COMMIT_MESSAGE" value="Added Julius van der Werf's example to test agains" />
</component>
<component name="XDebuggerManager">
<breakpoint-manager>
<breakpoints>
<line-breakpoint enabled="true" suspend="THREAD" type="python-line">
<url>file://$PROJECT_DIR$/Python/beefblup.py</url>
<line>44</line>
<option name="timeStamp" value="4" />
</line-breakpoint>
<line-breakpoint enabled="true" suspend="THREAD" type="python-line">
<url>file://$PROJECT_DIR$/Python/beefblup.py</url>
<line>59</line>
<option name="timeStamp" value="5" />
</line-breakpoint>
<line-breakpoint enabled="true" suspend="THREAD" type="python-line">
<url>file://$PROJECT_DIR$/Python/beefblup.py</url>
<line>71</line>
<option name="timeStamp" value="6" />
</line-breakpoint>
</breakpoints>
</breakpoint-manager>
</component>
</project>

BIN
Excel/van der Werf2.xlsx Normal file

Binary file not shown.

View file

@ -7,6 +7,9 @@
# Import the required packages # Import the required packages
using XLSX using XLSX
using LinearAlgebra
using Dates
using Gtk
# Display stuff # Display stuff
println("beefblup v 0.0.0.1") println("beefblup v 0.0.0.1")
@ -16,19 +19,23 @@ print("\n")
### Prompt User ### Prompt User
# Ask for an input spreadsheet # Ask for an input spreadsheet
# print("Enter the full filename of a beefblup worksheet> ") path = open_dialog_native(
# path = readline(stdin) "Select a beefblup worksheet",
path = "C:\\Users\\cclea\\source\\repos\\beefblup\\Excel\\Master BLUP Worksheet.xlsx" GtkNullContainer(),
("*.xlsx", GtkFileFilter("*.xlsx", name="beefblup worksheet"))
)
# Ask for an output text filename # Ask for an output text filename
# print("Enter the full filename of your desired results> ") savepath = save_dialog_native(
# savepath = readline(stdin) "Save your beefblup results",
savepath = "C:\\Users\\cclea\\source\\repos\\beefblup\\results.txt" GtkNullContainer(),
(GtkFileFilter("*.txt", name="Results file"),
"*.txt")
)
# Ask for heritability # Ask for heritability
# print("What is the heritability for this trait?> ") print("What is the heritability for this trait?> ")
# h2 = parse(Float64, readline(stdin)) h2 = parse(Float64, readline(stdin))
h2 = 0.4
### Import input filename ### Import input filename
print("[🐮]: Importing Excel file...") print("[🐮]: Importing Excel file...")
@ -39,7 +46,7 @@ data = XLSX.readxlsx(path)[1][:]
print("Done!\n") print("Done!\n")
### Process input file ### Process input file
print("[🐮]: Processing and formatting data ...") print("[🐮]: Processing and formatting data...")
# Extract the headers into a separate array # Extract the headers into a separate array
headers = data[1,:] headers = data[1,:]
@ -49,9 +56,9 @@ data = data[2:end,:]
data = sortslices(data, dims=1, lt=(x,y)->isless(x[2],y[2])) data = sortslices(data, dims=1, lt=(x,y)->isless(x[2],y[2]))
# Define fields to hold id values for animals and their parents # Define fields to hold id values for animals and their parents
ids = data[:,1] ids = string.(data[:,1])
damids = data[:,3] damids = string.(data[:,3])
sireids = data[:,4] sireids = string.(data[:,4])
numanimals = length(ids) numanimals = length(ids)
# Find the index values for animals and their parents # 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 i in 6:length(headers)
for j in numanimals:-1:1 for j in numanimals:-1:1
if !ismissing(data[j,i]) if !ismissing(data[j,i])
normal[i-5] = data[j,i] normal[i-5] = string(data[j,i])
break break
end end
end end
end end
# Print the results of the "normal" definition print("Done!\n")
println(" ")
println("For the purposes of this analysis, a 'normal' animal will be defined by") ### Create the fixed-effect matrix
println("the following traits:") print("[🐮]: Creating the fixed-effect matrix...")
for i in 6:length(headers)
print(headers[i]) # Form the fixed-effect matrix
print(": ") X = zeros(Int8, numanimals, floor(Int,sum(numgroups))-length(numgroups)+1)
print(normal[i-5]) X[:,1] = ones(Int8, 1, numanimals)
print("\n")
# 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 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") print("Done!\n")

View file

@ -9,4 +9,5 @@
using Pkg using Pkg
# Install requisite packages # Install requisite packages
Pkg.add("XLSX") Pkg.add("XLSX")
Pkg.add("Gtk")