lambdahatij branch test cod

lambdahatij_test.R

bamfiles <- "ERR188436.Aligned.sortedByCoord.out.bam"
names(bamfiles) <- "ERR188436"
load("fitpar_gc_str.rda")
load_all("~/proj/alpine/")
library(GenomicAlignments)
library(GenomicFeatures)
library(BSgenome.Hsapiens.UCSC.hg19)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

ebt <- exonsBy(txdb, by="tx")
models <- list("GC_str"=list(formula="count ~ ns(gc,knots=gc.knots,Boundary.knots=gc.bk) +           
                                      ns(relpos,knots=relpos.knots,Boundary.knots=relpos.bk) +       
                                      GC40.80 + GC40.90 + GC20.80 + GC20.90",
                 offset="fraglen"))
lib.sizes <- c("ERR188436"=1e6)

txdf <- select(txdb, keys(txdb, "GENEID")[1:400], "TXID", "GENEID")
table(txdf$GENEID)[table(txdf$GENEID) == 2]
two.iso <- names(table(txdf$GENEID)[table(txdf$GENEID) == 2])

ebt.count <- ebt[txdf$TXID[txdf$GENEID %in% two.iso]]
ebt.count <- keepStandardChromosomes(ebt.count)
count <- roughSummarizeOverlaps(ebt.count, bamfiles)
cbind(as.numeric(count), as.numeric(names(ebt.count)))

ebt.sub <- ebt[c("45207","45208")]
ebt.sub <- keepStandardChromosomes(ebt.sub)
plotGRL(ebt.sub)
roughSummarizeOverlaps(ebt.sub, bamfiles)

system.time({
              res <- estimateTheta(ebt.sub, bamfiles, fitpar, Hsapiens, models,
                     readlength=75, minsize=100, maxsize=250, subset=TRUE,
                     zerotopos=20, niter=100, lib.sizes=lib.sizes, optim=FALSE)
            })