#!/bioware/R-3.0.1/bin/Rscript library(baySeq) ########################## ### Use this for multi-threading # library(snow) cl <- makeCluster(4, type = "SOCK") # ### Otherwise use this #cl <- NULL ####### ######################### ### Load Data # Data should be target, length, sample_1, sample_2, ... sample_n # if no head change to FALSE # # check class of vectors with class() # BMAGING <- read.delim("5pairedExp.fpkm10.counts.txt", header = TRUE) ISOTIGS <- data.frame(BMAGING[,1]) ANNOTATION <- data.frame(BMAGING[,1:2]) LENGTHS <- as.matrix(BMAGING[,2]) # read all data # 1 = isotig, 2=length #COUNTSALL <- as.matrix(BMAGING[,-(1:2)]) # read some data # 3-4 = eggs, 5-6 = neonates, 7-8 = early, 9-10 mixed, 11-12 post COUNTS1 <-as.matrix(BMAGING[,3:6]) # eggs vs neonates COUNTS2 <-as.matrix(BMAGING[,5:8]) # neonates vs early COUNTS3 <-as.matrix(BMAGING[,7:10]) # early vs mixed COUNTS4 <-as.matrix(BMAGING[,9:12]) # mixed vs post also as.matrix(BMAGING[,-(1:4)]) COUNTS5 <-as.matrix(BMAGING[,c(3,4,7,8)]) # eggs vs early COUNTS6 <-as.matrix(BMAGING[,c(5,6,9,10)]) # neonates vs mixed COUNTS7 <-as.matrix(BMAGING[,c(7,8,11,12)]) # early vs post COUNTS8 <-as.matrix(BMAGING[,c(3,4,9,10)]) # eggs vs mixed COUNTS9 <-as.matrix(BMAGING[,c(5,6,11,12)]) # neonates vs post COUNTS10 <-as.matrix(BMAGING[,c(3,4,11,12)]) # eggs vs post #myLibsizeData <- getLibsizes(data = myCounts, estimationType = c("total")) SIZES1 <- getLibsizes(data = COUNTS1, estimationType = c("quantile")) SIZES2 <- getLibsizes(data = COUNTS2, estimationType = c("quantile")) SIZES3 <- getLibsizes(data = COUNTS3, estimationType = c("quantile")) SIZES4 <- getLibsizes(data = COUNTS4, estimationType = c("quantile")) SIZES5 <- getLibsizes(data = COUNTS5, estimationType = c("quantile")) SIZES6 <- getLibsizes(data = COUNTS6, estimationType = c("quantile")) SIZES7 <- getLibsizes(data = COUNTS7, estimationType = c("quantile")) SIZES8 <- getLibsizes(data = COUNTS8, estimationType = c("quantile")) SIZES9 <- getLibsizes(data = COUNTS9, estimationType = c("quantile")) SIZES10 <- getLibsizes(data = COUNTS10, estimationType = c("quantile")) ########################### ########################### REPS <- c(1,1,2,2) ########################### ######################## ## This is where you specifiy your hypotheses # 1,1,2,2,3,3,4,4,5,5 GROUPS <- list( NDE = c(1,1,1,1), DE = c(1,1,2,2), LIB = c(1,2,1,2) ) ######################## ######################## ### Build the countData objects # CD1 <- new("countData", seglens = LENGTHS, data = COUNTS1, replicates = REPS, libsizes = as.integer(SIZES1), groups = GROUPS, annotation = ANNOTATION) CD2 <- new("countData", seglens = LENGTHS, data = COUNTS2, replicates = REPS, libsizes = as.integer(SIZES2), groups = GROUPS, annotation = ANNOTATION) CD3 <- new("countData", seglens = LENGTHS, data = COUNTS3, replicates = REPS, libsizes = as.integer(SIZES3), groups = GROUPS, annotation = ANNOTATION) CD4 <- new("countData", seglens = LENGTHS, data = COUNTS4, replicates = REPS, libsizes = as.integer(SIZES4), groups = GROUPS, annotation = ANNOTATION) CD5 <- new("countData", seglens = LENGTHS, data = COUNTS5, replicates = REPS, libsizes = as.integer(SIZES5), groups = GROUPS, annotation = ANNOTATION) CD6 <- new("countData", seglens = LENGTHS, data = COUNTS6, replicates = REPS, libsizes = as.integer(SIZES6), groups = GROUPS, annotation = ANNOTATION) CD7 <- new("countData", seglens = LENGTHS, data = COUNTS7, replicates = REPS, libsizes = as.integer(SIZES7), groups = GROUPS, annotation = ANNOTATION) CD8 <- new("countData", seglens = LENGTHS, data = COUNTS8, replicates = REPS, libsizes = as.integer(SIZES8), groups = GROUPS, annotation = ANNOTATION) CD9 <- new("countData", seglens = LENGTHS, data = COUNTS9, replicates = REPS, libsizes = as.integer(SIZES9), groups = GROUPS, annotation = ANNOTATION) CD10 <- new("countData", seglens = LENGTHS, data = COUNTS10, replicates = REPS, libsizes = as.integer(SIZES10), groups = GROUPS, annotation = ANNOTATION) ######################## ######################## ### Negative Binomial approach ## calculate priors and put them in a new countData object # CDPriors.NB <- getPriors.NB(CD1, samplesize = 10000, estimation = "QL", cl = cl) ## calculate posterior likelihoods using CDPrios.NB # CDPost.NB <- getLikelihoods.NB(CDPriors.NB, pET = "BIC", cl = cl, nullData = TRUE) ######################### options(max.print=999999) options(width=800) sink("eggsVSneonates.bayseq.tab", append=FALSE, split=FALSE) topCounts(CDPost.NB, group = "DE", number = 40000) sink() save.image(file = "eggsVSneonates.bayseq.image") if(!is.null(cl)) stopCluster(cl) quit(save="no")