Using eisaR for Exon-Intron Split Analysis (EISA)
2024-10-29
1 Introduction
Exon-Intron Split Analysis has been described by Gaidatzis et al. (2015). It consists of separately quantifying exonic and intronic alignments in RNA-seq data, in order to measure changes in mature RNA and pre-mRNA reads across different experimental conditions. We have shown that this allows quantification of transcriptional and post-transcriptional regulation of gene expression.
The eisaR package contains convenience functions to facilitate the steps in an
exon-intron split analysis, which consists of:
1. preparing the annotation (exonic and gene body coordinate ranges, section 3)
2. quantifying RNA-seq alignments in exons and introns (sections 4.1 and 4.2)
3. calculating and comparing exonic and intronic changes across conditions (section 5)
4. visualizing the results (section 6)
For the steps 1. and 2. above, this vignette makes use of Bioconductor annotation and the QuasR package. It is also possible to obtain count tables for exons and introns using some other pipeline or approach, and directly start with step 3.
2 Installation
To install the eisaR package, start R and enter:
# BiocManager is needed to install Bioconductor packages
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
# Install eisaR
BiocManager::install("eisaR")3 Preparing the annotation
As mentioned, eisaR uses gene annotations from Bioconductor.
They are provided in the form of TxDb or EnsDb objects, e.g. via packages such as TxDb.Mmusculus.UCSC.mm10.knownGene or EnsDb.Hsapiens.v86.
You can see available annotations using the following code:
pkgs <- c(BiocManager::available("TxDb")
BiocManager::available("EnsDb"))If you would like to use an alternative source of gene annotations, you might
still be able to use eisaR by first converting your annotations into a TxDb
or an EnsDb (for creating a TxDb see makeTxDb in the txdbmaker
package, for creating an EnsDb see makeEnsembldbPackage in the ensembldb
package).
For this example, eisaR contains a small TxDb to illustrate how regions are extracted.
We will load it from a file. Alternatively, the object would be loaded using library(...),
for example using library(TxDb.Mmusculus.UCSC.mm10.knownGene).
# load package
library(eisaR)
# get TxDb object
txdbFile <- system.file("extdata", "hg19sub.sqlite", package = "eisaR")
txdb <- AnnotationDbi::loadDb(txdbFile)Exon and gene body regions are then extracted from the TxDb:
# extract filtered exonic and gene body regions
regS <- getRegionsFromTxDb(txdb = txdb, strandedData = TRUE)
#> extracting exon coordinates
#> total number of genes/exons: 12/32
#> removing overlapping/single-exon/ambiguous genes (8)
#> creating filtered regions for 4 genes (33.3%) with 20 exons (62.5%)
regU <- getRegionsFromTxDb(txdb = txdb, strandedData = FALSE)
#> extracting exon coordinates
#> total number of genes/exons: 12/32
#> removing overlapping/single-exon/ambiguous genes (9)
#> creating filtered regions for 3 genes (25%) with 17 exons (53.1%)
lengths(regS)
#> exons genebodies
#> 20 4
lengths(regU)
#> exons genebodies
#> 17 3
regS$exons
#> GRanges object with 20 ranges and 0 metadata columns:
#> seqnames ranges strand
#> <Rle> <IRanges> <Rle>
#> ENSG00000078808 chr1 17278-18194 -
#> ENSG00000078808 chr1 18828-21741 -
#> ENSG00000078808 chr1 23614-23747 -
#> ENSG00000078808 chr1 24202-24358 -
#> ENSG00000078808 chr1 27799-27854 -
#> ... ... ... ...
#> ENSG00000186891 chr1 5740-6070 -
#> ENSG00000186891 chr1 6755-7081 -
#> ENSG00000254999 chr3 2266-2513 +
#> ENSG00000254999 chr3 12300-12402 +
#> ENSG00000254999 chr3 12943-13884 +
#> -------
#> seqinfo: 3 sequences from an unspecified genomeAs you can see, the filtering procedure removes slightly more genes for unstranded data
(strandedData = FALSE), as overlapping genes cannot be discriminated even if
they reside on opposite strands.
You can also export the obtained regions into files. This may be useful if
you plan to align and/or quantify reads outside of R. For example, you can use
rtracklayer to export the regions in regS into .gtf files:
library(rtracklayer)
export(regS$exons, "hg19sub_exons_stranded.gtf")
export(regS$genebodies, "hg19sub_genebodies_stranded.gtf")4 Quantify RNA-seq alignments in exons and introns
For this example we will use the QuasR package for indexing and alignment of short reads, and a small RNA-seq dataset that is contained in that package. As mentioned, it is also possible to align or also quantify your reads using an alternative aligner/counter, and skip over these steps. For more details about the syntax, we refer to the documentation and vignette of the QuasR package.
4.1 Align reads
Let’s first copy the sample data from the QuasR package to the
current working directory, all contained in a folder named extdata:
library(QuasR)
#> Loading required package: parallel
#> Loading required package: Rbowtie
file.copy(system.file(package = "QuasR", "extdata"), ".", recursive = TRUE)
#> [1] TRUEWe next align the reads to a mini-genome (fasta file extdata/hg19sub.fa) using
qAlign. The sampleFile specifies the samples we want to use, and the paths
to the respective fastq files.
sampleFile <- "extdata/samples_rna_single.txt"
## Display the structure of the sampleFile
read.delim(sampleFile)
#> FileName SampleName
#> 1 rna_1_1.fq.bz2 Sample1
#> 2 rna_1_2.fq.bz2 Sample1
#> 3 rna_2_1.fq.bz2 Sample2
#> 4 rna_2_2.fq.bz2 Sample2
## Perform the alignment
proj <- qAlign(sampleFile = sampleFile,
genome = "extdata/hg19sub.fa",
aligner = "Rhisat2", splicedAlignment = TRUE)
#> Creating .fai file for: /tmp/Rtmpk4qUHb/Rbuildb2a7a5fbd25a6/eisaR/vignettes/extdata/hg19sub.fa
#> alignment files missing - need to:
#> create alignment index for the genome
#> create 4 genomic alignment(s)
#> Creating an Rhisat2 index for /tmp/Rtmpk4qUHb/Rbuildb2a7a5fbd25a6/eisaR/vignettes/extdata/hg19sub.fa
#> Finished creating index
#> Testing the compute nodes...OK
#> Loading QuasR on the compute nodes...preparing to run on 1 nodes...done
#> Available cores:
#> nebbiolo1: 1
#> Performing genomic alignments for 4 samples. See progress in the log file:
#> /tmp/Rtmpk4qUHb/Rbuildb2a7a5fbd25a6/eisaR/vignettes/QuasR_log_b320710f570d6.txt
#> Genomic alignments have been created successfully
alignmentStats(proj)
#> seqlength mapped unmapped
#> Sample1:genome 95000 5961 43
#> Sample2:genome 95000 5914 864.2 Count alignments in exons and gene bodies
Alignments in exons and gene bodies can now be counted using qCount and the
regU that we have generated earlier (assuming that the data is unstranded).
Intronic counts can then be obtained from the difference between gene bodies and
exons:
cntEx <- qCount(proj, regU$exons, orientation = "any")
#> counting alignments...done
#> collapsing counts by sample...done
#> collapsing counts by query name...done
cntGb <- qCount(proj, regU$genebodies, orientation = "any")
#> counting alignments...done
#> collapsing counts by sample...done
cntIn <- cntGb - cntEx
cntEx
#> width Sample1 Sample2
#> ENSG00000078808 4837 705 1065
#> ENSG00000186827 1821 37 8
#> ENSG00000186891 1470 26 2
cntIn
#> width Sample1 Sample2
#> ENSG00000078808 10307 5 15
#> ENSG00000186827 1012 3 0
#> ENSG00000186891 1734 3 0As mentioned, both alignments and counts can also be obtained using alternative approaches. It is required that the two resulting exon and intron count tables have identical structure (genes in rows, samples in columns, the same order of rows and columns in both tables).
4.3 Load full count tables
The above example only contains very few genes. For the rest of the vignette,
we will use count tables from a real RNA-seq experiment that are provided in the
eisaR package. The counts correspond to the raw data used in Figure 3a of Gaidatzis et al. (2015)
and are also available online from the supplementary material:
cntEx <- readRDS(system.file("extdata",
"Fig3abc_GSE33252_rawcounts_exonic.rds",
package = "eisaR"))
cntIn <- readRDS(system.file("extdata",
"Fig3abc_GSE33252_rawcounts_intronic.rds",
package = "eisaR"))5 Run EISA conveniently
All the further steps in exon-intron split analysis can now be performed using
a single function runEISA. If you prefer to perform the analysis step-by-step,
you can skip now to section 7.
# remove "width" column
Rex <- cntEx[, colnames(cntEx) != "width"]
Rin <- cntIn[, colnames(cntIn) != "width"]
# create condition factor (contrast will be TN - ES)
cond <- factor(c("ES", "ES", "TN", "TN"))
# run EISA
res <- runEISA(Rex, Rin, cond)
#> filtering quantifyable genes...keeping 11759 from 20288 (58%)
#> fitting statistical model...done
#> calculating log-fold changes...done5.1 Alternative implementations of EISA
There are six arguments in runEISA (modelSamples, geneSelection, effects,
statFramework, pscnt and sizeFactor) that control gene filtering,
calculation of contrasts and the statistical method used, summarized in the
bullet list below.
The default values of these arguments correspond to the currently recommended way
of running EISA. You can also run EISA exactly as it was described by Gaidatzis et al. (2015), by
setting method = "Gaidatzis2015". This will override the values of the six
other arguments and set them according to the published algorithm (as indicated
below).
modelSamples: Account for individual samples in statistical model? Possible values are:FALSE(method="Gaidatzis2015"): use a model of the form~ condition * region
TRUE(default): use a model adjusting for the baseline differences among samples, and with condition-specific region effects (similar to the model described in section 3.5 of the edgeR user guide)
geneSelection: How to select detected genes. Possible values are:"filterByExpr"(default): First, counts are normalized usingedgeR::calcNormFactors, treating intronic and exonic counts as individual samples. Then, theedgeR::filterByExprfunction is used with default parameters to select quantifiable genes."none": This will use all the genes provided in the count tables, assuming that an appropriate selection of quantifiable genes has already been done.
"Gaidatzis2015"(method="Gaidatzis2015"): First, intronic and exonic counts are linearly scaled to the mean library size (estimated as the sum of all intronic or exonic counts, respectively). Then, quantifiable genes are selected as the genes with countsxthat fulfilllog2(x + 8) > 5in both exons and introns.
statFramework: The framework withinedgeRthat is used for the statistical analysis. Possible values are:"QLF"(default): quasi-likelihood F-test usingedgeR::glmQLFitandedgeR::glmQLFTest. This framework is highly recommended as it gives stricter error rate control by accounting for the uncertainty in dispersion estimation.
"LRT"(method="Gaidatzis2015"): likelihood ratio test usingedgeR::glmFitandedgeR::glmLRT.
effects: How the effects (log2 fold-changes) are calculated. Possible values are:"predFC"(default): Fold-changes are calculated using the fitted model withedgeR::predFCand the value provided topscnt. Please note that if a sample factor is included in the statistical model (modelSamples=TRUE), effects cannot be obtained from that model. In that case, effects are obtained from a simpler model without sample effects."Gaidatzis2015"(method="Gaidatzis2015"): Fold-changes are calculated using the formulalog2((x + pscnt)/(y + pscnt)). Ifpscntis not set to 8,runEISAwill warn that this deviates from the method used in Gaidatzis et al., 2015.
pscnt: The pseudocount that is added to normalized counts before log transformation. ForgeneSelection="Gaidatzis2015",pscntis used both in gene selection as well as in the calculation of log2 fold-changes. Otherwise,pscntis only used in the calculation of log2 fold-changes inedgeR::predFC(, prior.count = pscnt).sizeFactor: How size factors (TMM normalization factors and library sizes) are calculated and used withineisaR:"exon"(default): Size factors are calculated for exonic counts and reused for the corresponding intronic counts."intron": Size factors are calculated for intronic counts and reused for the corresponding exonic counts."individual"(method="Gaidatzis2015"): Size factors are calculated independently for exonic and intronic counts.
While different values for these arguments typically yield similar results,
the defaults are often less stringent compared to method="Gaidatzis2015" when
selecting quantifiable genes, but more stringent when calling significant changes
(especially with low numbers of replicates).
Here is an illustration of how the results differ between method="Gaidatzis2015" and
the defaults:
res1 <- runEISA(Rex, Rin, cond, method = "Gaidatzis2015")
#> setting parameters according to Gaidatzis et al., 2015
#> filtering quantifyable genes...keeping 8481 from 20288 (41.8%)
#> fitting statistical model...done
#> calculating log-fold changes...done
res2 <- runEISA(Rex, Rin, cond)
#> filtering quantifyable genes...keeping 11759 from 20288 (58%)
#> fitting statistical model...done
#> calculating log-fold changes...done
# number of quantifiable genes
nrow(res1$DGEList)
#> [1] 8481
nrow(res2$DGEList)
#> [1] 11759
# number of genes with significant post-transcriptional regulation
sum(res1$tab.ExIn$FDR < 0.05)
#> [1] 469
sum(res2$tab.ExIn$FDR < 0.05)
#> [1] 139
# method="Gaidatzis2015" results contain most of default results
summary(rownames(res2$contrasts)[res2$tab.ExIn$FDR < 0.05] %in%
rownames(res1$contrasts)[res1$tab.ExIn$FDR < 0.05])
#> Mode FALSE TRUE
#> logical 46 93
# comparison of deltas
ids <- intersect(rownames(res1$DGEList), rownames(res2$DGEList))
cor(res1$contrasts[ids,"Dex"], res2$contrasts[ids,"Dex"])
#> [1] 0.989731
cor(res1$contrasts[ids,"Din"], res2$contrasts[ids,"Din"])
#> [1] 0.9893341
cor(res1$contrasts[ids,"Dex.Din"], res2$contrasts[ids,"Dex.Din"])
#> [1] 0.9673155
plot(res1$contrasts[ids,"Dex.Din"], res2$contrasts[ids,"Dex.Din"], pch = "*",
xlab = expression(paste(Delta, "exon", -Delta, "intron for method='Gaidatzis2015'")),
ylab = expression(paste(Delta, "exon", -Delta, "intron for default parameters")))
5.2 On the estimation of interactions in a split-plot design experiment
The calculation of the significance of interactions (here whether the fold-changes differ between exonic or intronic data) is well defined for experimental designs where all samples are independent from one another. Within EISA, this is not the case (each sample yields two data points, one for exons and one for introns). That results in a dependency between data points: If a sample is affected by a problem in the experiment, it might at the same time give rise to outlier values in both exonic and intronic counts.
In statistics, such an experimental design is often referred to as a split-plot
design, and a recommended way to analyze interactions in such experiments would
be to use a mixed effect model with the plot (in our case, the sample) as a random
effect. The disadvantage here however would be that available packages for mixed
effect models are not designed for count data, and we therefore use an alternative
approach to explicitly model the sample dependency, by introducing sample-specific
columns into the design matrix (for modelSamples=TRUE). That sample factor is
nested in the condition factor (no sample can belong to more than one condition).
Thus, we are in the situation described in section 3.5 (‘Comparisons both between and
within subjects’) of the edgeR user guide, and we use the approach
described there to define a design matrix with sample-specific baseline effects
as well as condition-specific region effects.
This has no impact on the effects (the log2 fold-changes of modelSamples=TRUE
and modelSamples=FALSE are nearly identical). However, in the presence of sample effects,
modelSamples=TRUE increases the sensitivity of detecting genes with significant
interactions. Here is a comparison of the EISA results with and without accounting
for the sample in the model:
res3 <- runEISA(Rex, Rin, cond, modelSamples = FALSE)
#> filtering quantifyable genes...keeping 11034 from 20288 (54.4%)
#> fitting statistical model...done
#> calculating log-fold changes...done
res4 <- runEISA(Rex, Rin, cond, modelSamples = TRUE)
#> filtering quantifyable genes...keeping 11759 from 20288 (58%)
#> fitting statistical model...done
#> calculating log-fold changes...done
ids <- intersect(rownames(res3$contrasts), rownames(res4$contrasts))
# number of genes with significant post-transcriptional regulation
sum(res3$tab.ExIn$FDR < 0.05)
#> [1] 5
sum(res4$tab.ExIn$FDR < 0.05)
#> [1] 139
# modelSamples=TRUE results are a super-set of
# modelSamples=FALSE results
summary(rownames(res3$contrasts)[res3$tab.ExIn$FDR < 0.05] %in%
rownames(res4$contrasts)[res4$tab.ExIn$FDR < 0.05])
#> Mode TRUE
#> logical 5
# comparison of contrasts
diag(cor(res3$contrasts[ids, ], res4$contrasts[ids, ]))
#> Dex Din Dex.Din
#> 0.9931259 0.9872635 0.9912837
plot(res3$contrasts[ids, 3], res4$contrasts[ids, 3], pch = "*",
xlab = "Interaction effects for modelSamples=FALSE",
ylab = "Interaction effects for modelSamples=TRUE")
# comparison of interaction significance
plot(-log10(res3$tab.ExIn[ids, "FDR"]), -log10(res4$tab.ExIn[ids, "FDR"]), pch = "*",
xlab = "-log10(FDR) for modelSamples=FALSE",
ylab = "-log10(FDR) for modelSamples=TRUE")
abline(a = 0, b = 1, col = "gray")
legend("bottomright", "y = x", bty = "n", lty = 1, col = "gray")
6 Visualize EISA results
We can now visualize the results by plotting intronic changes versus exonic changes (genes with significant interactions, which are likely to be post-transcriptionally regulated, are color coded):
plotEISA(res)
#> identified 139 genes to highlight
7 Run EISA step-by-step
As an alternative to runEISA (section 5) and plotEISA
(section 6) described above, you can also analyze the data step-by-step
as described in Gaidatzis et al. (2015). This may be preferable to understand each
individual step and be able to access intermediate results.
The results obtained in this way are identical to what you get with
runEISA(..., method = "Gaidatzis2015"), so if you are happy with runEISA you can
skip the rest of the vignette.
7.1 Normalization
Normalization is performed separately on exonic and intronic counts, assuming that varying exon over intron ratios between samples are of technical origin.
# remove column "width"
Rex <- cntEx[,colnames(cntEx) != "width"]
Rin <- cntIn[,colnames(cntIn) != "width"]
Rall <- Rex + Rin
fracIn <- colSums(Rin)/colSums(Rall)
summary(fracIn)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.2696 0.2977 0.3105 0.3459 0.3587 0.4929
# scale counts to the mean library size,
# separately for exons and introns
Nex <- t(t(Rex) / colSums(Rex) * mean(colSums(Rex)))
Nin <- t(t(Rin) / colSums(Rin) * mean(colSums(Rin)))
# log transform (add a pseudocount of 8)
NLex <- log2(Nex + 8)
NLin <- log2(Nin + 8)7.2 Identify quantifiable genes
Genes with very low counts in either exons or introns are better removed, as they will be too noisy to yield useful results. We use here a fixed cut-off on the normalized, log-transformed intron and exonic counts:
quantGenes <- rownames(Rex)[ rowMeans(NLex) > 5.0 & rowMeans(NLin) > 5.0 ]
length(quantGenes)
#> [1] 84817.3 Calculate \(\Delta I\), \(\Delta E\) and \(\Delta E - \Delta I\)
The count tables were obtained from a total RNA-seq experiments in mouse embryonic stem (MmES) cells and derived terminal neurons (MmTN), with two replicates for each condition.
We will now calculate the changes between neurons and ES cells in introns (\(\Delta I\)), in exons (\(\Delta E\)), and the difference between the two (\(\Delta E - \Delta I\)):
Dex <- NLex[,c("MmTN_RNA_total_a","MmTN_RNA_total_b")] - NLex[,c("MmES_RNA_total_a","MmES_RNA_total_b")]
Din <- NLin[,c("MmTN_RNA_total_a","MmTN_RNA_total_b")] - NLin[,c("MmES_RNA_total_a","MmES_RNA_total_b")]
Dex.Din <- Dex - Din
cor(Dex[quantGenes,1], Dex[quantGenes,2])
#> [1] 0.9379397
cor(Din[quantGenes,1], Din[quantGenes,2])
#> [1] 0.8449252
cor(Dex.Din[quantGenes,1], Dex.Din[quantGenes,2])
#> [1] 0.5518187Both exonic and intronic changes are correlated across replicates, and so are the differences, indicating a reproducible signal for post-transcriptional regulation.
7.4 Statistical analysis
Finally, we use the replicate measurement in the edgeR framework to calculate the significance of the changes:
# create DGEList object with exonic and intronic counts
library(edgeR)
#> Loading required package: limma
#>
#> Attaching package: 'limma'
#> The following object is masked from 'package:BiocGenerics':
#>
#> plotMA
cnt <- data.frame(Ex = Rex, In = Rin)
y <- DGEList(counts = cnt, genes = data.frame(ENTREZID = rownames(cnt)))
# select quantifiable genes and normalize
y <- y[quantGenes, ]
y <- calcNormFactors(y)
# design matrix with interaction term
region <- factor(c("ex","ex","ex","ex","in","in","in","in"), levels = c("in", "ex"))
cond <- rep(factor(c("ES","ES","TN","TN")), 2)
design <- model.matrix(~ region * cond)
rownames(design) <- colnames(cnt)
design
#> (Intercept) regionex condTN regionex:condTN
#> Ex.MmES_RNA_total_a 1 1 0 0
#> Ex.MmES_RNA_total_b 1 1 0 0
#> Ex.MmTN_RNA_total_a 1 1 1 1
#> Ex.MmTN_RNA_total_b 1 1 1 1
#> In.MmES_RNA_total_a 1 0 0 0
#> In.MmES_RNA_total_b 1 0 0 0
#> In.MmTN_RNA_total_a 1 0 1 0
#> In.MmTN_RNA_total_b 1 0 1 0
#> attr(,"assign")
#> [1] 0 1 2 3
#> attr(,"contrasts")
#> attr(,"contrasts")$region
#> [1] "contr.treatment"
#>
#> attr(,"contrasts")$cond
#> [1] "contr.treatment"
# estimate model parameters
y <- estimateDisp(y, design)
fit <- glmFit(y, design)
# calculate likelihood-ratio between full and reduced models
lrt <- glmLRT(fit)
# create results table
tt <- topTags(lrt, n = nrow(y), sort.by = "none")
head(tt$table[order(tt$table$FDR, decreasing = FALSE), ])
#> ENTREZID logFC logCPM LR PValue FDR
#> 14680 14680 6.374952 6.554051 98.12387 3.930119e-23 3.333134e-19
#> 75209 75209 5.339465 6.400361 89.61927 2.886985e-21 1.224226e-17
#> 93765 93765 3.849839 6.603142 52.47425 4.359257e-13 1.232362e-09
#> 17957 17957 4.342526 6.864176 51.81480 6.099022e-13 1.293145e-09
#> 268354 268354 9.855437 8.402066 50.71845 1.066128e-12 1.808366e-09
#> 19276 19276 5.164777 8.391296 47.65570 5.080440e-12 6.488859e-097.5 Visualize the results
Finally, we visualize the results by plotting intronic changes versus exonic changes (genes with significant interactions, which are likely to be post-transcriptionally regulated, are color coded):
sig <- tt$table$FDR < 0.05
sum(sig)
#> [1] 509
sig.dir <- sign(tt$table$logFC[sig])
cols <- ifelse(sig, ifelse(tt$table$logFC > 0, "#E41A1CFF", "#497AB3FF"), "#22222244")
# volcano plot
plot(tt$table$logFC, -log10(tt$table$FDR), col = cols, pch = 20,
xlab = expression(paste("RNA change (log2 ",Delta,"exon/",Delta,"intron)")),
ylab = "Significance (-log10 FDR)")
abline(h = -log10(0.05), lty = 2)
abline(v = 0, lty = 2)
text(x = par("usr")[1] + 3 * par("cxy")[1], y = par("usr")[4], adj = c(0,1),
labels = sprintf("n=%d", sum(sig.dir == -1)), col = "#497AB3FF")
text(x = par("usr")[2] - 3 * par("cxy")[1], y = par("usr")[4], adj = c(1,1),
labels = sprintf("n=%d", sum(sig.dir == 1)), col = "#E41A1CFF")
# Delta I vs. Delta E
plot(rowMeans(Din)[quantGenes], rowMeans(Dex)[quantGenes], pch = 20, col = cols,
xlab = expression(paste(Delta,"intron (log2 TN/ES)")),
ylab = expression(paste(Delta,"exon (log2 TN/ES)")))
legend(x = "bottomright", bty = "n", pch = 20, col = c("#E41A1CFF","#497AB3FF"),
legend = sprintf("%s (%d)", c("Up","Down"), c(sum(sig.dir == 1), sum(sig.dir == -1))))
8 Session information
The output in this vignette was produced under:
sessionInfo()
#> R Under development (unstable) (2024-10-21 r87258)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.21-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] parallel stats4 stats graphics grDevices utils datasets
#> [8] methods base
#>
#> other attached packages:
#> [1] edgeR_4.5.0 limma_3.63.0 QuasR_1.47.0
#> [4] Rbowtie_1.47.0 rtracklayer_1.67.0 GenomicFeatures_1.59.0
#> [7] AnnotationDbi_1.69.0 Biobase_2.67.0 GenomicRanges_1.59.0
#> [10] GenomeInfoDb_1.43.0 IRanges_2.41.0 S4Vectors_0.45.0
#> [13] BiocGenerics_0.53.0 eisaR_1.19.0 BiocStyle_2.35.0
#>
#> loaded via a namespace (and not attached):
#> [1] DBI_1.2.3 bitops_1.0-9
#> [3] deldir_2.0-4 httr2_1.0.5
#> [5] biomaRt_2.63.0 rlang_1.1.4
#> [7] magrittr_2.0.3 Rhisat2_1.23.0
#> [9] matrixStats_1.4.1 compiler_4.5.0
#> [11] RSQLite_2.3.7 png_0.1-8
#> [13] vctrs_0.6.5 txdbmaker_1.3.0
#> [15] stringr_1.5.1 pwalign_1.3.0
#> [17] pkgconfig_2.0.3 crayon_1.5.3
#> [19] fastmap_1.2.0 magick_2.8.5
#> [21] dbplyr_2.5.0 XVector_0.47.0
#> [23] utf8_1.2.4 Rsamtools_2.23.0
#> [25] rmarkdown_2.28 UCSC.utils_1.3.0
#> [27] tinytex_0.53 bit_4.5.0
#> [29] xfun_0.48 zlibbioc_1.53.0
#> [31] cachem_1.1.0 jsonlite_1.8.9
#> [33] progress_1.2.3 blob_1.2.4
#> [35] highr_0.11 DelayedArray_0.33.0
#> [37] BiocParallel_1.41.0 jpeg_0.1-10
#> [39] prettyunits_1.2.0 R6_2.5.1
#> [41] VariantAnnotation_1.53.0 bslib_0.8.0
#> [43] stringi_1.8.4 RColorBrewer_1.1-3
#> [45] jquerylib_0.1.4 Rcpp_1.0.13
#> [47] bookdown_0.41 SummarizedExperiment_1.37.0
#> [49] knitr_1.48 SGSeq_1.41.0
#> [51] igraph_2.1.1 tidyselect_1.2.1
#> [53] Matrix_1.7-1 abind_1.4-8
#> [55] yaml_2.3.10 codetools_0.2-20
#> [57] RUnit_0.4.33 hwriter_1.3.2.1
#> [59] curl_5.2.3 tibble_3.2.1
#> [61] lattice_0.22-6 ShortRead_1.65.0
#> [63] KEGGREST_1.47.0 evaluate_1.0.1
#> [65] BiocFileCache_2.15.0 xml2_1.3.6
#> [67] Biostrings_2.75.0 filelock_1.0.3
#> [69] pillar_1.9.0 BiocManager_1.30.25
#> [71] MatrixGenerics_1.19.0 generics_0.1.3
#> [73] RCurl_1.98-1.16 hms_1.1.3
#> [75] GenomicFiles_1.43.0 glue_1.8.0
#> [77] tools_4.5.0 interp_1.1-6
#> [79] BiocIO_1.17.0 BSgenome_1.75.0
#> [81] locfit_1.5-9.10 GenomicAlignments_1.43.0
#> [83] XML_3.99-0.17 grid_4.5.0
#> [85] latticeExtra_0.6-30 GenomeInfoDbData_1.2.13
#> [87] restfulr_0.0.15 cli_3.6.3
#> [89] rappdirs_0.3.3 fansi_1.0.6
#> [91] S4Arrays_1.7.0 dplyr_1.1.4
#> [93] sass_0.4.9 digest_0.6.37
#> [95] SparseArray_1.7.0 rjson_0.2.23
#> [97] memoise_2.0.1 htmltools_0.5.8.1
#> [99] lifecycle_1.0.4 httr_1.4.7
#> [101] statmod_1.5.0 bit64_4.5.2