Skip to content

Commit 6969f6f

Browse files
committed
Make changes suggested in JOSS review.
1 parent 78cd7b5 commit 6969f6f

File tree

12 files changed

+276
-13
lines changed

12 files changed

+276
-13
lines changed

.gitignore

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
.*
22
shiny_bookmarks/
33
paper.pdf
4+
varistran.sublime-*
45

56
# History files
67
.Rhistory

DESCRIPTION

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -10,19 +10,21 @@ Description: Transform RNA-Seq count data so that variance due to biological
1010
Authors@R: person("Paul", "Harrison", email = "[email protected]", role = c("aut", "cre"))
1111
Maintainer: Paul Harrison <[email protected]>
1212
URL: https://github.com/MonashBioinformaticsPlatform/varistran
13-
Version: 1.0.1
13+
Version: 1.0.2
1414
License: LGPL-2.1 | file LICENSE
1515
Depends:
1616
grid
1717
Imports:
1818
MASS,
1919
ggplot2,
20-
ggdendro,
2120
shiny,
2221
miniUI,
2322
seriation,
2423
gridBase
2524
Suggests:
2625
edgeR,
27-
limma
26+
limma,
27+
DESeq2,
28+
biomaRt,
29+
NBPSeq
2830
RoxygenNote: 6.0.1

NEWS

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
2+
1.0.2
3+
=====
4+
5+
Added n parameter to plot_heatmap, to show only the top n rows by span of expression levels.
6+
7+
Remove dependency on ggdenro, which is not available in R 3.4.1.
8+
9+
Added packages needed for testing to suggested dependencies:
10+
11+
* biomaRt
12+
* DESeq2
13+
* NBPSeq
14+

R/dendrogram.R

Lines changed: 182 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,182 @@
1+
2+
# This code for drawing dendrograms has been adapted from ggdendro package.
3+
4+
dendro_data <- function(model, type=c("rectangle","triangle"), ...) {
5+
dhc <- as.dendrogram(model)
6+
hcdata <- dendrogram_data(dhc, type=type, ...)
7+
list(
8+
segments = hcdata$segments,
9+
labels = hcdata$labels,
10+
class="hclust"
11+
)
12+
}
13+
14+
dendrogram_data <- function (x, type = c("rectangle", "triangle"), ...){
15+
16+
# Initialise variables that used to be in parameter list
17+
leaflab <- "perpendicular"
18+
center <- FALSE
19+
xlab <- ""
20+
ylab <- ""
21+
horiz <- FALSE
22+
#frame.plot <- FALSE
23+
xaxt <- "n"
24+
yaxt <- "s"
25+
dLeaf <- NULL
26+
edge.root <- is.leaf(x) || !is.null(attr(x, "edgetext"))
27+
28+
type <- match.arg(type)
29+
#leaflab <- match.arg(leaflab)
30+
hgt <- attr(x, "height")
31+
if (edge.root && is.logical(edge.root))
32+
edge.root <- 0.0625 * if (is.leaf(x)) 1 else hgt
33+
mem.x <- .memberDend(x)
34+
yTop <- hgt + edge.root
35+
if (center) {
36+
x1 <- 0.5
37+
x2 <- mem.x + 0.5
38+
}
39+
else {
40+
x1 <- 1
41+
x2 <- mem.x
42+
}
43+
xl. <- c(x1 - 1/2, x2 + 1/2)
44+
yl. <- c(0, yTop)
45+
if (edge.root) {
46+
if (!is.null(et <- attr(x, "edgetext"))) {
47+
my <- mean(hgt, yTop)
48+
}
49+
}
50+
51+
gg.plotNode <- function (x1, x2, subtree, type, center, leaflab, dLeaf, horiz=FALSE, ddsegments=NULL, ddlabels=NULL) {
52+
inner <- !is.leaf(subtree) && x1 != x2
53+
yTop <- attr(subtree, "height")
54+
bx <- plotNodeLimit(x1, x2, subtree, center)
55+
xTop <- bx$x
56+
hasP <- !is.null(nPar <- attr(subtree, "nodePar"))
57+
Xtract <- function(nam, L, default, indx) rep(if (nam %in%
58+
names(L)) L[[nam]] else default, length.out = indx)[indx]
59+
asTxt <- function(x){
60+
if (is.character(x) || is.expression(x))
61+
x else
62+
if (is.null(x)) "" else as.character(x)
63+
}
64+
i <- if (inner || hasP)
65+
1
66+
else 2
67+
lab.cex <- 1
68+
if (is.leaf(subtree)) {
69+
if (leaflab == "perpendicular") {
70+
Y <- yTop - dLeaf * lab.cex
71+
X <- xTop
72+
srt <- 90
73+
adj <- 1
74+
nodeText <- asTxt(attr(subtree, "label"))
75+
ddlabels <- rbind(ddlabels, data.frame(x=X, y=0, text=nodeText))
76+
}
77+
}
78+
else if (inner) {
79+
segmentsHV <- function(x0, y0, x1, y1) {
80+
# *************************
81+
data.frame(x0, y0, x1, y1) #AdV
82+
}
83+
for (k in seq_along(subtree)) {
84+
child <- subtree[[k]]
85+
yBot <- attr(child, "height")
86+
if (getOption("verbose"))
87+
cat("ch.", k, "@ h=", yBot, "; ")
88+
if (is.null(yBot))
89+
yBot <- 0
90+
xBot <- if (center)
91+
mean(bx$limit[k:(k + 1)])
92+
else bx$limit[k] + .midDend(child)
93+
if (type == "triangle") {
94+
ddsegments <- rbind(ddsegments, segmentsHV(xTop, yTop, xBot, yBot))
95+
}
96+
else {
97+
ddsegments <- rbind(ddsegments, segmentsHV(xTop, yTop, xBot, yTop))
98+
ddsegments <- rbind(ddsegments, segmentsHV(xBot, yTop, xBot, yBot))
99+
}
100+
vln <- NULL
101+
if (!is.null(attr(child, "edgetext"))) {
102+
edgeText <- asTxt(attr(child, "edgetext"))
103+
if (!is.null(vln)) {
104+
mx <- if (type == "triangle")
105+
(xTop + xBot + ((xTop - xBot)/(yTop - yBot)) * vln)/2
106+
else xBot
107+
my <- (yTop + yBot + 2 * vln)/2
108+
}
109+
else {
110+
mx <- if (type == "triangle")
111+
(xTop + xBot)/2
112+
else xBot
113+
my <- (yTop + yBot)/2
114+
}
115+
}
116+
plotNode_result <- gg.plotNode(bx$limit[k], bx$limit[k + 1], subtree = child,
117+
type, center, leaflab, dLeaf, horiz, ddsegments, ddlabels)
118+
ddsegments <- plotNode_result$segments
119+
ddlabels <- plotNode_result$labels
120+
}
121+
}
122+
return(list(segments=ddsegments, labels=ddlabels))
123+
}
124+
125+
ret <- gg.plotNode(x1, x2, x, type = type, center = center, leaflab = leaflab,
126+
dLeaf = dLeaf, horiz=FALSE,
127+
ddsegments=NULL, ddlabels=NULL)
128+
names(ret$segments) <- c("x", "y", "xend", "yend")
129+
names(ret$labels) <- c("x", "y", "label")
130+
ret
131+
}
132+
133+
134+
135+
.memberDend <- function (x)
136+
{
137+
r <- attr(x, "x.member")
138+
if (is.null(r)) {
139+
r <- attr(x, "members")
140+
if (is.null(r))
141+
r <- 1L
142+
}
143+
r
144+
}
145+
146+
147+
plotNodeLimit <- function (x1, x2, subtree, center)
148+
{
149+
inner <- !is.leaf(subtree) && x1 != x2
150+
if (inner) {
151+
K <- length(subtree)
152+
mTop <- .memberDend(subtree)
153+
limit <- integer(K)
154+
xx1 <- x1
155+
for (k in 1L:K) {
156+
m <- .memberDend(subtree[[k]])
157+
xx1 <- xx1 + (if (center)
158+
(x2 - x1) * m/mTop
159+
else m)
160+
limit[k] <- xx1
161+
}
162+
limit <- c(x1, limit)
163+
}
164+
else {
165+
limit <- c(x1, x2)
166+
}
167+
mid <- attr(subtree, "midpoint")
168+
center <- center || (inner && !is.numeric(mid))
169+
x <- if (center)
170+
mean(c(x1, x2))
171+
else x1 + (if (inner)
172+
mid
173+
else 0)
174+
list(x = x, limit = limit)
175+
}
176+
177+
178+
.midDend <- function (x)
179+
if (is.null(mp <- attr(x, "midpoint"))) 0 else mp
180+
181+
182+

R/grid_util.R

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -105,7 +105,7 @@ ordering_grob <- function(ordering, transpose=FALSE, mirror=FALSE, hint_size=uni
105105
if (is.null(ordering$dendrogram))
106106
return(nullGrob())
107107

108-
dd <- ggdendro::dendro_data(ordering$dendrogram)
108+
dd <- dendro_data(ordering$dendrogram)
109109
dds <- dd$segments
110110
x0 <- dd$segments$x
111111
x1 <- dd$segments$xend

R/heatmap.R

Lines changed: 17 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@
1313
#' @param baseline Baseline level for each row, to be subtracted when drawing the heatmap colors. If omitted, the row mean will be used.
1414
#' @param baseline_label Text description of what the baseline is.
1515
#' @param scale_label Text description of what the heatmap colors represent (after baseline is subtracted).
16+
#' @param n Show only this many rows. Rows are selected in order of greatest span of expression level.
1617
#'
1718
#' @return A grid grob. print()-ing this value will cause it to be displayed.
1819
#'
@@ -28,7 +29,8 @@ plot_heatmap <- function(
2829
feature_labels=NULL,
2930
baseline=NULL,
3031
baseline_label="row\nmean",
31-
scale_label="difference from\nrow mean") {
32+
scale_label="difference from\nrow mean",
33+
n=Inf) {
3234
y <- as.matrix(y)
3335

3436
if (is.null(sample_labels) && !is.null(colnames(y)))
@@ -52,7 +54,20 @@ plot_heatmap <- function(
5254
means <- baseline
5355
else
5456
means <- rowMeans(y, na.rm=TRUE)
55-
57+
58+
59+
# Show only a subset of rows, if desired
60+
if (n < nrow(y)) {
61+
y_span <- apply(y,1,max) - apply(y,1,min)
62+
selection <- rep(FALSE,nrow(y))
63+
selection[ order(-y_span)[ seq_len(n) ] ] <- TRUE
64+
65+
y <- y[selection,,drop=FALSE]
66+
feature_labels <- feature_labels[selection]
67+
means <- means[selection]
68+
}
69+
70+
5671
y_centered <- y - means
5772

5873
y_scaled <- y_centered / sqrt(rowMeans(y_centered*y_centered, na.rm=TRUE))

R/shiny_report.R

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -228,7 +228,7 @@ shiny_filter <- function(y, counts=NULL, sample_labels=NULL, feature_labels=NULL
228228
}
229229

230230

231-
#' Shiny report.
231+
#' Shiny report
232232
#'
233233
#' Produce an interactive Shiny report showing diagnostic plots of transformed counts.
234234
#'
@@ -253,7 +253,8 @@ shiny_filter <- function(y, counts=NULL, sample_labels=NULL, feature_labels=NULL
253253
#' @examples
254254
#'
255255
#' # Generate some random data.
256-
#' counts <- matrix(rnbinom(1000, size=1/0.01, mu=100), ncol=10)
256+
#' means <- runif(100,min=0,max=1000)
257+
#' counts <- matrix(rnbinom(1000, size=1/0.01, mu=rep(means,10)), ncol=10)
257258
#'
258259
#' y <- varistran::vst(counts)
259260
#' if (interactive())

R/varistran.R

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -153,6 +153,17 @@ vst_methods <- list(
153153
#' @references Anscombe, F.J. (1948) "The transformation of Poisson, binomial,
154154
#' and negative-binomial data", Biometrika 35 (3-4): 246-254
155155
#'
156+
#' @examples
157+
#'
158+
#' # Generate some random data.
159+
#' means <- runif(100,min=0,max=1000)
160+
#' counts <- matrix(rnbinom(1000, size=1/0.01, mu=rep(means,10)), ncol=10)
161+
#'
162+
#' y <- varistran::vst(counts)
163+
#'
164+
#' # Information about the transformation
165+
#' varistran::vst_advice(y)
166+
#'
156167
#' @export
157168
vst <- function(x, method="anscombe.nb", lib.size=NULL, cpm=FALSE, dispersion=NULL, design=NULL) {
158169
x <- as.matrix(x)

README.md

Lines changed: 23 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,23 @@ Varistran also uses various CRAN packages, which will be installed automatically
3636

3737
I intend to keep adding features to Varistran. Therefore, I recommend using its functions via `varistran::` rather than attaching its namespace with `library(varistran)`.
3838

39+
The examples below use the [Bottomly dataset](http://bowtie-bio.sourceforge.net/recount/ExpressionSets/bottomly_eset.RData) available from [ReCount](http://bowtie-bio.sourceforge.net/recount/).
40+
41+
```
42+
library(Biobase)
43+
44+
download.file("http://bowtie-bio.sourceforge.net/recount/ExpressionSets/bottomly_eset.RData",
45+
"bottomly_eset.RData")
46+
47+
load("bottomly_eset.RData")
48+
49+
counts <- exprs(bottomly.eset)
50+
51+
strain <- phenoData(bottomly.eset)$strain
52+
experiment.number <- factor( phenoData(bottomly.eset)$experiment.number )
53+
design <- model.matrix(~ strain + experiment.number)
54+
```
55+
3956
Say you have a count matrix `counts` and a design matrix `design`. To perform a variance stabilizing transformation:
4057

4158
```
@@ -48,7 +65,7 @@ An appropraite dispersion is estimated with the aid of the design matrix. If omi
4865

4966
### Diagnostic plots
5067

51-
`plot_stability` allows assessment of how well the variance has been stabilized:
68+
`plot_stability` allows assessment of how well the variance has been stabilized. Ideally this will produce a horizontal line, but counts below 5 will always show a drop off in variance.
5269

5370
```
5471
varistran::plot_stability(y, counts, design=design)
@@ -64,6 +81,10 @@ varistran::plot_biplot(y)
6481

6582
`plot_heatmap` draws a heatmap.
6683

84+
```
85+
varistran::plot_heatmap(y, n=50)
86+
```
87+
6788
<img src="doc/heatmap-example.png" height="300">
6889

6990

@@ -86,7 +107,7 @@ varistran::shiny_report(counts=counts)
86107

87108
## Test suite
88109

89-
After downloading the source code, a suite of tests can be run with:
110+
Download the source code, and ensure that the Bioconductor packages in the "Suggests:" field of the DESCRIPTION file are installed. Then a suite of tests can be run with:
90111

91112
```
92113
make test

man/plot_heatmap.Rd

Lines changed: 4 additions & 1 deletion
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

0 commit comments

Comments
 (0)