function [F, cv] = mainscan(y, cov, geno, nperm)
%scan a marker array for main-effect associations with a trait
%mainscan uses an unconstrained regression model
%   see f2scan for constrained models
%INPUT VARIABLES ARE:
%     y: nx1 vector of trait values
%   cov: nx1 vector of covariates,  use [] if none
%  geno: nxm genotype array. 
% nperm: number of permutations
%NOTES:
%this program handles F2 (0 1 2 coding), BC and RI (0 1 coding)
%stripped version...missing values are not allowed

F = genoscan(y, cov, geno);

%permutation test
if nargin > 3 & nperm > 0
   Fmax = permuteF(y, cov, geno, nperm);
  	cv=Fmax([ceil(.8*nperm),ceil(.9*nperm),ceil(.95*nperm),ceil(.99*nperm)]);
end
  
%%% SUBROUTINE GENOSCAN %%% 
function F = genoscan(y, cov, geno)

%constants and storage
[n,m] = size(geno);
F = zeros(m,1);

%set up F2 indicator
if max(geno(:))>1
   F2 = 1;
else
   F2 = 0;
end

%the null model regression
X0 = [ones(n,1),cov];
[rss0, df0] = RSS(y, X0);

%scan the genome doing regression at each marker
if F2
	for i = 1:m  
      X = [X0, geno(:,i)==1, geno(:,i)==2];
   	[rss1, df1] = RSS(y, X);
      F(i) = ((rss0-rss1)/(df0-df1))/(rss1/df1);
   end
else
 	for i = 1:m  
      X = [X0, geno(:,i)];
      [rss1, df1] = RSS(y, X);
      F(i) = ((rss0-rss1)/(df0-df1))/(rss1/df1);
   end
end
  

%%% SUBROUTINE RSS %%%
function [rss,df] = RSS(y,X)
[n,p] = size(X);
df = n - p;			
[Q, R]=qr(X,0);
b = R\(Q'*y);
yhat = X*b;		
rss = (y-yhat)'*(y-yhat);

%%% SUBROUTINE PERMUTE %%%
function Fmax = permuteF(y, cov, geno, nperm)
Fmax = zeros(nperm,1);
for j = 1:nperm 
   
   %permute the trait
   ord = randperm(length(y));
   ptrait = y(ord);
   
   %permute the covariates
   if ~(isempty(cov))
      pcov = cov(ord,:);
   else
      pcov = [];
   end

   Fmax(j) = max(genoscan(ptrait, pcov, geno));
end
Fmax = sort(Fmax);