inst/doc/tweeDEseq.R

In tweeDEseq: RNA-seq data analysis using the Poisson-Tweedie family of distributions

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

###################################################
### code chunk number 1: tweeDEseq.Rnw:57-59
###################################################
library(tweeDEseq)
library(tweeDEseqCountData)


###################################################
### code chunk number 2: tweeDEseq.Rnw:63-70
###################################################
data(pickrell)
countsNigerian <- exprs(pickrell.eset)
dim(countsNigerian)
countsNigerian[1:5, 1:5]
genderNigerian <- pData(pickrell.eset)[,"gender"]
head(genderNigerian)
table(genderNigerian)


###################################################
### code chunk number 3: tweeDEseq.Rnw:80-81 (eval = FALSE)
###################################################
## countsNigerianNorm <- normalizeCounts(countsNigerian, genderNigerian)


###################################################
### code chunk number 4: tweeDEseq.Rnw:83-85
###################################################
data(pickrellNorm)
countsNigerianNorm <- exprs(pickrellNorm.eset)


###################################################
### code chunk number 5: tweeDEseq.Rnw:87-88
###################################################
dim(countsNigerianNorm)


###################################################
### code chunk number 6: tweeDEseq.Rnw:92-94
###################################################
countsNigerianNorm <- filterCounts(countsNigerianNorm)
dim(countsNigerianNorm)


###################################################
### code chunk number 7: tweeDEseq.Rnw:114-118
###################################################
set.seed(123)
y <- rnbinom(1000, mu=8, size=1/0.2)
thetahat <- mlePoissonTweedie(y)
getParam(thetahat)


###################################################
### code chunk number 8: tweeDEseq.Rnw:126-127
###################################################
testShapePT(thetahat, a=0)


###################################################
### code chunk number 9: tweeDEseq.Rnw:131-137
###################################################
data(genderGenes)
data(hkGenes)

length(XiEgenes)
length(msYgenes)
length(hkGenes)


###################################################
### code chunk number 10: tweeDEseq.Rnw:141-147
###################################################
set.seed(123)
someHKgenes <- sample(hkGenes, size=25)
geneSubset <- unique(c("ENSG00000070031",
                       intersect(rownames(countsNigerianNorm),
                                 c(someHKgenes, msYgenes, XiEgenes))))
length(geneSubset)


###################################################
### code chunk number 11: tweeDEseq.Rnw:151-152
###################################################
chi2gof <- gofTest(countsNigerianNorm[geneSubset, ], a=0)


###################################################
### code chunk number 12: gof
###################################################
par(mfrow=c(1,2), mar=c(4, 5, 3, 4))
qq <- qqchisq(chi2gof, main="Chi2 Q-Q Plot", ylim = c(0, 60))
qq <- qqchisq(chi2gof, normal=TRUE, ylim = c(-5, 7))


###################################################
### code chunk number 13: secretin
###################################################
par(mfrow=c(1,2), mar=c(4, 5, 3, 2))
xf <- unlist(countsNigerianNorm["ENSG00000070031", genderNigerian=="female"])
compareCountDist(xf, main="Female samples")
xm <- unlist(countsNigerianNorm["ENSG00000070031", genderNigerian=="male"])
compareCountDist(xm, main="Male samples")


###################################################
### code chunk number 14: tweeDEseq.Rnw:192-194
###################################################
sort(xf)
sort(xm)


###################################################
### code chunk number 15: tweeDEseq.Rnw:198-201
###################################################
xf[which.max(xf)]
2^{log2(mean(xf))-log2(mean(xm))}
2^{log2(mean(xf[-which.max(xf)])) - log2(mean(xm))}


###################################################
### code chunk number 16: tweeDEseq.Rnw:213-214
###################################################
resPT <- tweeDE(countsNigerianNorm[geneSubset, ], group = genderNigerian)


###################################################
### code chunk number 17: tweeDEseq.Rnw:218-219
###################################################
print(resPT)


###################################################
### code chunk number 18: tweeDEseq.Rnw:225-227
###################################################
deGenes <- print(resPT, n=Inf, log2fc=log2(1.5), pval.adjust=0.05, print=FALSE)
dim(deGenes)


###################################################
### code chunk number 19: tweeDEseq.Rnw:231-234
###################################################
data(annotEnsembl63)
head(annotEnsembl63)
deGenes <- merge(deGenes, annotEnsembl63, by="row.names", sort=FALSE)


###################################################
### code chunk number 20: tabDEgenes
###################################################
library(xtable)
deGenes$Description <- gsub(" \\[.+$", "", deGenes$Description)
xtab <- xtable(deGenes[, c("Symbol", "Chr", "log2fc", "pval.adjust", "Description")],
               caption="Differentially expressed genes between female and male Nigerian individuals found by tweeDEseq.",
               label="tab:deGenes", align="|l|c|c|r|l|p{7cm}|", digits=c(0, 0, 0, 2, -2, 0))
print(xtab, floating=TRUE, sanitize.text.function=function(x) x, caption.placement="top", include.rownames=FALSE)


###################################################
### code chunk number 21: tweeDEseq.Rnw:253-255
###################################################
deGenes <- rownames(print(resPT, n=Inf, log2fc=log2(1.5), pval.adjust=0.05, print=FALSE))
length(deGenes)


###################################################
### code chunk number 22: MAnVplots
###################################################
hl <- list(list(genes=msYgenes, pch=21, col="blue", bg="blue", cex=0.7),
           list(genes=XiEgenes, pch=21, col="skyblue", bg="skyblue", cex=0.7),
           list(genes=deGenes, pch=1, col="red", lwd=2, cex=1.5)
          )
par(mfrow=c(1,2), mar=c(4, 5, 3, 2))
MAplot(resPT, cex=0.7, log2fc.cutoff=log2(1.5), highlight=hl, main="MA-plot")
Vplot(resPT, cex=0.7, highlight=hl, log2fc.cutoff=log2(1.5),
      pval.adjust.cutoff=0.05, main="Volcano plot")


###################################################
### code chunk number 23: tweeDEseq.Rnw:294-303
###################################################
genderGenes <- c(msYgenes[msYgenes %in% rownames(resPT)],
                 XiEgenes[XiEgenes %in% rownames(resPT)])
N <- nrow(resPT)
m <- length(genderGenes)
n <- length(deGenes)
k <- length(intersect(deGenes, genderGenes))
t <- array(c(k,n-k,m-k,N+k-n-m), dim=c(2,2), dimnames=list(SEX=c("in","out"),DE=c("yes","no")))
t
fisher.test(t, alternative="greater")


###################################################
### code chunk number 24: tweeDEseq.Rnw:310-313
###################################################
mod <- glmPT(countsNigerianNorm["ENSG00000070031",] ~ genderNigerian)
mod
summary(mod)


###################################################
### code chunk number 25: tweeDEseq.Rnw:318-319
###################################################
anova(mod)


###################################################
### code chunk number 26: tweeDEseq.Rnw:324-325
###################################################
resPTglm <- tweeDEglm( ~ genderNigerian, counts = countsNigerianNorm[geneSubset,])


###################################################
### code chunk number 27: tweeDEseq.Rnw:330-331
###################################################
head(resPTglm[sort(resPTglm$pval.adjust, index.return = TRUE)$ix,])


###################################################
### code chunk number 28: tweeDEseq.Rnw:342-344 (eval = FALSE)
###################################################
## tweeDEglm(~ genderNigerian, counts = countsNigerianNorm[geneSubset,],
##  offset = cqn.subset$offset)


###################################################
### code chunk number 29: tweeDEseq.Rnw:349-350
###################################################
sessionInfo()