knitr::opts_chunk$set(echo=TRUE)
Researchers may wish to superimpose a subset of the full dataset onto the full dataset. If a researcher is using the package to visualize RNA-seq data, then this subset of data is often differentially expressed genes (DEGs) returned from a model. In this case, the user may wish to use the dataMetrics
input parameter, which contains at least one quantitative variable returned from a model such as FDR, p-value, and log fold change.
As was shown in the article Data object, the data
object called soybean_ir_sub
contained 5,604 genes and two treatment groups, N and P [@soybeanIR]. We can examine the structure of its corresponding dataMetrics
object called soybean_ir_sub_metrics
as follows:
library(bigPint) data("soybean_ir_sub_metrics") str(soybean_ir_sub_metrics, strict.width = "wrap")
Similarly, as was shown in the data page, the data
object called soybean_cn_sub
contained 7,332 genes and three treatment groups, S1, S2, and S3 [@brown2015developmental]. We can examine the structure of its corresponding dataMetrics
object called soybean_cn_sub_metrics
as follows:
data("soybean_cn_sub_metrics") str(soybean_cn_sub_metrics, strict.width = "wrap")
As demonstrated in the two examples above, the dataMetrics
object must meet the following conditions:
list
data
object. For example, the soybean_ir_sub_metrics
object contains one list element ("N_P") and the soybean_cn_sub_metrics
object contains three list elements ("S1_S2", "S1_S3", "S2_S3").data.frame
^[a-zA-Z0-9]+_[a-zA-Z0-9]+
, wherecharacter
consisting of the unique names of the genesnumeric
or integer
consisting of a quantitative variable. This can be called anything. In the examples above, there are five of such columns called "logFC", "logCPM", "LR", "PValue", and "FDR".You can quickly double-check the names of the list elements in your dataMetrics
object as follows:
names(soybean_ir_sub_metrics) names(soybean_cn_sub_metrics)
If your dataMetrics
object does not fit this format, bigPint
will likely throw an informative error about why your format was not recognized.
If a researcher is using the bigPint
package to plot RNA-seq data, then many will create the dataMetrics
object using popular RNA-seq packages such as edgeR [@robinson2010edger], DESeq2 [@love2014moderated], and limma [@ritchie2015limma]. These R packages will output several interesting quantitative variables for each gene in the dataset that can be incorporated into the dataMetrics
object. bigPint
can then apply thresholds to these variables and create subsets of genes to superimpose. To create numerous bigPint
plots with the least effort, we recommend creating a dataMetrics
object that contains at least the following column types:
Many bigPint
plots use "FDR" to determine "significant" genes and subset them as overlay (FDR < 0.05). The bigPint
volcano plot uses "PValue" and "logFC". Naming these columns as above will save you time but the names and the default threshold values can be specified away from default when creating each bigPint
plot.
We now provide reproducible code for creating dataMetrics
objects with two or three treatment groups using both edgeR [@robinson2010edger] and DESeq2 [@love2014moderated].
The following example shows how to create the dataMetrics
object called soybean_ir_sub_metrics
, which was shown in the article Data metrics object [@soybeanIR]. The dataset from which it is created (soybean_ir_sub
) contains only two treatment groups, N and P. In this case, the edgeR [@robinson2010edger] package was primarily followed.
library(bigPint) library(edgeR) library(data.table) data(soybean_ir_sub) data = soybean_ir_sub rownames(data) = data[,1] y = DGEList(counts=data[,-1]) group = c(1,1,1,2,2,2) y = DGEList(counts=y, group=group) Group = factor(c(rep("N",3), rep("P",3))) design <- model.matrix(~0+Group, data=y$samples) colnames(design) <- levels(Group) y <- estimateDisp(y, design) fit <- glmFit(y, design) soybean_ir_sub_metrics <- list() for (i in 1:(ncol(fit)-1)){ for (j in (i+1):ncol(fit)){ contrast=rep(0,ncol(fit)) contrast[i]=1 contrast[j]=-1 lrt <- glmLRT(fit, contrast=contrast) lrt <- topTags(lrt, n = nrow(y[[1]]))[[1]] setDT(lrt, keep.rownames = TRUE)[] colnames(lrt)[1] = "ID" lrt <- as.data.frame(lrt) soybean_ir_sub_metrics[[paste0(colnames(fit)[i], "_", colnames(fit)[j])]] <- lrt } }
We can indeed examine the generated soybean_ir_sub_metrics
object as follows:
str(soybean_ir_sub_metrics, strict.width = "wrap")
And verify that it contains one list element:
names(soybean_ir_sub_metrics)
The following example shows how to create the dataMetrics
object called soybean_cn_sub_metrics
, which was shown in the article Data metrics object). The dataset from which it is created (soybean_cn_sub
) contains three treatment groups, S1, S2, and S3 [@brown2015developmental]. In this case, the edgeR [@robinson2010edger] package was primarily followed.
data(soybean_cn_sub) data = soybean_cn_sub rownames(data) = data[,1] y = DGEList(counts=data[,-1]) group = c(1,1,1,2,2,2,3,3,3) y = DGEList(counts=y, group=group) Group = factor(c(rep("S1",3), rep("S2",3), rep("S3",3))) design <- model.matrix(~0+Group, data=y$samples) colnames(design) <- levels(Group) y <- estimateDisp(y, design) fit <- glmFit(y, design) soybean_cn_sub_metrics <- list() for (i in 1:(ncol(fit)-1)){ for (j in (i+1):ncol(fit)){ contrast=rep(0,ncol(fit)) contrast[i]=1 contrast[j]=-1 lrt <- glmLRT(fit, contrast=contrast) lrt <- topTags(lrt, n = nrow(y[[1]]))[[1]] setDT(lrt, keep.rownames = TRUE)[] colnames(lrt)[1] = "ID" lrt <- as.data.frame(lrt) soybean_cn_sub_metrics[[paste0(colnames(fit)[i], "_", colnames(fit)[j])]] <- lrt } }
We can indeed examine the generated soybean_cn_sub_metrics
object as follows:
str(soybean_cn_sub_metrics, strict.width = "wrap")
And verify that it contains three list element:
names(soybean_cn_sub_metrics)
This example shows how to create a dataMetrics
object from (soybean_ir_sub
). In this case, the DESeq2 [@love2014moderated] package was used.
library(DESeq2) data(soybean_ir_sub) data = soybean_ir_sub rownames(data) = data[,1] data = as.matrix(data[,-1]) coldata = data.frame(row.names = colnames(data), treatment = unlist(lapply( colnames(data), function (x) unlist(strsplit(x, "[.]"))[1]))) dds = DESeqDataSetFromMatrix(countData = data, colData = coldata, design = ~ treatment) dds <- DESeq(dds) uTreat = unique(unlist(lapply(colnames(data), function (x) unlist(strsplit( x, "[.]"))[1]))) soybean_ir_sub_metrics <- list() for (i in 1:(length(uTreat)-1)){ for (j in (i+1):length(uTreat)){ res <- results(dds, contrast=c("treatment", uTreat[i], uTreat[j])) metrics = as.data.frame(res@listData) metrics = cbind(ID = res@rownames, metrics) metrics$ID = as.character(metrics$ID) metrics <- metrics[order(metrics$padj), ] soybean_ir_sub_metrics[[paste0(uTreat[i], "_", uTreat[j])]] <- metrics } }
By default, DESeq2
gives output with variables called "pvalue", "padj", and "log2FoldChange". Various functions in bigPint
expect column names like "FDR", "logFC", and "PValue" respectively in the dataMetrics
object. That can be modified manually using the threshVar
input parameter each time creating a plot. But it is easier to simply rename these parameters from the start in the dataMetrics
object.
for (df in seq_len(length(soybean_ir_sub_metrics))){ whichPadj = which(colnames(soybean_ir_sub_metrics[[df]])=="pvalue") colnames(soybean_ir_sub_metrics[[df]])[whichPadj] = "PValue" whichPadj = which(colnames(soybean_ir_sub_metrics[[df]])=="padj") colnames(soybean_ir_sub_metrics[[df]])[whichPadj] = "FDR" whichPadj = which(colnames(soybean_ir_sub_metrics[[df]])=="log2FoldChange") colnames(soybean_ir_sub_metrics[[df]])[whichPadj] = "logFC" }
We can indeed examine the generated soybean_ir_sub_metrics
object as follows:
str(soybean_ir_sub_metrics, strict.width = "wrap")
And verify that it contains one list element:
names(soybean_ir_sub_metrics)
This example shows how to create a dataMetrics
object from (soybean_cn_sub
). In this case, the DESeq2 [@love2014moderated] package was used. The DESeq
package expects a count table that contains integers and has gene-wise dispersion estimates larger than two orders of magnitude from the minimum value. To fit this requirement just for this didactic exercise, we multiply each value by ten and perform a ceiling()
function.
data(soybean_cn_sub) data = soybean_cn_sub rownames(data) = data[,1] data = as.matrix(ceiling(data[,-1] * 10)) coldata = data.frame(row.names = colnames(data), treatment = unlist(lapply( colnames(data), function (x) unlist(strsplit(x, "[.]"))[1]))) dds = DESeqDataSetFromMatrix(countData = data, colData = coldata, design = ~ treatment) dds <- DESeq(dds) uTreat <- unique(unlist(lapply(colnames(data), function (x) unlist(strsplit(x, "[.]"))[1]))) soybean_cn_sub_metrics <- list() for (i in 1:(length(uTreat)-1)){ for (j in (i+1):length(uTreat)){ res <- results(dds, contrast=c("treatment", uTreat[i], uTreat[j])) metrics = as.data.frame(res@listData) metrics = cbind(ID = res@rownames, metrics) metrics$ID = as.character(metrics$ID) metrics <- metrics[order(metrics$padj), ] soybean_cn_sub_metrics[[paste0(uTreat[i], "_", uTreat[j])]] <- metrics } }
By default, DESeq2
gives output with variables called "pvalue", "padj", and "log2FoldChange". Various functions in bigPint
expect column names like "FDR", "logFC", and "PValue" respectively in the dataMetrics
object. That can be modified manually using the threshVar
input parameter each time creating a plot. But it is easier to simply rename these parameters from the start in the dataMetrics
object.
for (df in seq_len(length(soybean_ir_sub_metrics))){ whichPadj = which(colnames(soybean_ir_sub_metrics[[df]])=="pvalue") colnames(soybean_ir_sub_metrics[[df]])[whichPadj] = "PValue" whichPadj = which(colnames(soybean_ir_sub_metrics[[df]])=="padj") colnames(soybean_ir_sub_metrics[[df]])[whichPadj] = "FDR" whichPadj = which(colnames(soybean_ir_sub_metrics[[df]])=="log2FoldChange") colnames(soybean_ir_sub_metrics[[df]])[whichPadj] = "logFC" }
We can indeed examine the generated soybean_cn_sub_metrics
object as follows:
str(soybean_cn_sub_metrics, strict.width = "wrap")
And verify that it contains three list element:
names(soybean_cn_sub_metrics)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.