/**********************************************************************************************/ /*********************** SYNCHRO1: Input data and checks ***************************************/ /**********************************************************************************************/ %macro SYNCHRO1(Message); options nonotes mautosource NOmerror nomprint noserror nosource nodate; %put SYNCHRO START; %let beta=0.02; %let ERR=;%let libname=;%let table=;%let block=;%let branch=; %do i=1 %to 8; %let m&i=0; %let f&i=0;%let ud=3; %end; %empieza: %window START color=gray irow=10 rows=23 icolumn=20 columns=32 #1 @2 "&Message" color=red attr=highlight #3 @2 "Table:" +1 libname 8 attr=underline color=blue "." table 8 attr=underline required=yes color=blue #5 @2 "Start Day:" +1 inicio 3 attr=underline required=yes color=blue #5 @18 "End Day:" +1 fin 3 attr=underline required=yes color=blue #7 @2 "MALE STAGES" attr=underline +5 "FEMALE STAGES" attr=underline #8 @5 "1" +1 m1 3 color=blue attr=underline +12 "1" +1 f1 3 color=blue attr=underline #9 @5 "2" +1 m2 3 color=blue attr=underline +12 "2" +1 f2 3 color=blue attr=underline #10 @5 "3" +1 m3 3 color=blue attr=underline +12 "3" +1 f3 3 color=blue attr=underline #11 @5 "4" +1 m4 3 color=blue attr=underline +12 "4" +1 f4 3 color=blue attr=underline #12 @5 "5" +1 m5 3 color=blue attr=underline +12 "5" +1 f5 3 color=blue attr=underline #13 @5 "6" +1 m6 3 color=blue attr=underline +12 "6" +1 f6 3 color=blue attr=underline #14 @5 "7" +1 m7 3 color=blue attr=underline +12 "7" +1 f7 3 color=blue attr=underline #15 @5 "8" +1 m8 3 color=blue attr=underline +12 "8" +1 f8 3 color=blue attr=underline #17 @7 "Press" +1 "F3" attr=underline +1 "when ready"; %display START; %if &libname= %then %let libname=work; %let message=; /**************** Data verifications *********************************************************/ PROC contents data=&libname.._ALL_ memtype=data out=verify noprint; data _null_; set verify end=eof; if memname=upcase("&table") then ERR=1; retain ERR; if eof and ERR^=1 then call symput('Message',"ERROR: Table does not exist"); run; %if %length(&Message)>0 %then %goto empieza; data _null_; do i=1 to 8 until(_error_); if symget('m'||left(i))>100 then call symput('Message',"ERROR: m"||substr(left(i),1,1)||" > 100"); if _ERROR_=1 then do;call symput('Message',"ERROR: m"||substr(left(i),1,1)||" not numeric");goto adeus;end; if symget('f'||left(i))>100 then call symput('Message',"ERROR: f"||substr(left(i),1,1)||" > 100"); if _ERROR_=1 then call symput('Message',"ERROR: f"||substr(left(i),1,1)||" not numeric"); end; if symget('inicio')>365 or symget('fin')>365 then call symput('Message',"ERROR: Start/End > 365"); if _ERROR_=1 then call symput('Message',"ERROR: Start/End not numeric"); adeus: run; %if %length(&Message)>0 %then %goto empieza; %if &inicio>=&fin %then %SYNCHRO1(ERROR: Start>=End); PROC contents data=&libname..&TABLE out=verify noprint; data verify; set verify; if substr(name,1,1)='D' then d=0+substr(name,2,3); proc sort; by d; data _null; set verify end=eof; if _n_=1 then do;day=0;id=0;block=0;sex=0; branch=0;inicio=0;fin=0;end; if name='ID' then ID=1; if name='RAMET' then ramet=1; if name='BLOCK' then block=1; if name='SEX' and type=2 then sex=1; if name='BRANCHE' then branch=1; if name="D&inicio" then inicio=1; if name="D&fin" then fin=1; if substr(name,1,1)='D' then do; if substr(name,2,3)>0 and substr(name,2,3)<365 then day+1; if type=2 then dayS=1; call symput('Day'||left(day),substr(name,2,3)); end; retain id ramet block sex dayS branch inicio fin; if eof then do; if ID^=1 then call symput('ERR',"ERROR: ID does not exist"); if block^=1 then call symput('BLOCK',1); if ramet^=1 then call symput('ERR',"ERROR: RAMET field does not exist"); if sex^=1 then call symput('ERR',"ERROR: Problems with SEX field"); if branch^=1 then call symput('BRANCH',1); if day<=2 then call symput('ERR',"ERROR: Too low day measurements"); if day=0 then call symput('ERR',"ERROR: No day measurements"); if dayS=1 then call symput('ERR',"ERROR: Some score fields are $"); if inicio^=1 then call symput('Message',"ERROR: D&inicio does not exist"); if fin^=1 then call symput('Message',"ERROR: D&fin does not exist"); call symput('n',day); end; run; %put &err; %if %length(&err)>0 %then %do; %window ERR color=gray irow=13 rows=3 icolumn=7 columns=40 #1 @2 "&ERR" color=red attr=highlight; %display ERR; %goto theend; %end; data pheno; set &libname..&table; ok1=0; if SEX in ('m','f','M','F') then ok1=1; if ok1=0 then call symput('ERR',"ERROR: SEX values ^= mf"); array d %original(D); do i=1 to &n; if d(i)=. or substr(left(D(i)),1,1) in ('1','2','3','4','5','6','7','8') then ok1=1; else call symput('ERR',"ERROR: Score values ^= 1-8"); end; if symget('BLOCK')=1 then BLOCK=1; if symget('BRANCH')=1 then BRANCH=_N_; drop ok1; run; %if %length(&Message)>0 %then %goto empieza; %if &BLOCK=1 %then %put NOTES: BLOCK field has been artificially created and set to 1; %if &BRANCH=1 %then %put NOTES: BRANCH field has been artificially created and set to _N_; %if %length(&err)>0 %then %do; %window ERR color=gray irow=13 rows=3 icolumn=7 columns=40 #1 @2 "&ERR" color=red attr=highlight; %display ERR; %goto theend; %end; %do i=2 %to &n; %let j=%eval(&i-1); %if &&day&i-&&day&j>=5 %then %put NOTES: Too many days between day &&day&j and day &&day&i; %end; %SYNCHRO2 %theend: %mend SYNCHRO1; /**********************************************************************************************/ /*********************** SYNCHRO2: indices computations ****************************************/ /**********************************************************************************************/ /***************** Relative pollen-sheddind and female receptivity ****************************/ %macro SYNCHRO2; data pheno; set pheno end=eof; array D(&n) %original(D); array PCT(&n) %original(PCT); SEX=upcase(SEX); if _N_=1 then miss=0; do i=1 to &n until(flag); if d(i)=. then miss+1; else if i>1 then do; if d(i-1) and D(i)-d(i-1)>2 then do; put 'ERROR: Too large stage change. Record has been deleted. ID=' id 'RAMET=' ramet 'BLOCK=' block 'branch=' branch; delete; flag=1; end; if d(i-1) and D(i)1 then do; if PCT(i)>0 and start=. then start=(symget('day'||left(i))+symget('day'||left(i-1)))/2; if PCT(i)=0 and start^=. and end=. then end=(symget('day'||left(i))+symget('day'||left(i-1)))/2; end; end; if pct(1)>0 then start=&inicio; if pct(&n)>0 then end=&fin; duration=end-start; if sum(of PCT(*))=0 then put 'NOTES: All receptivity/emission rates are null in ID=' id 'RAMET=' ramet 'BLOCK=' block 'BRANCH=' branch; if sum(of PCT(*))=. then put 'NOTES: All receptivity/emission rates are missing in ID=' id 'RAMET=' ramet 'BLOCK=' block 'BRANCH=' branch; if eof and miss>0 then put 'NOTES: ' miss 'score values are missing'; length IDS $ 8; IDS='ID'||substr(upcase(ID),1,6); keep IDS ramet block branch sex %original(PCT) start end duration; /***************** Means on clon basis (&ud=3) or ramet basis (&ud=15) *********************/ proc means data=PHENO noprint; class block ramet sex IDS; var %original(PCT) start end; output out=MEAN mean=%original(PCTM) min(start)=start max(end)=end; data MEAN; set MEAN; duration=end-start;run; /***************** Estimation of POij and RDPij ******************************************/ %let ud=15; %loop1: data M; set MEAN (where=(_type_=&ud)); if sex='M'; union='A'; rename idS=IDM RAMET=RAMETM block=blockM %rename(M) start=startM end=endM; drop sex _type_ _freq_ duration; data F; set MEAN (where=(_type_=&ud)); if sex='F'; union='A'; rename idS=IDF RAMET=RAMETF block=blockF %rename(F) start=startF end=endF; drop sex _type_ _freq_ duration; proc sql; create table MF (drop=union) as select * from M inner join F on M.union=F.union; %let keep=;%let classF=IDF;%let classM=IDM;%let result=MF; %if &ud=15 %then %do; %let classF=rametF;%let classM=rametM; %let result=RESULT2b; %let keep=%str(keep IDF RAMETF blockF IDM RAMETM blockM STARTM ENDM STARTF ENDF PO RD RDP %original(M) %original(F);); %end; data MF; set MF (where=(&classF^=&classM)); array m %original(M); array f %original(F); array den %original(DEN); array num %original(NUM); do i=1 to &n; num(i)=min(m(i),f(i)); den(i)=max(m(i),f(i)); end; DD=sum(of den(*)); NN=sum(of num(*)); PO=NN/DD; RD=min(endM,endF)-max(startF,startM); if RD<0 then RD=0; proc means data=MF noprint; class &classF &classM; var PO RD ; output out=MFM (where=(_type_=1 or _Type_=2)) mean(PO)=PO mean(RD)=RDmean min(PO)=POmin min(RD)=RDmin max(PO)=POmax max(RD)=RDmax sum(rd)=srd; proc sort data=mf; by &classF &classM;run; data &result; merge mf mfm (where=(_type_=2) keep=&classF srd _type_); by &classF; RDP=RD/SRD; &keep run; %if &ud=15 %then %do; data mfm; set mfm;length sex $ 1; if rametM='' then do; RAMET=RAMETF;sex='F';end; if rametF='' then do; RAMET=RAMETM;sex='M';end; proc sort; by sex ramet; proc sort data=mean; by sex ramet; data result2; merge MEAN (where=(_type_=15) keep=_type_ IDS RAMET BLOCK sex %original(PCTM) start end duration) MFM (keep=RAMET sex PO POmin POmax RDmean RDmin RDmax); by sex RAMET; drop _type_; run; %let ud=3; %goto loop1; %end; proc sort data=mean; by IDS; data _NULL_; set mean (where=(_type_=3)) end=eof; by IDS;if _n_=1 then n=0; if first.IDS then do; n+1;call symput('ID'||left(n),left(IDS));end; if eof then call symput('NoID',n); run; %let lista=%all; /***************** OUTPUT: Pollen-shedding and Female flowers receptivity *******************/ title1 'SYNCHRO.sas: A SAS program for analysing phenology synchronicity in seed orchards'; title2 "A total of &NOID clones assessed at &n days were included in calculations"; title3 'Start, duration and end of Pollen-shedding'; title4; proc print data=MEAN (where=(_type_=3 and sex='M') rename=(IDS=ID)); id ID; var start end duration; run; title1 'Start, duration and end of Female flowers receptivity'; title2;title3;title4; proc print data=MEAN (where=(_type_=3 and sex='F') rename=(IDS=ID)); id ID; var start end duration; run; /************************** OUTPUT: ASKEW & BLUSH 1990 **************************************/ proc sort data=mf; by IDF IDM; proc transpose data=MF out=PO; by IDF; var PO; id IDM; title1 'Cross Matrix: PO INDICES (Askew & Blush, 1990)'; title3 'Male'; proc print data=PO (drop=_name_ rename=(IDF=Female)) noobs round; var &lista; id Female; run; title1 'SUMMARY: PO values for all clones serving as female parents';title3; proc means data=MFM (where=(idf^='')) noprint; var PO POmin POmax; output out=POglobal mean(PO)=PO min(POMIN)=POmin max(POMAX)=POmax; data POGLOBAL; set mfm (where=(idf^='')) POglobal;if idf='' then idf='Overall'; proc print data=POGLOBAL (rename=(PO=POmean)) noobs round; id idf; var POmean POmin POmax; run; proc means data=MFM (where=(idm^='')) noprint; var PO POmin POmax; output out=POglobal mean(PO)=PO min(POMIN)=POmin max(POMAX)=POmax; title1 'SUMMARY: PO values for all clones serving as male parents'; data POGLOBAL; set mfm (where=(idm^='')) POglobal;if idm='' then idm='Overall'; proc print data=POGLOBAL (rename=(PO=POmean)) noobs round; id idm; var POmean POmin POmax; run; /***************************** OUTPUT: Xie et al 1994 ****************************************/ proc transpose data=MF out=RDP; by IDF; var RDP; id IDM; title1 'Cross Matrix: RDPij INDICES (Xie et al, 1990)'; title3 'Male'; proc print data=RDP (where=(Female^='Mean') drop=_name_ rename=(IDF=Female)) noobs round; var &lista; id Female; format _numeric_ 4.3; run; title1 'SUMMARY: RD values for all clones serving as female parents';title3; proc print data=MFM (where=(idf^='')) noobs round; id idf; var RDmean RDmin RDmax; run; title1 'SUMMARY: RD values for all clones serving as male parents'; proc print data=MFM (where=(idm^='')) noobs round; id idm; var RDmean RDmin RDmax; run; /***************************** ESTIMATION & OUTPUT: Gömöry et al 2003 *********************/ proc sort data=mean ; by ids sex; proc transpose data=mean (where=(_type_=3)) out=meanori; by ids sex; var %original(PCTM); data M; set meanori; if sex='M'; rename col1=Mvalue IDS=IDM; data F; set meanori; if sex='F'; rename col1=Fvalue IDS=IDF; proc sql; create table MF2 as select M.IDM,M.Mvalue,M._name_,F.IDF,F.Fvalue from M inner join F on M._name_=F._name_; proc sort data=MF2; by IDM _name_ Mvalue; proc means data=MF2 (where=(IDF^=IDM)) noprint; by IDM _name_ mvalue; var Fvalue; output out=corr mean=Fmean; proc corr data=corr outp=M noprint; by IDM; var Mvalue Fmean; proc sort data=MF2; by IDF _name_ Fvalue; proc means data=MF2 (where=(IDF^=IDM)) noprint; by IDF _name_ Fvalue; var Mvalue; output out=corr mean=Mmean; proc corr data=corr outp=F noprint; by IDF; var Fvalue Mmean; proc sort data=MF2; by IDF IDM _name_;run; proc corr data=MF2 outP=CORR noprint; by IDF IDM; var Fvalue Mvalue; proc transpose data=CORR (where=(_type_='CORR' and upcase(_name_)='FVALUE')) out=RPH; by IDF; var Mvalue; Id IDM; title1 'Pearson coefficients of phenological synchronization (rph) (Gomory et al, 2003)'; title3 'Male'; proc print data=RPH (drop=_name_ rename=(IDF=Female)) noobs round; id Female; run; data Rphmean; merge F (where=(_type_='CORR' and upcase(_name_)='FVALUE') rename=(mmean=rphF IDF=ID) drop=fvalue) M (where=(_type_='CORR' and upcase(_name_)='MVALUE') rename=(fmean=rphM IDM=ID) drop=Mvalue); by ID; drop _type_ _name_; title3; proc print data=rpHmean noobs round; var rphF rphM; id id; run; data mfm; set mfm; if IDM='' then do; IDS=IDF;sex='F';end; if IDF='' then do; IDS=IDM;sex='M';end; proc sort; by sex ids; data rphmean;set rphmean; rename ID=IDS; sex='M';rph=rphM;output; sex='F';rph=rphF;output; proc sort; by sex ids; proc sort data=mean; by sex ids; data result1; merge MEAN (where=(_type_=3) keep=_type_ IDS sex %original(PCTM) start end duration ) MFM (keep=IDS sex PO POmin POmax RDmean RDmin RDmax) RPHMEAN (keep=IDS sex RPH); by sex IDS; drop _type_; proc datasets library=work memtype=data nolist; delete m f mf mf2 mfm verify rphmean meanori corr;run;quit; %window grf color=gray irow=13 rows=7 icolumn=7 columns=45 #1 @2 "Do you want to draw a Phenogram (Y/N)" color=black +1 q1 1 attr=underline; %display grf; %put Phenological indices were succesfully calculated; %if %upcase(&q1)=Y %then %SYNCHRO3; %theend: %put SYNCHRO END; options notes merror mprint serror date source; %mend SYNCHRO2; /********************************************************************************************/ /*********************** SYNCHRO3: Drawing Phenograms ****************************************/ /********************************************************************************************/ %macro SYNCHRO3; %let lista=ALL; %let file=c:\MisDoc~1\; %let outgrf=D; %window grf color=gray irow=13 rows=14 icolumn=3 columns=100 #1 @2 "Select genotypes to be included" color=black +1 lista 100 color=blue #3 @2 "Use spaces as delimiter" color=black #4 @2 "Enter genotypes names preceded by ID" color=black #5 @2 "If you want to select all IDs, enter ALL" color=black #6 @2 "Do you want to print to a wmf File or to the Display? (F/D)" color=black +1 outgrf 1 color=blue attr=underline #6 @64 "Folder:" color=black +1 file 50 color=blue attr=underline #8 @10 "Press" +1 "F3" attr=underline +1 "when ready"; %display grf; %if %upcase(&outgrf)=F %then %do; %let device=wmf;%let display=display;%end; %else %do; %let device=win;%let display=nodisplay;%end; %if %upcase(&lista)=ALL %then %let lista=%all; %else %do; data _null_; array temp &lista; call symput('NoID',dim(temp));run; %end; %let listab=%listab; /******************************** DATA preparation ***************************************/ proc sort data=mean; by sex; proc transpose data=mean (where=(_type_=3)) out=graph; by sex; var %original(PCTM); id ids; data graph; set graph; day=substr(_name_,5,3); array IDS &lista; array IDSN &listab; do i=1 to &NoID; IDS(i)=IDS(i)/2; IDSN(i)=-IDS(i); IDS(i)=IDS(i)+i*120; IDSN(i)=IDSN(i)+i*120; end; d=day-1+1; run; /******************* Axis, symbols, patterns, annotate dataset *************************/ %let axisend1=%eval(120*(&NOID+1)); %let series=%eval(2*&Noid); %do k=4 %to &fin; %if %eval(((101-57)/&k)*&k=(101-57)) %then %do; %let step=&k; %let hminor=%eval(&k-1); %goto stepend; %end; %end; %stepend: %do k=1 %to &series; symbol&k color=black i=join value=none; %if %eval((&k/2)*2)=&k %then %str(pattern&k c=black v=solid;); %else %str(pattern&k c=black v=e;); %end; axis1 order= 0 to &axisend1 by 120 label=none value=('' %axis '') offset=(0,0); axis2 order=&inicio to &fin by &step offset=(0,0) length=45 pct label=('Day of the year'); %let listagrf=%listagrf; data anno; length function color style text $8; xsys='2'; ysys='2'; hsys='1'; function='poly';x=&fin-1;y=70; color='black'; style='solid';output; function='polycont';x=&fin;y=70; color=''; style='';output; x=&fin;y=170; output; x=&fin-1;y=170; output; x=&fin-1;y=70; output; function='label';text='100%';x=&fin; y=120;size=3;position='6';output; title1;title3; /****************** Phenograms (male, female and both together) *************************/ proc catalog catalog=work.gseg kill; run;quit; goptions &display DEVICE=&device gsfname=grf gsfmode=replace; filename grf "&file.Female.wmf"; proc gplot data=graph (where=(sex='F')) annotate=anno; title1 'Female Phenogram'; plot (&listagrf)*d / name='f' overlay areas=&series vaxis=axis1 haxis=axis2 frame hminor=&hminor vminor=0; run; filename grf "&file.Male.wmf"; proc gplot data=graph (where=(sex='M')) annotate=anno; title1 'Male Phenogram'; plot (&listagrf)*d / name='m' overlay areas=&series vaxis=axis1 haxis=axis2 frame hminor=&hminor vminor=0; run; goptions display device=&device gsfname=grf gsfmode=replace; filename grf "&file.MandF.wmf"; proc greplay nofs tc=work.gseg template=mf igout=gseg gout=gseg; tdef mf 1/color=white llx=-20 ulx=-20 urx=70 lrx=70 lly=0 uly=100 ury=100 lry=0 2/copy=1 xlatex=50; treplay 1:f 2:m; run;quit; /********************* Overall Phenology synchronity plot ***************************/ proc transpose data=MEAN (where=(_type_=2)) out=graph; var %original(PCTM); id sex; data graph; set graph; dia=substr(_name_,5,3)+0; rename m=Male F=Female; axis3 order=0 to 100 by 10 offset=(0,0) label=('%'); title1 'Overall Phenology Synchronity';title2; legend1 across=2 down=1 label=none value=('Pollen-shedding rate' 'Female receptivity rate'); symbol1 i=splines v=dot color=black h=0.7; symbol2 i=splines v=square color=black h=0.7; filename grf "&file.synchro.wmf"; proc gplot data=graph; plot (male female)*dia/ overlay frame legend=legend1 vaxis=axis3 haxis=axis2 vminor=0 hminor=&hminor; run; /**************** Frequency Histograms of the synchronity indices *******************/ pattern1 v=solid c=black; axis4 label=('N') minor=none length=70 pct; axis5 length=70 pct; axis6 value=('0.1' '0.2' '0.3' '0.4' '0.5' '0.6' '0.7' '0.8' '0.9' '1.0') length=70 pct; %let sex=M;%let titulo=Male; %gchart: proc gchart data=result1 (where=(sex="&sex")); title1 "Frecuency histograms for RPH indices (&titulo)"; filename grf "&file.RPH&sex..wmf"; vbar rph /midpoints=-0.9 -0.7 -0.5 -0.3 -0.1 0.1 0.3 0.5 0.7 0.9 axis=axis4 maxis=axis5 frame width=5 name="RPH&sex"; label rph="rph (&sex)"; proc gchart data=result1 (where=(sex="&sex")); title1 "Frecuency histograms for PO indices (&titulo)"; filename grf "&file.PO&sex..wmf"; vbar PO /midpoints=0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 axis=axis4 maxis=axis6 frame width=5 name="PO&sex"; label PO="PO (&sex)"; proc gchart data=result1 (where=(sex="&sex")); title1 "Frecuency histograms for RD indices (&titulo)"; filename grf "&file.RD&sex..wmf"; vbar RDMEAN /midpoints=4 to 24 by 2 axis=axis4 maxis=axis5 frame width=5 name="RD&sex"; label RDMEAN="RD (&sex)"; %if &sex=M %then %do; %let sex=F; %let titulo=Female; %goto gchart;%end; proc greplay nofs tc=work.gseg template=mf igout=gseg gout=gseg; filename grf "&file.allfreqs.wmf"; tdef mf 1/color=white llx=0 ulx=0 urx=33 lrx=33 lly=50 uly=100 ury=100 lry=50 2/copy=1 xlatex=33 3/copy=1 xlatex=66 4/copy=1 xlatey=-50 5/copy=2 xlatey=-50 6/copy=3 xlatey=-50; treplay 1:rphm 2:pom 3:rdm 4:rphf 5:pof 6:rdf; run;quit; proc datasets library=work memtype=data nolist; delete graph anno mean;run;quit; %put NOTES: All graphs were created succesfully; %MEND SYNCHRO3; /********************* Macros definition *************************/ %macro listab; %do j=1 %to &NoID; %scan(&lista,&j)b %end; %mend; %macro listagrf; %do j=1 %to &NoID; %scan(&listab,&j) %scan(&lista,&j) %end; %mend; %macro axis; %do j=1 %to &NoID; "%substr(%scan(&lista,&j),3,%eval(%length(%scan(&lista,&j))-2))" %end; %mend; %macro all; %do j=1 %to &NoID; &&ID&j %end; %mend; %macro rename(sex); %do i=1 %to &n; PCTM&&day&i=&sex&&day&i %end; %mend; %macro original(name); %do i=1 %to &n; &name&&day&i %end; %mend; /********************* Program execution *************************/ %SYNCHRO1()