inst/doc/TargetScore.R

In TargetScore: TargetScore: Infer microRNA targets using microRNA-overexpression data and sequence information

### R code from vignette source 'TargetScore.Rnw'

###################################################
### code chunk number 1: TargetScore
###################################################
library(TargetScore)


###################################################
### code chunk number 2: helppage (eval = FALSE)
###################################################
## ?TargetScore


###################################################
### code chunk number 3: toy
###################################################
trmt <- c(rnorm(10,mean=0.01), rnorm(1000,mean=1), rnorm(90,mean=2)) + 1e3

ctrl <- c(rnorm(1100,mean=1)) + 1e3

logFC <- log2(trmt) - log2(ctrl)

# 8 out of the 10 down-reg genes have prominent seq score A
seqScoreA <- c(rnorm(8,mean=-2), rnorm(1092,mean=0))

# 10 down-reg genes plus 10 more genes have prominent seq score B
seqScoreB <- c(rnorm(20,mean=-2), rnorm(1080,mean=0))

seqScores <- cbind(seqScoreA, seqScoreB)             

p.targetScore <- targetScore(logFC, seqScores, tol=1e-3)

# plot relation between targetScore and input variables
pairs(cbind(p.targetScore, logFC, seqScoreA, seqScoreB))


###################################################
### code chunk number 4: realtest
###################################################
extdata.dir <- system.file("extdata", package="TargetScore") 

# load test data
load(list.files(extdata.dir, "\\.RData$", full.names=TRUE))

myTargetScores <- lapply(mytestdata, function(x) targetScore(logFC=x$logFC, seqScores=x[,c(3,4)], tol=1e3, maxiter=200))

names(myTargetScores) <- names(mytestdata)

valid <- lapply(names(myTargetScores), function(x) table((myTargetScores[[x]] > 0.3), mytestdata[[x]]$validated))

names(valid) <- names(mytestdata)

# row pred and col validated targets
valid


###################################################
### code chunk number 5: getTargetScores demo (eval = FALSE)
###################################################
## library(TargetScoreData)
## library(gplots)
## 
## myTargetScores <- getTargetScores("hsa-miR-1", tol=1e-3, maxiter=200)
## 
## table((myTargetScores$targetScore > 0.1), myTargetScores$validated) # a very lenient cutoff
## 
## # obtain all of targetScore for all of the 112 miRNA (takes time)
## 
## logFC.imputed <- get_precomputed_logFC()
## 
## mirIDs <- unique(colnames(logFC.imputed))
## 
## # targetScoreMatrix <- mclapply(mirIDs, getTargetScores)
## 
## # names(targetScoreMatrix) <- mirIDs


###################################################
### code chunk number 6: getTargetScores limma demo (eval = FALSE)
###################################################
## # Demo using limma
## # download GEO data for hsa-miR-205
## library(limma)
## library(Biobase)
## library(GEOquery)
## library(gplots)
## 
## gset <- getGEO("GSE11701", GSEMatrix =TRUE, AnnotGPL=TRUE)
## 
## if (length(gset) > 1) idx <- grep("GPL6104", attr(gset, "names")) else idx <- 1
## 
## gset <- gset[[idx]]
## 
## geneInfo <- fData(gset)
## 
## x <- normalizeVSN(exprs(gset))
## 
## pData(gset)$title
## 
## design <- model.matrix(~0+factor(c(1,1,1,1,2,2,2,2)))
## 
## colnames(design) <- c("preNeg_control", "miR_205_transfect")
## 
## fit <- lmFit(x, design)
## 
## #contrast
## contrast.matrix <- makeContrasts(miR_205_transfect-preNeg_control,levels=design)
## 
## fit2 <- contrasts.fit(fit, contrast.matrix)
## 
## fit2 <- eBayes(fit2)
## 
## limma_stats <- topTable(fit2, coef=1, number=nrow(fit2))
## 
## limma_stats$gene_symbol <- geneInfo[match(limma_stats$ID, geneInfo$ID), "Gene symbol"]
## 
## logFC <- as.matrix(limma_stats$logFC)
## 
## rownames(logFC) <- limma_stats$gene_symbol
## 
## # targetScore
## myTargetScores <- getTargetScores("hsa-miR-205", logFC, tol=1e-3, maxiter=200)


###################################################
### code chunk number 7: ROC (eval = FALSE)
###################################################
## library(ggplot2)
## library(scales)
## library(ROCR)
## 
## # ROC
## roceval <- function(myscores, labels_true, method) {
##   	
## 	pred <- prediction(myscores, labels_true)
## 		
## 	perf <- performance( pred, "tpr", "fpr" )
## 
## 	auc <- unlist(slot(performance( pred, "auc" ), "y.values"))
## 
## 	fpr <- unlist(slot(perf, "x.values"))
## 	
## 	tpr <- unlist(slot(perf, "y.values"))
##   
## 	cutoffval <- unlist(slot(perf, "alpha.values"))
## 
## 	rocdf <- data.frame(x= fpr, y=tpr, cutoff=cutoffval, auc=auc, method=method,
## 		methodScore=sprintf("%s (%s)", method, percent(auc)), curveType="ROC")
## 
## 	return(rocdf)	
## }
## 
## limma_stats$p.val <- -log10(limma_stats$P.Value)
## 
## limma_stats$p.val[logFC > 0] <- 0
## 
## 
## myeval <- rbind(
##   roceval(myTargetScores$targetScore, myTargetScores$validated, "TargetScore"),
## 	roceval(limma_stats$p.val, myTargetScores$validated, "Limma"))
## 
## 
## ggroc <- ggplot(myeval, aes(x=x, y=y, color=methodScore)) +
## 
## 	geom_line(aes(linetype=methodScore), size=0.7) +
## 	
## 	scale_x_continuous("False positive rate", labels=percent) +
## 	
## 	scale_y_continuous("True positive rate", labels=percent) +
## 	
## 	theme(legend.title= element_blank())
## 
## print(ggroc)
## 


###################################################
### code chunk number 8: getTargetScores deseq demo (eval = FALSE)
###################################################
## library(DESeq)
## 
## cds <- makeExampleCountDataSet()
## 
## cds <- estimateSizeFactors( cds )
## 
## cds <- estimateDispersions( cds )
## 
## deseq_stats <- nbinomTest( cds, "A", "B" )
## 
## logFC <- deseq_stats$log2FoldChange[1:100]
## 
## # random sequence score
## seqScoreA <- rnorm(length(logFC))
## 
## seqScoreB <- rnorm(length(logFC))
## 
## seqScores <- cbind(seqScoreA, seqScoreB)
## 
## p.targetScore <- targetScore(logFC, seqScores, tol=1e-3)


###################################################
### code chunk number 9: sessi
###################################################
sessionInfo()