Triplicate Comparison for Microarrays

(e.g. 3X Wild type vs 3X Mutant)

Entering Data sets from WT and Modulated Experiments into R and combining

x46<-read.table("C://Documents and Settings/cucchecl/Desktop/46-48vs49-51.txt")

n<-length(x46[,2])

z0<-x46[,2]; z1<-x46[,3]; z2<-x46[,4]; z3<-x46[,5]; z4<-x46[,6]; z5<-x46[,7]; z6<-x46[,8]

zz<- cbind(z1,z2,z3,z4,z5,z6)

Normalizing Wild Type Data

qwz<-(1/3)*(sort(z1)+sort(z2)+sort(z3))

qw1<- qwz[rank(z1)]; qw2<- qwz[rank(z2)]; qw3<- qwz[rank(z3)]

Normalizing Modulated/Mutant Data

qmz<-(1/3)*(sort(z4)+sort(z5)+sort(z6))

qm1<- qmz[rank(z4)]; qm2<- qmz[rank(z5)]; qm3<- qmz[rank(z6)]

Combining Both Sets of Data and Displaying Boxplot

zz<- cbind(qw1,qw2,qw3,qm1,qm2,qm3)

boxplot(qw1,qw2,qw3,qm1,qm2,qm3, main="Quantile Normalized Box Plots")

Normalizing Data Together and Displaying Boxplot

 qnmz<- (1/6)*(sort(qw1)+sort(qw2) +sort(qw3) + sort(qm1) + sort(qm2) + sort(qm3))

 w1<- qnmz[rank(qw1)]; w2<- qnmz[rank(qw2)]; w3<- qnmz[rank(qw3)]; mm1<- qnmz[rank(qm1)]; mm2<- qnmz[rank(qm2)]; mm3<- qnmz[rank(qm3)]

 zzn<- cbind(w1,w2,w3,mm1,mm2,mm3)

 boxplot(w1,w2,w3,mm1,mm2,mm3, main="Quantile Normalized (wild-type and modulated) Box Plots")

Shapiro Test

m<-length(w1)

      stat1<- 1:m;  pv1<- 1:m

for(i in 1:m)

+ {

+ stat1[i] <- shapiro.test(zzn[i,])$ statistic

+ pv1[i] <- shapiro.test(zzn[i,])$ p.value

+ }

yn <- data.frame(z0,stat1,pv1)

yn1 <- yn[order(pv1, decreasing =F),]

yn2 <- length(which(yn1[,3] <= .05)) Removing Data that fails test conditions

yn3<- yn1[1:yn2,]

yn3[1:5,]

Wilcoxon Test on Data which failed Shapiro Test conditions

gsn <-as.numeric(rownames(yn3))

zzn1 <- zzn[gsn,]

m<- length(zzn1[,1])

stat<- 1:m

pv<- 1:m

for(i in 1:m)

+ {stat[i] <- wilcox.test(zzn1[i,], mu=0)$ statistic

+ pv[i] <- wilcox.test(zzn1[i,], mu=0)$ p.value

+ }

gn<-yn3[,1]

yn4<-data.frame(gsn, gn, stat, pv)

attach(yn4)

yn5 <- yn4[order(pv, decreasing =F),] Displays Data in order of decreasing pv

wn <- length(which(yn5[,4] <= .05)) Number genes that are significantly different using Wilcoxon (Mann-Whitney) non-parametric test

T-test conducted on Data which passed Normality Test

ynn2 <- length(which(yn1[,3] > .05)) Verification of total number of genes (yn2+ynn2=n)

m1<- yn2+1

ynn4<- yn1[m1:n,]

length(ynn4[,2])

gsn <-as.numeric(rownames(ynn4))

zzt <- zzn[gsn,]

m<- length(zzt[,1])

stat<-1:m

pv<-1:m

for(i in 1:m)

+ {

+ stat[i] <- t.test(zzt[i,], mu=0)$ statistic

+ pv[i] <- t.test(zzt[i,], mu=0)$ p.value

+ }

gn<-ynn4[,1]

yn6<-data.frame(gsn, gn, stat, pv)

attach(yn6)

yn7 <- yn6[order(pv, decreasing =F),] Displays Data in order of decreasing pv

tn <- length(which(yn7[,4] <= .05))

tn

Bonferroni

pvbf<-1:m

bfyn <- transform(yn7,pvbf =pv*m)

pvbf<-bfyn[,5]

{if (pvbf[i] >1) (pvbf[i] <-1)}

pvbf[1:5]

bfyn[1:5]

length(which(pvbf <=0.05))

bfyn[1:5,5]

length(which(pvbf >0.05))

length(which(bfyn[,5] <=0.05));length(which(bfyn[,5] <=0.05 & stat > 0))

length(which(bfyn[,5] <=0.05));length(which(bfyn[,5] <=0.05 & stat < 0))

bfyn[1:100,]

Principle Component Analysis

gs<-as.numeric(rownames(bfyn))

gs1<-gs[1:50]

vpc <- zzt[gs1,]

pc.cr<-princomp(vpc,cor=T,scores=T)

summary(pc.cr)

wpc1<-pc.cr$loadings

wpc1

biplot(pc.cr,main="Biplot with samples as variables",xlabs=as.character(1:50),ylabs=c(rep("W",3),rep("M",3)),lwd=2,cex=.7)

p1 <- prcomp(t(vpc), scale = T); summary(p1)

wpc <- p1$rotation

biplot(p1,main="Biplot with genes as variables",ylabs=as.character(1:50),xlabs=c(rep("W",3),rep("M",3)),lwd=2,cex=.7)

wpc2<- wpc[,1:3];  wpc2;  #Weights of the first 2 principal component

********************************************************************