R/ASCA2f.R

In SoniaTC/NOISeq: Exploratory analysis and differential expression for RNA-seq data

ASCA.2f<-function(X = X,Designa = Designa, Designb = Designb, Designc=NULL, Fac=c(1,2,2,2),Join = TRUE,Interaction=TRUE)

{
#--------------------------------------------------------------------------------------
#  Dimensions of the matrices: 
#  X (p x n) contains expression values of n genes (in columns) and p conditions (in rows)
#  Designa (p x I) contains 0's and 1's for the TIME-POINTS in the experiment
#  Designb (p x J) EXPERIMENTAL GROUP  
#  Designres (p x H) INDIVIDUALS

#  Join = TRUE if the analyses of the model b and ab is studied jointly
#  Interaction = TRUE to consider interaction "ab" in the separated model

n<-ncol(X)
p<-nrow(X)
I<-ncol(Designa)
J<-ncol(Designb)

Faca=Fac[1] # number components Model a (time)
Facb=Fac[2] # number components Model b  (second factor)
Facab=Fac[3] # number components Model ab (interaction) 
Facres=Fac[4] # number components Residues

#----------------------- Calculate Overall Mean --------------------------------------

offset<-apply(X,2,mean)
Xoff<-X-(cbind(matrix(1,nrow=p,ncol=1))%*%rbind(offset))


#-----------------------  PART I: Submodel a (TIME) -----------------------------------

	Model.a<-ASCAfun1(Xoff,Designa,Faca)
	Xres<-Xoff-Model.a$X

#-------------------------- PART II: Submodel b and ab -------------------------------

if(!Join) 
{
	Model.b<-ASCAfun1(Xoff,Designb,Facb)
	if (Interaction) {
		Model.ab<-ASCAfun2(Xoff,Designa,Designb,Facab)
	}
}

if(Join)
{
	Model.bab<-ASCAfun12(Xoff,Designa,Designb,Facab)
}

# ------------------------Collecting models ------------------------------------------

models <- ls(pattern="Model")
output <- vector(mode="list")
Xres <- Xoff
for (i in 1: length(models)) {
	mymodel <- get(models[i], envir=environment())
	output <- c(output, list(mymodel))
	Xres <- Xres - mymodel$X
	rm(mymodel)
	gc()
	
}
names(output) <- models

#------------------------- PART III: Submodel res -----------------------------------

Model.res<-ASCAfun.res(Xres,Facres)

Model.res<-list(Model.res)
names(Model.res)<-c("Model.res")
output<-c(output,Model.res)


#------------------------- Add Input data to the Output ----------------------------

Input<-list(X, Designa, Designb, Designc, Fac, Join,Interaction,Xoff)
names(Input)<-c("X", "Designa", "Designb", "Designc", "Fac", "Join","Interaction","Xoff")
Input<-list(Input)
names(Input)<-"Input"
output<-c(output,Input)

output 
}