Creating data metrics object

Data metrics object

Information about the dataMetrics object can be found in the article Data metrics object. 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 (Robinson, McCarthy, and Smyth 2010), DESeq2 (Love, Huber, and Anders 2014), and limma (Ritchie et al. 2015). 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:

Significance level (“PValue”)
Multiple comparison significance level (“FDR”)
Log fold change (“logFC”)

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 (Robinson, McCarthy, and Smyth 2010) and DESeq2 (Love, Huber, and Anders 2014).

Example: two treatments (edgeR)

The following example shows how to create the dataMetrics object called soybean_ir_sub_metrics, which was shown in the article Data metrics object (Lauter and Graham 2016). The dataset from which it is created (soybean_ir_sub) contains only two treatment groups, N and P. In this case, the edgeR (Robinson, McCarthy, and Smyth 2010) 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")

## List of 1
## $ N_P:'data.frame': 5604 obs. of 6 variables:
## ..$ ID : chr [1:5604] "Glyma.19G168700.Wm82.a2.v1" "Glyma.13G293500.Wm82.a2.v1"
##    "Glyma.05G188700.Wm82.a2.v1" "Glyma.13G173100.Wm82.a2.v1" ...
## ..$ logFC : num [1:5604] -5.92 2.99 -3.51 -3.91 -3.51 ...
## ..$ logCPM: num [1:5604] 7.52 8.08 8.83 8.27 10.19 ...
## ..$ LR : num [1:5604] 266 171 167 157 154 ...
## ..$ PValue: num [1:5604] 9.18e-60 3.65e-39 2.73e-38 6.04e-36 2.58e-35 ...
## ..$ FDR : num [1:5604] 5.14e-56 1.02e-35 5.09e-35 8.46e-33 2.89e-32 ...

And verify that it contains one list element:

names(soybean_ir_sub_metrics)

## [1] "N_P"

Example: three treatments (edgeR)

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 (Brown and Hudson 2015). In this case, the edgeR (Robinson, McCarthy, and Smyth 2010) 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")

## List of 3
## $ S1_S2:'data.frame': 7332 obs. of 6 variables:
## ..$ ID : chr [1:7332] "Glyma18g00690.1" "Glyma08g22380.1" "Glyma20g30460.1"
##    "Glyma07g09700.1" ...
## ..$ logFC : num [1:7332] 3.14 2.62 2.5 2.37 2.74 ...
## ..$ logCPM: num [1:7332] 8.04 8.05 8.14 8.19 7.93 ...
## ..$ LR : num [1:7332] 29.3 24.8 24 22.7 21 ...
## ..$ PValue: num [1:7332] 6.08e-08 6.50e-07 9.82e-07 1.94e-06 4.57e-06 ...
## ..$ FDR : num [1:7332] 0.000446 0.002384 0.0024 0.003557 0.006697 ...
## $ S1_S3:'data.frame': 7332 obs. of 6 variables:
## ..$ ID : chr [1:7332] "Glyma08g22380.1" "Glyma08g19290.1" "Glyma20g30460.1"
##    "Glyma18g00690.1" ...
## ..$ logFC : num [1:7332] 3.18 3.32 2.44 2.46 3.03 ...
## ..$ logCPM: num [1:7332] 8.05 8.14 8.14 8.04 7.85 ...
## ..$ LR : num [1:7332] 30.4 24.7 23.3 22.5 22.3 ...
## ..$ PValue: num [1:7332] 3.60e-08 6.53e-07 1.42e-06 2.09e-06 2.30e-06 ...
## ..$ FDR : num [1:7332] 0.000264 0.002393 0.003378 0.003378 0.003378 ...
## $ S2_S3:'data.frame': 7332 obs. of 6 variables:
## ..$ ID : chr [1:7332] "Glyma08g14670.3" "Glyma08g14670.2" "Glyma08g19290.1"
##    "Glyma06g46701.1" ...
## ..$ logFC : num [1:7332] 2.66 2.72 2.78 2.19 2.29 ...
## ..$ logCPM: num [1:7332] 7.86 7.82 8.14 7.63 7.84 ...
## ..$ LR : num [1:7332] 16.1 16 14.3 11.8 10.5 ...
## ..$ PValue: num [1:7332] 6.01e-05 6.30e-05 1.57e-04 5.81e-04 1.21e-03 ...
## ..$ FDR : num [1:7332] 0.231 0.231 0.385 1 1 ...

And verify that it contains three list element:

names(soybean_cn_sub_metrics)

## [1] "S1_S2" "S1_S3" "S2_S3"

Example: two treatments (DESeq2)

This example shows how to create a dataMetrics object from (soybean_ir_sub). In this case, the DESeq2 (Love, Huber, and Anders 2014) 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")

## List of 1
## $ N_P:'data.frame': 5604 obs. of 7 variables:
## ..$ ID : chr [1:5604] "Glyma.11G141800.Wm82.a2.v1" "Glyma.11G044800.Wm82.a2.v1"
##    "Glyma.05G188700.Wm82.a2.v1" "Glyma.13G173100.Wm82.a2.v1" ...
## ..$ baseMean: num [1:5604] 778 795 302 204 420 ...
## ..$ logFC : num [1:5604] -3.54 2.75 -3.53 -3.94 -3.1 ...
## ..$ lfcSE : num [1:5604] 0.249 0.207 0.267 0.31 0.244 ...
## ..$ stat : num [1:5604] -14.2 13.3 -13.2 -12.7 -12.7 ...
## ..$ PValue : num [1:5604] 5.49e-46 3.49e-40 5.34e-40 4.80e-37 5.76e-37 ...
## ..$ FDR : num [1:5604] 1.95e-42 6.20e-37 6.32e-37 4.09e-34 4.09e-34 ...

And verify that it contains one list element:

names(soybean_ir_sub_metrics)

## [1] "N_P"

Example: three treatments (DESeq2)

This example shows how to create a dataMetrics object from (soybean_cn_sub). In this case, the DESeq2 (Love, Huber, and Anders 2014) 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
    }
}

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")

## List of 3
## $ S1_S2:'data.frame': 7332 obs. of 7 variables:
## ..$ ID : chr [1:7332] "Glyma18g00690.1" "Glyma03g29150.1" "Glyma05g27450.2"
##    "Glyma10g31780.1" ...
## ..$ baseMean : num [1:7332] 50.4 33 46.3 45.3 29.8 ...
## ..$ log2FoldChange: num [1:7332] 3.22 3.14 -2.89 2.83 2.99 ...
## ..$ lfcSE : num [1:7332] 0.427 0.419 0.4 0.395 0.42 ...
## ..$ stat : num [1:7332] 7.55 7.5 -7.22 7.16 7.12 ...
## ..$ pvalue : num [1:7332] 4.30e-14 6.35e-14 5.03e-13 7.92e-13 1.06e-12 ...
## ..$ padj : num [1:7332] 1.92e-10 1.92e-10 1.01e-09 1.20e-09 1.28e-09 ...
## $ S1_S3:'data.frame': 7332 obs. of 7 variables:
## ..$ ID : chr [1:7332] "Glyma08g19290.1" "Glyma08g22380.1" "Glyma16g08810.1"
##    "Glyma04g37510.1" ...
## ..$ baseMean : num [1:7332] 55.4 50.8 34 39.1 45.7 ...
## ..$ log2FoldChange: num [1:7332] 3.39 3.26 3.17 -2.99 -2.9 ...
## ..$ lfcSE : num [1:7332] 0.391 0.404 0.424 0.405 0.4 ...
## ..$ stat : num [1:7332] 8.67 8.06 7.49 -7.4 -7.26 ...
## ..$ pvalue : num [1:7332] 4.29e-18 7.39e-16 7.15e-14 1.40e-13 3.90e-13 ...
## ..$ padj : num [1:7332] 2.66e-14 2.29e-12 1.48e-10 2.17e-10 4.83e-10 ...
## $ S2_S3:'data.frame': 7332 obs. of 7 variables:
## ..$ ID : chr [1:7332] "Glyma08g19290.1" "Glyma08g14670.2" "Glyma08g14670.3"
##    "Glyma02g06650.1" ...
## ..$ baseMean : num [1:7332] 55.4 40.7 42.3 41.4 45.7 ...
## ..$ log2FoldChange: num [1:7332] 2.81 2.77 2.68 2.31 1.92 ...
## ..$ lfcSE : num [1:7332] 0.396 0.395 0.387 0.4 0.342 ...
## ..$ stat : num [1:7332] 7.1 7.02 6.94 5.77 5.62 ...
## ..$ pvalue : num [1:7332] 1.28e-12 2.28e-12 3.88e-12 8.01e-09 1.94e-08 ...
## ..$ padj : num [1:7332] 6.91e-09 6.91e-09 7.83e-09 1.21e-05 2.35e-05 ...

And verify that it contains three list element:

names(soybean_cn_sub_metrics)

## [1] "S1_S2" "S1_S3" "S2_S3"

References

Brown, Anne V., and Karen A. Hudson. 2015. “Developmental Profiling of Gene Expression in Soybean Trifoliate Leaves and Cotyledons.” BMC Plant Biology 15 (1): 169.

Lauter, AN Moran, and MA Graham. 2016. “NCBI Sra Bioproject Accession: PRJNA318409.”

Love, Michael I., Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for Rna-Seq Data with Deseq2.” Genome Biology 15 (12): 550.

Ritchie, Matthew E., Belinda Phipson, Di Wu, Yifang Hu, Charity W. Law, Wei Shi, and Gordon K. Smyth. 2015. “Limma Powers Differential Expression Analyses for Rna-Sequencing and Microarray Studies.” Nucleic Acids Research 43 (7): e47–e47.

Robinson, Mark D., Davis J. McCarthy, and Gordon K. Smyth. 2010. “EdgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics 26 (1): 139–40.

Lindsay Rutter