2.1 DGEobj Nomenclature: Items, Types and Base Types

Multiple instances of a base type are accommodated by defining data “types” and “items”. Each data type is assigned a base type (e.g. geneData, GRanges, and Fit objects are all “types” of “baseType” = rowData).

Importantly, multiple instances of each type are allowed (exceptions described in Unique Items section) as long as each instance of a type is named uniquely. Each instance of a “type” is described as an “item” and each item must have a user-defined item name. The item name must be unique within a DGEobj.

The DGEobj package implements an S3 class data object, the DGEobj, that conceptually represents an extension of the capabilities of the RangedSummarizedExperiment (RSE) originally developed by Martin Morgan et al. Both the DGEobj and the RSE object capture the initial data for a differential gene expression analysis, namely a counts matrix along with associated gene and sample annotation. Additionally, the DGEobj extends this concept to support capture of various downstream data objects, the intermediate steps in an analysis, and thus capture the entire workflow of an analysis project.

The availability of a collection of results in a specific structured data object has multiple advantages. Sharing DGE results with other colleagues is simplified because the entire analysis is encapsulated and documented within the DGEobj. A recipient of data in this format can examine details of the analysis based on the annotation built into the DGE object. Downstream users of a dataset can drill back to the original data or tap any intermediate data items the original analyst captured. Importantly, assembling a collection of DGE projects in a structured format enables facile automation of integrative multi-project metanalyses.

The RSE object, which motivated building the DGEobj, can capture as many assays as desired. An “assay” is defined as any matrix with n genes (rows) and m samples (columns). A limitation of the RSE object however is that the RSE only allows 1 instance of row data (typically gene annotation) and 1 instance of column data (sample annotation with one row for every column of data in the assay slot). This limits the RSE in terms of its ability to hold downstream data objects because many of those objects meet the definition of row data also (e.g. DGElist, Fit objects, topTable output). Other types of data (e.g. design matrices, sample QC) meet the definition of column data. Thus, the DGEobj was modeled after the RSE object, but extended to accommodate multiple row and column data types. The DGEobj is thus uniquely suited to capturing the entire workflow of a DGE analysis.

3.1 Preserving the Original Data

In the course of a workflow an analyst will often need to subset the DGEobj based on low signal or other characteristics. However, a downstream analyst may need to start their analysis with the original unsubsetted data. For this reason, when the DGEobj is initialized, a copy of the starting data is also stored in a metadata slot. Metadata slots are carried along without subsetting so the original data may always be retrieved from these metadata items. The item names of the original data in the metadata slot have a “_orig” suffix.

4.1 Demonstration Dataset: Rat Liver Slice Compound Treatments

The following dataset was selected to demonstrate building and working with the DGEobj data structure.

Huang X, Cai H, Ammar R, Zhang Y et al. Molecular characterization of a precision-cut rat liver slice model for the evaluation of antifibrotic compounds. Am J Physiol Gastrointest Liver Physiol 2019 Jan 1;316(1):G15-G24. PMID: 30406699

Briefly, livers were removed from rats 4 weeks after bile duct ligation or sham operation. Rat liver slices were incubated in vitro with potential anti-fibrotic compounds. At the end of the incubation whole transcriptome RNA-Seq analysis was performed.

Files containing counts, sample annotations, and QC data associated with this project will be downloaded from the NCBI GEO resource.

The GEO data includes Ensembl gene IDs. Additional gene information such as chromosome positions, type of transcript, etc, will be downloaded from Ensembl using the biomaRt package.

4.2 Initializing a DGEobj from dataframes

Three properly formatted data frames are required to initialize a new DGEobj:

1. Counts matrix: For RNA-Seq, this is typically a genes x samples matrix or dataframe of numbers. Rownames must be present and contain a unique identifier (e.g. Ensembl geneid). Colnames should be a unique sample identifier.

2. Row Data: Additional annotations for the entities represented by each row. For the RNA-Seq case, this dataframe holds associated gene annotation. The rownames in the RowData must match the rownames in the counts matrix.

3. Col Data: Accordingly the ColData dataframe contains information about the samples (columns) of the counts matrix. We often refer to this as the “design” table because it typically contains the experimental details of the treatment of each sample. The rownames of ColData must match the colnames of the Counts Matrix.

In addition to these three items, the user must specify the “level” of the data (allowedLevels = c(“gene”, “isoform”, “exon”, “affy”)).

See the “Customizing the DGEobj” section below to define new levels.

# Get the raw counts, sample annotation ("design") and QC data from GEO
# Source: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE120804
getLocation <- "http://ftp.ncbi.nlm.nih.gov/geo/series/GSE120nnn/GSE120804/suppl"
countsFile  <- "GSE120804_counts.txt.gz"
designFile  <- "GSE120804_geo_sample_annotation_edit.csv.gz"
qcFile      <- "GSE120804_qc_metrics.txt.gz"

temp <- tempfile()
if (download.file(glue("{getLocation}/{countsFile}"), destfile = temp, mode = 'wb') > 0) {
    stop("Counts Download Failed")
}
counts <- read.delim(temp, stringsAsFactors = FALSE, row.names = 1)

temp <- tempfile()
if (download.file(glue("{getLocation}/{designFile}"), destfile = temp, mode = 'wb')) {
    stop("Design Download Failed")
}
design <- read.csv(temp, stringsAsFactors = FALSE)

temp <- tempfile()
if (download.file(glue("{getLocation}/{qcFile}"), destfile = temp, mode = 'wb')) {
    stop("Alignment QC Download Failed")
}
alignmentQC <- read.delim(temp, stringsAsFactors = FALSE)

rownames(design) <- str_sub(design$raw.file, start = 1, end = 21)

#correct the desired case/spelling of one column
design <- design %>%
    rename(ReplicateGroup = Replicate.group)

# Create a DiseaseStatus column by parsing ReplicateGroup
design$DiseaseStatus <- rep("Sham", nrow(design))
idx <- str_detect(design$ReplicateGroup, "BDL")
design$DiseaseStatus[idx] <- "BDL"

# Create an animal# column.  The animal number is encoded in the sample.name column
design$AnimalNum <- str_match(design$Sample.name, "r[0-9]{1,3}")

# Now get the gene annotation from Ensembl/biomaRt
ens.ds      <- "rnorvegicus_gene_ensembl"
ens.mart    <- useMart(biomart = "ensembl", dataset = ens.ds, host = "https://may2021.archive.ensembl.org")
ens.columns <- c("ensembl_gene_id", "rgd_symbol", "chromosome_name", "start_position",
                 "end_position", "strand", "gene_biotype", "description")
ens.data    <- getBM(attributes = ens.columns,
                     values     = rownames(counts),
                     mart       = ens.mart,
                     useCache   = F) %>%
               distinct(ensembl_gene_id, .keep_all = T)

# Filter the list to the genes used in the test dataset and properly format gene
# information for GenomicRanges use
gene.data <- left_join(data.frame(ensembl_gene_id = rownames(counts), stringsAsFactors = F),
                       ens.data,
                       by = "ensembl_gene_id") %>%
    dplyr::rename(start = start_position, end = end_position) %>%
    mutate(strand = case_when(strand == -1 ~ "-",
                              strand == 1  ~ "+",
                              TRUE         ~ "*"))
rownames(gene.data) <- gene.data$ensembl_gene_id

# Get transcript level data and keep max for each gene, or alternatively, use the cds length

ens.ds          <- "rnorvegicus_gene_ensembl"
ens.mart        <- useMart(biomart = "ensembl", dataset = ens.ds, host = "https://may2021.archive.ensembl.org")
ens.columns     <- c("ensembl_gene_id", "ensembl_transcript_id", "transcript_length")
transcript.data <- getBM(attributes = ens.columns,
                         values     = rownames(counts),
                         mart       = ens.mart,
                         useCache   = F) %>%
                   arrange(desc(transcript_length)) %>%
                   distinct(ensembl_gene_id, .keep_all = T)

#Add a transcript_length column to gene.data
gene.data <- left_join(gene.data,
                       select(transcript.data, ensembl_gene_id, ExonLength = transcript_length))
rownames(gene.data) <- gene.data$ensembl_gene_id
gene.data$ensembl_gene_id <- NULL

# enforce same order as counts
gene.data <- gene.data[rownames(counts),]

all(rownames(counts) == rownames(gene.data))

## [1] TRUE

all(colnames(counts) == rownames(design))

## [1] TRUE

dgeObj <- DGEobj::initDGEobj(primaryAssayData  = counts,
                             rowData = gene.data,
                             colData = design,
                             level = "gene",
                             customAttr = list(Genome    = "Rat.B6.0",
                                               GeneModel = "Ensembl.R89"))

# the annotation information is stored in a temporary .txt file
dgeObj <- annotateDGEobj(dgeObj, annotationFile)

5.1 Inventory function

To view a listing of the items in the DGEobj:

kable(inventory(dgeObj))

5.2 Examine DGEobj Metadata

Use the showMeta function to list the meta data associated with a project.

kable(showMeta(dgeObj))

5.3 DGEobj Length and Dimensions

The length of a DGEobj refers to the number of data items in the DGEobj.

length(dgeObj)

## [1] 8

The dimensions reported for a DGEobj are the dimensions of the assays contained in the DGEobj. That is, for the RNA-Seq application, the row dimension is the number of genes contained in the object and the column dimension is the number of samples contained in the object.

dim(dgeObj)

## [1] 32883    48

dim(dgeObj$counts)

## [1] 32883    48

5.4 Rownames and Colnames

The rownames and colnames of the DGEobj are defined by the primary “assay” matrix in a DGEobj, typically the counts matrix for RNA-Seq data.

Note, there is no facility provided to re-assign rownames or colnames of a DGEobj. The user should make any adjustments to row/colnames before instantiating the DGEobj.

5.5 Subsetting a DGEobj

A DGEobj may be square bracket subsetted similar the same way data frames and matrices may be subsetted. The subsetting function uses the base type to define how each data item is handled during subsetting. During subsetting, metadata items are carried along unchanged.

#subset to the first 100 genes
First100genes <- dgeObj[1:100,]
dim(First100genes)

## [1] 100  48

#subset to the first 10 samples
First10samples <- dgeObj[,1:10]
dim(First10samples)

## [1] 32883    10

#susbet genes and samples
AnotherSubset <- dgeObj[1:100, 1:10]
dim(AnotherSubset)

## [1] 100  10

Boolean vectors may also be used to subset dimensions of a DGEobj.

# select a subset of samples
idx <- dgeObj$design$DiseaseStatus == "BDL"
BDLonly <- dgeObj[,idx]
dim(BDLonly)

## [1] 32883    39

5.6 Inventory of a DGEobj

The inventory function prints a table of the data items present in a DGEobj. The output includes the item name, type, baseType, parent, class, and date created. If the verbose = TRUE argument is used, a funArgs column is also included.

kable(inventory(dgeObj))

If just the item names of data stored in the DGEobj are needed, use the base R names function.

names(dgeObj)

## [1] "counts_orig"   "counts"        "design_orig"   "design"        "geneData_orig" "geneData"      "granges_orig"  "granges"

5.7 Adding Items to a DGEobj

rownames(alignmentQC) <- alignmentQC$Metric

alignmentQC <- alignmentQC %>%
    select(-Metric) %>%
    t() %>%
    as.data.frame()

dgeObj <- addItem(dgeObj,
                  item = alignmentQC,
                  itemName = "AlignmentQC",
                  itemType = "alignQC")

# Perform TMM normalization on the raw counts
dgelist <- dgeObj$counts %>%
  DGEList() %>%
  calcNormFactors(method = "TMM")

# add the resulting edgeR DGEList item to the DGEobj
newdgeObj <- addItem(dgeObj,
                     item     = dgelist,
                     itemName = "normTMM",
                     itemType = "DGEList",
                     parent   = "counts",
                     funArgs  = "calcNormFactors; TMM"
)

kable(inventory(newdgeObj))

runNorm <- function(x, method){
  # Perform TMM normalization on the raw counts
  dgelist <- x$counts %>%
    DGEList() %>%
    calcNormFactors(method = method)

  # add the resulting edgeR DGEList item to the DGEobj
  x <- addItem(x,
               item     = dgelist,
               itemName = "normTMM",
               itemType = "DGEList",
               parent   = "counts",
               funArgs  = match.call()
  )
}

newdgeObj <- runNorm(dgeObj, method = "TMM")

#use verbose=TRUE to see the funArgs column
inv <- inventory(newdgeObj, verbose = TRUE)[ ,c(1:3,9)]
kable(inv)

# Perform TMM normalization on the raw counts
dgelist <- dgeObj$counts %>%
  DGEList() %>%
  calcNormFactors(method = "TMM")

# Formula must be composed of column names from the design table.
formula <- '~ 0 + ReplicateGroup'

# User-defined name for the designMatrix
designMatrixName <- "ReplicateGroupDesign"

# build the designMatrix
design <- getItem(dgeObj, "design")
designMatrix <- model.matrix (as.formula(formula), design)
#capture the formula as an attribute of the designMatrix
attr(designMatrix, "formula") <- formula

# Add both the dgelist and designMatrix to the DGEobj

NewDgeObj <- addItems(dgeObj,
                     itemList  = list(dgelist=dgelist, ReplicateGroupDesign=designMatrix),
                     itemTypes = list("DGEList", "designMatrix"),
                     parents   = list("counts", "design"))

5.8 Pre-defined Item Types

A number of pre-defined item types have been implemented to support the limma/voom DGE workflow for RNA-Seq and Affymetrix data types.

To see a complete list of pre-defined item types use the showTypes function.

kable(arrange(showTypes(dgeObj), BaseType, Type))

5.9 Adding a New Data Type

The capability to define new data types is core to the extensibility of the DGEobj data structure.

A set of data types has already been defined based on data types encountered in a typical limma/Voom workflow. The newTypes function is used to define data types to support a novel workflow.

# Add a new colData datatype called "sampleQC"
newDgeObj <- newType(dgeObj,
                    itemType = "sampleQC",
                    baseType = "col",
                    uniqueItem = FALSE)

After executing this call, the newDgeObj now contains the newly defined type definition and is ready to accept data items of that type.

6.1 Retrieve a single item

The getItem function allows the user to retrieve any item in a DGEobj by referencing its item name. Use the inventory function (described above) to check the contents of the DGEobj including the list of item names that are available.

counts <- getItem(dgeObj, "counts")

#peek at the first 4 columns
head(counts[ ,1:4])

##                    T_20170823MAN1_C05P01 T_20170823MAN1_H05P01 T_20170823MAN1_C04P01 T_20170823MAN1_D03P01
## ENSRNOG00000046319                     2                     0                 0.000                0.0000
## ENSRNOG00000047964                     0                     0                 0.000                0.0000
## ENSRNOG00000050370                     0                     0                 0.000                1.0000
## ENSRNOG00000032365                     0                     0                 0.000                0.0000
## ENSRNOG00000040300                     1                     6                 5.783                6.8135
## ENSRNOG00000058808                     0                     0                 0.000                0.0000

6.2 Retrieve multiple items

There are several ways to retrieve multiple items as a list or data objects.

The user can supply a list of item names to retrieve specific items. The items requested are returned in a list.

MyItems <- getItems(dgeObj, list("counts", "geneData"))

If all the items the user wants to retrieve are of the same type, the getType function may be used to retrieve them as a named list. For example, an analyst may want to retrieve a set of contrast results stored as topTable output.

# The below example object is a post-workflow object that has multiple contrasts
FullWorkflowDGEobj <- readRDS(system.file("exampleObj.RDS", package = "DGEobj"))

MyContrasts <- getType(FullWorkflowDGEobj, "topTable")
names(MyContrasts)

## [1] "BDL_vs_Sham"    "EXT1024_vs_BDL" "Nint_vs_BDL"    "Sora_vs_BDL"

Similarly, all items of the same baseType may be retrieved as a named list with the getBaseType function.

MyColData <- getBaseType(FullWorkflowDGEobj, "col")
names(MyColData)

## [1] "design"               "ReplicateGroupDesign"

Finally, the DGEobj can be recast as a named list.

MyList <- as.list(dgeObj)
class(MyList)

## [1] "list"

names(MyList)

## [1] "counts_orig"   "counts"        "design_orig"   "design"        "geneData_orig" "geneData"      "granges_orig"  "granges"       "AlignmentQC"

7.1 Function baseType

Returns the baseType for a given Type:

baseType(dgeObj, "DGEList")

## [1] "assay"

7.2 Querying Metadata Stored in Attributes

Certain metadata about a project can be stored in the attributes of the DGEobj. Similarly, each item deposited in a DGEobj can have it’s own set of attributes. Both of these types of attributes can be easily accessed with the base R attr function.

The showMeta function will return a named list of all attributes of a DGEobj.

To retrieve a specific attribute:

keywords <- attr(dgeObj, "Keywords")

To retrieve an attribute from an item, use dollar sign syntax to reference the item:

formula <- attr(FullWorkflowDGEobj$ReplicateGroupDesign, "formula")

7.3 Function rmItem

For common workflow usage, the DGEobj should be considered read-only. Items are added during the workflow and should not be modified further. Thus, under the normal intended usage, no items should be removed. Although, its use is discouraged, for completeness, a rmItem function exists.

To delete an item from a DGEobj:

MyDgeObj <- rmItem(dgeObj, "granges")

7.4 Function setAttributes

To add or update a single attribute, use the base Rattr function:

attr(dgeObj, "Tissue")  <- "Liver"

The setAttributes function allows attaching a named list of attributes. Unlike the base R attributes function, setAttributes will add or update the named attributes without deleting any existing attributes.

newDgeObj <- setAttributes(dgeObj, list(date=date(), analyst="JRT"))

7.5 Function showMeta

This function returns any project annotation attached as attributes to the DGEobj using function annotateDGEobj. It returns the attributes as a two column dataframe of Attribute/Value pairs. It only returns attributes with character values.

kable(showMeta(dgeObj))

8.1 Project Meta Data

User supplied meta data in the form of key/value pairs are attached to the DGEobj as attributes where the attribute name is the key and the attribute value is a character string. These metaData attributes can be loaded from a text file using function annotateDGEobj or can be manipulated individually with the base R attr function. The showMeta function displays these user supplied attributes.

8.2 Workflow attributes

Several attributes are used to capture information about each data item (“type”, “basetype”, “dateCreated” and “funArgs.”), as well as relationships between them (“parent”). The inventory function provides a convenient table view of these attributes. Operationally, the value of each of these attributes is a named list where the names are the item names and the value is the value for that item. For example, to query just the basetype attribute of the “design” itemName:

attr(dgeObj, "basetype")[["design"]]

## [1] "col"

8.3 Customizing the DGEobj Definition for New Data Types

A class DGEobjDef has been created to define the data types acceptable to a DGEobj and also allow for the definition of new datatypes as needed. A default DGEobjDef includes data types for gene expression (Affymetrix and RNA-Seq) applications. The user may define additional levels and data types when initializing a DGEobjDef object (initDGEobjDef function).

To return the default DGEobjDef:

defaultDGEobjDef <- initDGEobjDef()

The initDGEobj function uses a default DGEobj definition tailored for the limma/voom workflow. If new data types are needed to support a different workflow, the user should define the new data types during a call to initDGEobjDef, and before calling the DGEobj constructor (initDGEobj). The initDGEobjDef function returns a DGEobjDef class object that can then be passed explicitly to initDGEobj.

There are four optional arguments to initDGEobjDef:

• levels: A character string or vector of new level names.
• primaryAssayNames: A character string or vector of itemNames for the primary assay matrix for each level. Length must match length(levels).
• types: A named character vector of new data types. Names represent the type name and values define the basetype for this type.
• uniqueTypes: Adding a type name to this character vector restricts that data type to one instance per DGEobj.

The example below defines new types to accommodate metabolomics data.

# simulate some data
nr <- nrow(dgeObj)
nc <- ncol(dgeObj)

myMolecules <- dgeObj$geneData
design      <- dgeObj$design

myMSquant   <- rnorm(nr * nc) %>%
  exp() %>%
  matrix(nrow = nr, ncol = nc)

rownames(myMSquant) <- rownames(dgeObj)
colnames(myMSquant) <- colnames(dgeObj)

myDGEobjDef <- initDGEobjDef(levels            = "metabolite",
                             primaryAssayNames = "intensity",
                             types = c(normInt = "assay"))

MetabolomicObj <- initDGEobj(primaryAssayData = myMSquant,
                             rowData          = myMolecules,
                             colData          = design,
                             level            = "metabolite",
                             DGEobjDef        = myDGEobjDef)
kable(inventory(MetabolomicObj))

There are simple dimension-based criteria to determine the basetype of a new datatype.

• If both row and column dimensions match the primary assay data, the new object is an assay.
• If the row dimension of the new data matches the row dimension of the primary assay data, the new data type is a “row” basetype.
• If the row dimension of the new data matches the column dimension of the primary assay data, the new data type is “col” basetype.
• If the new data type does not match row or col dimensions, it is a “meta” basetype.

The DGEobj definitions created during initialization or added post-hoc by the newType function, are stored the “objDef” attribute of the DGEobj.

str(attr(dgeObj, "objDef"))

## List of 5
##  $ basetype         : Named chr [1:4] "row" "col" "assay" "meta"
##   ..- attr(*, "names")= chr [1:4] "row" "col" "assay" "meta"
##  $ type             : Named chr [1:42] "row" "col" "assay" "meta" ...
##   ..- attr(*, "names")= chr [1:42] "row" "col" "assay" "meta" ...
##  $ uniqueType       : chr [1:20] "counts" "counts_orig" "design" "design_orig" ...
##  $ allowedLevels    : chr [1:5] "gene" "exon" "isoform" "protein" ...
##  $ primaryAssayNames: Named chr [1:5] "counts" "counts" "intensities" "intensities" ...
##   ..- attr(*, "names")= chr [1:5] "gene" "exon" "isoform" "protein" ...
##  - attr(*, "class")= chr "DGEobjDef"