Source: https://github.com/markziemann/mesangial_rageko

Introduction

Looking at doing comparison of control and RAGE KO Mouse mesangial cells RNA-seq and miRNA-seq.

Methods

  • 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")
})

Data

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))

Samplesheet

#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

QC

TODO: rRNA carryover.

par(mar=c(5,8,3,1))
barplot(colSums(xx),horiz=TRUE,las=1,xlab="num reads",col=ss$cols)
grid()

MDS plot

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)

Differential expression

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)
Top gene expression differences
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")

Pathway analysis

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)
Top pathways
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")

MiR-214

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

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