%a example of a QTL analysis script
data_in

%setup phenotypes for genome scans
traits = pheno(:,2:6);
names = pnames(2:6);
ntraits = size(traits,2);

%initialize storage for genome scans
Fmn = zeros(nmarks, ntraits);
cvmn = zeros(4,ntraits);

%loop through the phenotypes
for j = 1:ntraits
   obsdx = find(~isnan(traits(:,j)));  %set missing data index
	[Fmn(:,j), cvmn(:,j)] = mainscan(traits(obsdx,j), [], genoi(obsdx,:), 1000);
end
uplim = max(Fmn(:));
save results01

%draw the main scan figures
figure
for j = 1:ntraits
   subplot(5,1,j)
   drawmainscan(Fmn(:,j), cvmn(:,j), chrid, [1:19], names(j)) 
   %axis([0 nmarks+1 0 uplim])
end

%initialize storage for pairwise genome scans
F1 = zeros(nmarks,nmarks,ntraits);
F2 = zeros(nmarks,nmarks,ntraits);
F3 = zeros(nmarks,nmarks,ntraits);
F4 = zeros(nmarks,nmarks,ntraits);
cvpair = zeros(4,ntraits);

%loop through the phenotypes
for j = 1:ntraits
   obsdx = find(~isnan(traits(:,j)));
   [F1(:,:,j), F2(:,:,j), F3(:,:,j), F4(:,:,j)] = ...
   	pairwise(traits(obsdx,j), [], genoi(obsdx,:), chrid);
   %cvpair(:,j) = pairpermute(traits(obsdx,j), [], genoi(obsdx,:), chrid, 100)
end
save results01

for j = 1:ntraits
   figure
   drawpairscan(F1(:,:,j), F2(:,:,j), chrid, [1:19], names(j))
   %figure
   %drawpairscan(F3(:,:,j), F4(:,:,j), chrid, [1:19], strcat('Beamer B6 x C3H',' ',names(j)))
end

fid = fopen('pair_results01.txt','w');
for j = 1:ntraits
   fprintf(fid,'%s\r\n',strcat('Beamer B6 x C3H',' ',names{j}))
   cv = [4.7 3.74 5.33 5.33];
   reportpairs(F1(:,:,j), F2(:,:,j), F3(:,:,j), F4(:,:,j), cv, mnames, chrid, fid)
   fprintf(fid,'\r\n')
end 
fclose(fid)
save results01