Example gene set analysis
Mark Ziemann & Kaumadi Wijesooriya
2022-01-18
Source: https://github.com/markziemann/SurveyEnrichmentMethods
Intro
Here we are performing an analysis of some gene expression data to demonstrate the difference between ORA and FCS methods and to highlight the differences caused by improper background gene set use.
The dataset being used is SRP128998 and we are comparing the cells grown in normal glucose condition (control) to the high glucose condition (case).
Data are obtained from http://dee2.io/
suppressPackageStartupMessages({
library("getDEE2")
library("DESeq2")
library("clusterProfiler")
library("mitch")
library("kableExtra")
library("eulerr")
})
dir.create("images")## Warning in dir.create("images"): 'images' already existsGet expression data and make an MDS plot
I’m using some RNA-seq data looking at the effect of hyperglycemia on hepatocytes.
name="SRP128998"
mdat<-getDEE2Metadata("hsapiens")
samplesheet <- mdat[grep("SRP128998",mdat$SRP_accession),]
samplesheet<-samplesheet[order(samplesheet$SRR_accession),]
samplesheet$trt<-as.factor(c(1,1,1,1,1,1,0,0,0,0,0,0))
samplesheet$VPA<-as.factor(c(0,0,0,1,1,1,0,0,0,1,1,1))
s1 <- subset(samplesheet,VPA==0)
s1 %>% kbl(caption = "sample sheet") %>% kable_paper("hover", full_width = F)| SRR_accession | QC_summary | SRX_accession | SRS_accession | SRP_accession | Sample_name | GEO_series | Library_name | trt | VPA | |
|---|---|---|---|---|---|---|---|---|---|---|
| 475895 | SRR6467479 | PASS | SRX3557428 | SRS2830728 | SRP128998 | GSM2932791 | GSE109140 | 1 | 0 | |
| 475896 | SRR6467480 | PASS | SRX3557429 | SRS2830730 | SRP128998 | GSM2932792 | GSE109140 | 1 | 0 | |
| 475897 | SRR6467481 | PASS | SRX3557430 | SRS2830729 | SRP128998 | GSM2932793 | GSE109140 | 1 | 0 | |
| 475901 | SRR6467485 | PASS | SRX3557434 | SRS2830733 | SRP128998 | GSM2932797 | GSE109140 | 0 | 0 | |
| 475902 | SRR6467486 | PASS | SRX3557435 | SRS2830734 | SRP128998 | GSM2932798 | GSE109140 | 0 | 0 | |
| 475903 | SRR6467487 | PASS | SRX3557436 | SRS2830735 | SRP128998 | GSM2932799 | GSE109140 | 0 | 0 |
w<-getDEE2("hsapiens",samplesheet$SRR_accession,metadata=mdat,legacy = TRUE)## For more information about DEE2 QC metrics, visit
## https://github.com/markziemann/dee2/blob/master/qc/qc_metrics.mdx<-Tx2Gene(w)
x<-x$Tx2Gene
# save the genetable for later
gt<-w$GeneInfo[,1,drop=FALSE]
gt$accession<-rownames(gt)
# counts
x1<-x[,which(colnames(x) %in% s1$SRR_accession)]Here show the number of genes in the annotation set, and those detected above the detection threshold.
# filter out lowly expressed genes
x1<-x1[which(rowSums(x1)/ncol(x1)>=(10)),]
nrow(x)## [1] 39297nrow(x1)## [1] 15635Now multidimensional scaling (MDS) plot to show the correlation between the datasets. If the control and case datasets are clustered separately, then it is likely that there will be many differentially expressed genes with FDR<0.05.
plot(cmdscale(dist(t(x1))), xlab="Coordinate 1", ylab="Coordinate 2", pch=19, col=s1$trt, main="MDS")
Differential expression
Now run DESeq2 for control vs case.
y <- DESeqDataSetFromMatrix(countData = round(x1), colData = s1, design = ~ trt)## converting counts to integer modey <- DESeq(y)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingde <- results(y)
de<-as.data.frame(de[order(de$pvalue),])
rownames(de)<-sapply(strsplit(rownames(de),"\\."),"[[",1)
head(de) %>% kbl() %>% kable_paper("hover", full_width = F)| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
|---|---|---|---|---|---|---|
| ENSG00000145050 | 5839.731 | -2.753692 | 0.1518338 | -18.13623 | 0 | 0 |
| ENSG00000149131 | 1346.633 | 2.161115 | 0.1427218 | 15.14215 | 0 | 0 |
| ENSG00000044574 | 124889.027 | -2.033391 | 0.1343040 | -15.14021 | 0 | 0 |
| ENSG00000128228 | 1676.368 | -2.836358 | 0.1895729 | -14.96183 | 0 | 0 |
| ENSG00000179218 | 78785.663 | -2.227516 | 0.1586383 | -14.04148 | 0 | 0 |
| ENSG00000090520 | 6751.044 | -2.138112 | 0.1538129 | -13.90073 | 0 | 0 |
Now let’s have a look at some of the charts showing differential expression. In particular, an MA plot and volcano plot.
maplot <- function(de,contrast_name) {
sig <-subset(de, padj < 0.05 )
up <-rownames(subset(de, padj < 0.05 & log2FoldChange > 0))
dn <-rownames(subset(de, padj < 0.05 & log2FoldChange < 0))
GENESUP <- length(up)
GENESDN <- length(dn)
DET=nrow(de)
SUBHEADER = paste(GENESUP, "up, ", GENESDN, "down", DET, "detected")
ns <-subset(de, padj > 0.05 )
plot(log2(de$baseMean),de$log2FoldChange,
xlab="log2 basemean", ylab="log2 foldchange",
pch=19, cex=0.5, col="dark gray",
main=contrast_name, cex.main=0.7)
points(log2(sig$baseMean),sig$log2FoldChange,
pch=19, cex=0.5, col="red")
mtext(SUBHEADER,cex = 0.7)
}
make_volcano <- function(de,name) {
sig <- subset(de,padj<0.05)
N_SIG=nrow(sig)
N_UP=nrow(subset(sig,log2FoldChange>0))
N_DN=nrow(subset(sig,log2FoldChange<0))
DET=nrow(de)
HEADER=paste(N_SIG,"@5%FDR,", N_UP, "up", N_DN, "dn", DET, "detected")
plot(de$log2FoldChange,-log10(de$padj),cex=0.5,pch=19,col="darkgray",
main=name, xlab="log2 FC", ylab="-log10 pval", xlim=c(-6,6))
mtext(HEADER)
grid()
points(sig$log2FoldChange,-log10(sig$padj),cex=0.5,pch=19,col="red")
}
maplot(de,name)
make_volcano(de,name)
Gene sets from Reactome
In order to perform gene set analysis, we need some gene sets.
if (! file.exists("ReactomePathways.gmt")) {
download.file("https://reactome.org/download/current/ReactomePathways.gmt.zip",
destfile="ReactomePathways.gmt.zip")
unzip("ReactomePathways.gmt.zip")
}
genesets<-gmt_import("ReactomePathways.gmt")FCS with Mitch
Mitch uses rank-ANOVA statistics for enrichment detection.
m <- mitch_import(de,DEtype = "DEseq2", geneTable = gt)## The input is a single dataframe; one contrast only. Converting
## it to a list for you.## Note: Mean no. genes in input = 15635## Note: no. genes in output = 14533## Note: estimated proportion of input genes in output = 0.93mres <- mitch_calc(m,genesets = genesets)## Note: When prioritising by significance (ie: small
## p-values), large effect sizes might be missed.m_up <- subset(mres$enrichment_result,p.adjustANOVA<0.05 & s.dist > 0)[,1]
m_dn <- subset(mres$enrichment_result,p.adjustANOVA<0.05 & s.dist < 0)[,1]
message(paste("Number of up-regulated pathways:",length(m_up) ))## Number of up-regulated pathways: 96message(paste("Number of down-regulated pathways:",length(m_dn) ))## Number of down-regulated pathways: 317head(mres$enrichment_result,10) %>% kbl() %>% kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 155 | Cellular responses to stimuli | 641 | 0 | -0.2223004 | 0 |
| 156 | Cellular responses to stress | 633 | 0 | -0.2167957 | 0 |
| 509 | Infectious disease | 637 | 0 | -0.2062708 | 0 |
| 1038 | SRP-dependent cotranslational protein targeting to membrane | 109 | 0 | -0.4728714 | 0 |
| 74 | Asparagine N-linked glycosylation | 252 | 0 | -0.3060943 | 0 |
| 87 | Axon guidance | 424 | 0 | -0.2324048 | 0 |
| 154 | Cellular response to starvation | 149 | 0 | -0.3633278 | 0 |
| 630 | Metabolism of proteins | 1550 | 0 | -0.1176239 | 0 |
| 708 | Nervous system development | 444 | 0 | -0.2086060 | 0 |
| 573 | L13a-mediated translational silencing of Ceruloplasmin expression | 108 | 0 | -0.4039887 | 0 |
m_up_nom <- subset(mres$enrichment_result,pANOVA<0.05 & s.dist > 0)[,1]
m_dn_nom <- subset(mres$enrichment_result,pANOVA<0.05 & s.dist < 0)[,1]ORA with clusterprofiler
Clusterprofiler uses a hypergeometric test. Firstly I will conduct the analysis separately for up and down regulated genes and with the correct backgound (as intended by the developers).
genesets2 <- read.gmt("ReactomePathways.gmt")
de_up <- rownames(subset(de,log2FoldChange>0,padj<0.05))
de_up <- unique(gt[which(rownames(gt) %in% de_up),1])
de_dn <- rownames(subset(de,log2FoldChange<0,padj<0.05))
de_dn <- unique(gt[which(rownames(gt) %in% de_dn),1])
de_bg <- rownames(de)
de_bg <- unique(gt[which(rownames(gt) %in% de_bg),1])
c_up <- as.data.frame(enricher(gene = de_up, universe = de_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="fdr"))
c_up <- rownames(subset(c_up, p.adjust < 0.05))
c_dn <- as.data.frame(enricher(gene = de_dn, universe = de_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="fdr"))
c_dn <- rownames(subset(c_dn, p.adjust < 0.05))
c_up_nom <- as.data.frame(enricher(gene = de_up, universe = de_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="none" ))
c_up_nom <- rownames(subset(c_up_nom, pvalue < 0.05))
c_dn_nom <- as.data.frame(enricher(gene = de_dn, universe = de_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="none"))
c_dn_nom <- rownames(subset(c_dn_nom, pvalue < 0.05))Now performing ORA with clusterprofiler with whole genome background list
wg_bg <- w$GeneInfo$GeneSymbol
f_up <- as.data.frame(enricher(gene = de_up, universe = wg_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="fdr"))
f_up <- rownames(subset(f_up, p.adjust < 0.05))
f_dn <- as.data.frame(enricher(gene = de_dn, universe = wg_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="fdr"))
f_dn <- rownames(subset(f_dn, p.adjust < 0.05))
f_up_nom <- as.data.frame(enricher(gene = de_up, universe = wg_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="none"))
f_up_nom <- rownames(subset(f_up_nom, pvalue < 0.05))
f_dn_nom <- as.data.frame(enricher(gene = de_dn, universe = wg_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="none"))
f_dn_nom <- rownames(subset(f_dn_nom, pvalue < 0.05))
f_de_nom <- union(f_up_nom,f_dn_nom)Bar chart of significance
Here the idea is to classify the pathways as 1 not significant, 2 nominally significant, or 3 FDR significant
c_up_df <- as.data.frame(enricher(gene = de_up, universe = de_bg, maxGSSize = 5000, TERM2GENE = genesets2,
pAdjustMethod="fdr" ,qvalueCutoff=1,pvalueCutoff=1))
n_c_up_ns <- nrow( subset(c_up_df,pvalue>0.05) )
n_c_up_nom <- nrow( subset(c_up_df,pvalue<0.05 ) )
n_c_up_fdr <- nrow( subset(c_up_df,p.adjust<0.05) )
c_dn_df <- as.data.frame(enricher(gene = de_dn, universe = de_bg, maxGSSize = 5000, TERM2GENE = genesets2,
pAdjustMethod="fdr" ,qvalueCutoff=1,pvalueCutoff=1))
n_c_dn_ns <- nrow( subset(c_dn_df,pvalue>0.05) )
n_c_dn_nom <- nrow( subset(c_dn_df,pvalue<0.05 ) )
n_c_dn_fdr <- nrow( subset(c_dn_df,p.adjust<0.05) )
n_m_up <- length(m_up)
n_m_dn <- length(m_dn)
n_m_up_nom <- length(m_up_nom)
n_m_dn_nom <- length(m_dn_nom)
par(mar=c(5,10,5,2))
ngenes <- c("ORA up fdr"=n_c_up_fdr,"ORA up nom"=n_c_up_nom,
"ORA dn fdr"=n_c_dn_fdr,"ORA dn nom"=n_c_dn_nom,
"FCS up fdr"=n_m_up,"FCS up nom"=n_m_up_nom,
"FCS dn fdr"=n_m_dn,"FCS dn nom"=n_m_dn_nom )
barplot(ngenes, horiz=TRUE,las=1)
c_up <- subset(c_up_df,p.adjust<0.05)$ID
c_dn <- subset(c_dn_df,p.adjust<0.05)$ID
c_up_nom <- subset(c_up_df,pvalue<0.05)$ID
c_dn_nom <- subset(c_dn_df,pvalue<0.05)$ID
n_c_fdr <- length(union(c_dn,c_up))
n_c_nom <- length(union(c_dn_nom,c_up_nom))
n_m_fdr <- length(subset(mres$enrichment_result,p.adjustANOVA<0.05 )[,1])
n_m_nom <- length(subset(mres$enrichment_result,pANOVA<0.05 )[,1])
par(mar=c(5,5,5,2))
ngenes <- c("ORA FDR<0.05"=n_c_fdr,"ORA p<0.05"=n_c_nom,"FCS FDR<0.05"=n_m_fdr,"FCS p<0.05"=n_m_nom)
barplot(ngenes,ylab="no. gene sets")
text((0:3*1.2)+0.7,ngenes-50,labels=ngenes,cex=1.1)
pdf("images/fdr1.pdf",width=4,height=3)
par(mar=c(5,7,2,2))
barplot(ngenes,xlab="no. gene sets",horiz=TRUE,las=1)
text( ngenes - 70, (0:3*1.2)+0.7,labels=ngenes,cex=1)
dev.off()## png
## 2Venn diagram comparison
Here I will omit combined analysis and include nominal p-value sets.
par(cex.main=0.5)
par(mar=c(2,2,2,2))
v1 <- list("FCS up"=m_up, "FCS dn"=m_dn,
"ORA up"=c_up,"ORA dn"=c_dn)
plot(euler(v1),quantities = TRUE, edges = "gray", main="FCS compared to ORA")
v0 <- list("FDR up"=m_up, "FDR dn"=m_dn,
"Nom up"=m_up_nom,"Nom dn"=m_dn_nom)
plot(euler(v0),quantities = TRUE, edges = "gray", main="Effect of FDR correction on FCS results")
v0 <- list("FDR up"=c_up, "FDR dn"=c_dn,
"Nom up"=c_up_nom,"Nom dn"=c_dn_nom)
plot(euler(v0),quantities = TRUE, edges = "gray", main="Effect of FDR correction on ORA results")
ora_nom <- union(c_up_nom,c_dn_nom)
ora_fdr <- union(c_up,c_dn)
fcs_nom <- union(m_up_nom,m_dn_nom)
fcs_fdr <- union(m_up,m_dn)
v3 <- list("ORA nom"=ora_nom, "ORA FDR"=ora_fdr,
"FCS nom"=fcs_nom,"FCS FDR"=fcs_fdr)
plot(euler(v3),quantities = TRUE, edges = "gray", main="Effect of FDR correction")
v2 <- list("ORA up"=c_up,"ORA dn"=c_dn,
"ORA* up"=f_up,"ORA* dn"=f_dn )
plot(euler(v2),quantities = TRUE, edges = "gray", main="Effect of inappropriate background* (whole genome)")
pdf("images/fcs_ora1.pdf",width=4,height=4)
plot(euler(v1),quantities = TRUE, edges = "gray", main="FCS vs ORA")
dev.off()## png
## 2pdf("images/orabg1.pdf",width=4,height=4)
plot(euler(v2),quantities = TRUE, edges = "gray", main="Effect of inappropriate background* (whole genome)")
dev.off()## png
## 2Jaccard calculation
# FCS vs ORA
cm <- length(intersect(c(c_up,c_dn), c(m_up,m_dn))) / length(union(c(c_up,c_dn), c(m_up,m_dn)))
#FCS fdr vs nom
fcs_fdr_nom <- length(intersect(c(fcs_nom), c(fcs_fdr))) / length(union(c(fcs_nom), c(fcs_fdr)))
#ORA fdr vs nom
ora_fdr_nom <- length(intersect(c(ora_nom), c(ora_fdr))) / length(union(c(ora_nom), c(ora_fdr)))
m_up <- gsub("^","up ",m_up)
m_dn <- gsub("^","dn ",m_dn)
m_de <- union(m_up,m_dn)
c_up <- gsub("^","up ",c_up)
c_dn <- gsub("^","dn ",c_dn)
c_de <- union(c_up,c_dn)
f_up <- gsub("^","up ",f_up)
f_dn <- gsub("^","dn ",f_dn)
f_de <- union(f_up,f_dn)
f_up_nom <- gsub("^","up ",f_up_nom)
f_dn_nom <- gsub("^","dn ",f_dn_nom)
f_de_nom <- union(f_up,f_dn_nom)
# ORA vs ORA*
cf <- length(intersect(c_de, f_de )) / length(union(c_de, f_de))
# ORA vs ORA*nom
cfn <- length(intersect(c_de, f_de_nom )) / length(union(c_de, f_de_nom))
dat <- c( "FCS vs ORA"=cm,
"FCS: FDR vs nominal"=fcs_fdr_nom,
"ORA: FDR vs nominal"= ora_fdr_nom,
"ORA vs ORA*"=cf,
"ORA vs ORA*nom"=cfn)
dat## FCS vs ORA FCS: FDR vs nominal ORA: FDR vs nominal ORA vs ORA*
## 0.6494118 0.7454874 0.5625000 0.4391785
## ORA vs ORA*nom
## 0.3781965saveRDS(dat,file = "ex1dat.rds")
par(mar=c(5,10,3,1))
barplot(rev(dat),xlab="Jaccard index",horiz = TRUE, las =1, xlim=c(0,.8) , main=name)
text( x=rev(dat)-0.05 , y= 1:length(rev(dat))*1.2-0.5, labels = signif(rev(dat),2))
pdf("images/jacbar1.pdf",width=4,height=3)
par(mar=c(5,9,3,1))
barplot(rev(dat),xlab="Jaccard index",horiz = TRUE, las =1, xlim=c(0,.8) , main=name)
text( x=rev(dat)-0.09 , y= 1:length(rev(dat))*1.2-0.5, labels = signif(rev(dat),2))
dev.off()## png
## 2png("images/jacbar1.png")
par(mar=c(5,10,3,1))
barplot(rev(dat),xlab="Jaccard index",horiz = TRUE, las =1, xlim=c(0,.75) , main=name)
text( x=rev(dat)+0.05 , y= 1:length(rev(dat))*1.2-0.5, labels = signif(rev(dat),2))
dev.off()## png
## 2Session information
sessionInfo()## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## locale:
## [1] LC_CTYPE=en_AU.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_AU.UTF-8 LC_COLLATE=en_AU.UTF-8
## [5] LC_MONETARY=en_AU.UTF-8 LC_MESSAGES=en_AU.UTF-8
## [7] LC_PAPER=en_AU.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] beeswarm_0.4.0 eulerr_6.1.1
## [3] kableExtra_1.3.4 mitch_1.4.1
## [5] clusterProfiler_4.0.5 DESeq2_1.32.0
## [7] SummarizedExperiment_1.22.0 Biobase_2.52.0
## [9] MatrixGenerics_1.4.3 matrixStats_0.61.0
## [11] GenomicRanges_1.44.0 GenomeInfoDb_1.28.4
## [13] IRanges_2.26.0 S4Vectors_0.30.2
## [15] BiocGenerics_0.38.0 getDEE2_1.2.0
##
## loaded via a namespace (and not attached):
## [1] shadowtext_0.1.1 fastmatch_1.1-3 systemfonts_1.0.3
## [4] plyr_1.8.6 igraph_1.2.11 lazyeval_0.2.2
## [7] polylabelr_0.2.0 splines_4.1.2 BiocParallel_1.26.2
## [10] ggplot2_3.3.5 digest_0.6.29 yulab.utils_0.0.4
## [13] htmltools_0.5.2 GOSemSim_2.18.1 viridis_0.6.2
## [16] GO.db_3.13.0 fansi_1.0.0 magrittr_2.0.1
## [19] memoise_2.0.1 Biostrings_2.60.2 annotate_1.70.0
## [22] graphlayouts_0.8.0 svglite_2.0.0 enrichplot_1.12.3
## [25] colorspace_2.0-2 rvest_1.0.2 blob_1.2.2
## [28] ggrepel_0.9.1 xfun_0.29 dplyr_1.0.7
## [31] crayon_1.4.2 RCurl_1.98-1.5 jsonlite_1.7.2
## [34] scatterpie_0.1.7 genefilter_1.74.1 survival_3.2-13
## [37] ape_5.6-1 glue_1.6.0 polyclip_1.10-0
## [40] gtable_0.3.0 zlibbioc_1.38.0 XVector_0.32.0
## [43] webshot_0.5.2 htm2txt_2.1.1 DelayedArray_0.18.0
## [46] scales_1.1.1 DOSE_3.18.3 DBI_1.1.2
## [49] GGally_2.1.2 Rcpp_1.0.7 viridisLite_0.4.0
## [52] xtable_1.8-4 gridGraphics_0.5-1 tidytree_0.3.7
## [55] bit_4.0.4 htmlwidgets_1.5.4 httr_1.4.2
## [58] fgsea_1.18.0 gplots_3.1.1 RColorBrewer_1.1-2
## [61] ellipsis_0.3.2 pkgconfig_2.0.3 reshape_0.8.8
## [64] XML_3.99-0.8 farver_2.1.0 sass_0.4.0
## [67] locfit_1.5-9.4 utf8_1.2.2 ggplotify_0.1.0
## [70] tidyselect_1.1.1 rlang_0.4.12 reshape2_1.4.4
## [73] later_1.3.0 AnnotationDbi_1.54.1 munsell_0.5.0
## [76] tools_4.1.2 cachem_1.0.6 downloader_0.4
## [79] generics_0.1.1 RSQLite_2.2.9 evaluate_0.14
## [82] stringr_1.4.0 fastmap_1.1.0 yaml_2.2.1
## [85] ggtree_3.0.4 knitr_1.37 bit64_4.0.5
## [88] tidygraph_1.2.0 caTools_1.18.2 purrr_0.3.4
## [91] KEGGREST_1.32.0 ggraph_2.0.5 nlme_3.1-153
## [94] mime_0.12 aplot_0.1.2 xml2_1.3.3
## [97] DO.db_2.9 rstudioapi_0.13 compiler_4.1.2
## [100] png_0.1-7 treeio_1.16.2 tibble_3.1.6
## [103] tweenr_1.0.2 geneplotter_1.70.0 bslib_0.3.1
## [106] stringi_1.7.6 highr_0.9 lattice_0.20-45
## [109] Matrix_1.4-0 vctrs_0.3.8 pillar_1.6.4
## [112] lifecycle_1.0.1 jquerylib_0.1.4 data.table_1.14.2
## [115] cowplot_1.1.1 bitops_1.0-7 httpuv_1.6.5
## [118] patchwork_1.1.1 qvalue_2.24.0 R6_2.5.1
## [121] promises_1.2.0.1 KernSmooth_2.23-20 echarts4r_0.4.3
## [124] gridExtra_2.3 gtools_3.9.2 MASS_7.3-54
## [127] assertthat_0.2.1 GenomeInfoDbData_1.2.6 grid_4.1.2
## [130] ggfun_0.0.4 tidyr_1.1.4 rmarkdown_2.11
## [133] ggforce_0.3.3 shiny_1.7.1