scripts/10c_Interaction-plots.Rnw

Version 1

% -*- TeX:Rnw:UTF-8 -*-
% ----------------------------------------------------------------
% .R knitr file  ************************************************
% ----------------------------------------------------------------
%%
<<echo=FALSE>>=
###############################################
##                                           ##
## (c) Andrej Blejec (andrej.blejec@nib.si)  ##
##                                           ##
###############################################
@
<<echo=FALSE,results='hide'>>=
options(width=70)
@
Input objects:
<<>>=
# y  gene expression data frame
# pd phenodata
# minsd minimal standard deviation to include the gene
@

<<>>=
y[1:5,1:5]
pd[1:5,1:5]
minsd
@

<<eval=FALSE,echo=FALSE>>=
# wide to long format
library(tidyr)
#y <- t18
sel <- which(apply(y,1,sd)>4)
y <- t(y[sel,])
y <- data.frame(y,id=rownames(y))
y <- gather(y, "id", "exprs")
y <- data.frame(y, pdata18[,c(1,4:8)])
y$nday <- as.numeric(levels(y$day)[as.numeric(y$day)])
head(y)
@


<<>>=
#minsd <- 5
minsd
sel <- which(apply(y,1,sd)> minsd)
y <- t(y[sel,])
head(y)
@

Plots are presented for genes with large standard deviation ( > \Sexpr{minsd} ).

\clearpage
\subsubsection{Interaction plots treatment x variety within days}

<<>>=
my.interaction.plot <- function(F1,F2,y,...){
    xlab <- deparse(substitute(F1))
    interaction.plot(F1, F2, y,
     col = 1:length(levels(F2)), 
     xlab = xlab,
     ...
     )
    points(as.numeric(F1)+(as.numeric(F2)-1.5)/15,y, col=as.numeric(F2))
    }
@
\clearpage
<<>>=
days <- levels(pd$day)
par(mfrow=c(1,length(days)),mar=c(5,4,4,3))
for(i in 1:ncol(y)){
    varname <- colnames(y)[i]
    ylim <- range(y[,i])
   for(dy in days){
    sel <- pd$day %in% dy
    with(pd[sel,],
       my.interaction.plot(variety,treat, y[sel,varname], ylim=ylim,type="b", pch=rep(16,2))
       )
    title(paste(varname, "\n", dy))
    }
}
@
\clearpage
\subsubsection{Interaction plots treatment x day within varieties}
<<>>=
days <- levels(pd$day)
par(mfrow=c(1,length(levels(pd$variety))),mar=c(5,4,4,3))
for(i in 1:ncol(y)){
    varname <- colnames(y)[i]
    ylim <- range(y[,i])
   for(var in levels(pd$variety)){
    sel <- pd$variety %in% var
    with(pd[sel,],
       my.interaction.plot(day,treat, y[sel,varname], 
       ylim=ylim,type="b", pch=rep(16,2),
       ))
    title(paste(varname,"\n", var))
    }
}
@