vignettes/vignette.R

In carinademel/RNAlife: Quantifying kinetic parameters to determine the RNA life cycle

### R code from vignette source '/Users/Carina.Demel/Library/R/3.1/library/RNAlife/doc/vignette.Rnw'

###################################################
### code chunk number 1: options
###################################################
options(digits=3, prompt=" ", continue=" ")


###################################################
### code chunk number 2: DTAplot
###################################################
library("RNAlife")
res=dta.plot(mu=100,lambda=0.05,unlabeledcol="orange")


###################################################
### code chunk number 3: LoadingLibrary
###################################################
library("RNAlife")


###################################################
### code chunk number 4: LoadingData
###################################################
data(spikein.counts)
data(gene.counts)
spikein.counts


###################################################
### code chunk number 5: LoadingAdditionalApikeInformation
###################################################
data(spikein.labeling)
data(spikein.lengths)
spikeins = rownames(spikein.counts)


###################################################
### code chunk number 6: LoadingAdditionalSampleInformation
###################################################
data(samples)
samples


###################################################
### code chunk number 7: DataframeGeneration
###################################################
sample.names = apply(samples,1,paste,collapse=" ")
mat = spikein.dataframe(counts=spikein.counts,
      spikeins=spikeins, 
  	  spikein.lengths=spikein.lengths, 
		  spikein.labeling=spikein.labeling, 
		  samples=sample.names,
		  conditions.labeling=samples$conditions.labeling)
summary(mat)
head(mat, n=10)


###################################################
### code chunk number 8: GLM (eval = FALSE)
###################################################
##   glm.nb(counts ~ offset(log.length)+sample+ccc+spikein, data=mat, link="log")


###################################################
### code chunk number 9: Normalization
###################################################
norm.res = spikein.normalization(spikein.counts,
		spikeins, 
		spikein.lengths, 
		spikein.labeling, 
		samples=colnames(spikein.counts), 
		samples$conditions.labeling)
summary(norm.res)
norm.res$sequencing.depth
norm.res$cross.contamination


###################################################
### code chunk number 10: SpikeinNormalization
###################################################
pchs = c(1,4,8,11, 21:24); names(pchs) = unique(mat$sample);
spike.colors = rainbow(n=length(unique(mat$spike))); names(spike.colors)=as.character(unique(mat$spike))#[order(spikes.labeling)] #col pal is ordered, spikes are now sorted according to L and U
mat$fitted.counts = norm.res$fitted.counts
plot(fitted.counts ~ counts, log="xy", col=spike.colors[as.character(mat$spike)], data=mat, pch=pchs[as.character(mat$sample)], bg=spike.colors[as.character(mat$spike)], xlab="Measured Counts", ylab="Fitted Counts", main="Correlation of fitted and measured Spike-in Counts")
names(spike.colors) = sapply(names(spike.colors), function(i){if(i %in% c("Spike2", "Spike4", "Spike8")){paste(i, " (4sU)")}else i})
legend('topleft', legend=names(spike.colors), col=spike.colors, pch=19) #fill=colors makes rectangles
legend('bottomright',legend=names(pchs), pch=pchs, pt.bg="black") #fill=colors makes rectangles
abline(0,1)


###################################################
### code chunk number 11: Dispersion
###################################################
data(gene.counts)
dispersion = estimateGeneDispersion(gene.counts, samples$conditions.labeling, samples$conditions)
head(dispersion)


###################################################
### code chunk number 12: Estimation
###################################################
seq.depths = norm.res$sequencing.depth
cross.cont = norm.res$cross.contamination
lengths = rnbinom(nrow(gene.counts), mu=1000, size=10)
fittingres = estimate.rates(gene.counts, dispersion, lengths, 
  samples$conditions, samples$conditions.labeling, samples$replicates, cross.cont, seq.depths,
  N=1, consider.replicates=TRUE, lab.time=c(5,5), alpha.initial=1, beta.initial=1, gene.indices=1:12)
names(fittingres)
# Synthesis rates
head(fittingres$mu)
# Degradation rates
head(fittingres$lambda)


###################################################
### code chunk number 13: GeneCounts
###################################################
plot(as.vector(gene.counts[,grep("L",colnames(gene.counts))]), as.vector(fittingres$exp.counts.L), xlab="Observed Genes Counts", ylab="Estimated Gene Counts",log="xy", pch=17, col="red")
points(as.vector(gene.counts[,grep("T",colnames(gene.counts))]), as.vector(fittingres$exp.counts.T),pch=15, col="blue")
legend("topleft", c("Labeled Gene Counts", "Total Gene Counts"), col=c("red", "blue"), pch=c(17, 15))
abline(0,1)


###################################################
### code chunk number 14: sessionInfo
###################################################
sessionInfo()