Difference between revisions of "Nunz-mik-12-12-05.R"

From Organic Design wiki
(additions)
m
 
(3 intermediate revisions by the same user not shown)
Line 1: Line 1:
−
auxillaryDir <- "/Volumes/HD2/Data/Nutrigenomics/MultipleScans/Auxillary"
+
# {{R}}
−
dataDir      <- "/Volumes/HD2/Data/Nutrigenomics/MultipleScans/GPR"
 
−
 
 
 
######################################################################################
 
######################################################################################
−
# Load required R libraries.
+
# Load required R libraries.
 
######################################################################################
 
######################################################################################
 
library(limma)
 
library(limma)
 +
datadir <- "/Volumes/HD2/Data/Nutrigenomics/MultipleScans/GPR"
 +
# Read targets (Miks file)
 +
#targets<-readTargets()
  
−
# Read targets
 
−
targets<-readTargets(path=auxillaryDir)
 
−
targets$Name <- paste(targets$Filetype, targets$Scan, sep="")
 
 
# Define weight function to read in the recorded flags in each GPR file
 
# Define weight function to read in the recorded flags in each GPR file
 
wt.obs<-function(gpr){
 
wt.obs<-function(gpr){
Line 16: Line 14:
 
}
 
}
  
−
files <- dir(dataDir, pattern="gpr")
 
 
# Read in GPR files
 
# Read in GPR files
−
RG<-read.maimages(files, names=substr(files,20,29), source="genepix",path=dataDir,columns=list(Rf = "F635 Median", Gf = "F532 Median",Rb = "B635 Median", Gb = "B532 Median")) # names=targets$Name,
+
# Currently not picking up any weight information
 +
# names=targets$Name,
 +
RG<-read.maimages(files=dir(datadir),source="genepix",path=datadir,columns=list(Rf = "F635 Median", Gf = "F532 Median",Rb = "B635 Median", Gb = "B532 Median"))
 
RG$printer <- getLayout(RG$genes)
 
RG$printer <- getLayout(RG$genes)
  
Line 24: Line 23:
 
# RG variables (otherwise weights are used in normalization).
 
# RG variables (otherwise weights are used in normalization).
 
RG.weights<-RG$weights
 
RG.weights<-RG$weights
−
 
+
−
#mwd weights were not loaded in
 
 
RG$weights=NULL
 
RG$weights=NULL
 
   
 
   
Line 36: Line 34:
 
library(limma)
 
library(limma)
  
−
#mwd There are none
 
 
RG$R[RG$R==0]<-0.5
 
RG$R[RG$R==0]<-0.5
 
RG$G[RG$G==0]<-0.5
 
RG$G[RG$G==0]<-0.5
Line 67: Line 64:
 
#########################################################################################
 
#########################################################################################
 
load('RGw.data')
 
load('RGw.data')
−
for(i in 1:(dim(RG$R)[2]/8)){
+
for(i in 1:(dim(RG$R)[2]/8)){ # BAD CODING ... Should be 1:(ncol(RG$R)/9)
−
  #mwd need to make the folders...
+
   fname<-paste("JPG/BGFG/BG_FG_",i,".png",sep="") ... pngs better
−
   fname<-paste("JPG/BGFG/BG_FG_",i,".jpg",sep="")
+
#  jpeg(fname, height = 900, width=1200, quality=80, bg="white")
−
  jpeg(fname, height = 900, width=1200, quality=80, bg="white")
+
  png(fname, height = 900, width=1200, pointsize=15) 
 
   par(mfrow=c(8,4))
 
   par(mfrow=c(8,4))
 
   for(j in 1:8){
 
   for(j in 1:8){
Line 113: Line 110:
 
# Boxplots for NBC data (Pre-norm and post-norm)
 
# Boxplots for NBC data (Pre-norm and post-norm)
 
#####################################################################################
 
#####################################################################################
−
lo<-grep("l",targets$Name)
+
lo<-grep("low", colnames(RG))
−
md<-grep("a",targets$Name)
+
md<-grep("average", colnames(RG))
−
hi<-grep("h",targets$Name)
+
hi<-grep("high", colnames(RG))
 
    
 
    
 
jpeg("JPG/Boxplots/boxplot-full-nbc-pre.jpg",height=900,width=1200,quality=80, bg="white")
 
jpeg("JPG/Boxplots/boxplot-full-nbc-pre.jpg",height=900,width=1200,quality=80, bg="white")
Line 188: Line 185:
 
load("normw.data")
 
load("normw.data")
 
#load("norm2w.data")
 
#load("norm2w.data")
−
sat.r<-RG$R==(2^16-1)
+
sat.r<-RG$R==(2^16-1) # Alot more saturated in R channel
 
sat.g<-RG$G==(2^16-1)
 
sat.g<-RG$G==(2^16-1)
 
sat<-sat.r | sat.g
 
sat<-sat.r | sat.g
  
−
# mwd minor name changes
+
#lo<-grep("l",targets$Name) Already done above
−
lo<-grep("L",targets$Name)
+
#md<-grep("a",targets$Name)
−
md<-grep("A",targets$Name)
+
#hi<-grep("h",targets$Name)
−
hi<-grep("H",targets$Name)
 
 
tsat<-sat[,lo] & sat[,md] & sat[,hi]
 
tsat<-sat[,lo] & sat[,md] & sat[,hi]
  
Line 234: Line 230:
  
 
targ<-targets[md,]
 
targ<-targets[md,]
−
design<-matrix(c(0, 0,1, 0, 0,0,-1,0,0, 0,0,0,  # C.25  
+
design<-matrix(c(0, 0,1, 0, 0,0,-1,0,0, 0,0,0,  # C.25
 
                 0, 0,0,-1, 0,0, 0,1,0, 0,0,0,  # M.I10.25
 
                 0, 0,0,-1, 0,0, 0,1,0, 0,0,0,  # M.I10.25
 
                 1, 0,0, 0,-1,0, 0,0,0, 0,0,0,  # C.12
 
                 1, 0,0, 0,-1,0, 0,0,0, 0,0,0,  # C.12
Line 241: Line 237:
 
                 ),12,5)
 
                 ),12,5)
 
colnames(design)<-c("C.25","M.I10.25","C.12","M.I0.12","M.I5.22")
 
colnames(design)<-c("C.25","M.I10.25","C.12","M.I0.12","M.I5.22")
−
rownames(design) <- c("P1","P10","P11","P12","P2","P3","P4","P5","P7","P8","P9","Ref")
+
rownames(design)<- substr(colnames(RG)[md], 20, 27) # targets$Name[md]
−
design
 
−
#rownames(design)<-targets$Name[md]
 
  
−
MM<-MCb.n[,c(1:5,7:11)]
+
MM<-MCb.n[,c(1:5,7:11)] # Pool 3 and ref-ref removed
 
design<-design[c(1:5,7:11),]
 
design<-design[c(1:5,7:11),]
  
 
fit<-lmFit(MM,design)
 
fit<-lmFit(MM,design)
−
cont.matrix<-makeContrasts(C.25vsM.I10.25=C.25-M.I10.25,       # Control versus MDR1 at 25wks
+
cont.matrix<-makeContrasts(C.25vsM.I10.25=C.25-M.I10.25,C.12vsM.I0.12=C.12-M.I0.12,
−
                          C.12vsM.I0.12=C.12-M.I0.12,         # Control versus MDR1 at 12 wks
+
                           C.25vsC.12=C.25-C.12,M.I10.25vsM.I0.12=M.I10.25-M.I0.12,
−
                           C.25vsC.12=C.25-C.12,               # Control at 12wks versus 25wks
 
−
                          M.I10.25vsM.I0.12=M.I10.25-M.I0.12, # MDR1 12wks versus 25wks
 
 
                           levels=design)
 
                           levels=design)
 
fit2<-contrasts.fit(fit,cont.matrix)
 
fit2<-contrasts.fit(fit,cont.matrix)
 
fit2<-eBayes(fit2)
 
fit2<-eBayes(fit2)
−
 
−
# ===============================================================================
 
−
# first data point
 
−
MM$M[1,]
 
−
# fit coefficients
 
−
fit$coef[1,]
 
−
# check
 
−
sum(MM$M[1,] * design[,1])/2
 
−
sum(MM$M[1,] * design[,3])/2
 
−
 
−
# Illustrates what contrasts matrix is doing
 
−
all.equal(fit$coef %*% cont.matrix, fit2$coef)
 
−
fit2$coef[1,]
 
−
# check for C.25vsC.12
 
−
sum(MM$M[1,] * design[,1])/2 - sum(MM$M[1,] * design[,3])/2
 
−
 
−
# ===============================================================================
 
  
 
par(mfrow=c(2,2))
 
par(mfrow=c(2,2))
−
p.n<-c("C.25vsM.I10.25","C.12vsM.I0.12","C.25vsC.12","M.I10.25vsM.I0.12")
+
p.n<-c("C.25vsM.I10.25","C.12vsM.I0.12","C.25vsC.12","M.I10.25vsM.I0.12") # Names on cont.matrix
 
tt<-kk<-list()
 
tt<-kk<-list()
 
for(i in 1:length(p.n)){
 
for(i in 1:length(p.n)){
Line 292: Line 267:
 
#
 
#
  
−
write.table(file="out.txt",RG$genes$Name[kk],quote=F,row.names=F,col.names=F)
+
#write.table(file="out.txt",RG$genes$Name[kk],quote=F,row.names=F,col.names=F) # this fails...
 +
write.table(file="out.txt",RG$genes$Name[sort(unique(unlist(kk)))],quote=F,row.names=F,col.names=F)
  
 
ord1<-order(as.numeric(rownames(tt[[1]])))
 
ord1<-order(as.numeric(rownames(tt[[1]])))
Line 301: Line 277:
 
               MCb.n$genes$Name,
 
               MCb.n$genes$Name,
 
               MCb.n$M[,1],-MCb.n$M[,5],
 
               MCb.n$M[,1],-MCb.n$M[,5],
−
              MCb.n$M[,3],MCb.n$M[,6],-MCb.n$M[,7],
+
#              MCb.n$M[,3],MCb.n$M[,6],-MCb.n$M[,7], # This should be only 4 + 11
 +
              MCb.n$M[,3],-MCb.n$M[,7], # This should be only 4 + 11
 
               -MCb.n$M[,2],MCb.n$M[,11],
 
               -MCb.n$M[,2],MCb.n$M[,11],
 
               -MCb.n$M[,4],MCb.n$M[,8],
 
               -MCb.n$M[,4],MCb.n$M[,8],
Line 315: Line 292:
 
               tt[[4]][ord4,1:3])
 
               tt[[4]][ord4,1:3])
 
colnames(out.mat)<-c("ProbeID","ProbeName",
 
colnames(out.mat)<-c("ProbeID","ProbeName",
−
                    paste("Pool",c("1.C12","10.M.I0.12","11.C25","12.M.I10.25","2.C.12","3.C.25","4.C25","5.M.I10.25","7.M.I5.22","8.M.I5.22","9.M.I0.12")[c(1,5,3,6,7,2,11,4,8,9,10)],sep=""),"Ref",
+
#                    paste("Pool",c("1.C12","10.M.I0.12","11.C25","12.M.I10.25","2.C.12","3.C.25","4.C25","5.M.I10.25","7.M.I5.22","8.M.I5.22","9.M.I0.12")[c(1,5,3,6,7,2,11,4,8,9,10)],sep=""),"Ref",
 +
                    paste("Pool",c("1.C12","10.M.I0.12","11.C25","12.M.I10.25","2.C.12","4.C25","5.M.I10.25","7.M.I5.22","8.M.I5.22","9.M.I0.12")[c(1,5,3,7,2,11,4,8,9,10)],sep=""),"Ref",
 
                     paste(c("log2FC","t.stat","p.value"),".",p.n[1],sep=""),
 
                     paste(c("log2FC","t.stat","p.value"),".",p.n[1],sep=""),
 
                     paste(c("log2FC","t.stat","p.value"),".",p.n[2],sep=""),
 
                     paste(c("log2FC","t.stat","p.value"),".",p.n[2],sep=""),
Line 323: Line 301:
 
                
 
                
 
write.table(out.mat,file='results.txt',sep='\t',row.names=F,col.names=T,quote=F)
 
write.table(out.mat,file='results.txt',sep='\t',row.names=F,col.names=T,quote=F)
 +
 +
 +
 +
# Checks
 +
# design maps the dye swaps
 +
# cont.matrix maps the contrasts
 +
 +
i <- 1
 +
p.n[i]
 +
colnames(design)
 +
toptable(fit2,coef=i,adjust="BH",number=nrow(RG$R))[ord1,][1:10,]
 +
 +
ind <- 1:10
 +
# Applying the cont.matrix[,1] to second group
 +
 +
apply(cbind(MM$M[ind,]%*%design[,1], -MM$M[ind,]%*%design[,2]), 1, mean)
 +
 +
# Checking with outputs
 +
tt[[1]][ord1,][1:10,]

Latest revision as of 00:43, 6 June 2007

Code snipits and programs written in R, S or S-PLUS

  2. Load required R libraries.

library(limma) datadir <- "/Volumes/HD2/Data/Nutrigenomics/MultipleScans/GPR"

  1. Read targets (Miks file)
  2. targets<-readTargets()
  1. Define weight function to read in the recorded flags in each GPR file

wt.obs<-function(gpr){ 

 flagged <- (gpr[, "Flags"])
 flagged

}

  1. Read in GPR files
  2. Currently not picking up any weight information
  3. names=targets$Name,

RG<-read.maimages(files=dir(datadir),source="genepix",path=datadir,columns=list(Rf = "F635 Median", Gf = "F532 Median",Rb = "B635 Median", Gb = "B532 Median")) RG$printer <- getLayout(RG$genes)

  1. Define object containing weights for each GPR file, and set weights to "NULL" in
  2. RG variables (otherwise weights are used in normalization).

RG.weights<-RG$weights

RG$weights=NULL

save(list=c("RG","RG.weights","targets"),file="RGw.data")

  2. Set zero values in foreground and background to 0.5

load("RGw.data") library(limma)

RG$R[RG$R==0]<-0.5 RG$G[RG$G==0]<-0.5 RG$Rb[RG$Rb==0]<-0.5 RG$Gb[RG$Gb==0]<-0.5

  2. Create objects containing background corrected intensity data, and non-background
  3. (RGh) corrected intensity data (RGb).

RGb<-backgroundCorrect(RG,method="none") RGh<-backgroundCorrect(RG,method="half")

  2. Create MA objects based on RGb and RGh, prior to normalization.

MAh<-normalizeWithinArrays(RGh,method="none") MAb<-normalizeWithinArrays(RGb,method="none")

  2. Perform loess normalization on MAh and MAb.

MAh.n<-normalizeWithinArrays(RGh,method='loess') MAb.n<-normalizeWithinArrays(RGb,method='loess')

save(list=c("RGb","RGh","MAh","MAb","MAh.n","MAb.n","RG.weights"),file="normw.data")

  2. Image plots

load('RGw.data') for(i in 1:(dim(RG$R)[2]/8)){ # BAD CODING ... Should be 1:(ncol(RG$R)/9) 

 fname<-paste("JPG/BGFG/BG_FG_",i,".png",sep="") ... pngs better
  1. jpeg(fname, height = 900, width=1200, quality=80, bg="white")
 png(fname, height = 900, width=1200, pointsize=15)  
 par(mfrow=c(8,4))
 for(j in 1:8){
   ar<-(i-1)*8+j
   imageplot(log2(RG$Rb[,ar]), RG$printer, low="white", high="red")
   text(1,80,colnames(RG$R),pos=4)
   imageplot(log2(RG$Gb[,ar]), RG$printer, low="white", high="green")
   text(1,80,colnames(RG$R),pos=4)
   imageplot(log2(RG$R[,ar]), RG$printer, low="white", high="red")
   imageplot(log2(RG$G[,ar]), RG$printer, low="white", high="green")
 }
 dev.off()

}

  2. Create histograms of red and green background and foreground intensities
  3. for each array. Four arrays per page.

load('RGw.data') load('normw.data')

for(i in 1:(ncol(RG$R)/4)){ 

 fname<-paste("JPG/Hist/BGFG_hist_",i,".jpg",sep="")
 jpeg(fname, height = 900, width=1200, quality=80, bg="white")
  1. par(mfrow=c(4,5))
 par(mfrow=c(4,4))
 for(j in 1:4){
   ar<-(i-1)*4+j
   hist(log2(RG$Rb[,ar]),50,main=colnames(RG$R)[ar],xlab="log2(Rb)",col='red')
   hist(log2(RG$Gb[,ar]),50,main=colnames(RG$R)[ar],xlab="log2(Gb)",col='green')
   hist(log2(RG$R[,ar]),50,main=colnames(RG$R)[ar],xlab="log2(R)",col='red')
   hist(log2(RG$G[,ar]),50,main=colnames(RG$R)[ar],xlab="log2(G)",col='green')
  1. plot(density(log2(RG$G[,ar])),col='green',lty=1)
  2. lines(density(log2(RG$Rb[,ar])),col='red',lty=2)
  3. lines(density(log2(RG$R[,ar])),col='red',lty=1)
  4. lines(density(log2(RG$Gb[,ar])),col='green',lty=2)
 }
 dev.off()

}

  2. Boxplots for NBC data (Pre-norm and post-norm)

lo<-grep("low", colnames(RG)) md<-grep("average", colnames(RG)) hi<-grep("high", colnames(RG))

jpeg("JPG/Boxplots/boxplot-full-nbc-pre.jpg",height=900,width=1200,quality=80, bg="white") boxplot(MAb$M~col(MAb$M),names=targets$Name,main="Pre-normalization (NBC): all arrays",cex.axis=0.7) dev.off() jpeg("JPG/Boxplots/boxplot-lo-nbc-pre.jpg",height=900,width=1200,quality=80,bg="white") boxplot(MAb$M[,lo]~col(MAb$M[,lo]),names=targets$Name[lo],main="Pre-normalization (NBC): low scans",cex.axis=0.7) dev.off() jpeg("JPG/Boxplots/boxplot-md-nbc-pre.jpg",height=900,width=1200,quality=80,bg="white") boxplot(MAb$M[,md]~col(MAb$M[,md]),names=targets$Name[md],main="Pre-normalization (NBC): medium scans",cex.axis=0.7) dev.off() jpeg("JPG/Boxplots/boxplot-hi-nbc-pre.jpg",height=900,width=1200,quality=80,bg="white") boxplot(MAb$M[,hi]~col(MAb$M[,hi]),names=targets$Name[hi],main="Pre-normalization (NBC): high arrays",cex.lab=0.7) dev.off()

jpeg("JPG/Boxplots/boxplot-full-nbc-post.jpg",height=900,width=1200,quality=80, bg="white") boxplot(MAb.n$M~col(MAb.n$M),names=targets$Name,main="Post-normalization (NBC): all arrays",cex.axis=0.7) dev.off() jpeg("JPG/Boxplots/boxplot-lo-nbc-post.jpg",height=900,width=1200,quality=80,bg="white") boxplot(MAb.n$M[,lo]~col(MAb.n$M[,lo]),names=targets$Name[lo],main="Post-normalization (NBC): low scans",cex.axis=0.7) dev.off() jpeg("JPG/Boxplots/boxplot-md-nbc-post.jpg",height=900,width=1200,quality=80,bg="white") boxplot(MAb.n$M[,md]~col(MAb.n$M[,md]),names=targets$Name[md],main="Post-normalization (NBC): medium scans",cex.axis=0.7) dev.off() jpeg("JPG/Boxplots/boxplot-hi-nbc-post.jpg",height=900,width=1200,quality=80,bg="white") boxplot(MAb.n$M[,hi]~col(MAb.n$M[,hi]),names=targets$Name[hi],main="Post-normalization (NBC): high arrays",cex.axis=0.7) dev.off()

  2. Boxplots for HBC data (Pre-norm and post-norm)

jpeg("JPG/Boxplots/boxplot-full-hbc-pre.jpg",height=900,width=1200,quality=80, bg="white") boxplot(MAh$M~col(MAh$M),names=targets$Name,main="Pre-normalization (HBC): all arrays",cex.axis=0.7) dev.off() jpeg("JPG/Boxplots/boxplot-lo-hbc-pre.jpg",height=900,width=1200,quality=80,bg="white") boxplot(MAh$M[,lo]~col(MAh$M[,lo]),names=targets$Name[lo],main="Pre-normalization (HBC): low scans",cex.axis=0.7) dev.off() jpeg("JPG/Boxplots/boxplot-md-hbc-pre.jpg",height=900,width=1200,quality=80,bg="white") boxplot(MAh$M[,md]~col(MAh$M[,md]),names=targets$Name[md],main="Pre-normalization (HBC): medium scans",cex.axis=0.7) dev.off() jpeg("JPG/Boxplots/boxplot-hi-hbc-pre.jpg",height=900,width=1200,quality=80,bg="white") boxplot(MAh$M[,hi]~col(MAh$M[,hi]),names=targets$Name[hi],main="Pre-normalization (HBC): high arrays",cex.lab=0.7) dev.off()

jpeg("JPG/Boxplots/boxplot-full-hbc-post.jpg",height=900,width=1200,quality=80, bg="white") boxplot(MAh.n$M~col(MAh.n$M),names=targets$Name,main="Post-normalization (HBC): all arrays",cex.axis=0.7) dev.off() jpeg("JPG/Boxplots/boxplot-lo-hbc-post.jpg",height=900,width=1200,quality=80,bg="white") boxplot(MAh.n$M[,lo]~col(MAh.n$M[,lo]),names=targets$Name[lo],main="Post-normalization (HBC): low scans",cex.axis=0.7) dev.off() jpeg("JPG/Boxplots/boxplot-md-hbc-post.jpg",height=900,width=1200,quality=80,bg="white") boxplot(MAh.n$M[,md]~col(MAh.n$M[,md]),names=targets$Name[md],main="Post-normalization (HBC): medium scans",cex.axis=0.7) dev.off() jpeg("JPG/Boxplots/boxplot-hi-hbc-post.jpg",height=900,width=1200,quality=80,bg="white") boxplot(MAh.n$M[,hi]~col(MAh.n$M[,hi]),names=targets$Name[hi],main="Post-normalization (HBC): high arrays",cex.axis=0.7) dev.off()

  1. Mb.nn.hi<-normalizeBetweenArrays(MAb.n$M[,odd])
  2. Mb.nn.med<-normalizeBetweenArrays(MAb.n$M[,even])
  3. Mh.nn.hi<-normalizeBetweenArrays(MAh.n$M[,odd])
  4. Mh.nn.med<-normalizeBetweenArrays(MAh.n$M[,even])
  5. save(list=c("Mh.nn.hi","Mh.nn.med","Mb.nn.hi","Mb.nn.med"),file="norm2w.data")
  2. Find saturated spots and create combined data sets

library(limma)

load("RGw.data") load("normw.data")

  1. load("norm2w.data")

sat.r<-RG$R==(2^16-1) # Alot more saturated in R channel sat.g<-RG$G==(2^16-1) sat<-sat.r | sat.g

  1. lo<-grep("l",targets$Name) Already done above
  2. md<-grep("a",targets$Name)
  3. hi<-grep("h",targets$Name)

tsat<-sat[,lo] & sat[,md] & sat[,hi]

sum(tsat) sum(apply(tsat,1,sum)>0)

MCb.n<-MAb.n MCb.n$M<-MAb.n$M[,hi] MCb.n$M[sat[,hi]]<-MAb.n$M[,md][sat[,hi]] MCb.n$M[sat[,md]]<-MAb.n$M[,lo][sat[,md]] MCb.n$A<-MAb.n$A[,hi] MCb.n$A[sat[,hi]]<-MAb.n$A[,md][sat[,hi]] MCb.n$A[sat[,md]]<-MAb.n$A[,lo][sat[,md]]

MCh.n<-MAh.n MCh.n$M<-MAh.n$M[,hi] MCh.n$M[sat[,hi]]<-MAh.n$M[,md][sat[,hi]] MCh.n$M[sat[,md]]<-MAh.n$M[,lo][sat[,md]] MCh.n$A<-MAh.n$A[,hi] MCh.n$A[sat[,hi]]<-MAh.n$A[,md][sat[,hi]] MCh.n$A[sat[,md]]<-MAh.n$A[,lo][sat[,md]]

  1. boxplot(MCh.n~col(MCh.n))
  2. boxplot(MCh.n~col(MCh.n))
  1. par(mfrow=c(3,4))
  2. for(i in 1:ncol(MCb.n$M)){
  3. plotMA(MCb.n,array=i)
  4. }
  5. par(mfrow=c(3,4))
  6. for(i in 1:ncol(MCh.n$M)){
  7. plotMA(MCh.n,array=i)
  8. }
  2. Linear model analysis

targ<-targets[md,] design<-matrix(c(0, 0,1, 0, 0,0,-1,0,0, 0,0,0, # C.25 

                0, 0,0,-1, 0,0, 0,1,0, 0,0,0,  # M.I10.25
                1, 0,0, 0,-1,0, 0,0,0, 0,0,0,  # C.12
                0,-1,0, 0, 0,0, 0,0,0, 0,1,0,  # M.I0.12
                0, 0,0, 0, 0,0, 0,0,1,-1,0,0,  # M.I5.22                 
                ),12,5)

colnames(design)<-c("C.25","M.I10.25","C.12","M.I0.12","M.I5.22") rownames(design)<- substr(colnames(RG)[md], 20, 27) # targets$Name[md]

MM<-MCb.n[,c(1:5,7:11)] # Pool 3 and ref-ref removed design<-design[c(1:5,7:11),]

fit<-lmFit(MM,design) cont.matrix<-makeContrasts(C.25vsM.I10.25=C.25-M.I10.25,C.12vsM.I0.12=C.12-M.I0.12, 

                          C.25vsC.12=C.25-C.12,M.I10.25vsM.I0.12=M.I10.25-M.I0.12,
                          levels=design)

fit2<-contrasts.fit(fit,cont.matrix) fit2<-eBayes(fit2)

par(mfrow=c(2,2)) p.n<-c("C.25vsM.I10.25","C.12vsM.I0.12","C.25vsC.12","M.I10.25vsM.I0.12") # Names on cont.matrix tt<-kk<-list() for(i in 1:length(p.n)){ 

 tti<-toptable(fit2,coef=i,adjust="BH",number=nrow(RG$R))
 kki<-as.numeric(rownames(tti)[1:20])

} for(i in 1:length(p.n)){ 

 MA<-list()
 MA$M<-fit2$coef[,i]
 MA$A<-fit2$Amean
 plotMA(MA,main=p.n[i],ylim=c(-4,4))
 points(fit2$Amean[kki],fit2$coef[kki,i],pch=16,col="red")
  1. points(fit2$Amean[kkabs(i-3)],fit2$coef[kkabs(i-3),i],pch=16,col="green")

}

  1. write.table(file="out.txt",RG$genes$Name[kk],quote=F,row.names=F,col.names=F) # this fails...

write.table(file="out.txt",RG$genes$Name[sort(unique(unlist(kk)))],quote=F,row.names=F,col.names=F)

ord1<-order(as.numeric(rownames(tt1))) ord2<-order(as.numeric(rownames(tt2))) ord3<-order(as.numeric(rownames(tt3))) ord4<-order(as.numeric(rownames(tt4))) out.mat<-cbind(MCb.n$genes$ID, 

              MCb.n$genes$Name,
              MCb.n$M[,1],-MCb.n$M[,5],
  1. MCb.n$M[,3],MCb.n$M[,6],-MCb.n$M[,7], # This should be only 4 + 11
              MCb.n$M[,3],-MCb.n$M[,7], # This should be only 4 + 11
              -MCb.n$M[,2],MCb.n$M[,11],
              -MCb.n$M[,4],MCb.n$M[,8],
              MCb.n$M[,9],-MCb.n$M[,10],
              MCb.n$M[,12],
  1. tt1[as.numeric(rownames(tt1)),1:3],
  2. tt2[as.numeric(rownames(tt2)),1:3],
  3. tt3[as.numeric(rownames(tt3)),1:3],
  4. tt4[as.numeric(rownames(tt4)),1:3])
              tt1[ord1,1:3],
              tt2[ord2,1:3],
              tt3[ord3,1:3],
              tt4[ord4,1:3])

colnames(out.mat)<-c("ProbeID","ProbeName",

  1. paste("Pool",c("1.C12","10.M.I0.12","11.C25","12.M.I10.25","2.C.12","3.C.25","4.C25","5.M.I10.25","7.M.I5.22","8.M.I5.22","9.M.I0.12")[c(1,5,3,6,7,2,11,4,8,9,10)],sep=""),"Ref",
                    paste("Pool",c("1.C12","10.M.I0.12","11.C25","12.M.I10.25","2.C.12","4.C25","5.M.I10.25","7.M.I5.22","8.M.I5.22","9.M.I0.12")[c(1,5,3,7,2,11,4,8,9,10)],sep=""),"Ref",
                    paste(c("log2FC","t.stat","p.value"),".",p.n[1],sep=""),
                    paste(c("log2FC","t.stat","p.value"),".",p.n[2],sep=""),
                    paste(c("log2FC","t.stat","p.value"),".",p.n[3],sep=""),
                    paste(c("log2FC","t.stat","p.value"),".",p.n[4],sep=""))


write.table(out.mat,file='results.txt',sep='\t',row.names=F,col.names=T,quote=F)


  1. Checks
  2. design maps the dye swaps
  3. cont.matrix maps the contrasts

i <- 1 p.n[i] colnames(design) toptable(fit2,coef=i,adjust="BH",number=nrow(RG$R))[ord1,][1:10,]

ind <- 1:10

  1. Applying the cont.matrix[,1] to second group

apply(cbind(MM$M[ind,]%*%design[,1], -MM$M[ind,]%*%design[,2]), 1, mean)

  1. Checking with outputs

tt1[ord1,][1:10,]

Retrieved from "https://organicdesign.nz/wiki/index.php?title=Nunz-mik-12-12-05.R&oldid=71174"