1
0
Fork 0
mirror of https://github.com/MillironX/beefblup.git synced 2024-11-10 18:23:08 +00:00
beefblup/MATLAB/beefblup.m

183 lines
No EOL
5.2 KiB
Matlab

% beefblup
% Main script for performing single-variate BLUP to find beef cattle
% breeding values
% Usage: beefblup
% (C) 2018 Thomas A. Christensen II
% Licensed under BSD-3-Clause License
% Prepare the workspace for computation
clear
clc
close all
% Import data from a suitable spreadsheet
[name, path] = uigetfile('*.xlsx','Select a beefblup worksheet');
fullname = [path name];
clear name path
[~, ~, data] = xlsread(fullname);
% Extract the headers into a separate array
headers = data(1,:);
data(1,:) = [];
% Convert the string dates to numbers
data(:,2) = num2cell(datenum(data(:,2)));
% Sort the array by date
data = sortrows(data,2);
% Coerce all id fields to string format
strings = [1 3 4];
data(:,strings) = cellfun(@num2str, data(:,strings), 'UniformOutput', false);
% Create a lookup lambda function to find the animal represented by a
% certain id
animalrow = @(id) find(strcmp(data(:,1), id));
ids = [data{:,1}];
numanimals = length(data(:,1));
% Store column numbers that need to be deleted
% Column 6 contains an intermediate Excel calculation and always needs to
% be deleted
colstodelete = 6;
% Coerce each group to string format
for i = 7:length(headers)
data(:,i) = cellfun(@num2str, data(:,i), 'UniformOutput', false);
end
% Find any columns that need to be deleted
for i = 7:length(headers)
if length(uniquecell(data(:,i))) <= 1
colname = headers{i};
disp(['Column "' colname '" does not have any unique animals and will be removed'])
disp('from this analysis');
colstodelete = [colstodelete i];
end
end
% Delete the appropriate columns from the datasheet and the headers
data(:,colstodelete) = [];
headers(colstodelete) = [];
% Determine how many contemporary groups there are
numgroups = ones(1, length(headers)-5);
for i = 6:length(headers)
numgroups(i-5) = length(uniquecell(data(:,i)));
end
% If there are more groups than animals, then the analysis cannot continue
if sum(numgroups) >= numanimals
disp('There are more contemporary groups than animals. The analysis will now abort.');
return
end
% Define a "normal" animal as one of the last in the groups, provided that
% all traits do not have null values
normal = cell([1 length(headers)-5]);
for i = 6:length(headers)
for j = numanimals:-1:1
if not(cellfun(@isempty, data(j,i)))
normal(i - 5) = data(j,i);
break
end
end
end
% Print the results of the "normal" definition
disp(' ')
disp('For the purposes of this analysis, a "normal" animal will be defined')
disp('by the following traits:')
for i = 6:length(headers)
disp([headers{i} ': ' normal{i-5}])
end
disp(' ')
disp('If no animal matching this description exists, the results may appear')
disp('outlandish, but are still as correct as the accuracy suggests')
disp(' ')
% Form the fixed effect matrix
X = zeros(numanimals, sum(numgroups)-length(numgroups)+1);
X(:,1) = ones(1, numanimals);
% Create an external counter that will increment through both loops
I = 2;
% Iterate through each group
for i = 1:length(normal)
% Find the traits that are present in this trait
traits = uniquecell(data(:,i+5));
% Remove the "normal" version from the analysis
normalindex = find(strcmp(traits, normal{i}));
traits(normalindex) = [];
% Iterate inside of the group
for j = 1:length(traits)
matchedindex = find(strcmp(data(:,i+5), traits{j}));
X(matchedindex, I) = 1;
I = I + 1;
end
end
% Create an empty matrix for the additive relationship matrix
A = sparse(numanimals, numanimals);
% Create lambdas to find sire and dam of each animal
dam = @(id) [data{animalrow(num2str(id)), 3}];
sire = @(id) [data{animalrow(num2str(id)), 4}];
% Create the additive relationship matrix by the FORTRAN method presented
% by Henderson
for i = 1:numanimals
if strcmp(sire(i), 'NaN') && strcmp(dam(i), 'NaN')
for j = 1:(i-1)
A(j,i) = 0;
end
A(i,i) = 1;
continue
end
if strcmp(sire(i), 'NaN') && not(strcmp(dam(i), 'NaN'))
for j = 1:(i-1)
A(j,i) = A(j,str2double(dam(i)));
A(i,j) = A(j,str2double(dam(i)));
end
A(i,i) = 1;
continue
end
if strcmp(dam(i), 'NaN') && not(strcmp(sire(i), 'NaN'))
for j=1:(i-1)
A(j,i) = A(j,str2double(sire(i)));
A(i,j) = A(j,str2double(sire(i)));
end
A(i,i) = 1;
continue
end
for j=1:(i-1)
A(j,i) = 0.5.*(A(j,str2double(sire(i))) + A(j,str2double(dam(i))));
A(i,j) = 0.5.*(A(j,str2double(sire(i))) + A(j,str2double(dam(i))));
end
A(i,i) = 1 + 0.5.*A(str2double(sire(i)), str2double(dam(i)));
continue
end
% Extract the observed data
Y = data{:, 5};
% The identity matrix for random effects
Z = speye(numanimals, numanimals);
% Prompt for heritablity
h2 = str2double(input('What is the heritablity for this trait? >> '));
lambda = (1-h2)/h2;
% Use the mixed-model equations
solutions = [X'*X X'*Z; Z'*X (Z'*Z)+(inv(A).*lambda)]\[X'*Y; Z'*Y];
% Find the accuracies
diaginv = diag(inv([X'*X X'*Z; Z'*X (Z'*Z)+(inv(A).*lambda)]));
reliability = 1 - diaginv.*lambda;