BLUPs.sas

************************************************************************************
GWAS BLUPs 2020 
Author: Mason Chizk
Date Last Modified: 04/14/2021

Description: This script takes a preformatted SAS dataset and estimates BLUPs,
outputting the results to a csv formatted for GWASpoly.  The code below provides a
basic framework, but each block of code should be inspected closely, as some of it
will need to be tailored to your specific scenerio.  Don't forget to check the log
for errors as you go!  Good luck!
***********************************************************************************;


************************************************************************************
STEP 1: This block assumes that you already have a phenotypic dataset called 
"GWAS_Pheno" read into SAS (since this will be very unique to each scenerio).  The 
model used in the mixed procedure below also assumes that the trait measured was 
observed over multiple years and reps. If this is not the case for your trait, simply 
remove the irrelevant terms from the model. Don't forget to use your own data and 
replace "trait1" with whatever trait you are using
***********************************************************************************;

proc mixed data=GWAS_Pheno method=Type3;
	class Year Rep Geno;
	model trait1=;
	* The "solution" option used below is what requests the BLUPs;
	random Year Rep(Year) Geno Year*Geno / solution;
	* ods output takes the SAS table called "SolutionR", which was 
	generated by "solution" and outputs it to the new "BLUPtrait1"
	dataset;
	ods output Mixed.SolutionR = BLUPtrait1;
run;

************************************************************************************
STEP 2: This block takes the raw BLUPs dataset from step 1 and filters out all of the
non-genotypic BLUPs (i.e. BLUPs for year, rep, and interactions).  Don't forget to 
replace "trait1" as needed
***********************************************************************************;

data Bluptrait1;
	*(where=()) is one of the most helpful snippets of code that I know of in SAS.
	It filters out data from the input dataset based on conditional logic;
	set Bluptrait1 (where=(Effect eq "Geno"));
	*keeping only the data that we want;
	keep Geno Estimate;
	rename Estimate=trait1;
run;

************************************************************************************
STEP 3: (ONLY IF USING MULTIPLE TRAITS!!!) If you are running BLUPs for multiple traits,
it will be convenient to go ahead and combine the data in SAS.  First, sort the genotypic
BLUPs by "Geno", or the name of the genotypes, since this is the order they will be combined
in.  And then, combine using a data merge step.  The "sortall" macro below is a nice way to 
perform multiple sorts in one run.
***********************************************************************************;

%macro sortall(tables,byvar);
%local i n table;
%let n=%sysfunc(countw(&tables));
%do i=1 %to &n;
	%let table=%scan(&tables,&i);
	proc sort data=&table;
	by &byvar;
	run;
%end;
%mend;

%sortall(BLUPtrait1 BLUPtrait2 BLUPtrait3, Geno);

data All_BLUPs;
	merge BLUPtrait1 BLUPtrait2 BLUPtrait3;
	by Geno;
run;

************************************************************************************
STEP 4: Write the All_BLUPs SAS dataset to a CSV file for GWASpoly
***********************************************************************************;

proc export data=All_BLUPs
	dbms="csv"
	outfile="whatever_your_path_is.csv";
run;