gistfile1.r

library(lumi)
library(arrayQualityMetrics)

setwd('D:/Storage/My Dropbox/Shared/Mike_LungCancer/')
#setwd('~/Dropbox/Shared/Mike_LungCancer/')
data <- lumiR.batch(c('Batch1_EA10095_20101202_FinalReportNonNormNoBack.txt',
                      'Batch2_EA10095_20110729_FinalReportNonNormNoBack_clean.txt'),
                    sampleInfoFile='SampleDesc.txt')

png('density.png', width=1200, height=800)
density(data)
dev.off()

# Basic normalization for QC
data.T <- lumiT(data, method='log2',
                ifPlot=T, verbose=T)
data.N <- lumiN(data.T)
arrayQualityMetrics(data.N,
                    outdir='log2_report')

png('before_after.png', width=1600, height=2200)
par(mfrow=c(2,1))
boxplot(exprs(data.T))
boxplot(exprs(data.N))
dev.off()

# Trying Levi's ffpe package with the unnormalized data
library(ffpe)

# 25-75th percentile of expression intensities, sorted by 
# smallest to largest inter-quartile range
png('sortedIQR.png')
sortedIqrPlot(data, dolog2=T)
dev.off()

# Use IQR as a QC metric vs all chips
png('iqrMetrics.png')
par(mfrow=c(2,1))
sampleQC(data,
         logtransform=F,
         xaxis='index',
         cor.to='pseudochip',
         QCmeasure='IQR')

# Again, but using +/- 5 arrays
qc <- sampleQC(data,
         logtransform=F,
         xaxis='index',
         cor.to='similar',
         goby=5,
         QCmeasure='IQR')
dev.off()

# Why am I getting different results for notindex?
png('iqrMetrics_notIndex.png')
sampleQC(data,
         logtransform=F,
         xaxis='notindex',
         cor.to='pseudochip',
         QCmeasure='IQR')
dev.off()

# Carry on with the 11 array sliding window; can filter more 
# stringently on the probe level to compensate
data.filtered <- data[, !qc$rejectQC]
data.filtered.T <- lumiT(data.filtered, method='log2')
data.filtered.N <- lumiN(data.filtered.T)

# No technical replicates; not much that can be done in checking
# probe variance. Just filter for basic variance
probe.var <- apply(exprs(data.filtered.N), 1, var)
data.variance <- data.filtered.N[probe.var > median(probe.var),]