# load data
X <- read.table('all.reads.num')
SL <- read.table('all.SL.num')

# Reads/seqLen
r1SL <- X[,1]/SL
r2SL <- X[,2]/SL

# total read number/sample
s1 <- sum(r1SL)
s2 <- sum(r2SL)

# total reads normalizing factors
f2 <- s1/s2

# normalize by read number factors
r1N <- r1SL
r2N <- r2SL*f2

# join lists into a matrix
XN  <- cbind(r1N,r2N)

# check for non linear bias
rat12t <- log2(XN[,1]/XN[,2])
int <- apply(log2(XN),1,mean)
plot(int,rat12t)

# normalize using splines
require(affy)
XN_spl <- normalize.qspline(XN)

# check again for non linear bias
rat12_spl <- log2(XN_spl[,1]/XN_spl[,2])
int_spl <- apply(log2(XN_spl),1,mean)
plot(int_spl,rat12_spl)

sd_spl <- sd(t(XN_spl))

# std.devs
sd_spl  <- sd(t(XN_spl))

# write(XN_spl, file="reads.normalized.num", ncolumns=2)
# write(sd_spl, file="all.stddev.num", ncolumns=1)