Source: https://github.com/markziemann/mesangial_rageko
Looking at doing comparison of control and RAGE KO Mouse mesangial cells RNA-seq and miRNA-seq.
RAGE KO data sequenced at Baker. 61 nt length.
WT data sequenced at BGI with 91 nt length.
Quality trimming with skewer 0.2.2 to ensure 3’ end base quality > 20.
GENECODE mouse transcriptome v30 was downloaded and indexed.
Kallisto 0.46.0 was used for mapping paired end reads to the transcriptome.
Count data was loaded into R, followed by DE analysis with DESeq2.
Enrichment analysis was performed with mitch, using REACTOME pathways.
OS, R and package versions are shown in the session info section at the end of this report.
suppressPackageStartupMessages({
library("reshape2")
library("DESeq2")
library("gplots")
library("mitch")
library("dplyr")
library("biomaRt")
library("kableExtra")
})
tmp <- read.table("3col.tsv.gz",header=F)
x <- as.matrix(acast(tmp, V2~V1, value.var="V3", fun.aggregate = sum))
x <- as.data.frame(x)
txinfo <- data.frame(do.call(rbind,strsplit(rownames(x),"\\|")))
colnames(txinfo) <- c("tx","GeneID","altGeneID","altTxID",
"txname","genename","length","biotype")
x$GeneID <- paste(txinfo$GeneID , txinfo$genename)
# write out transcript table
x2 <- x
x2$GeneID = NULL
xrpm <- apply(x2, 2, function(x){x/sum(x,na.rm=T)}) * 1000000
colnames(xrpm) <- gsub("$","_RPM",colnames(xrpm))
colnames(x2) <- gsub("$","_counts",colnames(x2))
x2 <- cbind(xrpm,x2)
write.table(x2,file="Tx_RPM_count_table_biotype.tsv")
xx <- aggregate(. ~ GeneID, x , sum)
colnames(xx) <- gsub("mRNA_","",colnames(xx))
rownames(xx) <- xx$GeneID
xx$GeneID <- NULL
xx <- round(xx)
dim(xx)
## [1] 56691 12
head(xx)
## rageko_ctrl_1 rageko_ctrl_2 rageko_ctrl_3
## ENSMUSG00000000001.5 Gnai3 2856 2448 2865
## ENSMUSG00000000003.16 Pbsn 0 0 0
## ENSMUSG00000000028.16 Cdc45 124 95 112
## ENSMUSG00000000031.18 H19 1 0 1
## ENSMUSG00000000037.18 Scml2 9 7 12
## ENSMUSG00000000049.12 Apoh 4 2 4
## rageko_TGF_1 rageko_TGF_2 rageko_TGF_3 wt_ctrl_1_R1
## ENSMUSG00000000001.5 Gnai3 2279 2047 2717 727
## ENSMUSG00000000003.16 Pbsn 0 0 0 0
## ENSMUSG00000000028.16 Cdc45 102 89 114 28
## ENSMUSG00000000031.18 H19 0 0 0 1452
## ENSMUSG00000000037.18 Scml2 7 6 9 0
## ENSMUSG00000000049.12 Apoh 0 1 0 4
## wt_ctrl_2_R1 wt_ctrl_3_R1 wt_TGF_1_R1 wt_TGF_2_R1
## ENSMUSG00000000001.5 Gnai3 780 621 1127 1150
## ENSMUSG00000000003.16 Pbsn 0 0 0 0
## ENSMUSG00000000028.16 Cdc45 29 31 49 52
## ENSMUSG00000000031.18 H19 1099 1065 11 6
## ENSMUSG00000000037.18 Scml2 6 1 1 3
## ENSMUSG00000000049.12 Apoh 4 5 0 0
## wt_TGF_3_R1
## ENSMUSG00000000001.5 Gnai3 1136
## ENSMUSG00000000003.16 Pbsn 0
## ENSMUSG00000000028.16 Cdc45 52
## ENSMUSG00000000031.18 H19 6
## ENSMUSG00000000037.18 Scml2 2
## ENSMUSG00000000049.12 Apoh 0
write.table(x=xx,file="mRNA_count_matrix.tsv",sep="\t")
Now quantify the gene biotypes. Starting with all samples.
txinfo2 <- txinfo
txinfo2$gene <- paste(txinfo2$GeneID , txinfo2$genename)
biotype <- unique(txinfo2[,c("gene","biotype")])
xxbiotype <- merge(xx,biotype,by.x=0,by.y="gene")
xxbiotype$Row.names=NULL
xxbiotype <- aggregate(. ~ biotype,xxbiotype,sum)
rownames(xxbiotype) <- xxbiotype$biotype
xxbiotype$biotype = NULL
# remove biotypes with very few reads
xxbiotype <- xxbiotype[which(rowMeans(xxbiotype)>10),]
# keep top 10 classes for chart
n=10
xxbiotype <- xxbiotype[tail(order(rowMeans(xxbiotype)),n),]
# make percent
xxbiotype2 <- apply(xxbiotype, 2, function(x){x*100/sum(x,na.rm=T)})
# create color palette:
library(RColorBrewer)
pal <- brewer.pal(n, "Paired")
par(mfrow=c(1,2))
par(mar=c(10,5,3,1))
barplot(as.matrix(xxbiotype2),col=pal,las=2)
plot.new()
legend("topleft", legend=rev(rownames(xxbiotype2)),fill=rev(pal),bg="white",cex=1)
#restore default setting
par(mar=c(5.1, 4.1, 4.1, 2.1))
Now with just the samples not treated with TGFb.
txinfo2 <- txinfo
txinfo2$gene <- paste(txinfo2$GeneID , txinfo2$genename)
biotype <- unique(txinfo2[,c("gene","biotype")])
xxbiotype <- merge(xx,biotype,by.x=0,by.y="gene")
xxbiotype$Row.names=NULL
xxbiotype <- aggregate(. ~ biotype,xxbiotype,sum)
rownames(xxbiotype) <- xxbiotype$biotype
xxbiotype$biotype = NULL
# remove TGFb samples
xxbiotype <- xxbiotype[,grep("TGF",colnames(xxbiotype),invert=TRUE)]
# remove biotypes with very few reads
xxbiotype <- xxbiotype[which(rowMeans(xxbiotype)>10),]
# keep top 10 classes for chart
n=10
xxbiotype <- xxbiotype[tail(order(rowMeans(xxbiotype)),n),]
# make percent
xxbiotype2 <- apply(xxbiotype, 2, function(x){x*100/sum(x,na.rm=T)})
# create color palette:
library(RColorBrewer)
pal <- brewer.pal(n, "Paired")
par(mfrow=c(1,2))
par(mar=c(10,5,3,1))
barplot(as.matrix(xxbiotype2),col=pal,las=2)
plot.new()
legend("topleft", legend=rev(rownames(xxbiotype2)),fill=rev(pal),bg="white",cex=1)
#restore default setting
par(mar=c(5.1, 4.1, 4.1, 2.1))
#colnames(xx) <- sapply(strsplit(colnames(xx),"_"),"[[",1)
mysamples <- colnames(xx)
ko <- as.numeric(grepl("ko",mysamples))
tgf <- as.numeric(grepl("TGF",mysamples))
ss <- data.frame(mysamples,ko,tgf)
rownames(ss) <- mysamples
ss
## mysamples ko tgf
## rageko_ctrl_1 rageko_ctrl_1 1 0
## rageko_ctrl_2 rageko_ctrl_2 1 0
## rageko_ctrl_3 rageko_ctrl_3 1 0
## rageko_TGF_1 rageko_TGF_1 1 1
## rageko_TGF_2 rageko_TGF_2 1 1
## rageko_TGF_3 rageko_TGF_3 1 1
## wt_ctrl_1_R1 wt_ctrl_1_R1 0 0
## wt_ctrl_2_R1 wt_ctrl_2_R1 0 0
## wt_ctrl_3_R1 wt_ctrl_3_R1 0 0
## wt_TGF_1_R1 wt_TGF_1_R1 0 1
## wt_TGF_2_R1 wt_TGF_2_R1 0 1
## wt_TGF_3_R1 wt_TGF_3_R1 0 1
TODO: rRNA carryover.
par(mar=c(5,8,3,1))
barplot(colSums(xx),horiz=TRUE,las=1,xlab="num reads",col=ss$cols)
grid()
And correlation heat map.
par(mar=c(5,5,3,3))
mds <- cmdscale(dist(t(xx)))
plot(mds, xlab="Coordinate 1", ylab="Coordinate 2",
type = "p",bty="n", cex=4 )
text(mds, labels=rownames(mds) ,col="black")
heatmap.2(cor(xx),trace="n",main="Pearson correlation heatmap",margin=c(8,8),cexRow=0.7,cexCol=0.7)
Compare wt control to RAGE KO control.
maplot <- function(de,contrast_name) {
de <- de[which(!is.na(de$padj)),]
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=1)
points(log2(sig$baseMean),sig$log2FoldChange,
pch=19, cex=0.5, col="red")
mtext(SUBHEADER,cex = 1)
}
make_volcano <- function(de,name) {
de <- de[which(!is.na(de$padj)),]
de$pvalue[which(de$pvalue==0)] <- 1e-320
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$pval),cex=0.5,pch=19,col="darkgray",
main=name, xlab="log2 FC", ylab="-log10 pval")
mtext(HEADER)
grid()
points(sig$log2FoldChange,-log10(sig$pval),cex=0.5,pch=19,col="red")
}
ss1 <- subset(ss,tgf==0)
xx1 <- xx[,which(colnames(xx) %in% rownames(ss1) )]
dim(xx1)
## [1] 56691 6
xx1 <- xx1[which(rowMeans(xx1)>10),]
dim(xx1)
## [1] 16920 6
dds <- DESeqDataSetFromMatrix(countData = xx1 , colData = ss1 , design = ~ ko )
## converting counts to integer mode
## the design formula contains one or more numeric variables with integer values,
## specifying a model with increasing fold change for higher values.
## did you mean for this to be a factor? if so, first convert
## this variable to a factor using the factor() function
res <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
dge[1:20,1:6] %>% kbl(caption = "Top gene expression differences") %>%
kable_paper("hover", full_width = F)
baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
---|---|---|---|---|---|---|
ENSMUSG00000001020.9 S100a4 | 4192.9386 | 5.170496 | 0.1377328 | 37.54005 | 0 | 0 |
ENSMUSG00000020044.14 Timp3 | 10285.4824 | 2.798266 | 0.0714929 | 39.14049 | 0 | 0 |
ENSMUSG00000020256.15 Aldh1l2 | 3821.4378 | 5.927939 | 0.1206065 | 49.15108 | 0 | 0 |
ENSMUSG00000022037.16 Clu | 4747.7435 | -7.309490 | 0.1278150 | -57.18804 | 0 | 0 |
ENSMUSG00000023036.15 Pcdhgc4 | 1683.2516 | -5.431429 | 0.1290685 | -42.08175 | 0 | 0 |
ENSMUSG00000027204.14 Fbn1 | 4626.3922 | 4.515813 | 0.0992572 | 45.49608 | 0 | 0 |
ENSMUSG00000029661.17 Col1a2 | 55636.8192 | 11.118047 | 0.1666730 | 66.70573 | 0 | 0 |
ENSMUSG00000029810.16 Tmem176b | 2962.7398 | -3.721489 | 0.0945140 | -39.37500 | 0 | 0 |
ENSMUSG00000035202.9 Lars2 | 440443.9343 | 7.995812 | 0.1385554 | 57.70843 | 0 | 0 |
ENSMUSG00000037071.4 Scd1 | 6218.4972 | 3.718488 | 0.0894171 | 41.58586 | 0 | 0 |
ENSMUSG00000040152.9 Thbs1 | 13191.7357 | 4.102414 | 0.0771886 | 53.14789 | 0 | 0 |
ENSMUSG00000069833.14 Ahnak | 36337.4655 | 3.850237 | 0.0825342 | 46.65023 | 0 | 0 |
ENSMUSG00000089774.3 Slc5a3 | 3102.6441 | 4.271342 | 0.1133767 | 37.67390 | 0 | 0 |
ENSMUSG00000101249.2 Gm29216 | 134092.2407 | -4.575209 | 0.0892471 | -51.26454 | 0 | 0 |
ENSMUSG00000119584.1 Rn18s-rs5 | 1624546.6620 | 8.192274 | 0.1498508 | 54.66955 | 0 | 0 |
ENSMUSG00000006205.14 Htra1 | 2118.8069 | 4.111353 | 0.1106220 | 37.16579 | 0 | 0 |
ENSMUSG00000002266.18 Zim1 | 4852.0107 | -3.413030 | 0.0930200 | -36.69137 | 0 | 0 |
ENSMUSG00000006519.12 Cyba | 1888.1110 | -3.138444 | 0.0904406 | -34.70173 | 0 | 0 |
ENSMUSG00000023904.11 Hcfc1r1 | 3021.9117 | -2.955816 | 0.0856391 | -34.51480 | 0 | 0 |
ENSMUSG00000031980.11 Agt | 922.7918 | -4.844777 | 0.1423726 | -34.02887 | 0 | 0 |
d1up <- rownames(subset(dge,padj <= 0.05 & log2FoldChange > 0))
d1dn <- rownames(subset(dge,padj <= 0.05 & log2FoldChange < 0))
txinfo$GeneID <- sapply(strsplit(txinfo$GeneID,"\\."),"[[",1)
g <- sapply(strsplit(rownames(dge),"\\."),"[[",1)
dge$biotype <- sapply(g, function(x) {
paste(unique(txinfo[txinfo$GeneID==x,"biotype"]),collapse=",") }
)
write.table(dge,file="mRNA_DE_ctrl_vs_rageko.tsv",quote=FALSE,sep="\t")
maplot(dge,"Cont1: Effect of RAGE KO")
make_volcano(dge,"Cont1: Effect of RAGE KO")
#agerdat <- assay(vsd)[grep("Ager",rownames(assay(vsd))),]
#par(mar=c(5,8,3,1))
#barplot(agerdat,horiz=TRUE,las=1,xlab="Normalised Ager expression",xlim=c(0,10))
Make a heatmap of the top 30 genes.
rpm1 <- apply(xx1, 2, function(x){x/sum(x,na.rm=T)}) * 1000000
agerdat <- rpm1[grep("Ager",rownames(rpm1)),]
par(mar=c(5,8,3,1))
barplot(agerdat,horiz=TRUE,las=1,xlab="Normalised Ager expression",xlim=c(0,16))
colfunc <- colorRampPalette(c("blue", "white", "red"))
rpm2 <- rpm1[which(rownames(rpm1) %in% rownames(head(dge,30))),]
heatmap.2(as.matrix(rpm2),margin=c(8, 22),cexRow=0.85,trace="none",
cexCol=0.9,col=colfunc(20),scale="row")
rownames(rpm2) <- sapply(strsplit(rownames(rpm2)," "),"[[",2)
heatmap.2(as.matrix(rpm2),margin=c(8, 22),cexRow=0.85,trace="none",
cexCol=0.9,col=colfunc(20),scale="row")
First need ortholog table.
Downloaded separately from the ensembl biomart website because the R package code is timing out.
orth <- read.table("../../ref/mart_export.txt",sep="\t",header=TRUE)
orth <- orth[which(orth$Human.gene.name != ""),]
orth <- orth[,c("Gene.stable.ID","Human.gene.name")]
orth <- unique(orth)
rownames(dge) <- sapply(strsplit(rownames(dge),"\\."),"[[",1)
Run mitch.
See the other report called “mrna_analysis_mitch.html” for the detailed analysis.
#download.file("https://reactome.org/download/current/ReactomePathways.gmt.zip", destfile="ReactomePathways.gmt.zip")
#unzip("ReactomePathways.gmt.zip")
# downloaded 19th August 2022
genesets <- gmt_import("ReactomePathways.gmt")
y <- mitch_import(dge, DEtype="deseq2",geneTable=orth)
## The input is a single dataframe; one contrast only. Converting
## it to a list for you.
## Note: Mean no. genes in input = 16920
## Note: no. genes in output = 12290
## Note: estimated proportion of input genes in output = 0.726
head(y)
## x
## A4GALT 0.5563104
## AAAS -9.3505565
## AACS 17.4368054
## AAGAB 0.6211840
## AAK1 10.8699941
## AAMDC -5.4631631
res <- mitch_calc(y, genesets, priority="effect",cores=16)
## Note: Enrichments with large effect sizes may not be
## statistically significant.
if (! file.exists("mrna_analysis_mitch.html")) {
mitch_report(res, "mrna_analysis_mitch.html")
}
head(res$enrichment_result,20) %>%
kbl(row.names=FALSE,caption = "Top pathways") %>%
kable_paper("hover", full_width = F)
set | setSize | pANOVA | s.dist | p.adjustANOVA |
---|---|---|---|---|
tRNA processing in the mitochondrion | 18 | 0.0000003 | -0.7014613 | 0.0000088 |
Cholesterol biosynthesis | 26 | 0.0000000 | 0.6917369 | 0.0000000 |
rRNA processing in the mitochondrion | 21 | 0.0000009 | -0.6199170 | 0.0000254 |
MET activates PTK2 signaling | 12 | 0.0006111 | 0.5713064 | 0.0081165 |
ALK mutants bind TKIs | 10 | 0.0021024 | 0.5616612 | 0.0206912 |
Laminin interactions | 18 | 0.0000371 | 0.5616037 | 0.0007146 |
Apoptosis induced DNA fragmentation | 10 | 0.0050416 | 0.5121824 | 0.0355514 |
Membrane binding and targetting of GAG proteins | 14 | 0.0012044 | -0.4998371 | 0.0142897 |
Synthesis And Processing Of GAG, GAGPOL Polyproteins | 14 | 0.0012044 | -0.4998371 | 0.0142897 |
Formation of a pool of free 40S subunits | 79 | 0.0000000 | -0.4903578 | 0.0000000 |
Peptide chain elongation | 67 | 0.0000000 | -0.4901281 | 0.0000000 |
Syndecan interactions | 19 | 0.0002684 | 0.4829573 | 0.0039479 |
Eukaryotic Translation Termination | 71 | 0.0000000 | -0.4766048 | 0.0000000 |
Formation of the ternary complex, and subsequently, the 43S complex | 47 | 0.0000000 | -0.4742580 | 0.0000008 |
Viral mRNA Translation | 67 | 0.0000000 | -0.4739731 | 0.0000000 |
Telomere C-strand synthesis initiation | 13 | 0.0033769 | -0.4695773 | 0.0281679 |
Mitochondrial translation initiation | 88 | 0.0000000 | -0.4695169 | 0.0000000 |
Mitochondrial translation termination | 88 | 0.0000000 | -0.4647096 | 0.0000000 |
Regulation of cholesterol biosynthesis by SREBP (SREBF) | 54 | 0.0000000 | 0.4634064 | 0.0000002 |
Eukaryotic Translation Elongation | 70 | 0.0000000 | -0.4553718 | 0.0000000 |
write.table(res$enrichment_result,file="mitch_pathways.tsv",sep="\t")
We have a list of mir-214 targets.
mir214_targ <- readLines("miR214_targets.txt")
mir214_targ <- list(mir214_targ)
names(mir214_targ) <- "mir214_targ"
res2 <- mitch_calc(y, mir214_targ, priority="effect",cores=16)
## Note: Enrichments with large effect sizes may not be
## statistically significant.
res2$enrichment_result
## set setSize pANOVA s.dist p.adjustANOVA
## 1 mir214_targ 3522 5.682073e-34 0.1395454 5.682073e-34
if (! file.exists("mrna_analysis_mitch_mir214.html")) {
mitch_report(res2, "mrna_analysis_mitch_mir214.html")
}
par(mfrow=c(2,1))
MIN=min(res2$ranked_profile[,1])
MAX=max(res2$ranked_profile[,1])
beeswarm(res2$detailed_sets$mir214_targ,pch=19,
lty="blank",cex=0.45,horiz=TRUE,ylim=c(MIN,MAX),
main="miR-214 targets",
xlab="position in rank (Ctrl vs RAGE KO)"
)
abline(v=0,lwd=2,lty=2,col="red")
abline(v=MIN,lwd=2,lty=2,col="black")
abline(v=MAX,lwd=2,lty=2,col="black")
grid()
par(mfrow=c(1,1))
It appears that they are upregulated in this contrast.
sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [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: Australia/Melbourne
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] pkgload_1.3.2 GGally_2.1.2
## [3] ggplot2_3.4.1 beeswarm_0.4.0
## [5] gtools_3.9.4 tibble_3.2.0
## [7] echarts4r_0.4.4 RColorBrewer_1.1-3
## [9] kableExtra_1.3.4 biomaRt_2.54.0
## [11] dplyr_1.1.0 mitch_1.10.0
## [13] gplots_3.1.3 DESeq2_1.38.3
## [15] SummarizedExperiment_1.28.0 Biobase_2.58.0
## [17] MatrixGenerics_1.10.0 matrixStats_0.63.0
## [19] GenomicRanges_1.50.2 GenomeInfoDb_1.34.6
## [21] IRanges_2.32.0 S4Vectors_0.36.1
## [23] BiocGenerics_0.44.0 reshape2_1.4.4
##
## loaded via a namespace (and not attached):
## [1] DBI_1.1.3 bitops_1.0-7 gridExtra_2.3
## [4] rlang_1.1.0 magrittr_2.0.3 compiler_4.3.1
## [7] RSQLite_2.3.0 systemfonts_1.0.4 png_0.1-8
## [10] vctrs_0.6.0 rvest_1.0.3 stringr_1.5.0
## [13] pkgconfig_2.0.3 crayon_1.5.2 fastmap_1.1.1
## [16] dbplyr_2.3.1 XVector_0.38.0 ellipsis_0.3.2
## [19] caTools_1.18.2 utf8_1.2.3 promises_1.2.0.1
## [22] rmarkdown_2.20 bit_4.0.5 xfun_0.37
## [25] zlibbioc_1.44.0 cachem_1.0.7 jsonlite_1.8.4
## [28] progress_1.2.2 blob_1.2.4 highr_0.10
## [31] later_1.3.0 DelayedArray_0.24.0 reshape_0.8.9
## [34] BiocParallel_1.32.5 prettyunits_1.1.1 parallel_4.3.1
## [37] R6_2.5.1 bslib_0.4.2 stringi_1.7.12
## [40] jquerylib_0.1.4 Rcpp_1.0.10 knitr_1.42
## [43] httpuv_1.6.9 Matrix_1.5-4.1 tidyselect_1.2.0
## [46] rstudioapi_0.14 yaml_2.3.7 codetools_0.2-19
## [49] curl_5.0.0 lattice_0.21-8 plyr_1.8.8
## [52] withr_2.5.0 shiny_1.7.4 KEGGREST_1.38.0
## [55] evaluate_0.20 BiocFileCache_2.6.0 xml2_1.3.3
## [58] Biostrings_2.66.0 filelock_1.0.2 pillar_1.8.1
## [61] KernSmooth_2.23-21 generics_0.1.3 RCurl_1.98-1.10
## [64] hms_1.1.2 munsell_0.5.0 scales_1.2.1
## [67] xtable_1.8-4 glue_1.6.2 tools_4.3.1
## [70] webshot_0.5.4 annotate_1.76.0 locfit_1.5-9.7
## [73] XML_3.99-0.13 grid_4.3.1 AnnotationDbi_1.60.0
## [76] colorspace_2.1-0 GenomeInfoDbData_1.2.9 cli_3.6.0
## [79] rappdirs_0.3.3 fansi_1.0.4 viridisLite_0.4.1
## [82] svglite_2.1.1 gtable_0.3.2 sass_0.4.5
## [85] digest_0.6.31 geneplotter_1.76.0 htmlwidgets_1.6.2
## [88] memoise_2.0.1 htmltools_0.5.4 lifecycle_1.0.3
## [91] httr_1.4.5 mime_0.12 bit64_4.0.5
## [94] MASS_7.3-58.3