Protoarray/Huprotarray common Analysis

library(limma)
library(ggplot2)
library(RColorBrewer)
targets.protoarray <- readTargets('../rawdata/protoarray_common.csv', sep=',')
targets.huprot <- readTargets('../rawdata_huprot/huprot_common.csv', sep=',')
targets.huprot$Condition <- as.factor(targets.huprot$Condition)
targets.protoarray$Condition <- as.factor(targets.protoarray$Condition)
targetColors.huprot <-
    with(targets.huprot,
         data.frame(Condition = levels(Condition),
                    color = I(brewer.pal(nlevels(Condition), name = 'Dark2'))))
targetColors.protoarray <-
    with(targets.protoarray,
         data.frame(Condition = levels(Condition),
                    color = I(brewer.pal(nlevels(Condition), name = 'Dark2'))))
x.protoarray <- read.maimages(targets.protoarray$FileName,
                   path=file.path('..', 'rawdata'),
                   source='genepix',
                   columns=list(G='F635 Mean', Gb='B635 Median'),
                   green.only=TRUE)

Read ../rawdata/control_67336_18052013.gpr 
Read ../rawdata/control_68051_21052013.gpr 
Read ../rawdata/control_68213_11052013.gpr 
Read ../rawdata/control_68473_05062013.gpr 
Read ../rawdata/Grade_L_68204_11052013.gpr 
Read ../rawdata/Grade_L_67341_18052013.gpr 
Read ../rawdata/Grade_L_68210_13052013.gpr 
Read ../rawdata/Grade_III_67346_18052013.gpr 
Read ../rawdata/Grade_III_68001_21052013.gpr 
Read ../rawdata/Grade_III_68052_21052013.gpr 
Read ../rawdata/Grade_III_68207_11052013.gpr 
Read ../rawdata/Grade_III_68211_11052013.gpr 
Read ../rawdata/Grade_III_68476_18052013_2.gpr 
Read ../rawdata/Grade_IV_68037_21052013.gpr 
Read ../rawdata/Grade_IV_68056_21052013_2.gpr 
Read ../rawdata/Grade_IV_68209_11052013.gpr 
Read ../rawdata/Grade_IV_68212_11052013.gpr 
Read ../rawdata/Grade_IV_68481_18052013.gpr 
Read ../rawdata/Grade_IV_68055_21052013_2.gpr 
Read ../rawdata/Grade_IV_67347_18092013_2.gpr 

x.huprot <- read.maimages(targets.huprot$FileName,
                   path=file.path('..', 'rawdata_huprot'),
                   source='genepix',
                   columns=list(G='F635 Mean', Gb='B635 Median'),
                   green.only=TRUE)

Read ../rawdata_huprot/Control_1256587_2000143978.gpr 
Read ../rawdata_huprot/Control_1311449_2000143966.gpr 
Read ../rawdata_huprot/Control_1312497_2000143964.gpr 
Read ../rawdata_huprot/Control_1313058_2000143976.gpr 
Read ../rawdata_huprot/Grade_L_CJ_14434_2000143958.gpr 
Read ../rawdata_huprot/Grade_L_CH_27096_2000155729.gpr 
Read ../rawdata_huprot/Grade_L_CJ_3659_2000143955.gpr 
Read ../rawdata_huprot/Grade_III_CH_15796_2000143692.gpr 
Read ../rawdata_huprot/Grade_III_CJ_11032_2000143957.gpr 
Read ../rawdata_huprot/Grade_III_CJ_27572_2000144442.gpr 
Read ../rawdata_huprot/Grade_III_CJ_9761_2000144443.gpr 
Read ../rawdata_huprot/Grade_III_CK_5121_2000144441.gpr 
Read ../rawdata_huprot/Grade_III_CK_6451_2000144440.gpr 
Read ../rawdata_huprot/Grade_IV_CH_12724_2000143977.gpr 
Read ../rawdata_huprot/Grade_IV_CH_29738_2000143975.gpr 
Read ../rawdata_huprot/Grade_IV_CH_31148_2000143961.gpr 
Read ../rawdata_huprot/Grade_IV_CH_32705_2000155727.gpr 
Read ../rawdata_huprot/Grade_IV_CJ_11291_2000143969.gpr 
Read ../rawdata_huprot/Grade_IV_CJ_12441_2000143979.gpr 
Read ../rawdata_huprot/Grade_IV_CK_6570_2000143967.gpr 

x.huprot$genes$Name <- gsub('\\n', '', x.huprot$genes$Name)
x.protoarray$genes$Name <- gsub('\\n', '', x.protoarray$genes$Name)

Extract common targets

controls.protoarray <- c('Buffer', 
                         'BSA~N/A', 
                         'GST1~RFU:322.74',
                         'GST2~RFU:578.54', 
                         'GST3~RFU:1770.2',
                         'GST4~RFU:4230.18',
                         'AlexaAntiMouseAb~N/A',
                         'HumanIgG1~N/A',
                         'HumanIgG2~N/A', 
                         'HumanIgG3~N/A',
                         'HumanIgG4~N/A', 
                         'Anti-HumanIgG1~N/A',
                         'Anti-HumanIgG2~N/A',
                         'Anti-HumanIgG3~N/A',
                         'Anti-HumanIgG4~N/A', 
                         'V5control5~N/A',
                         'V5control4~N/A', 
                         'V5control3~N/A', 
                         'V5control2~N/A',
                         'V5control1~N/A'
                         )
clean_name <- function(name){
  if (name %in% controls.protoarray){
    return(name)
  }
  name <- strsplit(name, ":")[[1]][2]
  name <- strsplit(name, "~")[[1]][1]
  return(name)
}

negative.controls.huprot <- x.huprot$genes$Name %in% c('BSA', 
                                         'HeLa cell lysates',
                                         'p300-BHC')
positive.controls.huprot <- x.huprot$genes$Name %in% c('H2A', 
                                         'H2B', 
                                          'H3', 
                                          'H4', 
                                          'GST10n',
                                          'GST50n',
                                          'GST100n', 
                                          'GST500n')

new_name <- sapply(as.vector(x.protoarray$genes$Name), clean_name)
x.protoarray$genes$Name <- new_name
x.protoarray <- x.protoarray[x.protoarray$genes$Name %in% c(x.huprot$genes$ID, controls.protoarray),]
negative.controls.protoarray <- x.protoarray$genes$Name %in% c('Buffer', 
                                         'BSA~N/A', 
                                         'GST1~RFU:322.74',
                                         'GST2~RFU:578.54', 
                                         'GST3~RFU:1770.2',
                                         'GST4~RFU:4230.18')
positive.controls.protoarray <- x.protoarray$genes$Name %in% c('AlexaAntiMouseAb~N/A',
                                         'HumanIgG1~N/A',
                                         'HumanIgG2~N/A', 
                                         'HumanIgG3~N/A',
                                         'HumanIgG4~N/A', 
                                         'Anti-HumanIgG1~N/A',
                                         'Anti-HumanIgG2~N/A',
                                         'Anti-HumanIgG3~N/A',
                                         'Anti-HumanIgG4~N/A', 
                                         'V5control5~N/A',
                                         'V5control4~N/A', 
                                         'V5control3~N/A', 
                                         'V5control2~N/A',
                                         'V5control1~N/A')

par(mar=c(8,4,4,4)+.1)
boxplot(log2(x.huprot$E),las=2, main="Pre-normalization Huprot",ylab="log2-intensity", col=targetColors.huprot$color[match(targets.huprot[colnames(x.huprot$E),]$Condition, targetColors.huprot$Condition)], pch=19, cex.axis=0.7, cex=0.7)
par(mar=c(8,4,4,4)+.1)

boxplot(log2(x.protoarray$E),las=2, main="Pre-normalization Protoarray",ylab="log2-intensity", col=targetColors.protoarray$color[match(targets.protoarray[colnames(x.protoarray$E),]$Condition, targetColors.protoarray$Condition)], pch=19, cex.axis=0.7, cex=0.7)

x.huprot$genes$Status  <- rep('regular', nrow(x.huprot))
x.huprot$genes$Status[positive.controls.huprot] <- 'POSITIVE'
x.huprot$genes$Status[negative.controls.huprot] <- 'NEGATIVE'
y.huprot <- neqc(x.huprot)
x.protoarray$genes$Status  <- rep('regular', nrow(x.protoarray))
x.protoarray$genes$Status[positive.controls.protoarray] <- 'POSITIVE'
x.protoarray$genes$Status[negative.controls.protoarray] <- 'NEGATIVE'
y.protoarray <- neqc(x.protoarray)

par(mar=c(8,4,4,4)+.1)
boxplot(log2(y.huprot$E),las=2, main="Post-normalization Huprot",ylab="log2-intensity", col=targetColors.huprot$color[match(targets.huprot[colnames(y.huprot$E),]$Condition, targetColors.huprot$Condition)], pch=19, cex.axis=0.7, cex=0.7)
par(mar=c(8,4,4,4)+.1)

boxplot(log2(y.protoarray$E),las=2, main="Post-normalization Protoarray",ylab="log2-intensity", col=targetColors.protoarray$color[match(targets.protoarray[colnames(y.protoarray$E),]$Condition, targetColors.protoarray$Condition)], pch=19, cex.axis=0.7, cex=0.7)

opar <- par(pch = 19)
pc = prcomp( t( y.huprot$E  ), scale=TRUE )
pc.summary <- summary(pc)
pc1.var <- pc.summary$importance[2,1]
pc2.var <- pc.summary$importance[2,2]
plot(pc$x[, 1:2], col=targetColors.huprot$color[match(targets.huprot$Condition, targetColors.huprot$Condition)], cex=1.5, xlab = paste('PC1: ', 100*pc1.var, '%', sep=''), ylab =  paste('PC2: ', 100*pc2.var, '%', sep=''),     main= "PCA normalized intensities huprot")
legend(x = 'topright', 
       legend = as.character(targetColors.huprot$Condition),
       col = targetColors.huprot$color, pch = par("pch"), bty = 'n', xjust = 1)

opar <- par(pch = 19)
pc = prcomp( t( y.protoarray$E  ), scale=TRUE )
pc.summary <- summary(pc)
pc1.var <- pc.summary$importance[2,1]
pc2.var <- pc.summary$importance[2,2]
plot(pc$x[, 1:2], col=targetColors.protoarray$color[match(targets.protoarray$Condition, targetColors.protoarray$Condition)], cex=1.5, xlab = paste('PC1: ', 100*pc1.var, '%', sep=''), ylab =  paste('PC2: ', 100*pc2.var, '%', sep=''),     main= "PCA normalized intensities protoarray")
legend(x = 'topright', 
       legend = as.character(targetColors.protoarray$Condition),
       col = targetColors.protoarray$color, pch = par("pch"), bty = 'n', xjust = 1)

MDS plot

opar <- par(pch = 19)
plotMDS(y.huprot,  col=targetColors.huprot$color[match(targets.huprot$Condition, targetColors.huprot$Condition)], cex=1.2, pch=19, main='MDS Huprot')
legend(x = 'topright', 
       legend = as.character(targetColors.huprot$Condition),
       col = targetColors.huprot$color, pch = par("pch"), bty = 'n', xjust = 1)

opar <- par(pch = 19)
plotMDS(y.protoarray,  col=targetColors.protoarray$color[match(targets.protoarray$Condition, targetColors.protoarray$Condition)], cex=1.2, pch=19, main='MDS Protoarray')
legend(x = 'topright', 
       legend = as.character(targetColors.protoarray$Condition),
       col = targetColors.protoarray$color, pch = par("pch"), bty = 'n', xjust = 1)

f.huprot <- factor(targets.huprot$Condition, levels = unique(targets.huprot$Condition))
design.huprot <- model.matrix(~0 + f.huprot)
colnames(design.huprot) <- levels(f.huprot)
f.protoarray <- factor(targets.protoarray$Condition, levels = unique(targets.protoarray$Condition))
design.protoarray <- model.matrix(~0 + f.protoarray)
colnames(design.protoarray) <- levels(f.protoarray)
corfit.protoarray <- duplicateCorrelation(y.protoarray,
                               design.protoarray,
                               ndups=2,
                               spacing=1
                               )
fit.protoarray <- lmFit(y.protoarray,
                        design.protoarray,
                        ndups=2,
                       correlation=corfit.protoarray$consensus.correlation)
contrast.matrix.protoarray <- makeContrasts(Grade2vsGrade1=Grade2-Grade1, 
                                 Grade3vsGrade2=Grade3-Grade2, 
                                 Grade4vsGrade3=Grade4-Grade3,
                                 Grade1vsControl=Grade1-Control,
                                 Grade2vsControl=Grade2-Control,
                                 Grade3vsControl=Grade3-Control,
                                 Grade4vsControl=Grade4-Control,                              GradevsControl=(Grade1+Grade2+Grade3+Grade4)/4-Control,
                              HighGradevsControl=(Grade3+Grade4)/2-Control,
                              LowGradevsControl=(Grade1+Grade2)/2-Control,
                  HighGradevsLowGrade=(Grade3+Grade4)/2-(Grade1+Grade2)/2,
                                 Grade3vsLowGrade=Grade3-(Grade1+Grade2)/2,
                                 Grade4vsLowGrade=Grade4-(Grade1+Grade2)/2,
                                 levels=design.protoarray)
fit2.protoarray <- contrasts.fit(fit.protoarray, contrast.matrix.protoarray)
fit2.protoarray <- eBayes(fit2.protoarray)
corfit.huprot <- duplicateCorrelation(y.huprot,
                               design.huprot,
                               ndups=2,
                               spacing=1
                               )
fit.huprot <- lmFit(y.huprot,
                        design.huprot,
                        ndups=2,
                       correlation=corfit.huprot$consensus.correlation)
contrast.matrix.huprot <- makeContrasts(Grade2vsGrade1=Grade2-Grade1, 
                                 Grade3vsGrade2=Grade3-Grade2, 
                                 Grade4vsGrade3=Grade4-Grade3,
                                 Grade1vsControl=Grade1-Control,
                                 Grade2vsControl=Grade2-Control,
                                 Grade3vsControl=Grade3-Control,
                                 Grade4vsControl=Grade4-Control,                              GradevsControl=(Grade1+Grade2+Grade3+Grade4)/4-Control,
                              HighGradevsControl=(Grade3+Grade4)/2-Control,
                              LowGradevsControl=(Grade1+Grade2)/2-Control,
                  HighGradevsLowGrade=(Grade3+Grade4)/2-(Grade1+Grade2)/2,
                                 Grade3vsLowGrade=Grade3-(Grade1+Grade2)/2,
                                 Grade4vsLowGrade=Grade4-(Grade1+Grade2)/2,
                                 levels=design.huprot)
fit2.huprot <- contrasts.fit(fit.huprot, contrast.matrix.huprot)
fit2.huprot <- eBayes(fit2.huprot)

coefs <- c('Grade2vsGrade1',
           'Grade3vsGrade2',
           'Grade4vsGrade3',
           'Grade1vsControl',
           'Grade2vsControl',
           'Grade3vsControl',
           'Grade4vsControl',
           'GradevsControl',
           'HighGradevsControl',
           'LowGradevsControl',
           'HighGradevsLowGrade',
           'Grade3vsLowGrade',
           'Grade4vsLowGrade')
for (coef in coefs){
  tt <- topTable(fit2.huprot, coef=coef, number = Inf)
  tt.sig <- subset(tt, tt$adj.P.Val<0.01)
  write.table(tt, file.path('..', 'results_common/huprot', paste(coef,'csv', sep='.') ), row.names = FALSE)
  write.table(tt.sig, file.path('..', 'results_common/huprot', 'significant_filtered', paste(coef,'csv', sep='.') ), row.names = FALSE)
}

for (coef in coefs){
  tt <- topTable(fit2.protoarray, coef=coef, number = Inf)
  tt.sig <- subset(tt, tt$adj.P.Val<0.01)
  write.table(tt, file.path('..', 'results_common/protoarray', paste(coef,'csv', sep='.') ), row.names = FALSE)
  write.table(tt.sig, file.path('..', 'results_common/protoarray', 'significant_filtered', paste(coef,'csv', sep='.') ), row.names = FALSE)
}

Similarity b/w Protoarray & Huprot array results

summary.df <- data.frame(coefficient=coefs, huprot=0, protoarray=0, common=0, jaccard=0, names=0, stringsAsFactors = FALSE)
row.names(summary.df) <- coefs
for (coef in coefs){
  huprot.results <- read.table(file.path('..', 'results_common/huprot', 'significant_filtered', paste(coef,'csv', sep='.') ), header = T)
  protoarray.results <- read.table(file.path('..', 'results_common/protoarray', 'significant_filtered', paste(coef,'csv', sep='.') ), header = T)
  #new_name <- sapply(as.vector(protoarray.results$Name), clean_name)
  #protoarray.results$Name <- new_name
  total <- length(huprot.results$ID)+ length(protoarray.results$Name)
  common <- length(intersect(huprot.results$ID, protoarray.results$Name))
  summary.df[coef,]$huprot <-   length(huprot.results$ID)
  summary.df[coef,]$protoarray <-   length(protoarray.results$Name)
  summary.df[coef,]$common <- common  
  summary.df[coef,]$jaccard <- common/total  
  summary.df[coef,]$names <- toString(intersect(huprot.results$ID, protoarray.results$Name))
}
summary.df