Differential analysis of omics data
The first step of the analysis pipeline is to perform differential analysis separately on each omic data. During the differential analysis, we can retrieve insights into the “changes” that happen during the biological process separately for each omic data. This exploratory analysis is essential to understand what alterations in the transcriptome can explain difference in the phenotype of interest, and what variations in any of the regulatory players could potentially contribute to these changes.
Generate MOList object using multi
omics data
To get started, we first need to generate a
MOList object. The
MOList object stands for
multi-omics list, which is the core of
entire analysis pipeline as it preserves our input data and the results
of the differential analysis in a centralized object. The
MOList object can be generated using the
MOList function.
The MOList constructor has a large
flexibility in terms of what input omics data are available. At the
minimum, when constructing a MOList object
for the first time, we need to provide the RNAseq data and the sample
grouping information. For any count-based omics data (RNAseq, small
RNAseq, and proteomics), the constructor takes a numeric matrix of
counts, where each row represents a gene or a protein, and each column
represents a sample. The grouping information can be any vector of the
same length as the number of samples, in the same order to the columns
of the count matrix. A continuous or a binary variable can be used.
However, we need to make sure we are supplying the same grouping
criteria for all omics data. Here we create a
myMOList using only RNAseq and small
RNAseq data, using the age of the samples
as the grouping criteria for downstream differential analysis.
myMOList <- MOList(
RNAseq = RNAseq_heart,
RNAGroupBy = RNAseq_heart_samples$Age,
smallRNAseq = smallRNAseq_heart,
smallRNAGroupBy = smallRNAseq_heart_samples$Age
)
#> RNAseq: RNAseq data has been added.
#> Small RNAseq: Small RNAseq data has been added.
#> Proteomics: Optional proteomics data has NOT been provided.
#> ATACseq: Optional ATAC peaks data has NOT been provided.As long as a MOList object is
constructed with RNAseq data, additional omics data can be added to the
existing object using the same syntax. Any existing data can also be
replaced by new data using the same syntax.
Specifically for ATACseq data, we need to provide the path to the
peak files instead of the count matrix. Since most peak callers output
peak files in BED
format, we can simply provide a path to a BED file, for each
condition. Our package has already provided two simulated peak files in
the extdata folder.
atacPeak1Path <- system.file("extdata", "peak1.bed", package = "IntegraTRN")
atacPeak2Path <- system.file("extdata", "peak2.bed", package = "IntegraTRN")A consensus peak file should be provided for each condition. When
biological replicates are used for the ATACseq experiment, the peak
files should be merged into one and filtered for blacklisted regions
before being provided to the constructor (Amemiya, Kundaje, and Boyle 2019). Here, we
simply use the previously created myMOList
object to append the ATACseq data.
Differential analysis of omics data
Once the MOList object is created, we
can perform differential analysis on each omics data separately.
However, since differential analysis on count-based data is different
that that of the epigenomic data, the process is slightly different.
Differential analysis of count-based omics data
For count-based omics data, we can use the
diffOmics function to perform differential
analysis. To simplify the process, the
diffOmics function essentially performs
differential expression analysis on all available count based omics data
in the MOList object. Since we can specify
to use either DESeq2 or
edgeR for the differential analysis, our
underlying assumption is that the expression data follows a negative
binomial distribution (Love, Huber, and Anders
2014) (Robinson, McCarthy, and Smyth
2010). Furthermore, the differential expression analysis utilizes
a generalized linear model (GLM) to allow adjustment on covariates. In
particular, the function allows adjustment on the batch effect, which is
a common issue in next-generation sequencing (NGS) data.
myMOList <- diffOmics(myMOList,
rnaseqBatch = RNAseq_heart_samples$Batch,
program = "DESeq2"
)
#> Performing differential analysis on RNAseq data
#>
#> Performing differential analysis on smallRNAseq data
#>
#> Processing peaks...The batch argument follow the same
criteria as the grouping information provided during the construction of
the MOList object. Therefore, this
argument is not limited to the batch effect, but can be any covariate
that we want to adjust for.
# NOT run, just illustrate the usage of the batch argument
myMOList <- diffOmics(myMOList,
smallRnaBatch = smallRNAseq_heart_samples$Sex,
program = "edgeR"
)Exploring differential gene expression analysis
To explore the results of the differential analysis, several
visualization methods are provided. First, to see the distribution of
differentially expressed mRNAs as well as exploring a optimal cutoff for
defining genes as differentially expressed, we can visualize the RNAseq
data using a volcano plot. We can specify the cutoffs for the adjusted
p-value and the log2 fold change (log2FC) to define differentially
expressed genes. Here, we define log2FC > 0.1 and
adjusted p-value < 0.05, which is the default, as the
cutoffs.
Figure 1. Volcano plot of RNAseq data
The plot also directly indicates the number of up- and down-regulated genes.
How do we efficiently explore the raw differential expression
results? The IntegraTRN package has
conveniently defined a set of data structures, primarily S4 classes and
methods, that allows users to easily explore the results. The previously
mentioned MOList object essentially works
as a list, which contains differential results stored in a
DETag object for each omics data. This
means that we can easily retrieve the differential expression results
via simple subsetting.
A simple snap shot of the DETag object
can be directly invoked.
myRNADETag
#> A DETag S4 object
#> Method for differential analysis: DESeq2
#>
#> Number of genes: 5000
#> Number of genes with adjusted p-value < 0.05: 3223
#>
#> A snap shot for differential analysis results:
#> baseMean log2FoldChange lfcSE stat pvalue
#> B4GALT3 958.13887 -0.0078640979 0.009959527 -0.78960552 4.297582e-01
#> SLC26A11 458.38254 -0.0019483551 0.009851067 -0.19778113 8.432163e-01
#> BEAN1 114.23147 -0.0006515378 0.020920562 -0.03114342 9.751552e-01
#> STAT3 2121.37501 -0.0343187743 0.008212107 -4.17904601 2.927345e-05
#> STEAP3 56.69544 0.0311278206 0.016709689 1.86286052 6.248189e-02
#> PIP5K1C 1552.90544 -0.0978374269 0.012828549 -7.62653887 2.411402e-14
#> padj
#> B4GALT3 4.955699e-01
#> SLC26A11 8.747057e-01
#> BEAN1 9.810414e-01
#> STAT3 8.402252e-05
#> STEAP3 9.057972e-02
#> PIP5K1C 3.367881e-13
#>
#> To access the full results, use the exportDE S4 methodTo further obtain the full results of the differential expression
analysis, we can use the exportDE
function. This function returns a data frame containing the outputs of
DESeq2 or
edgeR, with choice of using a consistent
interface.
exportDE(myRNADETag, original = TRUE) %>%
head(5)
#> baseMean log2FoldChange lfcSE stat pvalue
#> B4GALT3 958.13887 -0.0078640979 0.009959527 -0.78960552 4.297582e-01
#> SLC26A11 458.38254 -0.0019483551 0.009851067 -0.19778113 8.432163e-01
#> BEAN1 114.23147 -0.0006515378 0.020920562 -0.03114342 9.751552e-01
#> STAT3 2121.37501 -0.0343187743 0.008212107 -4.17904601 2.927345e-05
#> STEAP3 56.69544 0.0311278206 0.016709689 1.86286052 6.248189e-02
#> padj
#> B4GALT3 4.955699e-01
#> SLC26A11 8.747057e-01
#> BEAN1 9.810414e-01
#> STAT3 8.402252e-05
#> STEAP3 9.057972e-02By having original = TRUE, the data
frame remains exactly the same as the
DESeq2 output.
Some times, we may also want to rank the genes to perform gene set
enrichment analysis (GSEA). (Reimand et al.
2019) To this end, we can make use of the subclass of
DETag, which is the
TOPTag object. It automatically performs
piValue transformation (Reimand et al.
2019) and ranking based on a user-specified criteria. For
example, let’s retrieve genes that have at least 0.1 log2 fold change
and select the top 200 ranked genes.
myRNATOPTag <- TOPTag(myRNADETag,
pCutoff = 0.05,
logFCCutoff = 0.1,
topGenes = 200,
direction = "both"
)To see a snap shot:
myRNATOPTag
#> A TOPTag object.
#> The top 200 differential genes are selected.
#>
#> A snapshot of the top up-regulated genes:
#> logFC pvalue padj piValue rank
#> SLC35F1 0.3168227 6.830735e-61 1.707684e-57 56.76759 1
#> JAM2 0.1955125 2.053533e-51 2.566917e-48 47.59059 2
#> FHL2 0.2654241 6.220696e-49 6.220696e-46 45.20616 3
#> FBXW7 0.2408516 6.853848e-44 4.283655e-41 40.36819 4
#> CNOT11 0.1348388 4.780297e-43 2.655720e-40 39.57582 5
#>
#> A snapshot of the top down-regulated genes:
#> logFC pvalue padj piValue rank
#> RAI1 -0.1236128 4.680130e-10 3.175122e-09 -8.498240 -1
#> ATP8B2 -0.1022981 4.871540e-10 3.291581e-09 -8.482595 -2
#> SORCS2 -0.1413809 5.456246e-10 3.652106e-09 -8.437457 -3
#> PRRC2B -0.1897221 5.944744e-10 3.931709e-09 -8.405419 -4
#> SLC47A1 -0.1106573 6.664062e-10 4.384252e-09 -8.358105 -5
#>
#> Method for differential expression analysis: DESeq2
#> Log fold change cutoff: 0.1
#> Adjusted p-value cutoff: 0.05Due to the inheritance and object-oriented nature of the S4 class
system, we can handle the TOPTag object in
the same way as the DETag object.
# Export the desired results with rank information
myTopRNAResults <- exportDE(myRNATOPTag)
# Make numeric values easier to read
myTopRNAResults[, 1:4] <- format(myTopRNAResults[, 1:4], digits = 4, scientific = TRUE)
head(myTopRNAResults)
#> logFC pvalue padj piValue rank
#> SLC35F1 3.168e-01 6.831e-61 1.708e-57 5.677e+01 1
#> JAM2 1.955e-01 2.054e-51 2.567e-48 4.759e+01 2
#> FHL2 2.654e-01 6.221e-49 6.221e-46 4.521e+01 3
#> FBXW7 2.409e-01 6.854e-44 4.284e-41 4.037e+01 4
#> CNOT11 1.348e-01 4.780e-43 2.656e-40 3.958e+01 5
#> LDHD 2.558e-01 2.245e-41 9.356e-39 3.803e+01 6Here, a rank of 1 and -1 indicates the most up- and down-regulated gene.
For publication purposes, we can also label some genes that we are interested in on the volcano plot. Let’s first retrieve a few differentially expressed genes
genesToLabel <- exportDE(myRNADETag) %>%
dplyr::filter(padj < 0.05, abs(logFC) > 0.1) %>%
dplyr::arrange(desc(logFC)) %>%
rownames()
# Get an up- and a down- regulated gene
# For example use the 10th gene with highest fold change
genesToLabel <- genesToLabel[c(10, length(genesToLabel) - 10)]Then we again plot the volcano plot, with some genes labeled. We can
make use of the highlight argument to
specify the genes to label.
Figure 2. Volcano plot of RNAseq data with some genes labelled
Exploring differential small RNA expression analysis
Before moving on to small RNAseq data, we first need to annotate the small RNAs for their types. To do so, we need to provide a data frame of small RNA annotation, as in the name of the transcripts for each type of small RNA. The following is an example:
smallRNAAnnotation <- data.frame(
transcript = c("hsa-miR-1-3p", "hsa-miR-1-5p", "hsa-miR-2-5p"),
type = c("miRNA", "miRNA", "miRNA")
)To annotate the small RNAseq data, we can use the
annotateSmallRNA function.
# NOT run, just illustrate the usage of the smallRNAAnnotation argument
myMOList <- annotateSmallRNA(myMOList, anno = smallRNAAnnotation)We need to make sure that the names of the transcripts are consistent with the names in the count matrix, since different annotation databases may provide different names for the same transcript.
Alternatively, the package internally provides an annotation database for human small RNAs. This database is curated from several small RNA annotation databases, including miRBase, circBase, piRBase, piRNABank, GtRNAdb, GENCODE, and piRNACluster. Since we are dealing with human data, we can make use of this internal database.
The annotateSmallRNA function will
automatically checks for the coverage of each small RNA type and will
alert the user if the annotation does not cover all existing
transcripts.
Similar to the RNAseq data, we can visualize the distribution of differentially expressed small RNAs using a volcano plot.
Figure 3. Volcano plot of small RNAseq data
Here, we use the default P-value cutoff, which is 0.05.
In the mean time, the separately colors the up and down regulated small
RNAs, with the number of each annotated respectively.
On the other hand, small RNAs are usually a lot more complex and diverse than mRNAs. Each type of small RNA has its own function and mechanism of action. Therefore, we can also specifically color each type of small RNAs on the volcano plot, using a more specialized function.
By default, the plotVolcanoSmallRNA
function generates a volcano plot with each type of small RNA colored.
By default, it utilizes the BuPu palette
defined in the RColorBrewer package (Neuwirth 2022).
Figure 4. Volcano plot of small RNAseq data with each type of small RNA coloured
Alternatively, we can also specify a color scheme using the
colScheme argument.
# NOT run, just illustrate the usage of the colScheme argument
plotVolcanoSmallRNA(myMOList,
adjP = 0.05,
log2FC = 0,
colScheme = RColorBrewer::brewer.pal(6, "Set1"))Furthermore, although the function of microRNAs (miRNAs) is well studied, the understanding of other types of small RNAs such as piwi-interacting RNAs (piRNAs) and small nucleolar RNAs (snoRNAs) is still limited. Therefore, we can assess the contribution of each type of small RNA to the differences observed between samples. To do so, we can perform a principal component analysis (PCA) on each type of small RNA separately.
The plotSmallRNAPCAs function returns a
list of PCA plots, one for each type of small RNA. We can check the
length of the list to see how many PCA plots are generated.
In this case, we have 6 PCA plots, one for each small RNA type. We can plot the PCA plot for miRNA using the following command:
Figure 5. PCA plot of miRNA
Or, we can generate a combined PCA plot for all small RNA types.
pcaPlotList$miRNA + pcaPlotList$circRNA + pcaPlotList$piRNA + pcaPlotList$snoRNA +
pcaPlotList$snRNA + pcaPlotList$tRNA
Figure 6. Combined PCA plot of all small RNA types
Differential analysis of chromatin accessibility
Analysis of epigenomic data requires several more steps than that of the count-based omics data. In the case of ATACseq data, genomic regions are profiled for their accessibility to the hyperactive Tn5 transposase. This accessibility corresponds to the chromatin state of the region, which can be either open or closed for DNA binding proteins, particularly transcription factors (TFs). Therefore, to analyze differential chromatin accessibility, we need to annotate the peaks with the genes that are nearby. Furthermore, since we are interested in transcriptional regulation, we can parse through the open chromatin regions and identify the TF binding motifs that are enriched in these regions. Differentially enriched motifs can then be used to infer the TFs that are potentially involved in the regulation of the differentially expressed genes.
To do so, we first need additional annotation information. Here, we
make use of the annotation database from the
TxDb.Hsapiens.UCSC.hg38.knownGene package,
which is a TxDb object that contains the
genomic coordinates of all known genes in the hg38 genome
assembly (Team and Maintainer 2023). We
also need to provide a more general annotation database, which is the
org.Hs.eg.db package to extract the gene
information (Carlson 2023). Finally,
scanning the open chromatin regions for TF binding motifs requires
identifying the sequence of the open chromatin regions. Here, we use the
corresponding genome FASTA database from the
BSgenome.Hsapiens.UCSC.hg38 package (Pagès 2023).
# Annotate ATAC Peaks
library("TxDb.Hsapiens.UCSC.hg38.knownGene")
library("org.Hs.eg.db")
library("BSgenome.Hsapiens.UCSC.hg38")
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
annoDb <- "org.Hs.eg.db"
bsgenome <- BSgenome.Hsapiens.UCSC.hg38Motif enrichment analysis would also require a known set of position
weight matrices (PWMs) for all known curated TFs. The PWM is a matrix
that contains the probability of each nucleotide (A, T, C, G) at every
position of the binding motif. The
IntegraTRN package has already provided a
list of PWMs from the JASPAR2022 database
for vertebrates. We can load the PWMs using the following command:
This is a PWMatrixList object defined
in the TFBSTools package (Tan and Lenhard 2016). This package also
provides additional functionality to retrieve PWMs from other databases
or tools, such as HOMER etc. Here, we will use our package-provided PWMs
from JASPAR2022 (Castro-Mondragon et al. 2022).
The annotateATACPeaksMotif function
performs both peak annotation and motif enrichment analysis via a single
call. Here, we define the promoter region as the region
3kb upstream and 3kb downstream of the
transcription start site.
myMOList <- annotateATACPeaksMotif(myMOList,
tssRegion = c(-3000, 3000),
TxDb = txdb,
annoDb = annoDb,
bsgenome = bsgenome,
pwmL = jasparVertebratePWM
)Exploring differential chromatin accessibility
After performing peak annotation and motif enrichment analysis, we can explore the analysis via several methods.
First, let us take a look at the annotation on the peaks for their genomic features. This pie chart shows the percentage of peaks that are annotated as each genomic feature.
Figure 7. Annotation of ATAC peaks
To further access the quality of the ATACseq experiment, we can plot the coverage of the ATACseq peaks on the genome. This plot shows the where the peaks are located in each chromosome.
Figure 8. ATACseq coverage
On the left of the plots shows all differentially enriched motifs.
Beyond globally exploring the ATACseq data with the above plots, we
can also inspect the annotation on each individual peak. Once again, the
IntegraTRN package utilizes the S4
inheritance system to create a PEAKTag
object once motif enrichment analysis is performed. This object
implements the generic as.data.frame
method, conveniently allowing us to export the results as a data frame
to see where each individual peaks are located.
myMotifs <- as.data.frame(myMOList$DEATAC)
head(myMotifs[, 1:8]) # Many annotations, see the first 8 attributes
#> seqnames start end width strand Condition annotation geneChr
#> 1 chr1 1597754 1597883 130 * - Promoter (2-3kb) 1
#> 2 chr1 7783291 7783789 499 * - Promoter (<=1kb) 1
#> 3 chr1 7783791 7784289 499 * - Promoter (<=1kb) 1
#> 4 chr1 7784291 7784789 499 * - Promoter (<=1kb) 1
#> 5 chr1 9996206 9996704 499 * - Promoter (<=1kb) 1
#> 6 chr1 9996706 9997204 499 * - Promoter (<=1kb) 1Then we can explore the key part of the analysis, which is the
differentially enriched motifs. Typically, the differential motif
enrichment can be very stringent especially when the ATACseq libraries
are not sequenced very deeply (less than 100 million reads) (Yan et al. 2020). Therefore, we can use a more
relaxed cutoff for the p-value. The
plotATACMotifHeatmap function by default
utilizes adjusted p-value, but when a unadjusted p-value is provided, it
will use the unadjusted one and ignore the default. Since our ATACseq
data is simulated, here we use a relaxed cutoff.
Figure 9. Motif enrichment on ATACseq peaks
