Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- library(limma)
- library(vsn)
- # DATA PREPROCESSING
- # Download RAW data from GEO experiment
- setInternet2(use=FALSE)
- download.file('http://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE34571&format=file', 'GSE34571_RAW.tar', mode='wb')
- untar('GSE34571_RAW.tar', exdir='gpr')
- setwd('gpr')
- # Set function to assign weights
- f <- function(x) {ok.flags <- x$Flags > (-99)
- ok.snr <- x[,'SNR 532'] > 3
- as.numeric(ok.snr & ok.flags)
- }
- # Read single-channel .gpr files as double-channel data and perform bg correction
- RG <- read.maimages(dir(pattern='*.gpr.gz'), source="genepix",
- columns=list(R="F532 Mean",G="F532 Mean",Rb="B532",Gb="B532"),
- wt.fun=f)
- RG.b = backgroundCorrect(RG=RG, method='normexp')
- # Delete fake second channel
- rownames(RG.b$R) <- RG.b$genes$Name
- RG.b$G <- NULL
- RG.b$Gb <- NULL
- # Specify contrast and design matrices
- design <- model.matrix(~ 0+factor(c(1,1,1,1,2,2,2,2)))
- colnames(design) <- c("group1", "group2")
- contrast.matrix <- makeContrasts(group1-group2, levels=design)
- # Peform normalisation and stats
- norm.vsn <- normalizeVSN(RG.b$R)
- fit.vsn <- lmFit(norm.vsn, design)
- fit2.vsn <- contrasts.fit(fit.vsn, contrast.matrix)
- fit3.vsn <- eBayes(fit2.vsn)
- limma::volcanoplot(fit3.vsn)
- RG$printer$ngrid.r = 1
- imageplot(log2(RG.b$R[,2]), RG$printer)
- # DIAGNOSTIC PLOTS
- # maplot.1
- par(mfrow=c(3,2))
- limma::plotMA(log2(RG$R[,c(1:4)]))
- limma::plotMA(log2(RG$R[,c(5:8)]))
- limma::plotMA(log2(RG.b$R[,c(1:4)]))
- limma::plotMA(log2(RG.b$R[,c(5:8)]))
- limma::plotMA(norm.vsn[,c(1:4)])
- limma::plotMA(norm.vsn[,c(5:8)])
- # maplot.2.final
- par(mfrow=c(1,1))
- limma::plotMA(fit3.vsn)
- # maplot.3.filtered
- par(mfrow=c(2,1))
- limma::plotMA(log2(RG.b$R[((RG.b$weights[,1] == 1) & (RG.b$weights[,2] == 1)),c(1:2)]))
- limma::plotMA(log2(RG.b$R[,c(1:2)]))
- boxplot(RG$R)
- boxplot(norm.vsn)
- # plotdensities
- plotDensities(RG, singlechannels=c(1:8))
- plotDensities(RG.b, singlechannels=c(1:8))
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement