### R code from vignette source 'HTSFilter.Rnw' ################################################### ### code chunk number 1: style-Sweave ################################################### BiocStyle::latex() ################################################### ### code chunk number 2: HTSFilter.Rnw:26-27 ################################################### options(width=60) ################################################### ### code chunk number 3: loadingPackages ################################################### library(HTSFilter) library(edgeR) library(DESeq2) data("sultan") hist(log(exprs(sultan)+1), col="grey", breaks=25, main="", xlab="Log(counts+1)") pData(sultan) dim(sultan) ################################################### ### code chunk number 4: matrix ################################################### mat <- exprs(sultan) conds <- as.character(pData(sultan)$cell.line) ## Only 25 tested thresholds to reduce computation time filter <- HTSFilter(mat, conds, s.min=1, s.max=200, s.len=25) mat <- filter$filteredData dim(mat) dim(filter$removedData) ################################################### ### code chunk number 5: refilter ################################################### par(mfrow = c(1,2), mar = c(4,4,2,2)) filter.2 <- HTSFilter(mat, conds, s.len=25) dim(filter.2$removedData) hist(log(filter.2$filteredData+1), col="grey", breaks=25, main="", xlab="Log(counts+1)") ################################################### ### code chunk number 6: dge ################################################### dge <- DGEList(counts=exprs(sultan), group=conds) dge <- calcNormFactors(dge) dge <- estimateDisp(dge) et <- exactTest(dge) et <- HTSFilter(et, DGEList=dge, s.len=25, plot=FALSE)$filteredData dim(et) class(et) topTags(et) ################################################### ### code chunk number 7: dgetoptags ################################################### topTags(et) ################################################### ### code chunk number 8: dge2 ################################################### design <- model.matrix(~conds) dge <- DGEList(counts=exprs(sultan), group=conds) dge <- calcNormFactors(dge) dge <- estimateDisp(dge, design) fit <- glmFit(dge,design) lrt <- glmLRT(fit,coef=2) lrt <- HTSFilter(lrt, DGEGLM=fit, s.len=25, plot=FALSE)$filteredData dim(lrt) class(lrt) ################################################### ### code chunk number 9: dge2toptags ################################################### topTags(lrt) ################################################### ### code chunk number 10: cds2 (eval = FALSE) ################################################### ## ## dds <- DESeqDataSetFromMatrix(countData = exprs(sultan), ## colData = data.frame(cell.line = conds), ## design = ~ cell.line) ## dds <- DESeq(dds) ## filter <- HTSFilter(dds, s.len=25, plot=FALSE)$filteredData ## class(filter) ## dim(filter) ## res <- results(filter, independentFiltering=FALSE) ## head(res) ## ################################################### ### code chunk number 11: ses2 (eval = FALSE) ################################################### ## library(EDASeq) ## ses <- newSeqExpressionSet(exprs(sultan), ## phenoData=pData(sultan)) ## ses.norm <- betweenLaneNormalization(ses, which="full") ################################################### ### code chunk number 12: ses3 (eval = FALSE) ################################################### ## filter <- HTSFilter(counts(ses.norm), conds, s.len=25, norm="none", ## plot=FALSE) ## head(filter$on) ## table(filter$on) ################################################### ### code chunk number 13: sessionInfo ################################################### sessionInfo()