### R code from vignette source 'DTA.Rnw' ################################################### ### code chunk number 1: Loading library ################################################### library(DTA) ################################################### ### code chunk number 2: Loading Miller2011 ################################################### data(Miller2011) ################################################### ### code chunk number 3: Loading datamat ################################################### colnames(Sc.datamat)[1:6] rownames(Sc.datamat)[1:6] Sc.datamat[1:6,1:3] ################################################### ### code chunk number 4: Loading phenomat ################################################### head(Sc.phenomat) ################################################### ### code chunk number 5: Loading tnumber ################################################### head(Sc.tnumber) ################################################### ### code chunk number 6: Total expression histogram ################################################### Totals = Sc.datamat[,which(Sc.phenomat[,"fraction"]=="T")] Total = apply(log(Totals),1,median) plotsfkt = function(){ par(mar = c(5,4,4,2)+0.1+1) par(mai = c(1.1,1.1,1.3,0.7)) hist(Total,breaks = seq(0,ceiling(max(Total)),1/4), cex.main=1.5,cex.lab=1.25,cex.axis=1.125, main="Histogram of log(Total)", xlab="gene-wise median of total samples") hist(Total[Total >= 5],breaks = seq(0,ceiling(max(Total)),1/4), col = "#08306B",add = TRUE) } DTA.plot.it(filename = "histogram_cut_off",plotsfkt = plotsfkt, saveit = TRUE,notinR = TRUE) ################################################### ### code chunk number 7: Loading reliable ################################################### head(Sc.reliable) ################################################### ### code chunk number 8: Estimation procedure ################################################### res = DTA.estimate(Sc.phenomat,Sc.datamat,Sc.tnumber, ccl = 150,mRNAs = 60000,reliable = Sc.reliable, condition = "real_data",save.plots = TRUE,notinR = TRUE,folder = ".") ################################################### ### code chunk number 9: Entries of res ################################################### names(res) ################################################### ### code chunk number 10: Entries of res6 ################################################### names(res[["6"]]) ################################################### ### code chunk number 11: plot half-lives/sythesis rates ################################################### data(Sc.ribig.ensg) data(Sc.rpg.ensg) data(Sc.tf.ensg) data(Sc.stress.ensg) ################################################### ### code chunk number 12: plot half-lives/sythesis rates ################################################### plotsfkt = function(){ par(mar = c(5,4,4,2)+0.1+1) par(mai = c(1.1,1.1,1.3,0.7)) x=res[["6"]][["hl"]][c(Sc.reliable,Sc.rpg.ensg)] y=res[["6"]][["Rsr"]][c(Sc.reliable,Sc.rpg.ensg)] ellipsescatter(x,y, groups = list("Stress"=Sc.stress.ensg,"RiBi-genes"=Sc.ribig.ensg, "Rp-genes"=Sc.rpg.ensg,"TF"=Sc.tf.ensg), colors = c("darkred","darkgreen","darkblue","grey20"), xlim=c(0,150),ylim=c(0,600),xlab="half-lives",ylab="synthesis rates", cex.main=1.5,cex.lab=1.25,cex.axis=1.125,main="Ellipsescatter")} DTA.plot.it(filename = "ellipse_scatter",plotsfkt = plotsfkt, saveit = TRUE,notinR = TRUE) ################################################### ### code chunk number 13: Estimation procedure (no bias correction) ################################################### res.nobias = DTA.estimate(Sc.phenomat,Sc.datamat,Sc.tnumber, ccl = 150, mRNAs = 60000,reliable = Sc.reliable,save.plots = TRUE, notinR = TRUE,folder = ".",bicor = FALSE,condition="no_bias_correction") ################################################### ### code chunk number 14: Creating simulation object ################################################### sim.object = DTA.generate(timepoints=rep(c(6,12),2)) ################################################### ### code chunk number 15: Entries of sim.object ################################################### names(sim.object) ################################################### ### code chunk number 16: Estimation procedure for simulated data ################################################### res.sim = DTA.estimate(ratiomethod = "bias",save.plots = TRUE, notinR = TRUE,folder = ".",simulation = TRUE,sim.object = sim.object, condition = "simulation") ################################################### ### code chunk number 17: Estimation procedure for simulated data ################################################### res.sim = DTA.estimate(ratiomethod = "bias",save.plots = TRUE, notinR = TRUE,folder = ".",simulation = TRUE,sim.object = sim.object, bicor = FALSE,condition = "simulation_no_bias_correction") ################################################### ### code chunk number 18: Loading Miller2011dynamic ################################################### data(Miller2011dynamic) ################################################### ### code chunk number 19: Estimation procedure ################################################### timecourse.res = DTA.dynamic.estimate(Sc.phenomat.dynamic,Sc.datamat.dynamic, Sc.tnumber,ccl = 150,mRNAs = 60000,reliable = Sc.reliable.dynamic, LtoTratio = rep(0.1,7),save.plots = TRUE,notinR = TRUE,folder = ".", condition = "timecourse") ################################################### ### code chunk number 20: Entries of timecourse.res ################################################### names(timecourse.res) ################################################### ### code chunk number 21: Entries of timecourse.res1 ################################################### names(timecourse.res[["1"]]) ################################################### ### code chunk number 22: Creating timecourse simulation ################################################### nrgenes = 5000 truesynthesisrates = rf(nrgenes,5,5)*18 steady = rep(1,nrgenes) shock = 1/pmax(rnorm(nrgenes,mean = 8,sd = 4),1) induction = pmax(rnorm(nrgenes,mean = 8,sd = 4),1) changes.mat = cbind(steady,shock,shock*induction) mu.values.mat = changes.mat*truesynthesisrates mu.breaks.mat = cbind(rep(12,nrgenes),rep(18,nrgenes)) truehalflives = rf(nrgenes,15,15)*12 truelambdas = log(2)/truehalflives changes.mat = cbind(steady,shock,shock*induction,steady) lambda.values.mat = changes.mat*truelambdas lambda.breaks.mat = cbind(rep(12,nrgenes),rep(18,nrgenes),rep(27,nrgenes)) ################################################### ### code chunk number 23: Creating timecourse simulation object ################################################### timecourse.sim.object = DTA.dynamic.generate(duration = 36,lab.duration = 6, nrgenes = nrgenes,mu.values.mat = mu.values.mat,mu.breaks.mat = mu.breaks.mat, lambda.values.mat = lambda.values.mat,lambda.breaks.mat = lambda.breaks.mat) ################################################### ### code chunk number 24: Entries of timecourse.sim.object ################################################### names(timecourse.sim.object) ################################################### ### code chunk number 25: Estimation procedure for simulated timecourse data ################################################### timecourse.res.sim = DTA.dynamic.estimate(save.plots = TRUE,notinR = TRUE, simulation = TRUE,sim.object = timecourse.sim.object,ratiomethod = "tls", folder = ".",condition = "simulation_timecourse",check = FALSE) ################################################### ### code chunk number 26: Loading Doelken2008 (eval = FALSE) ################################################### ## data(Doelken2008) ################################################### ### code chunk number 27: Mm map tnumber (eval = FALSE) ################################################### ## Mm.tnumber = DTA.map.it(Mm.tnumber,Mm.enst2ensg) ################################################### ### code chunk number 28: Mm DTA estimation procedure (eval = FALSE) ################################################### ## res = DTA.estimate(Mm.phenomat,Mm.datamat,Mm.tnumber, ## ccl = 24*60,reliable = Mm.reliable,ratiomethod = "tls") ################################################### ### code chunk number 29: Hs map tnumber (eval = FALSE) ################################################### ## Hs.tnumber = DTA.map.it(Hs.tnumber,Hs.enst2ensg) ################################################### ### code chunk number 30: Hs DTA estimation procedure (eval = FALSE) ################################################### ## res = DTA.estimate(Hs.phenomat,Hs.datamat,Hs.tnumber, ## ccl = 50*60,reliable = Hs.reliable,ratiomethod = "tls", ## bicor = FALSE) ################################################### ### code chunk number 31: sessionInfo ################################################### toLatex(sessionInfo())