#Differential abundance analysis: All time points setwd("C:\\Princeton R script\\Publication timepoints -6h") library(limma) targets <- readTargets("Targets.txt") spottypes <- readSpotTypes("Spot_Types.txt") RG <- read.maimages(targets$FileName,source="genepix",wt.fun=wtflags(0)) RG$genes<-readGAL("august gal2.gal") RG$genes$Status<-controlStatus(spottypes, RG) RG$printer<- list(ngrid.r=8, ngrid.c=4, nspot.r=11, nspot.c=24, ndups=1, npins=32, start="topleft") library(convert) mraw<-as(RG, "marrayRaw") # Start of loop for diagnostic images for (i in 1:20){ image(mraw[,i],xvar="maGb", main=RG$targets[[1]][[i]]) x11() image(mraw[,i],xvar="maRb", main=RG$targets[[1]][[i]]) } # Background correction RGb<-backgroundCorrect(RG, method="normexp", offset=50) library(convert) mraw <- as(RGb,"marrayRaw") # Loop for print-tip boxplots for (i in 1:20){ x11() title = paste("Print-tip boxplots for array", RG$targets[[1]][[i]], ": prenormalization", separator="") boxplot(mraw[,i],xvar="maPrintTip", yvar="maM", main=title) } #Array boxplot to check for differences between arrays x11() boxplot(mraw, yvar="maM", main="Array boxplots: pre-normalization") # Loop for MA plots for (i in 1:20){ title = paste("MA plot Slide ", RG$targets[[1]][[i]], separator=" ") x11() plotMA(RG,array=i, main=title) } # NormalizeWithinArrays MA<-normalizeWithinArrays(RGb, method="robustspline") library(convert) mraw.norm<-as(MA, "marrayNorm") # Loop for post-normalization plots for (i in 1:20){ x11() title = paste("Print-tip boxplots for array", RG$targets[[1]][[i]], ": post robustspline normalization", separator=" ") boxplot(mraw.norm[,i],xvar="maPrintTip", yvar="maM", main=title) } #Array boxplot to check for differences between arrays x11() boxplot(mraw.norm, yvar="maM", main="Array boxplots time: post robustspline within normalization") MAbet<-normalizeBetweenArrays(MA,method="Gquantile") x11() plotDensities(MAbet) # Loop for MA plots for (i in 1:20){ x11() title = RG$targets[[1]][[i]] plotMA(MAbet,array=i, main=) } # Check between slide normalization library(convert) mraw.norm.bet <- as(MAbet,"marrayNorm") x11() boxplot(mraw.norm.bet, yvar="maA", main="Array intensity boxplots: post-Gquantile between normalization") x11() boxplot(mraw.norm.bet, yvar="maM", main="Array boxplots: post-Gquantile between normalization") #Linear model differential abundance analysis compared to relative time zero = UTt1 design<-modelMatrix(targets,ref="Reference") design fit<-lmFit(MAbet,design) contrast.matrix<-makeContrasts(Tt1-UTt1, levels=design) #Change for different comparisons e.g. Tt2 and Tt3 vs UTt1 contrast.matrix fit2<-contrasts.fit(fit, contrast.matrix) fit2<-eBayes(fit2) # extract all genes for T compared to UTt1 and apply FDR # adjustment to p-values to adjust for multiple testing TvUTt1<-topTable(fit2, adjust="fdr", number=8448) TvUTt1 write.table(TvUTt1,file="C:/Princeton R script/Publication timepoints -6h/top_table.txt",sep="\t",row.names=FALSE) x11() volcanoplot(fit2,names=fit$genes$Name)