scripts/50a_limma-fit-swp.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=60)
@

<<>>=
testvar <- "Vitvi01g01391"
@

\clearpage
\subsection{Overall Model}

Simple model with water potential for all data:
<<>>=
design1 <- with(pd, model.matrix(~ swp))
#design1[,2] <- (-design1[,2])
head(design1)
tail(design1)
@
Linear model
<<>>=
fit1 <- lmFit(exprs, design1)
fit1[testvar,]$coefficients
fit1 <- eBayes(fit1)
fit1[testvar,]$coefficients
head(fit1[1:6,]$coefficients)
fit <- fit1
fit[testvar,]$coef
@

Models for each variety

<<>>=
coefs <- fit1$coefficients
plot(coefs[,1],coefs[,2])
@

<<>>=
dim(coefs[abs(coefs[,2])>2,])

@

Top table

<<>>=
which <- 2
  tt <- topTable(fit1,coef=which, number=Inf)
  head(tt)
@

\clearpage
\subsection{Models for variety}

Model with water potential by variety. 
<<>>=
design2 <- with(pd, model.matrix(~ swp*variety))
head(design2)
tail(design2)

@
Linear model
<<>>=
fit2 <- lmFit(exprs, design2)
fit2[testvar,]$coefficients
@

<<>>=
varname <- testvar
par(mfrow=c(1,3))
p <- plot.ewp(varname,model="")
plot(pd$swp,exprs[testvar,],  xlim=c(0,-1.6),
   ,pch=16, col=6-as.numeric(pd$variety)*2,cex=pd$day/2)
c(p$C$coef[2],p$F$coef[2])
fit2$coef[testvar,]
p <- plot.ewp(varname,model="", exprs=exprs, pd=pd)
@

<<>>=
#
fit2 <- eBayes(fit2)
fit2[testvar,]$coefficients
head(fit2[1:6,]$coefficients)
@

For comparison, here is the beginning of coefficients for overall model.
<<>>=
head(fit1[1:6,]$coefficients)
@
<<>>=
fita[2]
fit1[testvar,]$coefficients[2]
fitc[2]
fit2[testvar,]$coefficients[2]
fitf[2]
fit2[testvar,]$coefficients[2]+fit2[testvar,]$coefficients[4]
@

Plots of coefficients

<<>>=
oldpar <- par(mfrow=c(2,2), mar=c(4,4,1,1))
coefs1 <- fit1$coefficients
coefs2 <- fit2$coefficients
plot(coefs1[,2],coefs2[,2], xlab="Slope All",ylab="Slope C")
abline(c(0,1),col=2)
cf24 <- coefs2[,2] + coefs2[,4]
plot(coefs2[,2], cf24, xlab="Slope C",ylab="Slope F")
abline(c(0,1),col=2)
plot(coefs1[,2], cf24, xlab="Slope All",ylab="Slope F")
abline(c(0,1),col=2)
par(oldpar)
@

<<warnings=FASLE,echo=FALSE>>=
source("../doc/Graph3d.r")
persp.3d(coefs2[,2],cf24, coefs1[,2],
phi=45, zlab="All",xlab="C",ylab="F")
@

\subsection{Mean expressions}

<<>>=
my.predict <- function(fit, x=1){
cf <- fit$coefficients
mc <- cf%*%c(1,x,0,0)
mf <- cf%*%c(1,x,1,x)
return(data.frame(mc=mc,mf=mf))
}
fit2$coefficients[1:2,]
my.predict(fit2)[1:2,]
my.predict(fit2)[testvar,]
ma <- fit1$coefficients%*%c(1,mean(swp))
macf <- data.frame(ma=ma,my.predict(fit2, x=mean(swp)))
head(macf)
@

<<>>=
macf[testvar,]
@

\subsubsection{ Plots of variety coefficients}

I will try to set the variety slopes compared to slope of complete data


For comparison, here is the beginning of coefficients for overall model.

<<>>=
head(fit1[1:6,]$coefficients)

@

Extract slopes

<<>>=
A <- fit1$coefficients[,2]
C <- fit2$coefficients[,2]
F <- fit2$coefficients[,2]+fit2$coefficients[,4]
ma <- macf[,"ma"]
mc <- macf[,"mc"]
mf <- macf[,"mf"]
coefs <- data.frame(A,C,F,ma,mc,mf)
head(coefs)
ca <- C-A
fa <- F-A
fc <- F-C
coefs <- data.frame(mc,mf,C,F)
@

<<>>=
par(mfrow=c(2,3))
hist(A,xlim=range(c(A,C,F)))
hist(C,xlim=range(c(A,C,F)))
hist(F,xlim=range(c(A,C,F)))
hist(ma)
hist(mc)
hist(mf)
@



<<>>=
delta <- 2
large_abs_A <- abs(A)> delta
small_fc <- abs(fc) < delta/5
small_ca <- abs(ca) < delta/5
oldpar <- par(mfrow=c(2,3), mar=c(4,4,1,1))
plot(A,ca, xlab="Slope All",ylab="Slope C - slope All")
points(A[large_abs_A],ca[large_abs_A],col=2,pch=16)
abline(h=0,v=0,col=2)
plot(A,fa, xlab="Slope All",ylab="Slope F - slope All")
points(A[large_abs_A],fa[large_abs_A],col=2,pch=16)
abline(h=0,v=0,col=2)
plot(ca,fa, xlab="Slope C - slope All",ylab="Slope F - slope All")
points(ca[large_abs_A],fa[large_abs_A],col=2,pch=16)
abline(h=0,v=0,col=2)
#
plot(A,fc, xlab="slope All",ylab="Slope F - slope C")
points(A[large_abs_A],fc[large_abs_A],col=2,pch=16)
points(A[large_abs_A&small_fc],fc[large_abs_A&small_fc],col=4,pch=16)
abline(h=0,v=0,col=2)
#
plot(C,fc, xlab="slope C",ylab="Slope F - slope C")
points(C[large_abs_A],fc[large_abs_A],col=2,pch=16)
points(C[large_abs_A&small_fc],fc[large_abs_A&small_fc],col=4,pch=16)
abline(h=0,v=0,col=2)
#
plot(C,F, xlab="slope C",ylab="Slope F ")
points(C[large_abs_A],F[large_abs_A],col=2,pch=16)
points(C[large_abs_A&small_fc],F[large_abs_A&small_fc],col=4,pch=16)
type1 <- large_abs_A&small_fc&small_ca
points(
   C[type1],
   F[type1],
   col=3,pch=16)

abline(h=0,v=0,col=2)
abline(c(0,1),col=3)
par(oldpar)
@

\subsubsection{Situation types and Type selection criteria}

Boundaries:
\begin{description}
  \item[min] for at least ... ( > )
  \item[max] for at most ... ( < )
  \item[delta] difference at most ... ( < )
\end{description}


<<>>=
#
delta <- 2
#
deltam <- 2
minA <- 2       # slope A larger
maxC <- 0.5     # slope C smaller
maxF <- maxC    # slope F smaller
maxfc <- 0.5    # difference F - C smaller
minfc <- 1      # differende F - C larger
minC <- minA    # slope C larger
minF <- minA    # slope F larger
maxca <- maxfc  # difference C - A smaller
maxfa <- maxfc  # difference F - A smaller
#
large_abs_A <- abs(A)> minA
large_abs_C <- abs(C)> minC
large_abs_F <- abs(F)> minF
#
small_fc <- abs(fc) < maxfc
small_ca <- abs(ca) < maxca
small_fa <- abs(fa) < maxfa
plot(C,F, xlab="slope C",ylab="Slope F ")
#points(C[large_abs_A],F[large_abs_A],col=2,pch=16)
#

type10 <- large_abs_A & small_fc & small_ca
type1 <- (large_abs_C|large_abs_F) & small_fc & (abs(mc-mf)<deltam)
table(type10,type1)
points(
   C[type1],
   F[type1],
   col=2, pch=16)
#
type2 <- (large_abs_C & large_abs_F) & (abs(fc) > minfc)
type2 <- type2 & ( (C*F >0) )
points(
   C[type2],
   F[type2],
   col=3,pch=16)
#
type30 <- ((abs(C)<maxC)|(abs(F)<maxF)) & (abs(fc)> minfc)
type3 <- ((abs(C)<maxC)|(abs(F)<maxF)) & (abs(fc)> minfc)&(abs(mc-mf)>deltam)
table(type3,type30)
points(
   C[type3],
   F[type3],
   col=4,pch=16)
#
type4 <- (large_abs_C|large_abs_F)&(abs(mc-mf)>deltam)&(C*F>0)
points(
   C[type4],
   F[type4],
   col=5,pch=16)
#
type5 <- ((abs(C)<maxC)|(abs(F)<maxF)) & (abs(fc)> minfc*2)&(abs(mc-mf)<deltam)
points(
   C[type5],
   F[type5],
   col=6,pch=16)

table(type3,type5)
points(
   C[type3],
   F[type3],
   col=4,pch=16)
# type 3 is a subset of type5, type3 & large mean difference

#
abline(h=0,v=0,col=2)
abline(c(0,1),col=3)
legend("bottomright",pch=rep(16,5),col=2:6,legend=paste("Type",1:5),bty="n")
types <- data.frame(type1,type2,type3,type4,type5)

@

<<>>=
head(types)
plot.types <- function(x=C, y=F, type=types,title,...){
xlab <- deparse(substitute(x))
ylab <- deparse(substitute(y))

plot(x,y, xlab=xlab,ylab=ylab)
m <- ncol(type)
sek <- c(1,2,4,5,3)
for( i in sek){
points(
   x[types[,i]],
   y[types[,i]],
   col=i+1,pch=16)
   }
   abline(h=0,v=0,col=2)
abline(c(0,1),col=3)
if(!missing(title)) title(title)
legend("bottomright",pch=rep(16,m),col=1+(1:m),legend=paste("Type",1:m),bty="n")
}
plot.types()

@


<<>>=
x <- coefs[type1,]
ylim <- range(x)
my.parallelplot <- function(x, ylim = range(x), title,  ...){
if(missing(title)) title <- deparse(substitute(x))
m <- ncol(x)
plot(0,0, type = "n"
    , xlim = c(1,m)
    , ylim = ylim
    , axes = FALSE
    , ann=FALSE
    )
title(title)
mtext(paste("n profiles:", nrow(x)),cex=0.75)
axis(2)
axis(1, at=1:m,labels=gsub("\\.\\.\\.","-",colnames(x)),cex=0.5)  
abline(v=1:m,col=8)
x$col <- 1:nrow(x)    
apply(x,1,function(y) lines(1:m,y[-(m+1)],col=y[m+1]))
invisible(nrow(x))
}
if(interactive()) my.parallelplot(coefs[type1,])
@

<<>>=
par(mfrow=c(2,3))
tpar <- data.frame(C,F,F-C,mc,mf,mf-mc)
ylim <- range(tpar)
my.parallelplot(tpar[type1,], ylim=ylim)
my.parallelplot(tpar[type2,], ylim=ylim)
my.parallelplot(tpar[type3,], ylim=ylim)
my.parallelplot(tpar[type4,], ylim=ylim)
my.parallelplot(tpar[type5,], ylim=ylim)
plot.types()

@

<<>>=
par(mfrow=c(2,3))
tpar <- data.frame(C,F,F-C,mc,mf,mf-mc)
abstpar <- abs(tpar)
ylim <- range(abstpar)
my.parallelplot(abstpar[type1,], ylim=ylim)
my.parallelplot(abstpar[type2,], ylim=ylim)
my.parallelplot(abstpar[type3,], ylim=ylim)
my.parallelplot(abstpar[type4,], ylim=ylim)
my.parallelplot(abstpar[type5,], ylim=ylim)
plot.types(abs(C),abs(F),title="Absolute slopes")

@

<<>>=
plot.types(C,F,title="Slopes")
plot.types(mc,mf,title="Mean expression")
@

<<>>=
varname <- testvar
plot.ewp(varname,model="", exprs=exprs, pd=pd)
@

<<>>=
apply(data.frame(type1,type2,type3,type4, type5),2,sum)
@


<<>>=
par(mfrow=c(2,3))
filt <- type1
print(sum(filt))
for(i in 1:min(sum(filt),6)){
varname <- rownames(exprs)[filt][i]
p <- plot.ewp(varname,model="", exprs=exprs, pd=pd, ylim=range(exprs)
     , cex= 1.5
     , pch = c(16,16)
     )
mtext(paste(round(coefs[varname,],1),collapse=" | "),cex=0.75)
mtext(paste(round(c(p$C$coefficients[2],p$F$coefficients[2]),1)
,collapse=" | "),1,cex=0.75)
}
@

<<>>=
par(mfrow=c(2,3))
filt <- type2
print(sum(filt))
for(i in 1:min(sum(filt),6)){
varname <- rownames(exprs)[filt][i]
p <- plot.ewp(varname,model="", exprs=exprs, pd=pd, ylim=range(exprs)
     , cex= 1.5
     , pch = c(16,16)
     )
mtext(paste(round(coefs[varname,],1),collapse=" | "),cex=0.75)
mtext(paste(round(c(p$C$coefficients[2],p$F$coefficients[2]),1)
,collapse=" | "),1,cex=0.75)
}
@

<<>>=
par(mfrow=c(2,3))
filt <- type3
print(sum(filt))
for(i in 1:min(sum(filt),6)){
varname <- rownames(exprs)[filt][i]
p <- plot.ewp(varname,model="", exprs=exprs, pd=pd, ylim=range(exprs)
     , cex= 1.5
     , pch = c(16,16)
     )
mtext(paste(round(coefs[varname,],1),collapse=" | "),cex=0.75)
mtext(paste(round(c(p$C$coefficients[2],p$F$coefficients[2]),1)
,collapse=" | "),1,cex=0.75)
}
@

<<>>=
par(mfrow=c(2,3))
filt <- type4
print(sum(filt))
for(i in 1:min(sum(filt),6)){
varname <- rownames(exprs)[filt][i]
p <- plot.ewp(varname,model="", exprs=exprs, pd=pd, ylim=range(exprs)
     , cex= 1.5
     , pch = c(16,16)
     )
mtext(paste(round(coefs[varname,],1),collapse=" | "),cex=0.75)
mtext(paste(round(c(p$C$coefficients[2],p$F$coefficients[2]),1)
,collapse=" | "),1,cex=0.75)
}
@

<<>>=
par(mfrow=c(2,3))
filt <- type5
print(sum(filt))
for(i in 1:min(sum(filt),6)){
varname <- rownames(exprs)[filt][i]
p <- plot.ewp(varname,model="", exprs=exprs, pd=pd, ylim=range(exprs)
     , cex= 1.5
     , pch = c(16,16)
     )
mtext(paste(round(coefs[varname,],1),collapse=" | "),cex=0.75)
mtext(paste(round(c(p$C$coefficients[2],p$F$coefficients[2]),1)
,collapse=" | "),1,cex=0.75)
}
@

<<>>=
which <- 2
  tt <- topTable(fit,coef=which, number=Inf)
  head(tt)
@

\clearpage
\subsubsection{Export output files}

Export fit in a form:

<<>>=
fit <- fit2
fitfn <- "lm-fit-swp*variety-statistics.txt"
write.fit(head(fit),file="", adjust="BH")
write.fit(fit,file=file.path(.oroot,fitfn), adjust="BH")
@

Export top table

<<results='asis'>>=
ttfn <- "TopTable.txt"
tt <- topTable(fit, number=Inf)
tt$type <- ""
types <- types[rownames(tt),]
all(rownames(tt)==rownames(types))
#ind <- sapply(rownames(tt),FUN=function(x) match(x,rownames(types)))
for(i in 1:ncol(types)) tt$type[types[,i]] <- colnames(types)[i]
ind <- sapply(rownames(tt),FUN=function(x) match(x,fdata$geneID))
ttext <- data.frame(tt,fdata[ind,2:5] )
my.write.table(ttext,file=file.path(.oroot,ttfn))
@


<<>>=
topn <- .topn
@

\clearpage
\subsection{Genes in types}

Plots of situation for genes in determined types

Figures caption:


<<echo=FALSE>>=
fig.cap <- "Mean expression (E) and mean water potential (WP). Upper panels show WP vs E and regression lines for three models. Left: all C (blue) and all F (blue). Middle: separate models for combinations of treatment and variety, same slope (+). Right: separate for combinations, different slopes (interaction, *). Black line is regression of all data, regardles of treatment and variety. Lower panels have swaped variables: E vs WP. Symbol sizes are related to time, large is later."
@

\Sexpr{fig.cap}

\clearpage

<<>>=
out <- ""
(cname <- "")
i <- 1
for(i in 1:ncol(types)){
out <- paste(out,"
\\clearpage
\\subsection{",colnames(types)[i],"}")
varnames <- rownames(tt)[types[,i]][1:5]
varname <- varnames[1]
for(varname in varnames)
 out <- paste(out,knit_child(file.path("../doc/","50aa_one-gene-plots.Rnw"),quiet=TRUE))
}
@
\Sexpr{out}



<<>>=

@