具有遗传性疾病和性状的遗传位点分析MATLAB代码

有代码，但运行不出结果，谁可以帮忙看看怎么改一下代码，运行出结果

% This is the program for problem2
% usage: this file contains many parts, you can run every part at once
% by removing or
% commenting other parts, pay attention not to run one part for many times!
% part1:load origin data file(genotypetype.dat) and convert it to cell

% initial
p_value=zeros(9445,1);
rs_str={'AA','CC','GG','TT','AC','AG','AT','CA','CG','CT',...
'GA','GC','GT','TA','TC','TG'};
% part2:encode the origin data to 0 1 2
for j=1:9445
temp1=17;
temp2=17;
for i=1:1000
for k=1:16
if strcmp(mat(i,j),rs_str(k)) == 1
if k<5
if temp1==17
temp1=k;
break;
elseif k~=temp1
temp2=k;
break;
end
end
end
end
if temp2 ~= 17
break;
end
end

if temp1>temp2
temp1=temp2;
end

for i=1:1000
for k=1:16
if strcmp(mat(i,j),rs_str(k)) == 1
if k==temp1
mat(i,j)=num2cell(0);
elseif k>4
mat(i,j)=num2cell(1);
else
mat(i,j)=num2cell(2);
end
end
end
end
end
% part2 end

% part3:remove the dirty data using WAF
for i=1:9445
[a,] = size(find(cell2mat(mat(:,i))==0));
[b,] = size(find(cell2mat(mat(:,i))==1));
[c,~] = size(find(cell2mat(mat(:,i))==2));
% if p_value is less than 0.05, it is dirty
if a<c
p_value(i,1)=(2a+b)/2/(a+b+c);
else
p_value(i,1)=(2c+b)/2/(a+b+c);
end
end
% part3 end

% part4:remove the dirty data using Hary-Weinberg
chi=zeros(9445,1);
for i=1:9445
[a,] = size(find(cell2mat(mat(:,i))==0));
[b,] = size(find(cell2mat(mat(:,i))==1));
[c,~] = size(find(cell2mat(mat(:,i))==2));
ex_a = ((2a+b)^2)/(4(a+b+c));
ex_b = (a+b+c)-((2a+b)^2)/4/(a+b+c)-((2c+b)^2)/4/(a+b+c);
ex_c = ((2*c+b)^2)/4/(a+b+c);
% if chi value is more than 6.64, it is dirty
chi(i,1) = (a-ex_a)^2 / ex_a + (b-ex_b)^2 / ex_b + (c-ex_c)^2 / ex_c;
end
% part4 end
% part6: plot using manhatten
manhattan=zeros(9445,1);
for i=1:9445
manhattan(i,1)=-log10(p_value(i,1));
end
plot(manhattan,'.');
%part6 end