Customizing the differential analysis of omics data
Key data structures explained
A key component of the package is the data structures it provides. The following are a list of the key S4 classes:
MOList: AnMOListobject is singly responsible for storing all the omics data and the results of each step of the analysis. While the raw data of each omic type are stored internally in the slots, all analysis results are kept as list elements that allow easy subsetting and extraction. By exploiting this feature, we can customize our analysis by extracting and manipulating the results of the analysis.DETag: ADETagobject is the core S4 class that stores the results of the differential expression analysis. It provides the users with a composite structure for exploring the results of the analysis.TOPTag: ATOPTagobject is a subclass of theDETagobject. Unlike theDETagobject, theTOPtagobject further limits the results to a subset of genes that we are interested in. It provides an user-friendly interface for filtering and ranking top differentially expressed genes/transcripts.PEAKTag: Similar toTOPTag, aPEAKTagobject is also a subclass of theDETagobject. However, it is specifically designed for handling peak annotation and motif enrichment analysis.TRNet: Finally, as the final component of the traditional workflow, aTRNetobject is the core S4 class that stores the final TRN. Although itself operates as a data storage object, it can be easily integrated with any third-party tools for further topological analysis or customized network visualization.
Let us begin with the first step of the analysis pipeline, which is to perform differential analysis on each omics data separately. We begin by exploring how these data structures facilitate flexible and detailed analysis of the results.
Exploring the results of the differential analysis
We will first construct a MOList object
using the data we have prepared. We consider the age of the samples as
the grouping criteria for the differential analysis.
myMOList <- MOList(
RNAseq = RNAseq_heart,
RNAGroupBy = RNAseq_heart_samples$Age,
smallRNAseq = smallRNAseq_heart,
smallRNAGroupBy = smallRNAseq_heart_samples$Age,
proteomics = protein_heart,
proteomicsGroupBy = protein_heart_samples$Age,
ATACpeak1 = atacPeak1Path,
ATACpeak2 = atacPeak2Path
)
#> RNAseq: RNAseq data has been added.
#> Small RNAseq: Small RNAseq data has been added.
#> Proteomics: Proteomics data has been added.
#> ATACseq: ATAC peaks data has been added.With the package’s one-step differential analysis function, we can easily explore the results of the differential analysis.
myMOList <- diffOmics(myMOList, program = "edgeR")
#> Performing differential analysis on RNAseq data
#>
#> Performing differential analysis on smallRNAseq data
#>
#> Performing differential analysis on proteomics data
#>
#> Processing peaks...In the IntegraTRN package, several
visualization methods are provided to explore the results of the
differential analysis. These visualization methods are highly efficient
in generating an overview of the results, allowing users to quickly
identify the key changes that happen during the biological process.
However, as biologists who generated the data, we may be interested in
some particular aspects of the results. Therefore, we can further
customize our experience by performing easy extraction of the
results.
To obtain the DETag object for a
particular omics data, we can simply subset the
MOList object using the name of the omics
data. For example, to retrieve differential analysis result of the
RNAseq data:
would work just fine. The DETag class
has implemented several methods for users to easily extract the results.
For example, we can directly export the results as a data frame using
the exportDE function.
exportDE(myRNADETag) %>%
head(5)
#> logFC pvalue padj
#> SLC35F1 0.3301368 1.435088e-30 7.175439e-27
#> HMGA2 -0.3862208 1.161249e-29 2.463844e-26
#> CDH11 -0.2866885 1.478307e-29 2.463844e-26
#> FHL2 0.2719520 1.650073e-28 2.062591e-25
#> JAM2 0.1827767 4.298455e-27 3.991485e-24Furthermore, in certain conditions, such as generating a heatmap, we
may want to explore only the differentially expressed genes as well as a
normalized count matrix for visualization purposes. Besides subsetting
the data frame returned by the exportDE
function, we can also further utilize the functionality of its subclass,
the TOPTag object. As the name suggests,
the TOP part implies that we could easily
perform filtering to obtain what we consider to be the “top”
differentially expressed genes. For example, we can retrieve the top 20
differentially expressed genes with any fold change, with an adjusted
p-value less than 0.05 considered to be significant.
myRNATOPTag <- TOPTag(myRNADETag,
pCutoff = 0.05,
logFCCutoff = 0,
topGenes = 20,
direction = "both"
)Similarly, we can also export the results as a data frame.
Furthermore, we can see that a ranking based on
piValue is also provided.
exportDE(myRNATOPTag) %>%
head(5)
#> logFC pvalue padj piValue rank
#> SLC35F1 0.3301368 1.435088e-30 7.175439e-27 26.14415 1
#> FHL2 0.2719520 1.650073e-28 2.062591e-25 24.68559 2
#> JAM2 0.1827767 4.298455e-27 3.991485e-24 23.39887 3
#> LDHD 0.2818013 1.218915e-26 7.475262e-24 23.12637 4
#> RD3L 0.4080196 1.495052e-26 7.475262e-24 23.12637 5In case where we want to obtain a normalized count matrix for visualization, such as when plotting a heatmap, we can also do so by a simple function call.
The piValue here ranks the genes based
on the adjusted p-value and the fold change, which is a standard ranking
for rank-based enrichment analysis. The
piValue is calculated as follows:
\[piValue = -log_{10}(padj) * sign(logFC)\]
where \(padj\) is the adjusted p-value and \(logFC\) is the log fold change.
We can then make a heatmap using the
ComplexHeatmap package.
heatmapMatrix <- topRNANormCounts %>%
as.matrix() %>%
t() %>%
scale() %>%
t()
exprRange <- range(heatmapMatrix)
col <- circlize::colorRamp2(
c(exprRange[1], 0, exprRange[2]),
c("blue", "white", "red")
)
ComplexHeatmap::Heatmap(heatmapMatrix,
col = col,
show_column_names = FALSE
)
Figure 1. Heatmap of top 20 differentially expressed genes
Similarly, we can use the list of selected genes to perform gene set
enrichment analysis (GSEA) (Reimand et al.
2019). While a rank based GSEA analysis utilizes all available
genes, we can again utilize the TOPTag
constructor to perform ranking while not filtering any genes. This can
be done by setting a non-selective cutoff. Here, a
topGenes = 1 means that we are selecting all genes that
satisfy the cutoffs.
rankedResult <- TOPTag(myRNADETag,
pCutoff = 1,
logFCCutoff = 0,
topGenes = 1,
direction = "both"
) %>% exportDE()This can then be used by the users as import to the GSEA analysis.
Alternative differential analysis methods
For most highly controlled experiments, the differential analysis
between two conditions is fairly straightforward. However, in certain
cases, especially for patient-derived samples, a much higher biological
variability may be observed. How could we handle such variability? With
the DETag class and the list behavior of
the MOList object, this has been made
possible.
In the following section, we will explore an alternative differential
analysis with multiple covariates. However, any alternative differential
expression analysis can be performed at the user’s discretion, and the
results can be easily integrated into the
MOList object.
In the case of a population-based cohort, many attributes could be contributing to the biological variability. In our example, we can see that the following attributes are available:
In this case, in the differential expression analysis, we can adjust
for both sex and the batch effect. Again, we use the
edgeR package for the differential
expression analysis.
# Create the edgeR DGEList object
dgeList <- edgeR::DGEList(
counts = RNAseq_heart,
group = RNAseq_heart_samples$Age
)
# Filter lowly expressed genes
keep <- edgeR::filterByExpr(dgeList)
dgeList <- dgeList[keep, ]
# Perform differential expression analysis
design <- model.matrix(~ RNAseq_heart_samples$Age +
RNAseq_heart_samples$Sex +
RNAseq_heart_samples$Batch)
dgeList <- edgeR::calcNormFactors(dgeList, method = "TMM") # or a compatible method
dgeList <- edgeR::estimateDisp(dgeList, design)
fit <- edgeR::glmQLFit(dgeList, design) # or other suitable models
qlf <- edgeR::glmQLFTest(fit, coef = ncol(RNAseq_heart_samples))With complex designs, the performance of user-defined differential
expression analysis may be better than the one-step differential
expression analysis provided by the
IntegraTRN package. As we have obtained
the results of the differential expression analysis, integrating the
results back into the MOList object is
fairly straightforward.
To do so, we need to construct a DETag
object. The core component of the DETag
object is a result data frame, as the original output of
edgeR or
DESeq2. Additionally is a normalized count
matrix.
# Obtaining the DE results and normalized count
deResult <- edgeR::topTags(qlf, n = nrow(dgeList)) %>%
as.data.frame()
normCounts <- edgeR::cpm(dgeList) %>%
as.matrix()
# Construct the DETag object
myNewRNADETag <- DETag(DEResult = deResult, method = "edgeR", normalizedCounts = normCounts)With the DETag object, we can then
integrate the results back into the MOList
object.
This would further allow us to perform downstream analysis using the customized differential expression results.
For the purpose of this specific example, we will not replace the original analysis results.
While we have demonstrated the flexibility of the
IntegraTRN package in case of the RNAseq
data, other omics data can be handled in the same way. With these
flexibility, we can not only explore the differential analysis results
in a more detailed manner, but also customize the analysis that are
specific to our data and biological question of interest.