Nothing
R/DE.plot.R
In NOISeq: Exploratory analysis and differential expression for RNA-seq data
Defines functions
mypretty
DE.plot <- function (output, q = NULL, graphic = c("MD","expr","chrom","distr"),
pch = 20, cex = 0.5, col = 1, pch.sel = 1, cex.sel = 0.6, col.sel = 2,
log.scale = TRUE, chromosomes = NULL, join = FALSE,...) {
mypar <- par(no.readonly = TRUE)
if (class(output) != "Output")
stop("Error. Output argument must contain an object generated by noiseq or noiseqbio functions.\n")
graphic <- match.arg(graphic)
noiseqbio = "theta" %in% colnames(output@results[[1]])[1:4]
## MD plot
if (graphic == "MD") {
if (noiseqbio) {
M = output@results[[1]][,"log2FC"]
D = abs(output@results[[1]][,1] - output@results[[1]][,2])+1
names(M) = names(D) = rownames(output@results[[1]])
plot(M, D, pch = pch, xlab = "M", ylab = "D", cex = cex, col = col, log = "y",...)
if(!is.null(q)) {
mySelection = rownames(degenes(output,q))
points(M[mySelection], D[mySelection], col = col.sel, pch = pch.sel, cex = cex.sel)
}
} else {
plot(output@results[[1]][,"M"], (1+output@results[[1]][,"D"]), pch = pch,
xlab = "M", ylab = "D", cex = cex, col = col, log = "y",...)
if(!is.null(q)) {
mySelection = rownames(degenes(output,q))
points(output@results[[1]][mySelection,"M"], output@results[[1]][mySelection,"D"]+1, col = col.sel, pch = pch.sel, cex = cex.sel)
}
}
}
## Expression plot: Condition1 vs Condition2
else if (graphic == "expr") {
data <- cbind(output@results[[1]][1],output@results[[1]][2])
rownames(data) <- rownames(output@results[[1]])
colnames(data) <- colnames(output@results[[1]][c(1,2)])
if (log.scale) {
k <- min(data, na.rm=TRUE)
escala = logscaling(data, base = 2, k = k)
plot(escala$data, pch = pch, cex = cex, col = col, yaxt = "n", xaxt = "n",...)
axis(side = 1, at = escala$at, labels = escala$labels)
axis(side = 2, at = escala$at, labels = escala$labels)
if(!is.null(q)) {
mySelection = rownames(degenes(output,q))
points(escala$data[mySelection,], col = col.sel, pch = pch.sel, cex = cex.sel)
}
} else {
plot(data, pch = pch, cex = cex, col = col,...)
if(!is.null(q)) {
mySelection = rownames(degenes(output,q))
points(data[mySelection,], col = col.sel, pch = pch.sel, cex = cex.sel)
}
}
}
## MANHATTAN PLOT
else if (graphic == "chrom") {
mydata <- data.frame(as.character(output@results[[1]][,"Chrom"]),
output@results[[1]][,c("GeneStart","GeneEnd")],output@results[[1]][,1:2])
mydata <- na.omit(mydata)
colnames(mydata) <- c("chr", "start", "end", colnames(mydata)[-c(1:3)])
if (is.null(chromosomes)) { # todos los cromosomas
chromosomes <- unique(mydata$chr)
chromosomes = sort(chromosomes)
}
print("REMEMBER. You are plotting these chromosomes and in this order:")
print(chromosomes)
# logarithmic scale
if (log.scale) {
if(min(mydata[,-c(1:3)], na.rm = TRUE) < 1) {
kk <- -min(mydata[,-c(1:3)], na.rm = TRUE)+1
} else { kk <- 0 }
mydata[,-c(1:3)] <- log2(mydata[,-c(1:3)]+kk)
}
# Selecting chromosomes and ordering positions
ordenat = NULL
for (cromo in chromosomes) {
myselec = mydata[mydata[,"chr"] == cromo,]
myselec = myselec[order(myselec[,"start"]),]
ordenat = rbind(ordenat, myselec)
}
sel.ord <- NULL
if (!is.null(q)) {
# up-regulated in first condition
cond1 <- rownames(degenes(output,q,M="up"))
# up-regulated in second condition
cond2 <- rownames(degenes(output,q,M="down"))
cond1 <- (rownames(ordenat) %in% cond1)*1
cond2 <- (rownames(ordenat) %in% cond2)*1
sel.ord <- cbind(cond1, cond2)
rownames(sel.ord) <- rownames(ordenat)
}
chr.long <- aggregate(ordenat[,"end"],
by = list(as.character(ordenat[,"chr"])), max)
chr.long <- chr.long[match(chromosomes, chr.long[,1]),]
if (join) { # si todos los cromosomas van en el mismo plot
total.long <- sum(as.numeric(chr.long$x))
chr.start <- cumsum(c(1,chr.long$x[-length(chr.long$x)]))
names(chr.start) <- chr.long[,1]
plot(c(1,total.long), c(-max(ordenat[,5]), max(ordenat[,4])),
type = "n", xlab = "", ylab = "Expression data", xaxt = "n")
axis(side = 1, at = chr.start, labels = chr.long[,1], font = 2)
abline(h = 0, lty = 2, lwd = 0.5)
for (ch in chromosomes) {
dat.chr <- ordenat[which(ordenat[,"chr"] == ch),]
dat.chr[,c("start","end")] <- dat.chr[,c("start","end")] +
chr.start[ch] - 1
rect(xleft = dat.chr[,"start"], ybottom = 0,
xright = dat.chr[,"end"], ytop = dat.chr[,4],
col = "grey", border = NA)
rect(xleft = dat.chr[,"start"], ybottom = -dat.chr[,5],
xright = dat.chr[,"end"], ytop = 0,
col = "grey", border = NA)
if (!is.null(q)) {
aux <- which(rownames(sel.ord) %in% rownames(dat.chr))
sel.chr1 <- dat.chr[,4]*sel.ord[aux,1]
sel.chr2 <- -dat.chr[,5]*sel.ord[aux,2]
rect(xleft = dat.chr[,"start"], ybottom = 0,
xright = dat.chr[,"end"], ytop = sel.chr1,
col = 2, border = NA)
rect(xleft = dat.chr[,"start"], ybottom = sel.chr2,
xright = dat.chr[,"end"], ytop = 0,
col = 3, border = NA)
}
segments(x0 = dat.chr[,"start"], y0 = 0, x1 = dat.chr[,"end"],
y1 = 0, col = 4, lwd = 0.5) # annotated genes
}
text(c(1,1), 0.9*c(max(ordenat[,4]), -max(ordenat[,5])),
colnames(mydata)[4:5], font = 3, adj = 0, col = 2:3)
# a plot for each chromosome
} else {
num.chr <- length(chromosomes)
k <- 20
long.prop <- round((chr.long$x / min(chr.long$x))*k, 0)
while (max(long.prop)*num.chr > 500 | max(long.prop) > 50) {
k <- k-1
long.prop <- round((chr.long$x / min(chr.long$x))*k, 0)
}
forlayout <- matrix(0, num.chr, max(long.prop))
#print(dim(forlayout))
for (i in 1:num.chr) {
forlayout[i, 1:long.prop[i]] <- i
}
layout(forlayout)
if (num.chr > 1) par(mar=c(2,4.1,0.1,0.1))
miylim <- c(-max(ordenat[,5], na.rm=TRUE), max(ordenat[,4],na.rm=TRUE))
for (i in 1:num.chr) {
apintar <- ordenat[which(ordenat[,"chr"] == chromosomes[i]), 2:5]
plot(c(1,chr.long[i,2]), miylim, type = "n",
xlab = "", ylab = chromosomes[i], xaxt = "n",
font.lab = 2, cex.lab = 1.3)
rect(xleft = apintar[,"start"], ybottom = 0, xright = apintar[,"end"],
ytop = apintar[,3], col = "grey", border = NA)
rect(xleft = apintar[,"start"], ybottom = -apintar[,4],
xright = apintar[,"end"], ytop = 0, col = "grey", border = NA)
if (!is.null(q)) {
aux <- which(rownames(sel.ord) %in% rownames(apintar))
asel <- sel.ord[aux,]*apintar[,3:4]
rect(xleft = apintar[,"start"], ybottom = 0, xright = apintar[,"end"],
ytop = asel[,1], col = 2, border = NA)
rect(xleft = apintar[,"start"], ybottom = -asel[,2],
xright = apintar[,"end"], ytop = 0, col = 3, border = NA)
}
abline(h = 0, lty = 2, lwd = 0.5)
segments(x0 = apintar[,"start"], y0 = 0, x1 = apintar[,"end"],
y1 = 0, col = 4, lwd = 0.5) # annotated genes
etiq <- mypretty(c(1,chr.long[i,2]), n = 10)
axis(side = 1, at = etiq, labels = etiq)
text(c(1,1), 0.9*miylim, colnames(mydata)[5:4],
font = 3, adj = 0, col = 3:2)
}
}
}
## DEG distribution across biotypes/chromosomes
else if (graphic == "distr") {
if(!is.null(q)) { # Computing DEG
mySelection = rownames(degenes(output,q))
detect = rownames(output@results[[1]]) %in% mySelection
} else {
stop("You must specify a valid value for q\n")
}
if ("Chrom" %in% colnames(output@results[[1]])) {
numplot = 1
infobio = output@results[[1]][,"Chrom"]
genome <- 100*table(infobio)/sum(table(infobio))
ordre <- order(genome, decreasing = TRUE)
perdet1 <- genome*table(infobio, detect)[names(genome),2] / table(infobio)[names(genome)]
perdet2 <- 100*table(infobio, detect)[names(genome),2] / sum(table(infobio, detect)[,2])
ceros <- rep(0, length(genome))
biotable1 <- as.matrix(rbind(genome[ordre], perdet1[ordre], perdet2[ordre], ceros))
rownames(biotable1) <- c("genome", "degVSgenome", "deg", "ceros")
if (!is.null(chromosomes)) biotable1 = biotable1[,chromosomes]
ymaxL1 <- ceiling(max(biotable1, na.rm = TRUE))
} else { numplot = 0 }
if ("Biotype" %in% colnames(output@results[[1]])) {
numplot = c(numplot, 1)
infobio = output@results[[1]][,"Biotype"]
genome <- 100*table(infobio)/sum(table(infobio))
ordre <- order(genome, decreasing = TRUE)
perdet1 <- genome*table(infobio, detect)[names(genome),2] / table(infobio)[names(genome)]
perdet2 <- 100*table(infobio, detect)[names(genome),2] / sum(table(infobio, detect)[,2])
ceros <- rep(0, length(genome))
biotable2 <- as.matrix(rbind(genome[ordre], perdet1[ordre], perdet2[ordre], ceros))
rownames(biotable2) <- c("genome", "degVSgenome", "deg", "ceros")
higher2 = which(biotable2[1,] > 2)
lower2 = which(biotable2[1,] <= 2)
if (length(higher2) > 0) {
ymaxL2 <- ceiling(max(biotable2[,higher2], na.rm = TRUE))
if (length(lower2) > 0) {
ymaxR2 <- ceiling(max(biotable2[,lower2], na.rm = TRUE))
biotable2[,lower2] <- biotable2[,lower2]*ymaxL2/ymaxR2
} else { ymaxR2 = ymaxL2 }
} else {
ymaxR2 <- ceiling(max(biotable2[,lower2], na.rm = TRUE))
ymaxL2 = ymaxR2
}
} else { numplot = c(numplot, 0) }
# Plot
if (sum(numplot) == 0) stop("Biotype or chromosome information is needed for this plot\n")
if (sum(numplot) == 1) { # 1 Plot
par(mar = c(10, 4, 2, 2))
if (numplot[1] == 1) {
barplot(biotable1[c(1,3),], main = "DEG distribution across chromosomes",
xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
beside = TRUE, col = c("grey", 2), las = 2,
ylim = c(0, ymaxL1), border = c("grey", 2))
barplot(biotable1[c(2,4),], main = "DEG distribution across chromosomes",
xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
beside = TRUE, col = c(2, 1), las = 2, density = 30,
ylim = c(0, ymaxL1), border = 2, add = TRUE)
legend(x = "topright", bty = "n", horiz = FALSE,
fill = c("grey", 2, 2), density = c(NA,30,NA),
border = c("grey", 2, 2),
legend = c("%Chrom in genome", "%DEG in Chrom", "%Chrom in DEG"))
}
if (numplot[2] == 1) {
barplot(biotable2[c(1,3),], main = "DEG distribution across biotypes",
xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
beside = TRUE, col = c("grey", 4), las = 2,
ylim = c(0, ymaxL2), border = c("grey", 4))
barplot(biotable2[c(2,4),], main = "DEG distribution across biotypes",
xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
beside = TRUE, col = c(4, 1), las = 2, density = 30,
ylim = c(0, ymaxL2), border = 4, add = TRUE)
axis(side=4, at = pretty(c(0,ymaxL2), n = 5),
labels = round(pretty(c(0,ymaxL2), n = 5)*ymaxR2/ymaxL2, 1))
if (ymaxR2 != ymaxL2) {
abline(v = 3*length(higher2) + 0.5, col = 3, lwd = 2, lty = 2)
}
legend(x = "topright", bty = "n", horiz = FALSE,
fill = c("grey", 4, 4), density = c(NA,30,NA),
border = c("grey", 4, 4),
legend = c("%Biotype in genome", "%DEG in Biotype", "%Biotype in DEG"))
}
}
if (sum(numplot) == 2) { # 2 Plots
par(mar = c(10, 4, 2, 2), mfrow = c(1,2))
# Chromosomes
barplot(biotable1[c(1,3),], main = "DEG distribution across chromosomes",
xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
beside = TRUE, col = c("grey", 2), las = 2,
ylim = c(0, ymaxL1), border = c("grey", 2))
barplot(biotable1[c(2,4),], main = "DEG distribution across chromosomes",
xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
beside = TRUE, col = c(2, 1), las = 2, density = 30,
ylim = c(0, ymaxL1), border = 2, add = TRUE)
legend(x = "topright", bty = "n", horiz = FALSE,
fill = c("grey", 2, 2), density = c(NA,30,NA), border = c("grey", 2, 2),
legend = c("%Chrom in genome", "%DEG in Chrom", "%Chrom in DEG"))
# Biotypes
barplot(biotable2[c(1,3),], main = "DEG distribution across biotypes",
xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
beside = TRUE, col = c("grey", 4), las = 2,
ylim = c(0, ymaxL2), border = c("grey", 4))
barplot(biotable2[c(2,4),], main = "DEG distribution across biotypes",
xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
beside = TRUE, col = c(4, 1), las = 2, density = 30,
ylim = c(0, ymaxL2), border = 4, add = TRUE)
axis(side=4, at = pretty(c(0,ymaxL2), n = 5),
labels = round(pretty(c(0,ymaxL2), n = 5)*ymaxR2/ymaxL2, 1))
if (ymaxR2 != ymaxL2) {
abline(v = 3*length(higher2) + 0.5, col = 3, lwd = 2, lty = 2)
}
legend(x = "topright", bty = "n", horiz = FALSE,
fill = c("grey", 4, 4), density = c(NA,30,NA),
border = c("grey", 4, 4),
legend = c("%Biotype in genome", "%DEG in Biotype", "%Biotype in DEG"))
}
}
par(mypar)
}
## Auxiliar function
mypretty <- function(x, n = 5) {
mywidth <- diff(x)
mybin <- ceiling(mywidth/(n-1))
mylabels0 <- x[1] + mybin*(0:(n-1))
ndig <- nchar(as.character(mylabels0))
#print(ndig)
mylabels <- mylabels0
k <- 0
for (i in 2:(n-1)) {
mylabels[i] <- (round(mylabels0[i]/10^(ndig[i]-1),k))*10^(ndig[i]-1)
while (mylabels[i] == mylabels[i-1]) {
k <- k+1
mylabels[i] <- (round(mylabels0[i]/10^(ndig[i]-1),k))*10^(ndig[i]-1)
}
}
mylabels
}
Try the NOISeq package in your browser
Any scripts or data that you put into this service are public.
NOISeq documentation built on Nov. 8, 2020, 5:10 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions mypretty
DE.plot <- function (output, q = NULL, graphic = c("MD","expr","chrom","distr"),
pch = 20, cex = 0.5, col = 1, pch.sel = 1, cex.sel = 0.6, col.sel = 2,
log.scale = TRUE, chromosomes = NULL, join = FALSE,...) {
mypar <- par(no.readonly = TRUE)
if (class(output) != "Output")
stop("Error. Output argument must contain an object generated by noiseq or noiseqbio functions.\n")
graphic <- match.arg(graphic)
noiseqbio = "theta" %in% colnames(output@results[[1]])[1:4]
## MD plot
if (graphic == "MD") {
if (noiseqbio) {
M = output@results[[1]][,"log2FC"]
D = abs(output@results[[1]][,1] - output@results[[1]][,2])+1
names(M) = names(D) = rownames(output@results[[1]])
plot(M, D, pch = pch, xlab = "M", ylab = "D", cex = cex, col = col, log = "y",...)
if(!is.null(q)) {
mySelection = rownames(degenes(output,q))
points(M[mySelection], D[mySelection], col = col.sel, pch = pch.sel, cex = cex.sel)
}
} else {
plot(output@results[[1]][,"M"], (1+output@results[[1]][,"D"]), pch = pch,
xlab = "M", ylab = "D", cex = cex, col = col, log = "y",...)
if(!is.null(q)) {
mySelection = rownames(degenes(output,q))
points(output@results[[1]][mySelection,"M"], output@results[[1]][mySelection,"D"]+1, col = col.sel, pch = pch.sel, cex = cex.sel)
}
}
}
## Expression plot: Condition1 vs Condition2
else if (graphic == "expr") {
data <- cbind(output@results[[1]][1],output@results[[1]][2])
rownames(data) <- rownames(output@results[[1]])
colnames(data) <- colnames(output@results[[1]][c(1,2)])
if (log.scale) {
k <- min(data, na.rm=TRUE)
escala = logscaling(data, base = 2, k = k)
plot(escala$data, pch = pch, cex = cex, col = col, yaxt = "n", xaxt = "n",...)
axis(side = 1, at = escala$at, labels = escala$labels)
axis(side = 2, at = escala$at, labels = escala$labels)
if(!is.null(q)) {
mySelection = rownames(degenes(output,q))
points(escala$data[mySelection,], col = col.sel, pch = pch.sel, cex = cex.sel)
}
} else {
plot(data, pch = pch, cex = cex, col = col,...)
if(!is.null(q)) {
mySelection = rownames(degenes(output,q))
points(data[mySelection,], col = col.sel, pch = pch.sel, cex = cex.sel)
}
}
}
## MANHATTAN PLOT
else if (graphic == "chrom") {
mydata <- data.frame(as.character(output@results[[1]][,"Chrom"]),
output@results[[1]][,c("GeneStart","GeneEnd")],output@results[[1]][,1:2])
mydata <- na.omit(mydata)
colnames(mydata) <- c("chr", "start", "end", colnames(mydata)[-c(1:3)])
if (is.null(chromosomes)) { # todos los cromosomas
chromosomes <- unique(mydata$chr)
chromosomes = sort(chromosomes)
}
print("REMEMBER. You are plotting these chromosomes and in this order:")
print(chromosomes)
# logarithmic scale
if (log.scale) {
if(min(mydata[,-c(1:3)], na.rm = TRUE) < 1) {
kk <- -min(mydata[,-c(1:3)], na.rm = TRUE)+1
} else { kk <- 0 }
mydata[,-c(1:3)] <- log2(mydata[,-c(1:3)]+kk)
}
# Selecting chromosomes and ordering positions
ordenat = NULL
for (cromo in chromosomes) {
myselec = mydata[mydata[,"chr"] == cromo,]
myselec = myselec[order(myselec[,"start"]),]
ordenat = rbind(ordenat, myselec)
}
sel.ord <- NULL
if (!is.null(q)) {
# up-regulated in first condition
cond1 <- rownames(degenes(output,q,M="up"))
# up-regulated in second condition
cond2 <- rownames(degenes(output,q,M="down"))
cond1 <- (rownames(ordenat) %in% cond1)*1
cond2 <- (rownames(ordenat) %in% cond2)*1
sel.ord <- cbind(cond1, cond2)
rownames(sel.ord) <- rownames(ordenat)
}
chr.long <- aggregate(ordenat[,"end"],
by = list(as.character(ordenat[,"chr"])), max)
chr.long <- chr.long[match(chromosomes, chr.long[,1]),]
if (join) { # si todos los cromosomas van en el mismo plot
total.long <- sum(as.numeric(chr.long$x))
chr.start <- cumsum(c(1,chr.long$x[-length(chr.long$x)]))
names(chr.start) <- chr.long[,1]
plot(c(1,total.long), c(-max(ordenat[,5]), max(ordenat[,4])),
type = "n", xlab = "", ylab = "Expression data", xaxt = "n")
axis(side = 1, at = chr.start, labels = chr.long[,1], font = 2)
abline(h = 0, lty = 2, lwd = 0.5)
for (ch in chromosomes) {
dat.chr <- ordenat[which(ordenat[,"chr"] == ch),]
dat.chr[,c("start","end")] <- dat.chr[,c("start","end")] +
chr.start[ch] - 1
rect(xleft = dat.chr[,"start"], ybottom = 0,
xright = dat.chr[,"end"], ytop = dat.chr[,4],
col = "grey", border = NA)
rect(xleft = dat.chr[,"start"], ybottom = -dat.chr[,5],
xright = dat.chr[,"end"], ytop = 0,
col = "grey", border = NA)
if (!is.null(q)) {
aux <- which(rownames(sel.ord) %in% rownames(dat.chr))
sel.chr1 <- dat.chr[,4]*sel.ord[aux,1]
sel.chr2 <- -dat.chr[,5]*sel.ord[aux,2]
rect(xleft = dat.chr[,"start"], ybottom = 0,
xright = dat.chr[,"end"], ytop = sel.chr1,
col = 2, border = NA)
rect(xleft = dat.chr[,"start"], ybottom = sel.chr2,
xright = dat.chr[,"end"], ytop = 0,
col = 3, border = NA)
}
segments(x0 = dat.chr[,"start"], y0 = 0, x1 = dat.chr[,"end"],
y1 = 0, col = 4, lwd = 0.5) # annotated genes
}
text(c(1,1), 0.9*c(max(ordenat[,4]), -max(ordenat[,5])),
colnames(mydata)[4:5], font = 3, adj = 0, col = 2:3)
# a plot for each chromosome
} else {
num.chr <- length(chromosomes)
k <- 20
long.prop <- round((chr.long$x / min(chr.long$x))*k, 0)
while (max(long.prop)*num.chr > 500 | max(long.prop) > 50) {
k <- k-1
long.prop <- round((chr.long$x / min(chr.long$x))*k, 0)
}
forlayout <- matrix(0, num.chr, max(long.prop))
#print(dim(forlayout))
for (i in 1:num.chr) {
forlayout[i, 1:long.prop[i]] <- i
}
layout(forlayout)
if (num.chr > 1) par(mar=c(2,4.1,0.1,0.1))
miylim <- c(-max(ordenat[,5], na.rm=TRUE), max(ordenat[,4],na.rm=TRUE))
for (i in 1:num.chr) {
apintar <- ordenat[which(ordenat[,"chr"] == chromosomes[i]), 2:5]
plot(c(1,chr.long[i,2]), miylim, type = "n",
xlab = "", ylab = chromosomes[i], xaxt = "n",
font.lab = 2, cex.lab = 1.3)
rect(xleft = apintar[,"start"], ybottom = 0, xright = apintar[,"end"],
ytop = apintar[,3], col = "grey", border = NA)
rect(xleft = apintar[,"start"], ybottom = -apintar[,4],
xright = apintar[,"end"], ytop = 0, col = "grey", border = NA)
if (!is.null(q)) {
aux <- which(rownames(sel.ord) %in% rownames(apintar))
asel <- sel.ord[aux,]*apintar[,3:4]
rect(xleft = apintar[,"start"], ybottom = 0, xright = apintar[,"end"],
ytop = asel[,1], col = 2, border = NA)
rect(xleft = apintar[,"start"], ybottom = -asel[,2],
xright = apintar[,"end"], ytop = 0, col = 3, border = NA)
}
abline(h = 0, lty = 2, lwd = 0.5)
segments(x0 = apintar[,"start"], y0 = 0, x1 = apintar[,"end"],
y1 = 0, col = 4, lwd = 0.5) # annotated genes
etiq <- mypretty(c(1,chr.long[i,2]), n = 10)
axis(side = 1, at = etiq, labels = etiq)
text(c(1,1), 0.9*miylim, colnames(mydata)[5:4],
font = 3, adj = 0, col = 3:2)
}
}
}
## DEG distribution across biotypes/chromosomes
else if (graphic == "distr") {
if(!is.null(q)) { # Computing DEG
mySelection = rownames(degenes(output,q))
detect = rownames(output@results[[1]]) %in% mySelection
} else {
stop("You must specify a valid value for q\n")
}
if ("Chrom" %in% colnames(output@results[[1]])) {
numplot = 1
infobio = output@results[[1]][,"Chrom"]
genome <- 100*table(infobio)/sum(table(infobio))
ordre <- order(genome, decreasing = TRUE)
perdet1 <- genome*table(infobio, detect)[names(genome),2] / table(infobio)[names(genome)]
perdet2 <- 100*table(infobio, detect)[names(genome),2] / sum(table(infobio, detect)[,2])
ceros <- rep(0, length(genome))
biotable1 <- as.matrix(rbind(genome[ordre], perdet1[ordre], perdet2[ordre], ceros))
rownames(biotable1) <- c("genome", "degVSgenome", "deg", "ceros")
if (!is.null(chromosomes)) biotable1 = biotable1[,chromosomes]
ymaxL1 <- ceiling(max(biotable1, na.rm = TRUE))
} else { numplot = 0 }
if ("Biotype" %in% colnames(output@results[[1]])) {
numplot = c(numplot, 1)
infobio = output@results[[1]][,"Biotype"]
genome <- 100*table(infobio)/sum(table(infobio))
ordre <- order(genome, decreasing = TRUE)
perdet1 <- genome*table(infobio, detect)[names(genome),2] / table(infobio)[names(genome)]
perdet2 <- 100*table(infobio, detect)[names(genome),2] / sum(table(infobio, detect)[,2])
ceros <- rep(0, length(genome))
biotable2 <- as.matrix(rbind(genome[ordre], perdet1[ordre], perdet2[ordre], ceros))
rownames(biotable2) <- c("genome", "degVSgenome", "deg", "ceros")
higher2 = which(biotable2[1,] > 2)
lower2 = which(biotable2[1,] <= 2)
if (length(higher2) > 0) {
ymaxL2 <- ceiling(max(biotable2[,higher2], na.rm = TRUE))
if (length(lower2) > 0) {
ymaxR2 <- ceiling(max(biotable2[,lower2], na.rm = TRUE))
biotable2[,lower2] <- biotable2[,lower2]*ymaxL2/ymaxR2
} else { ymaxR2 = ymaxL2 }
} else {
ymaxR2 <- ceiling(max(biotable2[,lower2], na.rm = TRUE))
ymaxL2 = ymaxR2
}
} else { numplot = c(numplot, 0) }
# Plot
if (sum(numplot) == 0) stop("Biotype or chromosome information is needed for this plot\n")
if (sum(numplot) == 1) { # 1 Plot
par(mar = c(10, 4, 2, 2))
if (numplot[1] == 1) {
barplot(biotable1[c(1,3),], main = "DEG distribution across chromosomes",
xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
beside = TRUE, col = c("grey", 2), las = 2,
ylim = c(0, ymaxL1), border = c("grey", 2))
barplot(biotable1[c(2,4),], main = "DEG distribution across chromosomes",
xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
beside = TRUE, col = c(2, 1), las = 2, density = 30,
ylim = c(0, ymaxL1), border = 2, add = TRUE)
legend(x = "topright", bty = "n", horiz = FALSE,
fill = c("grey", 2, 2), density = c(NA,30,NA),
border = c("grey", 2, 2),
legend = c("%Chrom in genome", "%DEG in Chrom", "%Chrom in DEG"))
}
if (numplot[2] == 1) {
barplot(biotable2[c(1,3),], main = "DEG distribution across biotypes",
xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
beside = TRUE, col = c("grey", 4), las = 2,
ylim = c(0, ymaxL2), border = c("grey", 4))
barplot(biotable2[c(2,4),], main = "DEG distribution across biotypes",
xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
beside = TRUE, col = c(4, 1), las = 2, density = 30,
ylim = c(0, ymaxL2), border = 4, add = TRUE)
axis(side=4, at = pretty(c(0,ymaxL2), n = 5),
labels = round(pretty(c(0,ymaxL2), n = 5)*ymaxR2/ymaxL2, 1))
if (ymaxR2 != ymaxL2) {
abline(v = 3*length(higher2) + 0.5, col = 3, lwd = 2, lty = 2)
}
legend(x = "topright", bty = "n", horiz = FALSE,
fill = c("grey", 4, 4), density = c(NA,30,NA),
border = c("grey", 4, 4),
legend = c("%Biotype in genome", "%DEG in Biotype", "%Biotype in DEG"))
}
}
if (sum(numplot) == 2) { # 2 Plots
par(mar = c(10, 4, 2, 2), mfrow = c(1,2))
# Chromosomes
barplot(biotable1[c(1,3),], main = "DEG distribution across chromosomes",
xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
beside = TRUE, col = c("grey", 2), las = 2,
ylim = c(0, ymaxL1), border = c("grey", 2))
barplot(biotable1[c(2,4),], main = "DEG distribution across chromosomes",
xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
beside = TRUE, col = c(2, 1), las = 2, density = 30,
ylim = c(0, ymaxL1), border = 2, add = TRUE)
legend(x = "topright", bty = "n", horiz = FALSE,
fill = c("grey", 2, 2), density = c(NA,30,NA), border = c("grey", 2, 2),
legend = c("%Chrom in genome", "%DEG in Chrom", "%Chrom in DEG"))
# Biotypes
barplot(biotable2[c(1,3),], main = "DEG distribution across biotypes",
xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
beside = TRUE, col = c("grey", 4), las = 2,
ylim = c(0, ymaxL2), border = c("grey", 4))
barplot(biotable2[c(2,4),], main = "DEG distribution across biotypes",
xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
beside = TRUE, col = c(4, 1), las = 2, density = 30,
ylim = c(0, ymaxL2), border = 4, add = TRUE)
axis(side=4, at = pretty(c(0,ymaxL2), n = 5),
labels = round(pretty(c(0,ymaxL2), n = 5)*ymaxR2/ymaxL2, 1))
if (ymaxR2 != ymaxL2) {
abline(v = 3*length(higher2) + 0.5, col = 3, lwd = 2, lty = 2)
}
legend(x = "topright", bty = "n", horiz = FALSE,
fill = c("grey", 4, 4), density = c(NA,30,NA),
border = c("grey", 4, 4),
legend = c("%Biotype in genome", "%DEG in Biotype", "%Biotype in DEG"))
}
}
par(mypar)
}
## Auxiliar function
mypretty <- function(x, n = 5) {
mywidth <- diff(x)
mybin <- ceiling(mywidth/(n-1))
mylabels0 <- x[1] + mybin*(0:(n-1))
ndig <- nchar(as.character(mylabels0))
#print(ndig)
mylabels <- mylabels0
k <- 0
for (i in 2:(n-1)) {
mylabels[i] <- (round(mylabels0[i]/10^(ndig[i]-1),k))*10^(ndig[i]-1)
while (mylabels[i] == mylabels[i-1]) {
k <- k+1
mylabels[i] <- (round(mylabels0[i]/10^(ndig[i]-1),k))*10^(ndig[i]-1)
}
}
mylabels
}
Try the NOISeq package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.