mjktfw

microarray.analysis

Nov 10th, 2013
126
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. library(limma)
  2. library(vsn)
  3.  
  4. # DATA PREPROCESSING
  5. # Download RAW data from GEO experiment
  6. setInternet2(use=FALSE)
  7. download.file('http://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE34571&format=file', 'GSE34571_RAW.tar', mode='wb')
  8. untar('GSE34571_RAW.tar', exdir='gpr')
  9. setwd('gpr')
  10.  
  11. # Set function to assign weights
  12. f <- function(x) {ok.flags <- x$Flags > (-99)
  13.                   ok.snr <- x[,'SNR 532'] > 3
  14.                   as.numeric(ok.snr & ok.flags)
  15. }
  16.  
  17. # Read single-channel .gpr files as double-channel data and perform bg correction  
  18. RG <- read.maimages(dir(pattern='*.gpr.gz'), source="genepix",
  19.                     columns=list(R="F532 Mean",G="F532 Mean",Rb="B532",Gb="B532"),
  20.                     wt.fun=f)
  21. RG.b = backgroundCorrect(RG=RG, method='normexp')
  22.  
  23. # Delete fake second channel
  24. rownames(RG.b$R) <- RG.b$genes$Name
  25. RG.b$G <- NULL
  26. RG.b$Gb <- NULL
  27.  
  28. # Specify contrast and design matrices
  29. design <- model.matrix(~ 0+factor(c(1,1,1,1,2,2,2,2)))
  30. colnames(design) <- c("group1", "group2")
  31. contrast.matrix <- makeContrasts(group1-group2, levels=design)
  32.  
  33. # Peform normalisation and stats
  34. norm.vsn <- normalizeVSN(RG.b$R)
  35. fit.vsn <- lmFit(norm.vsn, design)
  36. fit2.vsn <- contrasts.fit(fit.vsn, contrast.matrix)
  37. fit3.vsn <- eBayes(fit2.vsn)
  38. limma::volcanoplot(fit3.vsn)
  39.  
  40.  
  41. RG$printer$ngrid.r = 1
  42. imageplot(log2(RG.b$R[,2]), RG$printer)
  43.  
  44. # DIAGNOSTIC PLOTS
  45.  
  46. # maplot.1
  47. par(mfrow=c(3,2))
  48. limma::plotMA(log2(RG$R[,c(1:4)]))
  49. limma::plotMA(log2(RG$R[,c(5:8)]))
  50.  
  51. limma::plotMA(log2(RG.b$R[,c(1:4)]))
  52. limma::plotMA(log2(RG.b$R[,c(5:8)]))
  53.  
  54. limma::plotMA(norm.vsn[,c(1:4)])
  55. limma::plotMA(norm.vsn[,c(5:8)])
  56.  
  57. # maplot.2.final
  58. par(mfrow=c(1,1))
  59. limma::plotMA(fit3.vsn)
  60.  
  61. # maplot.3.filtered
  62. par(mfrow=c(2,1))
  63. limma::plotMA(log2(RG.b$R[((RG.b$weights[,1] == 1) & (RG.b$weights[,2] == 1)),c(1:2)]))
  64. limma::plotMA(log2(RG.b$R[,c(1:2)]))
  65.  
  66. boxplot(RG$R)
  67. boxplot(norm.vsn)
  68.  
  69. # plotdensities
  70. plotDensities(RG, singlechannels=c(1:8))
  71. plotDensities(RG.b, singlechannels=c(1:8))
RAW Paste Data