## Load the package library("FAMT") ## Load the datasets data("expression") # gene expressions data("covariates") # explanatory variables data("annotations") # variables for interpretation ## Create a FAMTdata object containing the expression, covariates ## and annotations datasets ## the array identification is given in the second column of covariates chicken = as.FAMTdata(expression, covariates, annotations, idcovar=2) ## Summary of the FAMTdata summaryFAMT(chicken) ## Classical method: test the significance of the relationship between ## each gene expression and the abdominal fatness variable (6th column of covariates) ## taking into account the effect of the dam (3rd column of covariates) chicken.raw = modelFAMT(chicken, x=c(3,6), test=6, nbf=0) hist(chicken.raw$pval, main="Histogram of p-values", xlab="p-values") ## Summary of the FAMTmodel summaryFAMT(chicken.raw, alpha=0.05, info=c("ID","Name")) ## FAMT complete multiple testing procedure ## when the optimal number of factors is unknown ## (the function includes an algorithm to determine it) modelfinal=modelFAMT(chicken, x=c(3,6), test=6) modelfinal$nbf ## Summary of the FAMTmodel summaryFAMT(modelfinal, alpha=seq(0, 0.3, 0.05)) ## Plot the histogram of both the unadjusted and factors adjusted p-values par(mfrow=c(1,2)) hist(modelfinal$pval, main="Histogram of p-values",xlab="Unadjusted p-values") hist(modelfinal$adjpval, main="Histogram of adjusted p-values",xlab="Adjusted p-values") ## Determination of the optimal number of factors nbfactors(chicken, x=c(3,6), test=6, diagnostic.plot=TRUE) ## Estimation of the proportion of true null hypotheses based on the p-values density's estimation pi0FAMT(modelfinal, method="density", diagnostic.plot=TRUE) ## FAMT factors description chicken.defacto = defacto(modelfinal, axes=1:2, select.covar=4:5,select.annot=3:6, cex=0.6) ## Matrix of p-values chicken.defacto$covariates chicken.defacto$annotations ########################### # Second illustrative example ########################### # Data Expr=read.table("Expr.txt",sep="\t",header=TRUE) Covar=read.table("covar.txt",sep="\t",header=TRUE) # Create the FAMTdata Poulets=as.FAMTdata(Expr,Covar,idcovar=1) # Summary of the FAMTdata summaryFAMT(Poulets) # FAMT complete multiple testing procedure # with 5 factors model=modelFAMT(Poulets,x=2,nbf=5) # Extraction of the number of rejected genes rejections =summaryFAMT(model, alpha=seq(0.001,0.01,0.001))$nbreject # Plot the number of significant genes for various FDR control levels plot(rejections[,1],rejections[,3], type="l", lty=1, ylab="Number of significant genes", xlab="False Discovery Rate control level", ylim=c(0,7000)) lines(rejections[,1],rejections[,2], type="l", lty=2) legend("topleft",c("FAMT","Classical method"), lty=c(1,2))